#### Abstract

An efficient numerical scheme for solving delay differential equations with a piecewise constant delay function is developed in this paper. The proposed approach is based on a hybrid of block-pulse functions and Taylor’s polynomials. The operational matrix of delay corresponding to the proposed hybrid functions is introduced. The sparsity of this matrix significantly reduces the computation time and memory requirement. The operational matrices of integration, delay, and product are employed to transform the problem under consideration into a system of algebraic equations. It is shown that the developed approach is also applicable to a special class of nonlinear piecewise constant delay differential equations. Several numerical experiments are examined to verify the validity and applicability of the presented technique.

#### 1. Introduction

Many problems arising in diverse areas of science and engineering can be described by delay differential equations (DDEs). Typical examples are transmission lines, population dynamics, chemical processes, mathematical biology, physics, electrodynamics, quantum mechanics, robotics, and communication networks [1–3]. Many research works have been devoted to the numerical treatments of DDEs involving constant delay [4–16]. To our knowledge, a few papers have been paid to the numerical investigation of delay differential equations with a piecewise constant delay function [17–21]. It should be mentioned that, except for some special cases, it is either difficult or impossible to obtain a closed-form solution to a constant delay system. Obviously, the situation becomes more complicated when time-delay is a piecewise constant function.

It is generally supposed that the delay function is either constant or continuous. However, in some actual situations, time-delay is a piecewise constant function. Most of the existing computational techniques are capable of providing satisfactory approximation for problems whose solutions are infinitely smooth or well-behaved. To deal with a problem whose solution is a nonsmooth function such as time-delay systems, they encounter some challenges. The computational procedures developed in [17, 18] are based on block-pulse functions. Among piecewise constant basis functions, block-pulse functions have a simple structure. Simplicity in operations is a distinctive feature of these functions [22]. As a result, they can be implemented without too much effort. Although these functions possess some nice properties such as disjointedness, orthogonality, and completeness, they do not provide accurate solutions for DDEs. The numerical results obtained by the methods implemented in [17, 18] are not in good agreement with the exact solution of the mentioned systems. It should be pointed out that the analytical solution of these systems cannot be produced solely either by continuous basis functions or by piecewise constant basis functions. The main difficulty arises because of the lack of smoothness of the solution associated with DDEs. To overcome these drawbacks, it motivates one to design a flexible framework to accurately model the mixed features of continuity and jumps that occur in the solution of the problem under consideration.

Recently, an effective computational technique for solving DDEs with a piecewise constant delay function has been introduced by Marzban and Shahsiah [19]. The developed method is based on a hybrid of block-pulse functions and Chebyshev polynomials. More recently, an efficient numerical scheme has been presented for solving linear piecewise constant delay systems [20]. The proposed approach is based on a hybrid of block-pulse functions and Legendre polynomials. The purpose of this work is to present a simple framework to numerically solve piecewise constant delay systems. The proposed procedure is based on a hybrid of block-pulse functions and Taylor’s polynomials. The structure of this kind of hybrid functions is much simpler than that of those available in the literature, particularly, the hybrid of block-pulse functions and Chebyshev polynomials. This type of hybrid functions allows one to significantly reduce the computation time as well as the memory usage. The proposed procedure is also applicable to nonlinear piecewise constant delay systems, but some modifications are required.

The paper is organized as follows. In Section 2, we explain the fundamental properties of the hybrid of block-pulse functions and Taylor’s polynomials. In Section 3, the operational matrix of delay associated with the proposed hybrid functions is constructed. Section 4 is dedicated to the problem statement and its description. Finally, in Section 5, various types of piecewise constant delay systems are investigated to evaluate the performance and computational efficiency of the suggested numerical scheme.

#### 2. Hybrid Functions

##### 2.1. Hybrid of Block-Pulse Functions and Taylor’s Polynomials

Hybrid functions , , , are defined on the interval as [14]where is the order of block-pulse functions. Here , , are the well-known Taylor’s polynomials of degree . A detailed explanation of the block-pulse functions can be found in [22]. It is worth noting that the mentioned hybrid functions constitute a semiorthogonal set on the interval .

##### 2.2. Function Approximation

A function can be represented in terms of the hybrid functions aswhere the coefficients , , , are calculated by the following formula: If the series given by (2) is truncated, then it can be written in the formwhereIn order to illustrate the usefulness and advantages of the presented framework over the conventional Taylor’s expansion, we investigate two cases as follows.

*Case 1. *Assume thatClearly, the function has two corner points at the points and . Accordingly, Taylor’s expansion of about each of the mentioned points does not exist. This function can be exactly approximated by the hybrid functions as follows: in whichThe elements of the vector are defined by

*Case 2. *Suppose thatObviously, this function is discontinuous at the points and . However, it can be efficiently approximated by the proposed hybrid functions. To express this function in terms of hybrid functions, we choose . In addition, Assume that is an arbitrary positive integer number greater than . With the use of (2), we have where the hybrid functions , , , are defined by Also, the associated coefficients are given byThe hybrid expansion of over the interval is presented by Similarly, the hybrid approximation of on the interval is described by Moreover, the hybrid expansion of over the interval is expressed by which is precisely Taylor’s expansion of the function about , up to degree . Define It is easily verified that The obtained result indicates that the maximum error converges to zero as the value of increases. Some important aspects of the proposed hybrid functions are clarified in [23].

##### 2.3. Operational Matrices of Integration and Product

The integration of the vector presented by (6) can also be expanded in terms of hybrid functions asin which is the operational matrix of integration corresponding to the hybrid functions. The structure of this matrix is given by [14]where and is defined by It should be noted that is the operational matrix of integration associated with the Taylor’s polynomials on the interval . Furthermore, if is a vector of order , then the expression can be represented in terms of hybrid functions asThe matrix is called the operational matrix of product [14].

#### 3. Derivation of the Operational Matrix of Delay

In this section, we shall derive the operational matrix of delay associated with the hybrid of block-pulse functions and Taylor’s polynomials. For this purpose, we employ a procedure analogous to the one developed in [19] for the hybrid of block-pulse functions and Chebyshev polynomials. The delayed vector can be represented in terms of as follows:in which is called the operational matrix of delay. In addition, is a piecewise constant function defined by where . It is further supposed that , , are arbitrary rational numbers in the interval and , . To construct , we first obtain the value of . Let and take as the smallest positive integer number such that the following conditions are satisfied: Assume that . Next, we define as the greatest common divisor of the integers and , and ; that is,Definewhere denotes the greatest integer value less than or equal to . It is worth noting that is selected in such a way that the number of subintervals is minimized. As a consequence, we obtain the following subintervals: For convenience, we use the notations defined by Clearly . Consequently, we have in which As a result, the problem is reduced to find the delay operational matrix for the delay function defined by where To construct the matrix , we first derive the matrix for , in such a way that the following relation is satisfied: By the properties of the hybrid functions, it is easy to verify that, for the case , the only functions with nonzero values are , for . Notice that Accordingly, can be written in terms of as where is the -dimensional identity matrix and denotes the Kronecker product [24]. In addition, is an matrix in which the only nonzero entry is equal to one and located at the th row and th column. Finally, if we express in terms of , we find

*Remark 1. *If , then is a zero matrix of order . It is worth noting that the structure of the operational matrix of delay associated with the proposed hybrid functions is exactly the same as the structure of the hybrid of block-pulse functions and Chebyshev polynomials.

#### 4. Problem Statement and Its Approximation

Consider the time-varying piecewise constant delay system described bywhere , , , , and are matrices of appropriate dimensions, is a constant specified vector, and , , are arbitrary rational numbers in . The objective is to find , , satisfying (41)–(43).

##### 4.1. Description of the Method

Let Assume that each and each of , , can be written in terms of hybrid functions as Consequentlywhere and are identity matrices of order and , respectively. Furthermore, and are vectors of order and , respectively, given by Alsowhere is a vector of order , defined by Similarly, , , and can be represented in terms of hybrid functions aswhere , , and are of dimensions , , and , respectively. The delayed vector can also be expressed in terms of hybrid functions asin which is the delay operational matrix presented by (25). As a resultwhere , , and can be determined in a manner analogous to the derivation process of the matrix described by (24). Furthermorein which is the operational matrix of integration described by (20). By integrating both sides of (41) with respect to and using (47), (49), (51), (52), and (54), we obtainFrom the preceding equation, we deduce that

#### 5. Simulation Results

In this section, five examples of varying complexity are investigated to evaluate the usefulness and effectiveness of the proposed computational technique. The correct choice of would greatly improve the efficiency and accuracy of the presented procedure.

##### 5.1. Example 1

Consider the following piecewise constant delay system [17]:The exact solution to this problem is given by

To apply the method developed in the current paper, we first choose the value of , with the use of (30). Because and , hence, we select . We also take . LetwhereBy expanding in terms of hybrid functions, we getWe also havewhere is the delay operational matrix given by (25). Integrating both sides of (57) with respect to and using (62)–(65) yieldin which is the -dimensional identity matrix and is the operational matrix of integration given by (20). With the use of (66) and (62), we would obtain the same value as the exact value of .

##### 5.2. Example 2

Consider the following two-dimensional time-varying piecewise constant delay system:The exact solutions to this problem are given as follows: To solve this problem by the proposed approach, we first determine the value of with the use of (30). Because and , consequently, we select . Also, we choose . LetExpanding , , , and in terms of hybrid functions impliesAlso, using (70), we obtainMoreoverThereforein which , , and can be obtained in a way similar to the construction method of the matrix defined by (24). By integrating both sides of (67) with respect to and using (68), (70), (71), (73), (74), (75), (76), and (77), we getBy solving the resulting system described by (78) and then using (73), we would obtain the same values as the exact values of and .

##### 5.3. Example 3

Consider the following time-delay system [3]:The exact solution to this problem is given by [3]To employ the proposed method, we first determine the suitable value of . Using (30), it follows that Assume that is an arbitrary positive integer number. LetSimilarlyBy integrating both sides of (79) with respect to and using (85)-(86), we conclude thatwhere is the -dimensional identity matrix. Define in which and denote the approximate solution determined by the present approach and the exact solution, respectively. In Table 1, the numerical results of the maximum error denoted by , for , and various values of are reported. The simulation results obtained by the current approach show a significant improvement in the accuracy of the solution compared to the method developed in [19].

##### 5.4. Example 4

As a more complicated problem, consider the following problem:The analytical solution to this problem is described byAlthough the above-mentioned problem is a nonlinear delay differential equation, the method developed in the present paper can easily be implemented. Because and , we take . Furthermore, we take . LetSimilarly where , , , and are vectors of order . Also, we have in which Using (99) givesIntegrating both sides of (89) with respect to and using (95), (97), (98), (99), (100), and (102) imply thatBy solving (103) and using (95), the exact value of will be obtained.

##### 5.5. Example 5

As the final example, we consider the time-varying multidelay system defined bywhere

The analytical solution to this problem is described by

Although the above-mentioned system involves multiple piecewise constant delays, the proposed procedure can be applied. To employ our approach, we select . By choosing , the exact solution of will be obtained.

#### 6. Conclusion

An efficient and flexible framework has been successfully developed for solving piecewise constant delay systems. The foundation of the proposed method is based on a hybrid of block-pulse functions and Taylor’s polynomials. The operational matrix of delay has been constructed. The operational matrices of integration, delay, and product associated with the hybrid functions were utilized to transform the problem under consideration into a system of algebraic equations. It is worth emphasizing that the exact solutions of Examples , , , and cannot be produced solely either by block-pulse functions or by Taylor’s polynomials. After determining the suitable value of , the number of subintervals, small values for , the degree of the Taylor’s polynomials, are needed to achieve a specific accuracy. The simulation results demonstrate the effectiveness of the suggested procedure. In addition, the proposed numerical scheme can be extended to a class of nonlinear piecewise constant delay systems, but some modifications are required.

#### Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this manuscript.