# Numerical Solution of Piecewise Constant Delay Systems Based on a Hybrid Framework

Authors: H. R. Marzban; S. Hajiabdolrahmani
Journal: International Journal of Differential Equations
 (2016)

---

## 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.


---

## Body

## 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 Section2, 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ϕnm(t),  n=1,2,…,N,  m=0,1,…,M-1, are defined on the interval 0,tf as [14](1)
ϕ
n
m
t
=
T
m
N
t
-
n
-
1
t
f
,
t
∈
n
-
1
N
t
f
,
n
N
t
f
,
0
,
o
t
h
e
r
w
i
s
e
,where n is the order of block-pulse functions. Here Tm(t)=tm, m=0,1,…,M-1, are the well-known Taylor’s polynomials of degree m. 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 0,tf.


### 2.2. Function Approximation

A functionf(t) can be represented in terms of the hybrid functions as(2)ft=∑n=1∞∑m=0∞cnmϕnmt,where the coefficients cnm,  n=1,2,…,N,  m=0,1,…,M-1, are calculated by the following formula: (3)cnm=1Nmm!dmftdtmt=n-1/Ntf.If the series given by (2) is truncated, then it can be written in the form(4)ft≃∑n=1N∑m=0M-1cnmϕnmt=CTΦt,where(5)C=c10,…,c1,M-1,c20,…,c2,M-1,…,cN0,…,cN,M-1T,(6)Φt=ϕ10t,…,ϕ1,M-1t,ϕ20t,…,ϕ2,M-1t,…,ϕN0t,…,ϕN,M-1tT.In 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 that(7)ft=t,0≤t<1,1,1≤t<2,3-t,2≤t≤3.
Clearly, the functionf(t) has two corner points at the points t=1 and t=2. Accordingly, Taylor’s expansion of f(t) about each of the mentioned points does not exist. This function can be exactly approximated by the hybrid functions as follows: (8)ft=∑n=13∑m=01cnmϕnmt=CTΦt,in which(9)C=0,13,1,0,1,-13T,Φt=ϕ10t,ϕ11t,ϕ20t,ϕ21t,ϕ30t,ϕ31tT.The elements of the vector Φ(t) are defined by (10)ϕ10t=1,t∈0,1,0,otherwise,ϕ11t=3t,t∈0,1,0,otherwise,ϕ20t=1,t∈1,2,0,otherwise,ϕ21t=3t-3,t∈1,2,0,otherwise,ϕ30t=1,t∈2,3,0,otherwise,ϕ31t=3t-6,t∈2,3,0,otherwise.Case 2.
Suppose that(11)gt=expt,0≤t<1,t2+t,1≤t<2,t+4,2≤t≤3.Obviously, this function is discontinuous at the points t=1 and t=2. However, it can be efficiently approximated by the proposed hybrid functions. To express this function in terms of hybrid functions, we choose N=3. In addition, Assume that M is an arbitrary positive integer number greater than 2. With the use of (2), we have (12)gt=∑n=13∑m=0M-1cnmϕnmt,where the hybrid functions ϕnm(t), n=1,2,3, m=0,1,…,M-1, are defined by (13)ϕnmt=Tm3t-3n+3=3mt-n+1m.Also, the associated coefficients cnm are given by(14)c1m=13mm!,m=0,1,…,M-1,c20=2,c21=1,c22=19,c2m=0,m=3,4,…,M-1,c30=6,c31=13,c3m=0,m=2,3,…,M-1.
The hybrid expansion ofg(t) over the interval [1,2] is presented by (15)c20ϕ20t+c21ϕ21t+c22ϕ22t=2+3t-3+193t-32=t2+t.Similarly, the hybrid approximation of g(t) on the interval 2,3 is described by (16)c30ϕ30t+c31ϕ31t=6+133t-6=t+4.Moreover, the hybrid expansion of g(t) over the interval [0,1) is expressed by (17)gt=∑m=0M-1tmm!,which is precisely Taylor’s expansion of the function exp⁡t about t=0, up to degree M-1. Define (18)E=gt-CTΦt∞=max⁡gt-CTΦt:0≤t≤1.It is easily verified that (19)E=eM!.The obtained result indicates that the maximum error converges to zero as the value of M 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Φ(t) presented by (6) can also be expanded in terms of hybrid functions as(20)∫0tΦτdτ≃PΦt,in which P is the MN×MN operational matrix of integration corresponding to the hybrid functions. The structure of this matrix is given by [14](21)P=E∗HH⋯H0E∗H⋯H00E∗⋯H⋮⋮⋮⋱⋮000⋯E∗,where (22)H=1Ntf0⋯0tf220⋯0⋮⋮⋱⋮tfMM0⋯0,and E∗ is defined by (23)
E
∗
=
1
N
0
1
0
⋯
0
0
0
1
2
⋯
0
⋮
⋮
⋮
⋱
⋮
0
0
0
⋯
1
M
-
1
0
0
0
⋯
0
.It should be noted that E∗ is the operational matrix of integration associated with the Taylor’s polynomials on the interval [((n-1)/N)tf,(n/N)tf). Furthermore, if C is a vector of order MN×1, then the expression Φ(t)ΦT(t)C can be represented in terms of hybrid functions as(24)ΦtΦTtC≃C~Φt.The MN×MN matrix C~ is called the operational matrix of product [14].


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

Hybrid functionsϕnm(t),  n=1,2,…,N,  m=0,1,…,M-1, are defined on the interval 0,tf as [14](1)
ϕ
n
m
t
=
T
m
N
t
-
n
-
1
t
f
,
t
∈
n
-
1
N
t
f
,
n
N
t
f
,
0
,
o
t
h
e
r
w
i
s
e
,where n is the order of block-pulse functions. Here Tm(t)=tm, m=0,1,…,M-1, are the well-known Taylor’s polynomials of degree m. 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 0,tf.


## 2.2. Function Approximation

A functionf(t) can be represented in terms of the hybrid functions as(2)ft=∑n=1∞∑m=0∞cnmϕnmt,where the coefficients cnm,  n=1,2,…,N,  m=0,1,…,M-1, are calculated by the following formula: (3)cnm=1Nmm!dmftdtmt=n-1/Ntf.If the series given by (2) is truncated, then it can be written in the form(4)ft≃∑n=1N∑m=0M-1cnmϕnmt=CTΦt,where(5)C=c10,…,c1,M-1,c20,…,c2,M-1,…,cN0,…,cN,M-1T,(6)Φt=ϕ10t,…,ϕ1,M-1t,ϕ20t,…,ϕ2,M-1t,…,ϕN0t,…,ϕN,M-1tT.In 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 that(7)ft=t,0≤t<1,1,1≤t<2,3-t,2≤t≤3.
Clearly, the functionf(t) has two corner points at the points t=1 and t=2. Accordingly, Taylor’s expansion of f(t) about each of the mentioned points does not exist. This function can be exactly approximated by the hybrid functions as follows: (8)ft=∑n=13∑m=01cnmϕnmt=CTΦt,in which(9)C=0,13,1,0,1,-13T,Φt=ϕ10t,ϕ11t,ϕ20t,ϕ21t,ϕ30t,ϕ31tT.The elements of the vector Φ(t) are defined by (10)ϕ10t=1,t∈0,1,0,otherwise,ϕ11t=3t,t∈0,1,0,otherwise,ϕ20t=1,t∈1,2,0,otherwise,ϕ21t=3t-3,t∈1,2,0,otherwise,ϕ30t=1,t∈2,3,0,otherwise,ϕ31t=3t-6,t∈2,3,0,otherwise.Case 2.
Suppose that(11)gt=expt,0≤t<1,t2+t,1≤t<2,t+4,2≤t≤3.Obviously, this function is discontinuous at the points t=1 and t=2. However, it can be efficiently approximated by the proposed hybrid functions. To express this function in terms of hybrid functions, we choose N=3. In addition, Assume that M is an arbitrary positive integer number greater than 2. With the use of (2), we have (12)gt=∑n=13∑m=0M-1cnmϕnmt,where the hybrid functions ϕnm(t), n=1,2,3, m=0,1,…,M-1, are defined by (13)ϕnmt=Tm3t-3n+3=3mt-n+1m.Also, the associated coefficients cnm are given by(14)c1m=13mm!,m=0,1,…,M-1,c20=2,c21=1,c22=19,c2m=0,m=3,4,…,M-1,c30=6,c31=13,c3m=0,m=2,3,…,M-1.
The hybrid expansion ofg(t) over the interval [1,2] is presented by (15)c20ϕ20t+c21ϕ21t+c22ϕ22t=2+3t-3+193t-32=t2+t.Similarly, the hybrid approximation of g(t) on the interval 2,3 is described by (16)c30ϕ30t+c31ϕ31t=6+133t-6=t+4.Moreover, the hybrid expansion of g(t) over the interval [0,1) is expressed by (17)gt=∑m=0M-1tmm!,which is precisely Taylor’s expansion of the function exp⁡t about t=0, up to degree M-1. Define (18)E=gt-CTΦt∞=max⁡gt-CTΦt:0≤t≤1.It is easily verified that (19)E=eM!.The obtained result indicates that the maximum error converges to zero as the value of M 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Φ(t) presented by (6) can also be expanded in terms of hybrid functions as(20)∫0tΦτdτ≃PΦt,in which P is the MN×MN operational matrix of integration corresponding to the hybrid functions. The structure of this matrix is given by [14](21)P=E∗HH⋯H0E∗H⋯H00E∗⋯H⋮⋮⋮⋱⋮000⋯E∗,where (22)H=1Ntf0⋯0tf220⋯0⋮⋮⋱⋮tfMM0⋯0,and E∗ is defined by (23)
E
∗
=
1
N
0
1
0
⋯
0
0
0
1
2
⋯
0
⋮
⋮
⋮
⋱
⋮
0
0
0
⋯
1
M
-
1
0
0
0
⋯
0
.It should be noted that E∗ is the operational matrix of integration associated with the Taylor’s polynomials on the interval [((n-1)/N)tf,(n/N)tf). Furthermore, if C is a vector of order MN×1, then the expression Φ(t)ΦT(t)C can be represented in terms of hybrid functions as(24)ΦtΦTtC≃C~Φt.The MN×MN matrix C~ 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 Φ(t-a(t)) can be represented in terms of Φ(t) as follows:(25)Φt-at=DΦt,in which D is called the operational matrix of delay. In addition, a(t) is a piecewise constant function defined by (26)at=τ1,T0≤t<T1,τ2,T1≤t<T2,⋮⋮τr,Tr-1≤t<Tr,where 0=T0<T1<⋯<Tr-1<Tr=tf. It is further supposed that τi,  i=1,2,…,r, are arbitrary rational numbers in the interval [0,tf) and Ti∈Q,  i=0,1,…,r. To construct D, we first obtain the value of N. Let (27)A=i:τi≠0,and take γ as the smallest positive integer number such that the following conditions are satisfied: (28)γτi∈Z,γTj∈Z,i∈A,j=1,2,…,r.Assume that A=k. Next, we define λ as the greatest common divisor of the integers γτi and γTj, i∈A and j=1,2,…,r; that is,(29)λ=g.c.dγτ1,γτ2,…,γτk,γT1,γT2,…,γTr.Define(30)
N
=
γ
λ
t
f
,
if
γ
λ
t
f
∈
Z
,
γ
λ
t
f
+
1
,
otherwise
,where [(γ/λ)tf] denotes the greatest integer value less than or equal to (γ/λ)tf. It is worth noting that N is selected in such a way that the number of subintervals is minimized. As a consequence, we obtain the following subintervals: (31)0,λγ,λγ,2λγ,…,N-1λγ,Nλγ.For convenience, we use the notations defined by (32)h=λγ,Nj=γλTj-Tj-1,lj=γλτj,j=1,2,…,r.Clearly ∑j=1rNj=N. Consequently, we have (33)at=kih,i-1h≤t<ih,i=1,2,…,N,in which (34)ki=l1,1≤i≤N1,l2,N1+1≤i≤N1+N2,⋮⋮lj,∑k=1j-1Nk+1≤i≤∑k=1jNk,⋮⋮lr,∑k=1r-1Nk+1≤i≤∑k=1rNk.As a result, the problem is reduced to find the delay operational matrix for the delay function defined by (35)Φt-at=Φt-k1h,0≤t<t1,Φt-k2h,t1≤t<t2,⋮⋮Φt-kNh,tn-1≤t<tN,where (36)ti=ih,i=1,2,…,N.To construct the matrix D, we first derive the matrix Di for i=1,2,…,N, in such a way that the following relation is satisfied: (37)Φt-kih=DiΦt,ti-1≤t<ti.By the properties of the hybrid functions, it is easy to verify that, for the case ti-1≤t<ti, the only functions with nonzero values are ϕ(t-ki)m(t-kih), for m=0,1,…,M-1. Notice that (38)ϕt-kimt-kih=ϕimt,m=0,1,…,M-1.Accordingly, ϕ(t-ki)m(t-kih) can be written in terms of ϕim(t) as (39)Di=Si⊗IM,i=1,2,…,N,where IM is the M-dimensional identity matrix and ⊗ denotes the Kronecker product [24]. In addition, Si is an N×N matrix in which the only nonzero entry is equal to one and located at the (i-ki)th row and ith column. Finally, if we express Φ(t-a(t)) in terms of Φ(t), we find(40)D=D1+D2+⋯+DN.Remark 1.
Ifi-ki≤0, then Si is a zero matrix of order N×N. 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 by(41)X˙t=EtXt+FtXt-at+GtUt,0≤t≤tf,(42)X0=X0,(43)Xt=0,t<0,(44)at=τ1,T0≤t<T1,τ2,T1≤t<T2,⋮⋮τr,Tr-1≤t<Tr,where X(t)∈Rl,  U(t)∈Rq,  E(t),  F(t), and G(t) are matrices of appropriate dimensions, X0 is a constant specified vector, and τi,  i=1,2,…,r, are arbitrary rational numbers in [0,tf). The objective is to find X(t),  0≤t≤tf, satisfying (41)–(43).


### 4.1. Description of the Method

Let(45)Xt=x1t,x2t,…,xltT,Ut=u1t,u2t,…,uqtT.Assume that each xi(t) and each of uj(t),  i=1,2,…,q, can be written in terms of hybrid functions as (46)xit=ΦTtXi,ujt=ΦTtUj.Consequently(47)Xt=Il⊗ΦTtX,Ut=Iq⊗ΦTtU,where Il and Iq are identity matrices of order l×l and q×q, respectively. Furthermore, X and U are vectors of order lNM×1 and qNM×1, respectively, given by (48)X=X1T,X2T,…,XlTT,U=U1T,U2T,…,UqTT.Also(49)X0=Il⊗ΦTtd,where d is a vector of order lNM×1, defined by (50)d=d1T,d2T,…,dlTT.Similarly, E(t), F(t), and G(t) can be represented in terms of hybrid functions as(51)Et=ETIl⊗Φt,Ft=FTIl⊗Φt,Gt=GTIq⊗Φt,where ET, FT, and GT are of dimensions l×lNM, l×lNM, and l×qNM, respectively. The delayed vector X(t-a(t)) can also be expressed in terms of hybrid functions as(52)Xt-at=Il⊗ΦTtIl⊗DTX,in which D is the delay operational matrix presented by (25). As a result(53)EtXt=ETIl⊗ΦtIl⊗ΦTtX=Il⊗ΦTtE~TX,FtXt-at=FTIl⊗ΦtIl⊗ΦTtIl⊗DTX=Il⊗ΦTtF~TIl⊗DTX,GtUt=GTIq⊗ΦtIq⊗ΦTtU=Il⊗ΦTtG~TU,where E~, F~, and G~ can be determined in a manner analogous to the derivation process of the matrix C~ described by (24). Furthermore(54)∫0tIl⊗ΦTsds=Il⊗ΦTtIl⊗PT,∫0tFsXs-asds=Il⊗ΦTtIl⊗PTF~TIl⊗DTX,in which P is the operational matrix of integration described by (20). By integrating both sides of (41) with respect to t and using (47), (49), (51), (52), and (54), we obtain(55)Il⊗ΦTtX-Il⊗ΦTtd=Il⊗ΦTtIl⊗PTE~TX+Il⊗ΦTtIl⊗PTF~TIl⊗DTX+Il⊗ΦTtIl⊗PTG~TU.From the preceding equation, we deduce that(56)X=IMN-Il⊗PTE~T-Il⊗PTF~TIl⊗DT-1d+Il⊗PTG~TU.


## 4.1. Description of the Method

Let(45)Xt=x1t,x2t,…,xltT,Ut=u1t,u2t,…,uqtT.Assume that each xi(t) and each of uj(t),  i=1,2,…,q, can be written in terms of hybrid functions as (46)xit=ΦTtXi,ujt=ΦTtUj.Consequently(47)Xt=Il⊗ΦTtX,Ut=Iq⊗ΦTtU,where Il and Iq are identity matrices of order l×l and q×q, respectively. Furthermore, X and U are vectors of order lNM×1 and qNM×1, respectively, given by (48)X=X1T,X2T,…,XlTT,U=U1T,U2T,…,UqTT.Also(49)X0=Il⊗ΦTtd,where d is a vector of order lNM×1, defined by (50)d=d1T,d2T,…,dlTT.Similarly, E(t), F(t), and G(t) can be represented in terms of hybrid functions as(51)Et=ETIl⊗Φt,Ft=FTIl⊗Φt,Gt=GTIq⊗Φt,where ET, FT, and GT are of dimensions l×lNM, l×lNM, and l×qNM, respectively. The delayed vector X(t-a(t)) can also be expressed in terms of hybrid functions as(52)Xt-at=Il⊗ΦTtIl⊗DTX,in which D is the delay operational matrix presented by (25). As a result(53)EtXt=ETIl⊗ΦtIl⊗ΦTtX=Il⊗ΦTtE~TX,FtXt-at=FTIl⊗ΦtIl⊗ΦTtIl⊗DTX=Il⊗ΦTtF~TIl⊗DTX,GtUt=GTIq⊗ΦtIq⊗ΦTtU=Il⊗ΦTtG~TU,where E~, F~, and G~ can be determined in a manner analogous to the derivation process of the matrix C~ described by (24). Furthermore(54)∫0tIl⊗ΦTsds=Il⊗ΦTtIl⊗PT,∫0tFsXs-asds=Il⊗ΦTtIl⊗PTF~TIl⊗DTX,in which P is the operational matrix of integration described by (20). By integrating both sides of (41) with respect to t and using (47), (49), (51), (52), and (54), we obtain(55)Il⊗ΦTtX-Il⊗ΦTtd=Il⊗ΦTtIl⊗PTE~TX+Il⊗ΦTtIl⊗PTF~TIl⊗DTX+Il⊗ΦTtIl⊗PTG~TU.From the preceding equation, we deduce that(56)X=IMN-Il⊗PTE~T-Il⊗PTF~TIl⊗DT-1d+Il⊗PTG~TU.


## 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 ofN would greatly improve the efficiency and accuracy of the presented procedure.


### 5.1. Example  1

Consider the following piecewise constant delay system [17]:(57)x˙t=xt-at+ut,0≤t≤1,(58)xt=0,t≤0,(59)ut=1,t≥0,(60)at=0.1,0≤t<0.35,0.3,0.35≤t<0.7,0.5,0.7≤t≤1.The exact solution to this problem is given by(61)
x
t
=
t
,
0
≤
t
<
0.1
,
1
200
+
9
10
t
+
1
2
t
2
,
0.1
≤
t
<
0.2
,
11
3000
+
23
25
t
+
2
5
t
2
+
1
6
t
3
,
0.2
≤
t
<
0.3
,
961
240000
+
1831
2000
t
+
169
400
t
2
+
7
60
t
3
+
1
24
t
4
,
0.3
≤
t
<
0.35
,
290161
3840000
+
7
10
t
+
1
2
t
2
,
0.35
≤
t
<
0.4
,
83067
1280000
+
39
50
t
+
3
10
t
2
+
1
6
t
3
,
0.4
≤
t
<
0.5
,
259201
3840000
+
911
1200
t
+
29
80
t
2
+
1
12
t
3
+
1
24
t
4
,
0.5
≤
t
<
0.6
,
6417817
96000000
+
22937
30000
t
+
689
2000
t
2
+
17
150
t
3
+
1
60
t
4
+
1
120
t
5
,
0.6
≤
t
<
0.65
,
1070467
48000000
+
3496561
3840000
t
+
1
5
t
2
+
1
6
t
3
,
0.65
≤
t
<
0.7
,
37190303
192000000
+
3737
6000
t
+
129
400
t
2
+
1
20
t
3
+
1
24
t
4
,
0.7
≤
t
<
0.8
,
2444401
12800000
+
6399
10000
t
+
1679
6000
t
2
+
31
300
t
3
+
1
120
t
4
+
1
120
t
5
,
0.8
≤
t
<
0.85
,
7103237
64000000
+
3266161
3840000
t
+
1
10
t
2
+
1
6
t
3
,
0.85
≤
t
<
0.9
,
8852837
64000000
+
2799601
3840000
t
+
121
400
t
2
+
1
60
t
3
+
1
24
t
4
,
0.9
≤
t
≤
1
.To apply the method developed in the current paper, we first choose the value ofN, with the use of (30). Because γ=20 and λ=1, hence, we select N=20. We also take M=6. Let(62)xt=CTΦt,where(63)C=c10,…,c1,M-1,…,cN0,…,cN,M-1T,Φt=ϕ10t,…,ϕ1,M-1t,…,ϕN0t,…,ϕN,M-1tT.By expanding u(t) in terms of hybrid functions, we get(64)ut=d1TΦt.We also have(65)xt-at=CTDΦt,where D is the delay operational matrix given by (25). Integrating both sides of (57) with respect to t and using (62)–(65) yield(66)CT=d1TPIMN-DP-1,in which IMN is the MN-dimensional identity matrix and P 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 x(t).


### 5.2. Example  2

Consider the following two-dimensional time-varying piecewise constant delay system:(67)x˙1tx˙2t=11+t2tt2x1t-atx2t-at+12ut,0≤t≤1,(68)x1tx2t=00,t<0,(69)x10x20=11,(70)ut=t,0≤t<0.5,1-t,0.5≤t<1,(71)at=0.2,0≤t<0.3,0.7,0.3≤t≤1.The exact solutions to this problem are given as follows: (72)
x
1
t
=
1
+
1
2
t
2
,
0
≤
t
<
0.2
,
4319
7500
+
103
50
t
+
18
25
t
2
+
11
30
t
3
+
1
4
t
4
,
0.2
≤
t
<
0.3
,
147071
120000
+
1
2
t
2
,
0.3
≤
t
<
0.5
,
117071
120000
+
t
-
1
2
t
2
,
0.5
≤
t
<
0.7
,
-
10331
12000
+
747
200
t
-
161
200
t
2
+
1
30
t
3
+
1
4
t
4
,
0.7
≤
t
<
0.9
,
-
424558459
4200000000
+
1148571
500000
t
-
355883
2000000
t
2
-
4681
30000
t
3
+
1767
4000
t
4
+
71
750
t
5
-
7
120
t
6
+
1
35
t
7
0.9
≤
t
≤
1
,
x
2
t
=
1
+
t
2
,
0
≤
t
<
0.2
,
179473
187500
+
101
50
t
2
+
16
75
t
3
+
3
20
t
4
+
1
5
t
5
,
0.2
≤
t
<
0.3
,
3169351
3000000
+
t
2
,
0.3
≤
t
<
0.5
,
1669351
3000000
+
2
t
-
t
2
,
0.5
≤
t
<
0.7
,
-
220481
3000000
+
2
t
+
49
200
t
2
+
3
100
t
3
-
1
10
t
4
+
1
5
t
5
,
0.7
≤
t
<
0.9
,
36460339381
84000000000
+
2
t
-
63163
40000
t
2
+
4372217
3000000
t
3
-
11101
40000
t
4
+
1981
15000
t
5
+
191
900
t
6
-
11
140
t
7
+
1
40
t
8
,
0.9
≤
t
≤
1
.To solve this problem by the proposed approach, we first determine the value of N with the use of (30). Because γ=10 and λ=1, consequently, we select N=10. Also, we choose M=9. Let(73)x1t=C1TΦt,x2t=C2TΦt.Expanding 1, 1+t, 2t, and t2 in terms of hybrid functions implies(74)1=f1TΦt,1+t=f2TΦt,2t=f3TΦt,t2=f4TΦt.Also, using (70), we obtain(75)
∫
0
t
u
τ
d
τ
=
1
2
t
2
=
f
5
T
Φ
t
,
0
≤
t
<
0.5
,
-
1
4
+
t
-
1
2
t
2
=
f
6
T
Φ
t
,
0.5
≤
t
<
1
.Moreover(76)x1t-at=C1TDΦt,x2t-at=C2TDΦt.Therefore(77)1+tx2t-at=C2TDΦtΦTtf2=C2TDf~2Φt,2tx1t-at=C1TDΦtΦTtf3=C1TDf~3Φt,t2x2t-at=C2TDΦtΦTtf4=C2TDf~4Φt,in which f~2, f~3, and f~4 can be obtained in a way similar to the construction method of the matrix C~ defined by (24). By integrating both sides of (67) with respect to t and using (68), (70), (71), (73), (74), (75), (76), and (77), we get(78)C1TIMN-DP-C2TDf~2P=f1T+f5T+f6T,C2TIMN-Df~4P-C1TDf~3P=f1T+2f5T+2f6T.By solving the resulting system described by (78) and then using (73), we would obtain the same values as the exact values of x1(t) and x2(t).


### 5.3. Example  3

Consider the following time-delay system [3]:(79)x˙t=-5xt-5xt-at+2ut,0≤t≤2,(80)x0=1,(81)ut=1,t≥0,(82)at=0,0≤t<0.8,0.3,0.8≤t<1.4,0.6,1.4≤t≤1.7,0.9,1.7≤t≤2.The exact solution to this problem is given by [3](83)xt=0.2+0.8e-10t,0≤t<0.8,0.2+0.8e3-10t+0.8e-41-e3e-5t,0.8≤t<1.1,0.2+0.8e5-10t-4e-2.51-e3te-5t+0.81-e36.5e-2.5+e-4e-5t,1.1≤t<1.4,0.2+0.8e9-10t-4e-11-e3te-5t+0.81-e38e-1-0.5e-2.5+e-4e-5t,1.4≤t<1.7,0.2+0.8e12-10t-4e0.51-e3te-5t+0.81-e39.5e0.5-0.5e-1-0.5e-2.5+e-4e-5t,1.7≤t≤2.To employ the proposed method, we first determine the suitable value of N. Using (30), it follows that (84)N=γλtf=20.Assume that M is an arbitrary positive integer number. Let(85)xt=CTΦt.Similarly(86)x0=ut=dTΦt,xt-at=CTDΦt.By integrating both sides of (79) with respect to t and using (85)-(86), we conclude that(87)CT=dT+2dTPIMN+5P+5DP-1,where IMN is the MN-dimensional identity matrix. Define (88)E=max⁡xet-xt:0≤t≤1,in which x(t) and xe(t) 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 E, for N=20, and various values of M 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].Table 1
Computational results of maximum error for Example3.
Present method CPU time (s) E N = 20, M=1 0.112 1.165 e - 02 N = 20, M=2 0.202 3.231 e - 03 N = 20, M=4 0.398 2.573 e - 06 N = 20, M=6 0.769 4.836 e - 08 N = 20, M=8 0.853 3.942 e - 11 N = 20, M=10 1.282 2.376 e - 13 N = 20, M=12 1.841 5.875 e - 16 N = 20, M=14 2.167 3.624 e - 17 N = 20, M=16 2.536 5.117 e - 20


### 5.4. Example  4

As a more complicated problem, consider the following problem:(89)x˙t=tx2t-at+x3t-at+ut,0≤t≤1,(90)xt=0,t<0,(91)x0=1,(92)ut=t,0≤t<0.4,0.4,0.4≤t≤1,(93)at=0.4,0≤t<0.8,0.8,0.8≤t≤1.The analytical solution to this problem is described by(94)
x
t
=
1
+
1
2
t
2
,
0
≤
t
<
0.4
,
663251
1640625
+
25933
15625
t
-
729
6250
t
2
+
117
250
t
3
-
3
100
t
4
+
13
100
t
5
-
1
120
t
6
+
1
56
t
7
,
0.4
≤
t
<
0.8
,
235691
1640625
+
42187
15625
t
-
7623
6250
t
2
+
253
250
t
3
-
43
100
t
4
+
23
100
t
5
-
7
120
t
6
+
1
56
t
7
,
0.8
≤
t
≤
1
.Although the above-mentioned problem is a nonlinear delay differential equation, the method developed in the present paper can easily be implemented. Because γ=5 and λ=2, we take N=3. Furthermore, we take M=8. Let(95)xt=CTΦt.Similarly (96)x0=v1TΦt,t=v2TΦt,(97)∫0tuτdτ=12t2=v3TΦt,0≤t<0.4,-225+25t=v4TΦt,0.4≤t<1,where v1, v2, v3, and v4 are vectors of order MN×1. Also, we have(98)xt-at=CTDΦt,(99)x2t-at=CTDΦtΦTtDTC=CTDW~Φt,(100)x3t-at=CTDW~ΦtΦTtW=CTDW~2Φt, in which (101)W=DTC.Using (99) gives(102)tx2t-at=CTDW~ΦtΦTtv2=CTDW~v~2Φt.Integrating both sides of (89) with respect to t and using (95), (97), (98), (99), (100), and (102) imply that(103)CTIMN-DW~v~2P-DW~2P-v1T+v3T+v4T=0.By solving (103) and using (95), the exact value of xt will be obtained.


### 5.5. Example  5

As the final example, we consider the time-varying multidelay system defined by(104)x˙t=txt-a1t+t2x2t-a2t+ut,0≤t≤1,x0=1,ut=t,0≤t<12,12,12≤t≤1,where(105)
a
1
t
=
1
2
,
0
≤
t
<
3
4
,
1
4
,
3
4
≤
t
≤
1
,
a
2
t
=
1
4
,
0
≤
t
<
1
2
,
3
4
,
1
2
≤
t
≤
1
.The analytical solution to this problem is described by(106)
x
t
=
1
+
1
2
t
2
,
0
≤
t
<
1
4
,
20373
20480
+
1
2
t
2
+
11
32
t
3
-
1
16
t
4
+
1
10
t
5
,
1
4
≤
t
<
1
2
,
48191
61440
+
1
2
t
+
9
16
t
2
-
1
6
t
3
+
1
8
t
4
,
1
2
≤
t
<
3
4
,
1289743
1966080
+
1
2
t
+
14287
40960
t
2
+
187
384
t
3
-
1
256
t
4
+
1
24
t
5
+
1
48
t
6
,
3
4
≤
t
≤
1
.Although the above-mentioned system involves multiple piecewise constant delays, the proposed procedure can be applied. To employ our approach, we selectN=4. By choosing M=7, the exact solution of x(t) will be obtained.


## 5.1. Example  1

Consider the following piecewise constant delay system [17]:(57)x˙t=xt-at+ut,0≤t≤1,(58)xt=0,t≤0,(59)ut=1,t≥0,(60)at=0.1,0≤t<0.35,0.3,0.35≤t<0.7,0.5,0.7≤t≤1.The exact solution to this problem is given by(61)
x
t
=
t
,
0
≤
t
<
0.1
,
1
200
+
9
10
t
+
1
2
t
2
,
0.1
≤
t
<
0.2
,
11
3000
+
23
25
t
+
2
5
t
2
+
1
6
t
3
,
0.2
≤
t
<
0.3
,
961
240000
+
1831
2000
t
+
169
400
t
2
+
7
60
t
3
+
1
24
t
4
,
0.3
≤
t
<
0.35
,
290161
3840000
+
7
10
t
+
1
2
t
2
,
0.35
≤
t
<
0.4
,
83067
1280000
+
39
50
t
+
3
10
t
2
+
1
6
t
3
,
0.4
≤
t
<
0.5
,
259201
3840000
+
911
1200
t
+
29
80
t
2
+
1
12
t
3
+
1
24
t
4
,
0.5
≤
t
<
0.6
,
6417817
96000000
+
22937
30000
t
+
689
2000
t
2
+
17
150
t
3
+
1
60
t
4
+
1
120
t
5
,
0.6
≤
t
<
0.65
,
1070467
48000000
+
3496561
3840000
t
+
1
5
t
2
+
1
6
t
3
,
0.65
≤
t
<
0.7
,
37190303
192000000
+
3737
6000
t
+
129
400
t
2
+
1
20
t
3
+
1
24
t
4
,
0.7
≤
t
<
0.8
,
2444401
12800000
+
6399
10000
t
+
1679
6000
t
2
+
31
300
t
3
+
1
120
t
4
+
1
120
t
5
,
0.8
≤
t
<
0.85
,
7103237
64000000
+
3266161
3840000
t
+
1
10
t
2
+
1
6
t
3
,
0.85
≤
t
<
0.9
,
8852837
64000000
+
2799601
3840000
t
+
121
400
t
2
+
1
60
t
3
+
1
24
t
4
,
0.9
≤
t
≤
1
.To apply the method developed in the current paper, we first choose the value ofN, with the use of (30). Because γ=20 and λ=1, hence, we select N=20. We also take M=6. Let(62)xt=CTΦt,where(63)C=c10,…,c1,M-1,…,cN0,…,cN,M-1T,Φt=ϕ10t,…,ϕ1,M-1t,…,ϕN0t,…,ϕN,M-1tT.By expanding u(t) in terms of hybrid functions, we get(64)ut=d1TΦt.We also have(65)xt-at=CTDΦt,where D is the delay operational matrix given by (25). Integrating both sides of (57) with respect to t and using (62)–(65) yield(66)CT=d1TPIMN-DP-1,in which IMN is the MN-dimensional identity matrix and P 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 x(t).


## 5.2. Example  2

Consider the following two-dimensional time-varying piecewise constant delay system:(67)x˙1tx˙2t=11+t2tt2x1t-atx2t-at+12ut,0≤t≤1,(68)x1tx2t=00,t<0,(69)x10x20=11,(70)ut=t,0≤t<0.5,1-t,0.5≤t<1,(71)at=0.2,0≤t<0.3,0.7,0.3≤t≤1.The exact solutions to this problem are given as follows: (72)
x
1
t
=
1
+
1
2
t
2
,
0
≤
t
<
0.2
,
4319
7500
+
103
50
t
+
18
25
t
2
+
11
30
t
3
+
1
4
t
4
,
0.2
≤
t
<
0.3
,
147071
120000
+
1
2
t
2
,
0.3
≤
t
<
0.5
,
117071
120000
+
t
-
1
2
t
2
,
0.5
≤
t
<
0.7
,
-
10331
12000
+
747
200
t
-
161
200
t
2
+
1
30
t
3
+
1
4
t
4
,
0.7
≤
t
<
0.9
,
-
424558459
4200000000
+
1148571
500000
t
-
355883
2000000
t
2
-
4681
30000
t
3
+
1767
4000
t
4
+
71
750
t
5
-
7
120
t
6
+
1
35
t
7
0.9
≤
t
≤
1
,
x
2
t
=
1
+
t
2
,
0
≤
t
<
0.2
,
179473
187500
+
101
50
t
2
+
16
75
t
3
+
3
20
t
4
+
1
5
t
5
,
0.2
≤
t
<
0.3
,
3169351
3000000
+
t
2
,
0.3
≤
t
<
0.5
,
1669351
3000000
+
2
t
-
t
2
,
0.5
≤
t
<
0.7
,
-
220481
3000000
+
2
t
+
49
200
t
2
+
3
100
t
3
-
1
10
t
4
+
1
5
t
5
,
0.7
≤
t
<
0.9
,
36460339381
84000000000
+
2
t
-
63163
40000
t
2
+
4372217
3000000
t
3
-
11101
40000
t
4
+
1981
15000
t
5
+
191
900
t
6
-
11
140
t
7
+
1
40
t
8
,
0.9
≤
t
≤
1
.To solve this problem by the proposed approach, we first determine the value of N with the use of (30). Because γ=10 and λ=1, consequently, we select N=10. Also, we choose M=9. Let(73)x1t=C1TΦt,x2t=C2TΦt.Expanding 1, 1+t, 2t, and t2 in terms of hybrid functions implies(74)1=f1TΦt,1+t=f2TΦt,2t=f3TΦt,t2=f4TΦt.Also, using (70), we obtain(75)
∫
0
t
u
τ
d
τ
=
1
2
t
2
=
f
5
T
Φ
t
,
0
≤
t
<
0.5
,
-
1
4
+
t
-
1
2
t
2
=
f
6
T
Φ
t
,
0.5
≤
t
<
1
.Moreover(76)x1t-at=C1TDΦt,x2t-at=C2TDΦt.Therefore(77)1+tx2t-at=C2TDΦtΦTtf2=C2TDf~2Φt,2tx1t-at=C1TDΦtΦTtf3=C1TDf~3Φt,t2x2t-at=C2TDΦtΦTtf4=C2TDf~4Φt,in which f~2, f~3, and f~4 can be obtained in a way similar to the construction method of the matrix C~ defined by (24). By integrating both sides of (67) with respect to t and using (68), (70), (71), (73), (74), (75), (76), and (77), we get(78)C1TIMN-DP-C2TDf~2P=f1T+f5T+f6T,C2TIMN-Df~4P-C1TDf~3P=f1T+2f5T+2f6T.By solving the resulting system described by (78) and then using (73), we would obtain the same values as the exact values of x1(t) and x2(t).


## 5.3. Example  3

Consider the following time-delay system [3]:(79)x˙t=-5xt-5xt-at+2ut,0≤t≤2,(80)x0=1,(81)ut=1,t≥0,(82)at=0,0≤t<0.8,0.3,0.8≤t<1.4,0.6,1.4≤t≤1.7,0.9,1.7≤t≤2.The exact solution to this problem is given by [3](83)xt=0.2+0.8e-10t,0≤t<0.8,0.2+0.8e3-10t+0.8e-41-e3e-5t,0.8≤t<1.1,0.2+0.8e5-10t-4e-2.51-e3te-5t+0.81-e36.5e-2.5+e-4e-5t,1.1≤t<1.4,0.2+0.8e9-10t-4e-11-e3te-5t+0.81-e38e-1-0.5e-2.5+e-4e-5t,1.4≤t<1.7,0.2+0.8e12-10t-4e0.51-e3te-5t+0.81-e39.5e0.5-0.5e-1-0.5e-2.5+e-4e-5t,1.7≤t≤2.To employ the proposed method, we first determine the suitable value of N. Using (30), it follows that (84)N=γλtf=20.Assume that M is an arbitrary positive integer number. Let(85)xt=CTΦt.Similarly(86)x0=ut=dTΦt,xt-at=CTDΦt.By integrating both sides of (79) with respect to t and using (85)-(86), we conclude that(87)CT=dT+2dTPIMN+5P+5DP-1,where IMN is the MN-dimensional identity matrix. Define (88)E=max⁡xet-xt:0≤t≤1,in which x(t) and xe(t) 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 E, for N=20, and various values of M 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].Table 1
Computational results of maximum error for Example3.
Present method CPU time (s) E N = 20, M=1 0.112 1.165 e - 02 N = 20, M=2 0.202 3.231 e - 03 N = 20, M=4 0.398 2.573 e - 06 N = 20, M=6 0.769 4.836 e - 08 N = 20, M=8 0.853 3.942 e - 11 N = 20, M=10 1.282 2.376 e - 13 N = 20, M=12 1.841 5.875 e - 16 N = 20, M=14 2.167 3.624 e - 17 N = 20, M=16 2.536 5.117 e - 20


## 5.4. Example  4

As a more complicated problem, consider the following problem:(89)x˙t=tx2t-at+x3t-at+ut,0≤t≤1,(90)xt=0,t<0,(91)x0=1,(92)ut=t,0≤t<0.4,0.4,0.4≤t≤1,(93)at=0.4,0≤t<0.8,0.8,0.8≤t≤1.The analytical solution to this problem is described by(94)
x
t
=
1
+
1
2
t
2
,
0
≤
t
<
0.4
,
663251
1640625
+
25933
15625
t
-
729
6250
t
2
+
117
250
t
3
-
3
100
t
4
+
13
100
t
5
-
1
120
t
6
+
1
56
t
7
,
0.4
≤
t
<
0.8
,
235691
1640625
+
42187
15625
t
-
7623
6250
t
2
+
253
250
t
3
-
43
100
t
4
+
23
100
t
5
-
7
120
t
6
+
1
56
t
7
,
0.8
≤
t
≤
1
.Although the above-mentioned problem is a nonlinear delay differential equation, the method developed in the present paper can easily be implemented. Because γ=5 and λ=2, we take N=3. Furthermore, we take M=8. Let(95)xt=CTΦt.Similarly (96)x0=v1TΦt,t=v2TΦt,(97)∫0tuτdτ=12t2=v3TΦt,0≤t<0.4,-225+25t=v4TΦt,0.4≤t<1,where v1, v2, v3, and v4 are vectors of order MN×1. Also, we have(98)xt-at=CTDΦt,(99)x2t-at=CTDΦtΦTtDTC=CTDW~Φt,(100)x3t-at=CTDW~ΦtΦTtW=CTDW~2Φt, in which (101)W=DTC.Using (99) gives(102)tx2t-at=CTDW~ΦtΦTtv2=CTDW~v~2Φt.Integrating both sides of (89) with respect to t and using (95), (97), (98), (99), (100), and (102) imply that(103)CTIMN-DW~v~2P-DW~2P-v1T+v3T+v4T=0.By solving (103) and using (95), the exact value of xt will be obtained.


## 5.5. Example  5

As the final example, we consider the time-varying multidelay system defined by(104)x˙t=txt-a1t+t2x2t-a2t+ut,0≤t≤1,x0=1,ut=t,0≤t<12,12,12≤t≤1,where(105)
a
1
t
=
1
2
,
0
≤
t
<
3
4
,
1
4
,
3
4
≤
t
≤
1
,
a
2
t
=
1
4
,
0
≤
t
<
1
2
,
3
4
,
1
2
≤
t
≤
1
.The analytical solution to this problem is described by(106)
x
t
=
1
+
1
2
t
2
,
0
≤
t
<
1
4
,
20373
20480
+
1
2
t
2
+
11
32
t
3
-
1
16
t
4
+
1
10
t
5
,
1
4
≤
t
<
1
2
,
48191
61440
+
1
2
t
+
9
16
t
2
-
1
6
t
3
+
1
8
t
4
,
1
2
≤
t
<
3
4
,
1289743
1966080
+
1
2
t
+
14287
40960
t
2
+
187
384
t
3
-
1
256
t
4
+
1
24
t
5
+
1
48
t
6
,
3
4
≤
t
≤
1
.Although the above-mentioned system involves multiple piecewise constant delays, the proposed procedure can be applied. To employ our approach, we selectN=4. By choosing M=7, the exact solution of x(t) 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 Examples1, 2, 4, and 5 cannot be produced solely either by block-pulse functions or by Taylor’s polynomials. After determining the suitable value of N, the number of subintervals, small values for M, 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.



---

*Source: 9754906-2016-12-29.xml*