Analysis of Linear Piecewise Constant Delay Systems Using a Hybrid Numerical Scheme

An efficient computational technique for solving linear delay differential equations with a piecewise constant delay function is presented.Thenew approach is based on a hybrid of block-pulse functions and Legendre polynomials. A key feature of the proposed framework is the excellent representation of smooth and especially piecewise smooth functions. The operational matrices of delay, derivative, and product corresponding to the mentioned hybrid functions are implemented to transform the original problem into a system of algebraic equations. Illustrative examples are included to demonstrate the validity and applicability of the proposed numerical scheme.


Introduction
Delay differential equations (DDEs) naturally arise in diverse areas of science and engineering such as transmission lines, communication networks, biological models, population dynamics, and transportation systems [1,2]. So far, a large body of literature has been devoted to the theoretical aspects and numerical treatments of DDEs with constant delays [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. It is known that, except for some simple cases, it is either extremely difficult or impossible to analytically solve DDEs. Accordingly, a numerical algorithm has to be adopted in most cases. The situation becomes more complicated when the time-delay is a piecewise constant function. Owing to the nature of DDEs, none of the smooth basis functions is able to properly model the inherent behavior of this class of systems. This is due to the lack of smoothness of analytical solution associated with DDEs. It should be pointed out that the approximation of a piecewise smooth function by a finite number of smooth functions often fails to converge because of the existence of the well-known Gibbs phenomenon. Consequently, a suitable basis is required to accurately model the true locations of the switching points that occur in the exact solution of a delay differential equation. It is generally assumed that the delay function is either constant or continuous. However, in some situations, time-delay is a piecewise constant function. To the best of our knowledge, few research works have been dedicated to the development of computational algorithms for solving DDEs involving piecewise constant delay [4,19]. Recently, Marzban and Shahsiah proposed an efficient numerical scheme for solving DDEs containing piecewise constant delay. Their method is based on a hybrid of block-pulse functions and Chebyshev polynomials. It has been demonstrated that the method implemented in [19] is effective and produces very accurate results.
In what follows, we describe some similarities and differences between our method and the procedures developed in [15,16]. First, the foundations of the proposed framework and those used in the previous works are based on a hybrid of block-pulse functions and Legendre polynomials. Second, the current paper is an extension of our previous ones. More specifically, the time-delay systems considered in [15,16] involve constant delay, while here we investigate linear piecewise constant delay systems. Obviously, the latter systems are a general class of constant DDEs. Third, the approach employed here is based on the derivative matrix corresponding to the mentioned hybrid functions while the method implemented in our earlier works is based on the operational matrix of integration. Fourth, the operational matrix of delay associated with the piecewise constant delay systems is 2 Advances in Numerical Analysis constructed. To derive this matrix, we use the same approach as that of [19]. The purpose of this paper is to introduce an efficient numerical technique for solving DDEs with a piecewise constant delay. It should be emphasized that the analytical solutions of DDEs cannot be obtained solely either by blockpulse functions or by Legendre polynomials. Combining block-pulse functions and Legendre polynomials allows one to simultaneously make use of the best properties of the two mentioned bases. It is worth noting that the value of , the order of block-pulse functions, plays an essential role in modelling of the problem under consideration. Indeed, by selecting the suitable value of , we are able to correctly determine the exact locations of the switching points that occur in the solution associated with a piecewise constant delay system. The excellent properties of hybrid functions together with the operational matrices of delay, derivative, and product are then utilized to transform the delay differential equation under investigation into a system of algebraic equations whose solution is much easier than the original one.
The rest of the paper is organized as follows. In Section 2, the basic properties of hybrid of block-pulse functions and Legendre polynomials are presented. In Section 3, the operational matrices of derivative, product, and delay are presented. Section 4 is devoted to the problem statement and its approximation. In Section 5, three examples are tested to show the efficiency and accuracy of the proposed numerical scheme.

Hybrid Functions
where and are the orders of block-pulse functions and Legendre polynomials, respectively. Here, ( ) are the wellknown Legendre polynomials of order which are orthogonal on the interval [−1, 1] and satisfy the following recursive formula [20]: Since ( ) consists of block-pulse functions and Legendre polynomials, which are both complete and orthogonal, the set of the hybrid of block-pulse functions and Legendre polynomials is a complete orthogonal set in the Hilbert space L 2 [0, ).

Operational Matrices
In this section, we first present the operational matrix of derivative based on the weak representation of the derivative operator. We then state the operational matrices of delay and product corresponding to the proposed hybrid functions.

The Operational Matrix of Derivative.
We approximate the derivative of by the derivative of (3); that is, The right hand side of the preceding equation can be represented in terms of hybrid functions as where The relationship between the two vectors and is expressed by where ( +1)× ( +1) matrix D is the operational matrix of derivative. The structure of the operational matrix of derivative D is given by [21].
To obtain the operational matrix of delay corresponding to the proposed hybrid functions, we apply an approach analogous to the one devised in [19]. To do this, we first divide the time interval [0, ] into subintervals of the same length where the value of is obtained in the following manner.
Define as the smallest positive integer number in such a way that Suppose that | | = . Next, we choose as the greatest common divisor of the integers and , ∈ , = 0, 1, . . . , ; that is, where [( / ) ] denotes the greatest integer value less than or equal to ( / ) . It should be noted that is chosen in such a way that the number of subintervals is minimized. As a consequence, we obtain the following subintervals: where ℎ = / = / . Now, from (20) and (22) we can find the integer numbers and , such that = ℎ, ∈ , = ℎ, = 0, 1, . . . , . 4

Advances in Numerical Analysis
Note that = 0, for ∉ , 1 ≤ ≤ . Therefore, we can rewrite ( ) as follows: For convenience, we define ( ) as Therefore, the problem is reduced to find the operational matrix of delay for the following delay function: where = ℎ, = 1, 2, . . . , .
In order to find the matrix D, we first derive the matrix for = 1, 2, . . . , , so that the following relation holds: With the use of (1), it is worth noting that, in the case of −1 ≤ < , the only functions with nonzero values are ( − ( )) ( − ( )ℎ), for = 0, 1, . . . , . Since by expanding ( − ( )) ( − ( )ℎ) in terms of ( ), we get the ( + 1) × ( + 1) identity matrix as the corresponding coefficients. Therefore, we have where +1 is the ( + 1)-dimensional identity matrix and is an × matrix in which the only nonzero entry is equal to one and located at the ( − ( ))th row and th column.
As a consequence, if we expand Φ( − ( )) in terms of hybrid functions Φ( ), we deduce that To illustrate the derivation process of the operational matrix of delay associated with the mentioned hybrid functions, we present an example. Define It is easily verified that = 6. As a result, we get For = 1 and = 3, we conclude that − ( ) = 0. Therefore, 1 and 3 are zero matrices of order 6 × 6. If = 2, it follows that − ( ) = 1. Consequently, 2 is a 6 × 6 matrix in which the only nonzero entry is equal to one and located at the 1st row and 2nd column. If = 4, then we deduce that − ( ) = 1. Hence, 4 is a 6 × 6 matrix in which the only nonzero entry is equal to one and located at the 1st row and 4th column. A similar argument can be used for = 5 and = 6 to specify matrices 5 , and 6 . Using the above comments, we obtain ] .

Approximation Using Hybrid
where and are the × and × identity matrices, respectively, and ⊗ denotes the Kronecker product [22]. It should be noted that and are vectors of orders ( + 1) × 1 and ( + 1) × 1, respectively, given by Also where is a vector of order ( + 1) × 1 defined by Similarly ( ) = ( ⊗ Φ ( )) , where , , and are matrices of orders ( + 1) × , ( + 1) × , and ( + 1) × , respectively. The delay vector ( − ( )) can also be expanded in terms of hybrid functions as ( − ( )) = ( ⊗ Φ ( )) where D is the operational matrix of delay described by (18). Now, using (16), we have Advances in Numerical Analysis wherẽand̃can be calculated in a way similar to the construction method of matrix̃given by (17). Moreover, Consequently Now, using the operational matrix of derivative and substituting (44)-(50) in (37), we get Because the elements of Φ( ) are linearly independent functions over the interval [0, ], it follows that Therefore, the original problem is transformed into a system of linear algebraic equations. The aforementioned system with the associated initial condition can be easily solved by the well-known Tau method [20,23] for the unknown vector .

Illustrative Examples
In this section, three examples are investigated to evaluate the performance of the method.
The exact solution to this problem is given by [19] To solve this problem by the method developed in the present paper, we first determine the value of , the required number of subintervals, with the use of (22). For this problem, we choose = 10. Also, we select = 8. Let With the use of (58), we would obtain the same value as the exact value of ( ).
Example 2. Consider the following piecewise constant delay system [13]: The analytical solution to this problem is presented by [13] ( To solve this problem by the proposed approach, we first determine the value of . With the use of (22), we take = 20. In addition, let be an arbitrary positive integer number. Define where ( ) and ( ) denote the approximate solution obtained by the present method and the exact solution, respectively. In Table 1, the maximum absolute error of ( ) corresponding to = 20 and different values of are summarized. This table shows that there is an excellent agreement between the analytical and approximate solutions. After choosing , small values for are needed to achieve a satisfactory approximation.
The exact solution to this problem is described by Although the above problem is a nonlinear delay differential equation, the method developed in the current paper is applicable. To employ the procedure described in Section 4, we first select = 5. Moreover, we choose = 2. Let with the following initial condition: in which and D is the operational matrix of delay given by (18). Moreover,̃and̃can be determined in a way similar to the construction method of the matrix̃presented by (17). By solving (77) ] . (80) Using (73), we would obtain the same value as the exact value of ( ).

Conclusion
An efficient procedure has been successfully developed for solving delay differential equations with a piecewise constant delay function. The method is based upon a hybrid of blockpulse functions and Legendre polynomials. The nice properties of the hybrid functions together with the associated operational matrices were used to convert the original problem into a system of algebraic equations. The proposed framework allows one to simultaneously make use of the best advantages of the two mentioned bases. After determining the appropriate value of , small values for are required to obtain an admissible approximation. It is worth noting that the correct choice of has a fundamental effect on the solution accuracy. It should be emphasized that the analytical solutions of Examples 1 and 3 cannot be derived solely either by block-pulse functions or by Legendre polynomials. The simulation results demonstrate the reliability and effectiveness of the proposed approximation scheme.