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

Authors: H. R. Marzban; S. M. Hoseini
Journal: Advances in Numerical Analysis
 (2016)

---

## Abstract

An efficient computational technique for solving linear delay differential equations with a piecewise constant delay function is presented. The new 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.


---

## Body

## 1. 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–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 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 block-pulse 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 N, the order of block-pulse functions, plays an essential role in modelling of the problem under consideration. Indeed, by selecting the suitable value of N, 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 Section2, 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.


## 2. Hybrid Functions

Hybrid functionsϕ
n
m
t, n
=
1,2
,
…
,
N, m
=
0,1
,
…
,
M, are defined on the interval 0
,
t
f as [15](1)
ϕ
n
m
t
=
P
m
2
N
t
f
t
-
2
n
+
1
,
t
∈
n
-
1
N
t
f
,
n
N
t
f
,
0
,
otherwise,where n and m are the orders of block-pulse functions and Legendre polynomials, respectively. Here, P
m
(
t
) are the well-known Legendre polynomials of order m which are orthogonal on the interval -
1,1 and satisfy the following recursive formula [20]:(2)
P
0
t
=
1
,
P
1
t
=
t
,
P
m
+
1
t
=
2
m
+
1
m
+
1
t
P
m
t
-
m
m
+
1
P
m
-
1
t
,
m
=
1,2
,
3
,
…
.Since ϕ
n
m
t 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
,
t
f.


## 3. 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.A functionu
∈
L
2
0
,
t
f can be approximated by the hybrid functions as follows:(3)
u
t
≃
∑
n
=
1
N
∑
m
=
0
M
u
^
n
m
ϕ
n
m
t
=
Φ
T
t
U
,where(4)
U
=
u
^
10
,
…
,
u
^
1
M
,
u
^
20
,
…
,
u
^
2
M
,
…
,
u
^
N
0
,
…
,
u
^
N
M
Tis the vector of coefficients and(5)
Φ
t
=
ϕ
10
t
,
…
,
ϕ
1
M
t
,
ϕ
20
t
,
…
,
ϕ
2
M
t
,
…
,
ϕ
N
0
t
,
…
,
ϕ
N
M
t
T
.The coefficients u
^
n
m, n
=
1,2
,
…
,
N, m
=
0,1
,
…
,
M, are obtained by the following formula:(6)
u
^
n
m
=
N
2
m
+
1
t
f
∫
n
-
1
/
N
t
f
n
/
N
t
f
u
t
ϕ
n
m
t
d
t
.


### 3.1. The Operational Matrix of Derivative

We approximate the derivative ofu by the derivative of (3); that is,(7)
d
d
t
u
t
≃
d
d
t
∑
n
=
1
N
∑
m
=
0
M
u
^
n
m
ϕ
n
m
t
.The right hand side of the preceding equation can be represented in terms of hybrid functions as(8)
d
d
t
∑
n
=
1
N
∑
m
=
0
M
u
^
n
m
ϕ
n
m
t
=
∑
n
=
1
N
∑
m
=
0
M
s
^
n
m
ϕ
n
m
t
=
Φ
T
t
S
,where(9)
S
=
s
^
10
,
…
,
s
^
1
M
,
s
^
20
,
…
,
s
^
2
M
,
…
,
s
^
N
0
,
…
,
s
^
N
M
T
.The relationship between the two vectorsU and S is expressed by(10)
S
=
D
U
,where N
M
+
1
×
N
M
+
1 matrix D is the operational matrix of derivative.The structure of the operational matrix of derivativeD is given by [21].In the caseN
=
1,(11)
D
=
A
,in the case N
=
2,(12)
D
=
A
1
B
1
C
N
A
N
,and, in the case N
≥
3,(13)
D
=
A
1
B
1
C
2
A
2
B
2
⋱
⋱
⋱
C
N
-
1
A
N
-
1
B
N
-
1
C
N
A
N
,where, for m
,
j
=
0,1
,
…
,
M, we set(14)
A
m
j
=
N
2
m
+
1
t
f
1
-
-
1
m
+
j
-
q
m
j
,
A
1
m
j
=
N
2
m
+
1
t
f
1
2
-
-
1
m
+
j
-
q
m
j
,
A
n
m
j
=
N
2
m
+
1
t
f
1
2
-
-
1
m
+
j
2
-
q
m
j
,
n
=
2,3
,
…
,
N
-
1
,
A
N
m
j
=
N
2
m
+
1
t
f
1
-
-
1
m
+
j
2
-
q
m
j
,
B
n
m
j
=
N
2
m
+
1
t
f
-
1
j
2
,
n
=
1,2
,
…
,
N
-
1
,
C
n
m
j
=
N
2
m
+
1
t
f
-
-
1
m
2
,
n
=
2,3
,
…
,
N
,
(15)
q
m
j
=
2
,
i
f
∃
k
∈
Z
;
m
-
j
=
2
k
+
1
>
0
,
0
,
o
t
h
e
r
w
i
s
e
.


### 3.2. The Operational Matrix of Product

LetW be an arbitrary vector of order N
M
+
1
×
1. Then, the expression Φ
t
Φ
T
t
W can be expanded in terms of hybrid functions as follows:(16)
Φ
t
Φ
T
t
W
≃
W
~
Φ
t
,in which W
~ is a matrix of order N
M
+
1
×
N
M
+
1. This matrix is called the operational matrix of product and has the following structure:(17)
W
~
=
W
~
1
O
O
⋯
O
O
W
~
2
O
⋯
O
O
O
W
~
3
⋯
O
⋮
⋮
⋮
⋱
⋮
O
O
O
⋯
W
~
N
,where W
~
i, i
=
1,2
,
3
,
…
,
N, are M
+
1
×
M
+
1 matrices given in [15].


### 3.3. The Operational Matrix of Delay

The goal of this subsection is to determine the operational matrix of delay associated with the developed hybrid functions. For this purpose, let(18)
Φ
t
-
τ
t
=
D
Φ
t
,in which τ
t is defined by(19)
τ
t
=
τ
1
,
T
0
≤
t
<
T
1
,
τ
2
,
T
1
≤
t
<
T
2
,
⋮
⋮
τ
r
,
T
r
-
1
≤
t
≤
T
r
,where 0
=
T
0
<
T
1
<
T
2
<
⋯
<
T
r
-
1
<
T
r
=
t
f and τ
i, i
=
1
,
…
,
r, are known constants. Furthermore, assume that I
=
i
∶
τ
i
≠
0 is a nonempty set, τ
i, i
=
1,2
,
…
,
r, are nonnegative rational numbers, and T
i
∈
Q for i
=
0,1
,
…
,
r.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
,
t
f into N subintervals of the same length where the value of N is obtained in the following manner.Defineγ as the smallest positive integer number in such a way that(20)
γ
τ
i
∈
Z
,
i
∈
I
,
γ
T
j
∈
Z
,
j
=
0,1
,
…
,
r
.Suppose that I
=
l. Next, we choose λ as the greatest common divisor of the integers γ
τ
i and γ
T
j, i
∈
I, j
=
0,1
,
…
,
r; that is,(21)
λ
=
g.c.d
γ
τ
1
,
γ
τ
2
,
…
,
γ
τ
l
,
γ
T
1
,
…
,
γ
T
r
. Let(22)
N
=
γ
λ
t
f
,
i
f
γ
λ
t
f
∈
Z
,
γ
λ
t
f
+
1
,
o
t
h
e
r
w
i
s
e
,where [
γ
/
λ
t
f
] denotes the greatest integer value less than or equal to γ
/
λ
t
f. It should be noted that N is chosen in such a way that the number of subintervals is minimized. As a consequence, we obtain the following subintervals:(23)
0
,
h
,
h
,
2
h
,
…
,
N
-
1
h
,
N
h
,where h
=
λ
/
γ
=
t
f
/
N. Now, from (20) and (22) we can find the integer numbers k
i and l
j, such that(24)
τ
i
=
k
i
h
,
i
∈
I
,
T
j
=
l
j
h
,
j
=
0,1
,
…
,
r
.Note that k
i
=
0, for i
∉
I, 1
≤
i
≤
r. Therefore, we can rewrite τ
t as follows:(25)
τ
t
=
k
1
h
,
l
0
h
≤
t
<
l
1
h
,
k
2
h
,
l
1
h
≤
t
<
l
2
h
,
⋮
⋮
k
r
h
,
l
r
-
1
h
≤
t
≤
l
r
h
.For convenience, we define k
i as(26)
k
i
=
k
1
,
0
≤
i
≤
l
1
,
k
2
,
l
1
<
i
≤
l
2
,
⋮
⋮
k
r
,
l
r
-
1
<
i
≤
l
r
.Therefore, the problem is reduced to find the operational matrix of delay for the following delay function:(27)
Φ
t
-
τ
t
=
Φ
t
-
k
1
h
,
0
≤
t
<
t
1
,
Φ
t
-
k
2
h
,
t
1
≤
t
<
t
2
,
⋮
⋮
Φ
t
-
k
N
h
,
t
N
-
1
≤
t
<
t
N
,where (28)
t
i
=
l
i
h
,
i
=
1,2
,
…
,
N
.In order to find the matrixD, we first derive the matrix D
i for i
=
1,2
,
…
,
N, so that the following relation holds: (29)
Φ
t
-
k
i
h
=
D
i
Φ
t
,
t
i
-
1
≤
t
<
t
i
.With the use of (1), it is worth noting that, in the case of t
i
-
1
≤
t
<
t
i, the only functions with nonzero values are ϕ
i
-
k
i
m
t
-
k
i
h, for m
=
0,1
,
…
,
M. Since (30)
ϕ
i
-
k
i
m
t
-
k
i
h
=
ϕ
i
m
t
,
m
=
0,1
,
…
,
M
,by expanding ϕ
i
-
k
i
m
t
-
k
i
h in terms of ϕ
i
m
(
t
), we get the M
+
1
×
M
+
1 identity matrix as the corresponding coefficients. Therefore, we have (31)
D
i
=
S
i
⊗
I
M
+
1
,
i
=
1,2
,
…
,
N
,where I
M
+
1 is the M
+
1-dimensional identity matrix and S
i is an N
×
N matrix in which the only nonzero entry is equal to one and located at the i
-
k
ith row and ith column.Remark 1.
Ifi
-
k
i
≤
0, then S
i is a zero matrix of order N
×
N.As a consequence, if we expandΦ
t
-
τ
t in terms of hybrid functions Φ
t, we deduce that (32)
D
=
D
1
+
D
2
+
⋯
+
D
N
.To illustrate the derivation process of the operational matrix of delay associated with the mentioned hybrid functions, we present an example. Define(33)
τ
t
=
1
6
,
0
≤
t
<
1
3
,
1
2
,
1
3
≤
t
≤
1
.It is easily verified that N
=
6. As a result, we get (34)
k
i
=
1
,
1
≤
i
≤
2
,
3
,
3
≤
i
≤
6
.For i
=
1 and i
=
3, we conclude that i
-
k
i
=
0. Therefore, S
1 and S
3 are zero matrices of order 6
×
6. If i
=
2, it follows that i
-
k
i
=
1. Consequently, S
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 i
=
4, then we deduce that i
-
k
i
=
1. Hence, S
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 i
=
5 and i
=
6 to specify matrices S
5, and S
6. Using the above comments, we obtain(35)
D
1
=
S
1
⊗
I
M
=
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
D
2
=
S
2
⊗
I
M
=
0
I
M
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
D
3
=
S
3
⊗
I
M
=
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
D
4
=
S
4
⊗
I
M
=
0
0
0
I
M
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
D
5
=
S
5
⊗
I
M
=
0
0
0
0
0
0
0
0
0
0
I
M
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
D
6
=
S
6
⊗
I
M
=
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
I
M
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
.Consequently(36)
D
=
∑
i
=
1
6
D
i
=
0
I
M
0
I
M
0
0
0
0
0
0
I
M
0
0
0
0
0
0
I
M
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
.


## 3.1. The Operational Matrix of Derivative

We approximate the derivative ofu by the derivative of (3); that is,(7)
d
d
t
u
t
≃
d
d
t
∑
n
=
1
N
∑
m
=
0
M
u
^
n
m
ϕ
n
m
t
.The right hand side of the preceding equation can be represented in terms of hybrid functions as(8)
d
d
t
∑
n
=
1
N
∑
m
=
0
M
u
^
n
m
ϕ
n
m
t
=
∑
n
=
1
N
∑
m
=
0
M
s
^
n
m
ϕ
n
m
t
=
Φ
T
t
S
,where(9)
S
=
s
^
10
,
…
,
s
^
1
M
,
s
^
20
,
…
,
s
^
2
M
,
…
,
s
^
N
0
,
…
,
s
^
N
M
T
.The relationship between the two vectorsU and S is expressed by(10)
S
=
D
U
,where N
M
+
1
×
N
M
+
1 matrix D is the operational matrix of derivative.The structure of the operational matrix of derivativeD is given by [21].In the caseN
=
1,(11)
D
=
A
,in the case N
=
2,(12)
D
=
A
1
B
1
C
N
A
N
,and, in the case N
≥
3,(13)
D
=
A
1
B
1
C
2
A
2
B
2
⋱
⋱
⋱
C
N
-
1
A
N
-
1
B
N
-
1
C
N
A
N
,where, for m
,
j
=
0,1
,
…
,
M, we set(14)
A
m
j
=
N
2
m
+
1
t
f
1
-
-
1
m
+
j
-
q
m
j
,
A
1
m
j
=
N
2
m
+
1
t
f
1
2
-
-
1
m
+
j
-
q
m
j
,
A
n
m
j
=
N
2
m
+
1
t
f
1
2
-
-
1
m
+
j
2
-
q
m
j
,
n
=
2,3
,
…
,
N
-
1
,
A
N
m
j
=
N
2
m
+
1
t
f
1
-
-
1
m
+
j
2
-
q
m
j
,
B
n
m
j
=
N
2
m
+
1
t
f
-
1
j
2
,
n
=
1,2
,
…
,
N
-
1
,
C
n
m
j
=
N
2
m
+
1
t
f
-
-
1
m
2
,
n
=
2,3
,
…
,
N
,
(15)
q
m
j
=
2
,
i
f
∃
k
∈
Z
;
m
-
j
=
2
k
+
1
>
0
,
0
,
o
t
h
e
r
w
i
s
e
.


## 3.2. The Operational Matrix of Product

LetW be an arbitrary vector of order N
M
+
1
×
1. Then, the expression Φ
t
Φ
T
t
W can be expanded in terms of hybrid functions as follows:(16)
Φ
t
Φ
T
t
W
≃
W
~
Φ
t
,in which W
~ is a matrix of order N
M
+
1
×
N
M
+
1. This matrix is called the operational matrix of product and has the following structure:(17)
W
~
=
W
~
1
O
O
⋯
O
O
W
~
2
O
⋯
O
O
O
W
~
3
⋯
O
⋮
⋮
⋮
⋱
⋮
O
O
O
⋯
W
~
N
,where W
~
i, i
=
1,2
,
3
,
…
,
N, are M
+
1
×
M
+
1 matrices given in [15].


## 3.3. The Operational Matrix of Delay

The goal of this subsection is to determine the operational matrix of delay associated with the developed hybrid functions. For this purpose, let(18)
Φ
t
-
τ
t
=
D
Φ
t
,in which τ
t is defined by(19)
τ
t
=
τ
1
,
T
0
≤
t
<
T
1
,
τ
2
,
T
1
≤
t
<
T
2
,
⋮
⋮
τ
r
,
T
r
-
1
≤
t
≤
T
r
,where 0
=
T
0
<
T
1
<
T
2
<
⋯
<
T
r
-
1
<
T
r
=
t
f and τ
i, i
=
1
,
…
,
r, are known constants. Furthermore, assume that I
=
i
∶
τ
i
≠
0 is a nonempty set, τ
i, i
=
1,2
,
…
,
r, are nonnegative rational numbers, and T
i
∈
Q for i
=
0,1
,
…
,
r.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
,
t
f into N subintervals of the same length where the value of N is obtained in the following manner.Defineγ as the smallest positive integer number in such a way that(20)
γ
τ
i
∈
Z
,
i
∈
I
,
γ
T
j
∈
Z
,
j
=
0,1
,
…
,
r
.Suppose that I
=
l. Next, we choose λ as the greatest common divisor of the integers γ
τ
i and γ
T
j, i
∈
I, j
=
0,1
,
…
,
r; that is,(21)
λ
=
g.c.d
γ
τ
1
,
γ
τ
2
,
…
,
γ
τ
l
,
γ
T
1
,
…
,
γ
T
r
. Let(22)
N
=
γ
λ
t
f
,
i
f
γ
λ
t
f
∈
Z
,
γ
λ
t
f
+
1
,
o
t
h
e
r
w
i
s
e
,where [
γ
/
λ
t
f
] denotes the greatest integer value less than or equal to γ
/
λ
t
f. It should be noted that N is chosen in such a way that the number of subintervals is minimized. As a consequence, we obtain the following subintervals:(23)
0
,
h
,
h
,
2
h
,
…
,
N
-
1
h
,
N
h
,where h
=
λ
/
γ
=
t
f
/
N. Now, from (20) and (22) we can find the integer numbers k
i and l
j, such that(24)
τ
i
=
k
i
h
,
i
∈
I
,
T
j
=
l
j
h
,
j
=
0,1
,
…
,
r
.Note that k
i
=
0, for i
∉
I, 1
≤
i
≤
r. Therefore, we can rewrite τ
t as follows:(25)
τ
t
=
k
1
h
,
l
0
h
≤
t
<
l
1
h
,
k
2
h
,
l
1
h
≤
t
<
l
2
h
,
⋮
⋮
k
r
h
,
l
r
-
1
h
≤
t
≤
l
r
h
.For convenience, we define k
i as(26)
k
i
=
k
1
,
0
≤
i
≤
l
1
,
k
2
,
l
1
<
i
≤
l
2
,
⋮
⋮
k
r
,
l
r
-
1
<
i
≤
l
r
.Therefore, the problem is reduced to find the operational matrix of delay for the following delay function:(27)
Φ
t
-
τ
t
=
Φ
t
-
k
1
h
,
0
≤
t
<
t
1
,
Φ
t
-
k
2
h
,
t
1
≤
t
<
t
2
,
⋮
⋮
Φ
t
-
k
N
h
,
t
N
-
1
≤
t
<
t
N
,where (28)
t
i
=
l
i
h
,
i
=
1,2
,
…
,
N
.In order to find the matrixD, we first derive the matrix D
i for i
=
1,2
,
…
,
N, so that the following relation holds: (29)
Φ
t
-
k
i
h
=
D
i
Φ
t
,
t
i
-
1
≤
t
<
t
i
.With the use of (1), it is worth noting that, in the case of t
i
-
1
≤
t
<
t
i, the only functions with nonzero values are ϕ
i
-
k
i
m
t
-
k
i
h, for m
=
0,1
,
…
,
M. Since (30)
ϕ
i
-
k
i
m
t
-
k
i
h
=
ϕ
i
m
t
,
m
=
0,1
,
…
,
M
,by expanding ϕ
i
-
k
i
m
t
-
k
i
h in terms of ϕ
i
m
(
t
), we get the M
+
1
×
M
+
1 identity matrix as the corresponding coefficients. Therefore, we have (31)
D
i
=
S
i
⊗
I
M
+
1
,
i
=
1,2
,
…
,
N
,where I
M
+
1 is the M
+
1-dimensional identity matrix and S
i is an N
×
N matrix in which the only nonzero entry is equal to one and located at the i
-
k
ith row and ith column.Remark 1.
Ifi
-
k
i
≤
0, then S
i is a zero matrix of order N
×
N.As a consequence, if we expandΦ
t
-
τ
t in terms of hybrid functions Φ
t, we deduce that (32)
D
=
D
1
+
D
2
+
⋯
+
D
N
.To illustrate the derivation process of the operational matrix of delay associated with the mentioned hybrid functions, we present an example. Define(33)
τ
t
=
1
6
,
0
≤
t
<
1
3
,
1
2
,
1
3
≤
t
≤
1
.It is easily verified that N
=
6. As a result, we get (34)
k
i
=
1
,
1
≤
i
≤
2
,
3
,
3
≤
i
≤
6
.For i
=
1 and i
=
3, we conclude that i
-
k
i
=
0. Therefore, S
1 and S
3 are zero matrices of order 6
×
6. If i
=
2, it follows that i
-
k
i
=
1. Consequently, S
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 i
=
4, then we deduce that i
-
k
i
=
1. Hence, S
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 i
=
5 and i
=
6 to specify matrices S
5, and S
6. Using the above comments, we obtain(35)
D
1
=
S
1
⊗
I
M
=
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
D
2
=
S
2
⊗
I
M
=
0
I
M
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
D
3
=
S
3
⊗
I
M
=
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
D
4
=
S
4
⊗
I
M
=
0
0
0
I
M
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
D
5
=
S
5
⊗
I
M
=
0
0
0
0
0
0
0
0
0
0
I
M
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
D
6
=
S
6
⊗
I
M
=
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
I
M
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
.Consequently(36)
D
=
∑
i
=
1
6
D
i
=
0
I
M
0
I
M
0
0
0
0
0
0
I
M
0
0
0
0
0
0
I
M
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
.


## 4. Problem Statement and Its Approximation

### 4.1. Piecewise Constant Delay Systems

Consider the linear time-varying piecewise constant delay system described by(37)
x
˙
t
=
A
t
x
t
+
B
t
x
t
-
τ
t
+
C
t
u
t
,
0
≤
t
≤
t
f
,
(38)
x
t
=
φ
t
,
t
≤
0
,
(39)
τ
t
=
τ
1
,
T
0
≤
t
<
T
1
,
τ
2
,
T
1
≤
t
<
T
2
,
⋮
⋮
τ
r
,
T
r
-
1
≤
t
≤
T
r
,where x
t
∈
R
p, u
t
∈
R
q, and A
t, B
t, and C
t are matrices with appropriate dimensions. The problem is to find x
t, t
∈
0
,
t
f, satisfying (37)–(39).


### 4.2. Approximation Using Hybrid Functions

We approximate the system dynamics by hybrid functions as follows. Let(40)
x
t
=
x
1
t
,
x
2
t
,
…
,
x
p
t
T
,
u
t
=
u
1
t
,
u
2
t
,
…
,
u
q
t
T
.Using (3), each x
i
t and each u
j
t, i
=
1,2
,
…
,
p, j
=
1,2
,
…
,
q, can be written in terms of hybrid functions as follows:(41)
P
M
N
x
i
=
∑
n
=
1
N
∑
m
=
0
M
x
^
i
n
m
ϕ
n
m
t
=
Φ
T
t
X
i
,
P
M
N
u
j
=
∑
n
=
1
N
∑
m
=
0
M
u
^
j
n
m
ϕ
n
m
t
=
Φ
T
t
U
j
.From (40), it follows that(42)
P
M
N
x
t
=
I
p
⊗
Φ
T
t
X
,
P
M
N
u
t
=
I
q
⊗
Φ
T
t
U
,where I
p and I
q are the p
×
p and q
×
q identity matrices, respectively, and ⊗ denotes the Kronecker product [22]. It should be noted that X and U are vectors of orders p
N
M
+
1
×
1 and q
N
M
+
1
×
1, respectively, given by (43)
X
=
X
1
T
,
X
2
T
,
…
,
X
p
T
T
,
U
=
U
1
T
,
U
2
T
,
…
,
U
q
T
T
.Also(44)
φ
t
-
τ
t
=
I
p
⊗
Φ
T
t
R
,where R is a vector of order p
N
M
+
1
×
1 defined by (45)
R
=
R
1
T
,
R
2
T
,
…
,
R
p
T
T
.Similarly (46)
A
t
=
A
T
I
p
⊗
Φ
t
,
B
t
=
B
T
I
p
⊗
Φ
t
,
C
t
=
C
T
I
q
⊗
Φ
t
,where A
,
B, and C are matrices of orders p
N
M
+
1
×
p, p
N
M
+
1
×
p, and q
N
M
+
1
×
p, respectively. The delay vector x
t
-
τ
t can also be expanded in terms of hybrid functions as(47)
P
M
N
x
t
-
τ
t
=
I
p
⊗
Φ
T
t
R
+
I
p
⊗
Φ
T
t
I
p
⊗
D
T
X
,
0
≤
t
≤
t
f
,where D is the operational matrix of delay described by (18).Now, using (16), we have(48)
A
t
x
t
=
A
T
I
p
⊗
Φ
t
I
p
⊗
Φ
T
t
X
=
I
p
⊗
Φ
T
t
A
~
T
X
,
C
t
u
t
=
C
T
I
q
⊗
Φ
t
I
q
⊗
Φ
T
t
U
=
I
p
⊗
Φ
T
t
C
~
T
U
,where A
~ and C
~ can be calculated in a way similar to the construction method of matrix U
~ given by (17). Moreover, (49)
B
t
x
t
-
τ
t
=
B
T
I
p
⊗
Φ
t
I
p
⊗
Φ
T
t
R
+
I
p
⊗
D
T
X
.Consequently(50)
B
t
x
t
-
τ
t
=
I
p
⊗
Φ
T
t
B
~
T
R
+
I
p
⊗
D
T
X
.Now, using the operational matrix of derivative and substituting (44)–(50) in (37), we get(51)
I
p
⊗
Φ
T
t
I
p
⊗
D
X
=
I
p
⊗
Φ
T
t
A
~
T
X
+
I
p
⊗
Φ
T
t
C
~
T
U
+
I
p
⊗
Φ
T
t
B
~
T
R
+
I
p
⊗
D
T
X
. Because the elements of Φ
(
t
) are linearly independent functions over the interval 0
,
t
f, it follows that(52)
I
p
⊗
D
X
=
A
~
T
X
+
C
~
T
U
+
B
~
T
R
+
I
p
⊗
D
T
X
.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 X.


## 4.1. Piecewise Constant Delay Systems

Consider the linear time-varying piecewise constant delay system described by(37)
x
˙
t
=
A
t
x
t
+
B
t
x
t
-
τ
t
+
C
t
u
t
,
0
≤
t
≤
t
f
,
(38)
x
t
=
φ
t
,
t
≤
0
,
(39)
τ
t
=
τ
1
,
T
0
≤
t
<
T
1
,
τ
2
,
T
1
≤
t
<
T
2
,
⋮
⋮
τ
r
,
T
r
-
1
≤
t
≤
T
r
,where x
t
∈
R
p, u
t
∈
R
q, and A
t, B
t, and C
t are matrices with appropriate dimensions. The problem is to find x
t, t
∈
0
,
t
f, satisfying (37)–(39).


## 4.2. Approximation Using Hybrid Functions

We approximate the system dynamics by hybrid functions as follows. Let(40)
x
t
=
x
1
t
,
x
2
t
,
…
,
x
p
t
T
,
u
t
=
u
1
t
,
u
2
t
,
…
,
u
q
t
T
.Using (3), each x
i
t and each u
j
t, i
=
1,2
,
…
,
p, j
=
1,2
,
…
,
q, can be written in terms of hybrid functions as follows:(41)
P
M
N
x
i
=
∑
n
=
1
N
∑
m
=
0
M
x
^
i
n
m
ϕ
n
m
t
=
Φ
T
t
X
i
,
P
M
N
u
j
=
∑
n
=
1
N
∑
m
=
0
M
u
^
j
n
m
ϕ
n
m
t
=
Φ
T
t
U
j
.From (40), it follows that(42)
P
M
N
x
t
=
I
p
⊗
Φ
T
t
X
,
P
M
N
u
t
=
I
q
⊗
Φ
T
t
U
,where I
p and I
q are the p
×
p and q
×
q identity matrices, respectively, and ⊗ denotes the Kronecker product [22]. It should be noted that X and U are vectors of orders p
N
M
+
1
×
1 and q
N
M
+
1
×
1, respectively, given by (43)
X
=
X
1
T
,
X
2
T
,
…
,
X
p
T
T
,
U
=
U
1
T
,
U
2
T
,
…
,
U
q
T
T
.Also(44)
φ
t
-
τ
t
=
I
p
⊗
Φ
T
t
R
,where R is a vector of order p
N
M
+
1
×
1 defined by (45)
R
=
R
1
T
,
R
2
T
,
…
,
R
p
T
T
.Similarly (46)
A
t
=
A
T
I
p
⊗
Φ
t
,
B
t
=
B
T
I
p
⊗
Φ
t
,
C
t
=
C
T
I
q
⊗
Φ
t
,where A
,
B, and C are matrices of orders p
N
M
+
1
×
p, p
N
M
+
1
×
p, and q
N
M
+
1
×
p, respectively. The delay vector x
t
-
τ
t can also be expanded in terms of hybrid functions as(47)
P
M
N
x
t
-
τ
t
=
I
p
⊗
Φ
T
t
R
+
I
p
⊗
Φ
T
t
I
p
⊗
D
T
X
,
0
≤
t
≤
t
f
,where D is the operational matrix of delay described by (18).Now, using (16), we have(48)
A
t
x
t
=
A
T
I
p
⊗
Φ
t
I
p
⊗
Φ
T
t
X
=
I
p
⊗
Φ
T
t
A
~
T
X
,
C
t
u
t
=
C
T
I
q
⊗
Φ
t
I
q
⊗
Φ
T
t
U
=
I
p
⊗
Φ
T
t
C
~
T
U
,where A
~ and C
~ can be calculated in a way similar to the construction method of matrix U
~ given by (17). Moreover, (49)
B
t
x
t
-
τ
t
=
B
T
I
p
⊗
Φ
t
I
p
⊗
Φ
T
t
R
+
I
p
⊗
D
T
X
.Consequently(50)
B
t
x
t
-
τ
t
=
I
p
⊗
Φ
T
t
B
~
T
R
+
I
p
⊗
D
T
X
.Now, using the operational matrix of derivative and substituting (44)–(50) in (37), we get(51)
I
p
⊗
Φ
T
t
I
p
⊗
D
X
=
I
p
⊗
Φ
T
t
A
~
T
X
+
I
p
⊗
Φ
T
t
C
~
T
U
+
I
p
⊗
Φ
T
t
B
~
T
R
+
I
p
⊗
D
T
X
. Because the elements of Φ
(
t
) are linearly independent functions over the interval 0
,
t
f, it follows that(52)
I
p
⊗
D
X
=
A
~
T
X
+
C
~
T
U
+
B
~
T
R
+
I
p
⊗
D
T
X
.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 X.


## 5. Illustrative Examples

In this section, three examples are investigated to evaluate the performance of the method.Example 1.
Consider the time-varying piecewise constant delay system(53)
x
˙
t
=
t
2
x
t
-
τ
t
+
u
t
,
0
≤
t
≤
1
,
(54)
x
0
=
1
,
(55)
u
t
=
2
t
+
1
,
t
≥
0
,
(56)
τ
t
=
0.2
,
0
≤
t
<
0.3
,
0.7
,
0.3
≤
t
<
1
.The exact solution to this problem is given by [19](57)
x
t
=
1
+
t
+
t
2
,
0
≤
t
<
0.2
,
62341
62500
+
t
+
t
2
+
7
25
t
3
+
3
20
t
4
+
1
5
t
5
,
0.2
≤
t
<
0.3
,
1006717
1000000
+
t
+
t
2
,
0.3
≤
t
<
0.7
,
2720369
3000000
+
t
+
t
2
+
79
300
t
3
-
1
10
t
4
+
1
5
t
5
,
0.7
≤
t
<
0.9
,
76256550101
84000000000
+
t
+
t
2
+
693817
3000000
t
3
+
459
40000
t
4
+
167
5000
t
5
+
7
50
t
6
-
11
140
t
7
+
1
40
t
8
,
0.9
≤
t
≤
1
.To solve this problem by the method developed in the present paper, we first determine the value of N, the required number of subintervals, with the use of (22). For this problem, we choose N
=
10. Also, we select M
=
8. Let(58)
x
t
=
X
T
Φ
t
. By expanding u
(
t
) and t
2 in terms of hybrid functions, we get(59)
u
t
=
11
10
,
1
10
,
0,0
,
0,0
,
0,0
,
0
,
13
10
,
1
10
,
0,0
,
0,0
,
0,0
,
0
,
3
2
,
1
10
,
0,0
,
0,0
,
0,0
,
0
,
17
10
,
1
10
,
0,0
,
0,0
,
0,0
,
0
,
19
10
,
1
10
,
0,0
,
0,0
,
0,0
,
0
,
21
10
,
1
10
,
0,0
,
0,0
,
0,0
,
0
,
23
10
,
1
10
,
0,0
,
0,0
,
0,0
,
0
,
5
2
,
1
10
,
0,0
,
0,0
,
0,0
,
0
,
27
10
,
1
10
,
0,0
,
0,0
,
0,0
,
0,0
,
29
10
,
1
10
,
0,0
,
0,0
,
0,0
,
0
Φ
t
=
U
T
Φ
t
,
(60)
t
2
=
1
300
,
1
200
,
1
600
,
0,0
,
0,0
,
0,0
,
7
300
,
3
200
,
1
600
,
0,0
,
0,0
,
0,0
,
19
300
,
1
40
,
1
600
,
0,0
,
0,0
,
0,0
,
37
300
,
7
200
,
1
600
,
0,0
,
0,0
,
0,0
,
61
300
,
9
200
,
1
600
,
0,0
,
0,0
,
0,0
,
91
300
,
11
200
,
1
600
,
0,0
,
0,0
,
0,0
,
127
300
,
13
200
,
1
600
,
0,0
,
0,0
,
0,0
,
169
300
,
3
40
,
1
600
,
0,0
,
0,0
,
0,0
,
217
300
,
17
200
,
1
600
,
0,0
,
0,0
,
0,0
,
271
300
,
19
200
,
1
600
,
0,0
,
0,0
,
0,0
Φ
t
=
e
T
Φ
t
. Substituting (58)–(60) in (53) yields(61)
D
X
=
e
~
T
D
T
X
+
U
,with the following initial condition:(62)
1
=
∑
m
=
0
M
-
1
m
X
m
,where D is the delay operational matrix described by (18). By solving the resulting algebraic system, we can find the vector X as(63)
X
=
79
75
,
11
200
,
1
600
,
0,0
,
0,0
,
0,0
,
0
,
88
75
,
13
200
,
1
600
,
0,0
,
0,0
,
0,0
,
0
,
658097
500000
,
109663
1400000
,
6059
2800000
,
1
36000
,
1
1750000
,
1
126000000
,
0,0
,
0,0
,
4440151
3000000
,
17
200
,
1
600
,
0,0
,
0,0
,
0,0
,
0
,
4980151
3000000
,
19
200
,
1
600
,
0,0
,
0,0
,
0,0
,
0
,
5580151
3000000
,
21
200
,
1
600
,
0,0
,
0,0
,
0,0
,
0
,
6240151
3000000
,
23
200
,
1
600
,
0,0
,
0,0
,
0,0
,
0
,
4695907
2000000
,
433113
2800000
,
14701
4200000
,
49
900000
,
13
14000000
,
1
126000000
,
0,0
,
0,0
,
3215527
1200000
,
2484429
14000000
,
691
168000
,
77
1125000
,
3
2800000
,
1
126000000
,
0,0
,
0,0
,
1542687364553
504000000000
,
57254691267
280000000000
,
27335101823
5544000000000
,
360730541
3960000000000
,
12971011
8008000000000
,
10861
504000000000
,
7481
27720000000000
,
1
308000000000
,
1
51480000000000
,
0
T
.With the use of (58), we would obtain the same value as the exact value of x
(
t
).Example 2.
Consider the following piecewise constant delay system [13]:(64)
x
˙
t
=
-
5
x
t
-
5
x
t
-
τ
t
+
2
u
t
,
0
≤
t
≤
2
,
x
0
=
1
,
u
t
=
1
,
t
≥
0
,
τ
t
=
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 analytical solution to this problem is presented by [13](65)
x
t
=
0.2
+
0.8
e
-
10
t
,
0
≤
t
<
0.8
,
0.2
+
0.8
e
3
-
10
t
+
0.8
e
-
4
1
-
e
3
e
-
5
t
,
0.8
≤
t
<
1.1
,
0.2
+
0.8
e
5
-
10
t
-
4
e
-
2.5
1
-
e
3
t
e
-
5
t
+
0.8
1
-
e
3
6.5
e
-
2.5
+
e
-
4
e
-
5
t
,
1.1
≤
t
<
1.4
,
0.2
+
0.8
e
9
-
10
t
-
4
e
-
1
1
-
e
3
t
e
-
5
t
+
0.8
1
-
e
3
8
e
-
1
-
0.5
e
-
2.5
+
e
-
4
e
-
5
t
,
1.4
≤
t
≤
1.7
,
0.2
+
0.8
e
12
-
10
t
-
4
e
0.5
1
-
e
3
t
e
-
5
t
+
0.8
1
-
e
3
9.5
e
0.5
-
0.5
e
-
1
-
0.5
e
-
2.5
+
e
-
4
e
-
5
t
,
1.7
≤
t
≤
2
.To solve this problem by the proposed approach, we first determine the value of N. With the use of (22), we take N
=
20. In addition, let M be an arbitrary positive integer number. Define(66)
E
=
max
⁡
x
e
t
-
x
t
∶
0
≤
t
≤
2
,where x
(
t
) and x
e
(
t
) denote the approximate solution obtained by the present method and the exact solution, respectively. In Table 1, the maximum absolute error of x
(
t
) corresponding to N
=
20 and different values of M are summarized. This table shows that there is an excellent agreement between the analytical and approximate solutions. After choosing N, small values for M are needed to achieve a satisfactory approximation.Table 1
Results of maximum error for Example2.
Present method
E
M
=
2
2.05
e
-
02
M
=
4
2.64
e
-
05
M
=
6
2.87
e
-
08
M
=
8
2.12
e
-
11
M
=
10
1.09
e
-
14
M
=
12
4.11
e
-
18
M
=
14
1.17
e
-
21Example 3.
As a more complicated problem, consider the following nonlinear piecewise constant delay system:(67)
x
˙
t
=
t
x
t
-
τ
t
+
t
x
2
t
-
τ
t
+
x
3
t
-
τ
t
+
u
t
,
0
≤
t
≤
1
,
(68)
x
0
=
1
,
(69)
x
t
=
-
1
,
t
<
0
,
(70)
u
t
=
1
,
0
≤
t
<
0.8
,
-
1
,
0.8
≤
t
≤
1
,
(71)
τ
t
=
0.4
,
0
≤
t
<
0.8
,
0.8
,
0.8
≤
t
<
1
.The exact solution to this problem is described by(72)
x
t
=
1
,
0
≤
t
<
0.4
,
1
25
+
2
t
+
t
2
,
0.4
≤
t
<
0.8
,
41
25
+
t
2
,
0.8
≤
t
≤
1
.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 N
=
5. Moreover, we choose M
=
2. Let(73)
x
t
=
X
T
Φ
t
.By expressing u
(
t
), t, and φ
(
t
-
τ
(
t
)
) in terms of hybrid functions, we obtain(74)
u
t
=
1,0
,
0,1
,
0,0
,
1,0
,
0,1
,
0,0
,
-
1,0
,
0
Φ
t
=
U
T
Φ
t
,
(75)
t
=
1
10
,
1
10
,
0
,
3
10
,
1
10
,
0
,
1
2
,
1
10
,
0
,
7
10
,
1
10
,
0
,
9
10
,
1
10
,
0
Φ
t
=
A
T
Φ
t
,
(76)
φ
t
-
τ
t
=
-
1,0
,
0
,
-
1,0
,
0,0
,
0,0
,
0,0
,
0,0
,
0,0
Φ
t
=
R
T
Φ
t
.Substituting (73)–(75) in (67) implies(77)
D
X
=
A
~
T
K
+
A
~
T
K
~
T
K
+
K
~
2
T
K
+
U
,with the following initial condition:(78)
1
=
∑
m
=
0
M
-
1
m
X
m
,in which (79)
K
=
R
+
D
T
Xand D is the operational matrix of delay given by (18). Moreover, A
~ and K
~ can be determined in a way similar to the construction method of the matrix U
~ presented by (17). By solving (77)-(78), it follows that (80)
X
=
1,0
,
0,1
,
0,0
,
97
75
,
3
10
,
1
150
,
29
15
,
17
50
,
1
150
,
184
75
,
9
50
,
1
150
T
.Using (73), we would obtain the same value as the exact value of x
(
t
).


## 6. 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 block-pulse 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 ofN, small values for M are required to obtain an admissible approximation. It is worth noting that the correct choice of N 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.



---

*Source: 8903184-2016-08-07.xml*