High order Runge–Kutta methods for impulsive delay differential equations

Applied Mathematics and Computation 313 (2017) 12–23

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

High order Runge–Kutta methods for impulsive delay differential equationsR Gui-Lai Zhang School of Mathematics and Statistics, Northeastern University at Qinhuangdao, Hebei 066004, China

a r t i c l e

i n f o

a b s t r a c t

Keywords: Impulsive delay differential equation Runge–Kutta method Convergence

The purpose of this paper is to obtain high order Runge–Kutta methods for nonlinear impulsive delay differential equations (IDDEs). In order to achieve the purpose, continuity and smoothness of the exact solutions of IDDEs are analyzed. And then the discontinuous points and non-smooth points are chosen as a part of the nodes of the Runge–Kutta methods. Furthermore, it is proved that the pth order Runge–Kutta method for IDDEs is also convergent of order p. And some simple numerical examples are given to demonstrate the theoretical results. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Impulsive delay differential equations (IDDEs) arise naturally from a wide variety of applications such as threshold phenomena in biology, bursting rhythm models in medicine, optimal control models in economics, circuit networks and frequency modulated systems, etc. Recently, the corresponding theory of the exact solutions of IDDEs has been significantly developed, see [1–3,5,7–9,11,12] and references therein. However, as we all know, it is difficult, sometime maybe impossible, to get the explicit solutions for IDDEs. Therefore, it is necessary to investigate the numerical methods for IDDEs in order to use these works for mathematical simulations. In this paper, we study Runge–Kutta methods for the following equation:



x (t ) = f (t , x(t ), x(t − τ )), x(τw ) = Iw (x(τw− )), x(t ) = (t ),

t ∈ [t0 , T ], t = τw , w = 1, 2, . . . , N, w = 1, 2, . . . , N , t ∈ [t0 − τ , t0 ]

(1.1)

where x(τw ) = x(τw+ ) − x(τw− ), the positive constant τ is the delay,  is the initial function, and t0 < τ1 < τ2 < · · · < τN = T . Assume that the function f : [t0 , T ] × Cd × Cd → Cd is continuous in t and Lipschitz continuous with respect to the second and third variable in the following sense: there are two positive real constants α and β such that

 f (t, x1 , y1 ) − f (t, x2 , y2 ) ≤ αx1 − x2  + βy1 − y2 

(1.2)

for arbitrary t ∈ [t0 , T], x1 , x2 , y1 , y2 ∈ where  ·  is any convenient norm on And also assume that all the functions Iw , w = 1, 2, . . . , N satisfy Lipschitz condition, which implies that there is a positive constant γ such that Cd ,

Iw (x ) − Iw (y ) ≤ γ x − y, for ∀x, y ∈ Cd .

Cd .

(1.3)

R This work is supported by the Natural Science Foundation of Hebei Province, A2015501130, the Research Project of Higher School Science and Technology in Hebei Province, ZD2015211 and the Fundamental Research Funds for Central Universities, N152304007. E-mail address: [email protected]

http://dx.doi.org/10.1016/j.amc.2017.05.054 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

13

Ding et al. in [4] proved that the explicit Euler for linear IDDEs is convergent of order 1, and pionted out that high order Runge–Kutta methods for (1.1) is a difficult unsolved problem which needs new strategies to solve. My present paper is to solve this problem, and obtains a perfect result: the pth order Runge–Kutta for IDDE (1.1) is also convergent of order p when special stepsizes are chosen. The paper is organized as follows. In Section 2, convergence of the Runge–Kutta methods for a special case of Eq. (1.1) is studied. In Section 3, continuity and smoothness of the exact solutions of IDDEs are analyzed. In Section 4, convergence of the Runge–Kutta methods for (1.1) is studied. In Section 5, some numerical examples are given to confirm the theoretical results. In order to look at the main results in the right theoretical perspective, we give the final section, that is, the proof of Theorem 4.1. 2. A special case High order numerical methods for some special cases of (1.1) have been studied in [13,14]. In this section, we also study a special case of (1.1) which has a wide range of applications. But the convergence of the numerical methods for this equation is a problem still remaining unsolved. This equation is the special case of (1.1) as τw = t0 + nw τ , nw ∈ Sˆ ⊆ Z+ , w = 1, 2, · · · and ni < nj if i < j, that is, the following equation:



x (t ) = f (t , x(t ), x(t − τ )), x(t0 + nw τ ) = Iw (x(t0 + nw τ − )), x(t ) = (t ),

t ≥ t0 , t = t0 + nw τ , nw ∈ Sˆ , w = 1, 2, · · · , t ∈ [t0 − τ , t0 ],

(2.1)

where the positive constant τ is the delay. Obviously, all the discontinuous points and non-smooth points of the exact solution of (2.1) are contained in the set {t : t = t0 + nτ , n ∈ N}. Denote the approximations to the exact solution x(tk, l ) and x(t0 + (k + 1 )τ − ) by yk,l (0 ≤ l ≤ m − 1 ) and yk, m , respecτ , m is a positive integer. tively, where tk,l = t0 + kτ + lh, k ∈ N, h = m The numerical method for (2.1) can be constructed as follows:

⎧  j j j i ⎪ Yk,l = yk,l + h sj=1 ai j f (tk,l , Yk,l , Zk,l ), k ∈ N, i = 1, 2, . . . , s, ⎪ ⎪ s ⎨y i i i = y + h b f ( t , Y , Z ), l = 0, 1, 2, . . . , m − 1, k,l+1 k,l i=1 i k,l k,l k,l  ⎪ yk,m + Iw (yk,m ), if ∃w ∈ Sˆ such that k + 1 = nw, ⎪ ⎪ ⎩yk+1,0 = yk,m ,

(2.2)

otherwise,

where y0,0 = (t0 ), i Zk,l

=

 (t0 − τ + lh + ci h ), if k = 0, Yki−1,l ,

if k > 0,

i = t + (km + l + c )h = t + kτ + lh + c h, (k ∈ N, l = 0, 1, . . . , m − 1, i = 1, 2, . . . , v). and tk,l 0 0 i i

Theorem 2.1. Assume that f and  are sufficient smooth and all the functions Ik are Lipschitz continuous. If ci , aij , bj are the coefficients of a pth order Runge–Kutta method, then the same method (2.2) for (2.1) is convergent of order p. The theorem can be proved the same as Theorem 4.1 in this paper. 3. Continuity and smoothness In this section, we analyze two kinds of points of the exact solution of (1.1): discontinuous points (impulsive points); non-smooth points (the propagation of discontinuous points for derivatives of the solution of IDDE along the integration interval). Two simple examples are provided to observe the discontinuous points and non-smooth points of the solution of IDDE (1.1) as follows. Example 3.1. Consider the following equation:

⎧ ⎪  ⎪ ⎨x (t ) = −x(t − 1 ),

t ≥ 0, t =

x(t ) = x(t − ), ⎪ ⎪ ⎩ x(t ) = (t ) ≡ 1,

3k , k ∈ Z+ , 2 t ∈ [−1, 0],

3k , k ∈ Z+ , 2

t=

where x(t ) = x(t + ) − x(t − ) and Z+ = {1, 2, · · · }. Proceeding as described above, we obtain

x(t ) = 1 − t for t ∈ [0, 1], x(t ) = 1 − t +

(t − 1 )2 2



for t ∈ 1,



3 , 2

(3.1)

14

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

1 0.8 0.6

x’ disc

x’ discontinous

0.4

x’’ disc

x

0.2

x’’ disc

0

x’ disc

−0.2 −0.4

x disc

x’’’ disc

−0.6

x disc

−0.8 −1 −1

−0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

t Fig. 1. The exact solution of (3.1).

x

3

2

= 2x

3− 2

3 =− , 4





5 (t − 1 )2 3 −t + for t ∈ ,2 , 8 2 2 5 5 (t − 1 )2 (t − 2 )3 x(t ) = − t + − for t ∈ 2, , 8 2 3! 2 5

5 5t (t − 1 )2 (t − 2 )3 x(t ) = − − + − for t ∈ ,3 , 16 8 2 3! 2 17 x ( 3 ) = 2x ( 3− ) = − , 24 2 5t (t − 1 )2 (t − 2 )3 (t − 3 )4 7 x(t ) = − − + − + for t ∈ [3, ], 3 8 2 3! 4! 2 7  3 4 2 135 5t 5 ( t − 2 ) ( t − 3 ) x(t ) = − − + + (t − 1 )2 − + for t ∈ , 4 , etc. 3 64 16 16 3! 4! 2 x(t ) =

The solution is displayed in Fig. 1. We observe that despite the fact that the differential equation f and initial function  are C∞ , the solution of (3.1) has discontinuous points and non-smooth points. In fact, every point of the solution of (3.1) is infinitely differentiable except the points where t belongs to the following set

N

  3 t : t = + n, n ∈ N}, where N = {0, 1, 2, · · · . 2

Example 3.2. Consider the following equation:

⎧ ⎪  ⎪ ⎨x (t ) = −x(t − 1 ),

t ≥ 0, t =

x(t ) = 2x(t − ), ⎪ ⎪ ⎩ x(t ) = (t ) = −t,

2k , k ∈ Z+ , 3 t ∈ [−1, 0].

2k , k ∈ Z+ , 3

t=

Proceeding as described above, we obtain

x(t ) = x

2

3





1 1 2 (t − 1 )2 − for t ∈ 0, , 2 2 3



2− = 3x 3

4 =− , 3





25 1 2 (t − 1 )2 − for t ∈ ,1 , 2 18 3 4

37 t (t − 2 )3 x(t ) = − + − for t ∈ 1, , 18 2 3! 3 x(t ) =

(3.2)

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

15

1

x’’’ disc 0

x’’ discontinuous

x

−1

x’ disc

−2

x disc −3

x disc

−4 −5

−1

−0.5

0

0.5

1

1.5

2

t Fig. 2. The exact solution of (3.2).

x

4

3



4− = 3x 3

=−

217 , 54





767 t (t − 2 )3 4 5 + − for t ∈ , , 162 2 3! 3 3

1007 25t (t − 2 )3 5 x(t ) = − + − for t ∈ , 2 , etc. 162 18 3! 3 x(t ) = −

The solution is displayed in Fig. 2. We observe that despite the fact that the initial function satisfies  (0− ) = −(−1 ), the solution of (3.2) has discontinuous points and non-smooth points. In fact, every point of the solution of (3.2) is infinitely differentiable except the points where t belongs to the following set

   2 4 t : t = + n, n ∈ N} {t : t = + n, n ∈ N . 3 3

N

Remark 3.3. From the above two examples, we can obtain that the time t of discontinuous points and non-smooth points of the exact solution of (1.1) should be in the following set:





{t0 + nτ }





n∈N

 



{ τk + n τ } .

k∈Z + n∈N

4. Runge–Kutta methods In this section, we study convergence of Runge–Kutta methods for (1.1). The Eq. (1.1) is as follows:



x (t ) = f (t , x(t ), x(t − τ )), x(τw ) = Iw (x(τw− )), x(t ) = (t ),

t ∈ [t0 , T ], t = τw , w = 1, 2, . . . , N, w = 1, 2, . . . , N , t ∈ [t0 − τ , t0 ].

For arbitrary τ i , there is a nonnegative integer ni such that τi ∈ [t0 + ni τ , t0 + (ni + 1 )τ ), i = 1, 2, . . . , N. Assume that

S1 = {t0 }

N 

{τi − ni τ } = {t0 , ξ1 , ξ2 , . . . , ξv } ⊂ [t0 , t0 + τ )

i=1

where ξ i < ξ j when i < j. Hence the discontinuous points and non-smooth points should be in the following set n

S=

N 

{t0 + kτ , ξ1 + kτ , ξ2 + kτ , . . . , ξv + kτ }



[0, T ].

k=0

In order to ensure the convergence of the method, we take the discontinuous points and non-smooth points as nodes. For convenience, define ξv+1 := t0 + τ and the set S˜ := {1, 2, . . . , N}. • S1 = {t0 } which implies that for ∀τw , w ∈ S˜ , there is a positive integer nw such that τw = t0 + nw τ . This case has been solved in section 2.

16

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

Table 1 The errors of Runge–Kutta methods for (5.1) at t = 3. m

2-Lobatto IIIA

Three-order Heun method

Classical four-order method

100 200 400 800 1600 Ratio

3.884493996769711e−06 9.710893991154990e−07 2.427702203849913e−07 6.069242387829466e−08 1.517308899356973e−08 4.0 0 0 047167351188

1.463641086107970e−08 1.856239106867719e−09 2.336651150347535e−10 2.930923559407717e−11 3.669883841261878e−12 7.946955564490793

1.082192113699421e−10 6.585940126591083e−12 4.061473379834979e−13 2.522981823460668e−14 1.873501354054952e−15 15.553018476324334

• S1 = {t0 } which implies that {t0 } ⊂ S1 and ∃ξ1 ∈ S1 , t0 < ξ1 < t0 + τ . So this case will be studied in the following. − Denote the quantities yk,l (0 ≤ l ≤ m − 1 ) and yk, m are the approximations to x(tk, l ) and x(tk,m ), respectively, where



t0 + lh0 , if k = 0,  −1 t0 + ( ki=0 mhi ) + lhk , if k > 0,

tk,l =

k and hk = h¯ u , u = k −  v+1 (v + 1 ),  ·  denotes the floor function and h¯ 0 = ξ1m−t0 , h¯ i = ξi+1m−ξi , i = 1, 2, . . . , v, m is an positive integer. The numerical methods for (1.1) can be constructed as follows:

⎧ i  j j j Yk,l = yk,l + hk sj=1 ai j f (tk,l , Yk,l , Zk,l ), k ∈ N, i = 1, 2, . . . , s, ⎪ ⎪  ⎨y s i i i b f ( t , Y , Z ), l = 0, 1, 2, . . . , m − 1, i k,l+1 = yk,l + hk i=1 k,l k,l k,l  ˜ ⎪ y + I (y ), if ∃w ∈ S such that tk,m = τw , ⎪ ⎩yk+1,0 = k,m w k,m yk,m ,

(4.1)

otherwise,

where y0,0 = y(t0,0 ) = y(t0 ) = (t0 ), i Zk,l =

 i (tk,l − τ ), if k ≤ v, Yki−v−1,l , if k > v,

i =t and tk,l k,l + ci hk , k ∈ N, l = 0, 1, . . . , m − 1, i = 1, 2, . . . , s.

Theorem 4.1. Assume that f and  are sufficient smooth and all the functions Iw (y ), w = 1, 2, . . . , N are Lipschitz continuous. If ci , aij , bj are the coefficients of a pth order Runge–Kutta method, then the same method (4.1) for (1.1) is convergent of order p. The theorem is proved in the final section. 5. Numerical experiments In this section, we will consider two simple numerical examples. Example 5.1. Consider the example of paper [4], which is the following equation:



x (t ) = −2x(t ) − x(t − 1 ), x(2k ) = − 12 x(2k− ), x(t ) = 1 + t,

t ≥ 0, t = 2k, k ∈ Z+, (5.1) t ∈ [−1, 0].

The solution of (5.1) for t ∈ [−1, 3] has the following form:

x(t ) =

⎧ 1 + t, ⎪ ⎪ ⎪ t 1 3 −2t ⎪ ⎪ ⎨− 2 + 4 + 4 e ,

t ∈ [−1, 0] t ∈ [0, 1 )

(

)

− − 34 e2 te−2t + 34 + 34 e2 e−2t , 7 − 8t + 16 + 38 e4 t 2 e−2t − 32 e4 + 34 e2 + 38 + 98 e2 + 21 e4 e−2t , 16

⎪ ⎪ ⎪ ⎪ ⎪ ⎩

t 4

1 2

(

(

t ∈ [1, 2 )

)te−2t

)

t ∈ [2, 3 ).

Ding et al. in [4] only illustrate that explicit Euler method for (5.1) is convergent of order 1. In fact, the Runge–Kutta methods (2.2) for (5.1) conserve their orders of convergence for ordinary differential equations by Theorem 2.1 of the present paper. Table 1 also roughly illustrates this point. Example 5.2. Consider the following equation:

⎧  ⎨x (t ) = −x(t ) + 2x(t − 1 ), x(t ) = x(t − ), ⎩ x(t ) = 1,

t ≥ 0, t = t=

3w ,w 2

3w ,w 2

= 1, 2, · · · ,

= 1, 2, · · · ,

t ∈ [−1, 0].

(5.2)

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

17

Table 2 The errors of Runge–Kutta methods for (5.2) at t = 2. m

2-Lobatto IIIA

Three-order Heun method

Classical four-order method

10 20 40 80 160 Ratio

2.934298408607816e−04 7.331901067564317e−05 1.832735089335458e−05 4.581687636395770e−06 1.145412400926915e−06 4.0 0 0696522123113

1.385257435781995e−05 1.682553590232061e−06 2.073170191074780e−07 2.572888746499302e−08 3.204559728686718e−09 8.108876122170052

1.987836086847494e−07 1.209103750454688e−08 7.454872275047819e−10 4.627498384479623e−11 2.884359417976157e−12 16.203226672820126

Obviously, in this example, S1 = {0, 12 }. So all the step sizes hk = h¯ 0 = h¯ 1 = The Runge–Kutta method (4.1) for (5.2) is changed as follows:

1 2m ,

m is a positive integer.

⎧ i  j j Yk,l = yk,l + hk sj=1 ai j (−Yk,l + 2Zk,l ), k ∈ N, i = 1, 2, . . . , s, ⎪ ⎪  ⎨ s i i yk,l+1 = yk,l + hk i=1 bi (−Yk,l + 2Zk,l ), l = 0, 1, 2, . . . , m − 1,  ⎪ 2 y , if k = 3 w − 1 , w = 1, 2, · · · , k,m ⎪ ⎩yk+1,0 = yk,m ,

(5.3)

otherwise,

where y0,0 = 1,



i Zk,l =

1,

if k ≤ 1,

Yki−2,l ,

if k > 1.

− Obviously, yk,l (0 ≤ l ≤ m − 1 ) and yk, m are the approximations to x(tk, l ) and x(tk,m ), respectively, where tk,l = k 2

k 2

+ lhk =

l 2m , 0

+ ≤ l ≤ m, k = 0, 1, 2, · · · . By Theorem 4.1, if ci , aij , bj are the coefficients of a pth order Runge–Kutta method, then the same methods (5.3) for (5.2) is also convergent of order p. Table 2 also roughly illustrates this point. 6. Proof of Theorem 4.1 The aim of this section is to prove Theorem 4.1 by the technique of [6, pp. 341, Proof of Theorem 17.1]. First, the relationship between the exact solution of (1.1) and the solutions of a series of initial value problems of ordinary differential equations (ODEs) on different intervals is studied. Second, the relationship between the Runge–Kutta method (4.1) for (1.1) and the same method for the series of ODEs without impulsive perturbations is studied. Finally, we will prove that there is a constant C such that the global errors between the numerical solution obtained from (4.1) for (1.1) and the exact solution of (1.1) satisfy

x(tk,l ) − yk,l  ≤ Ch p , k = 0, 1, . . . , nN · (v + 1 ), l = 1, 2, . . . , m,

(6.1)

where h := max{h¯ 1 , h¯ 2 , . . . , h¯ v } (h → 0 when m → ∞). 6.1. The exact solutions’ relationship The relationship between the exact solution of (1.1) and the solutions of a series of initial value problems of ODEs is studied as follows. (I) Consider the exact solution of (1.1) on the v + 1 subintervals of [t0 , t0 + τ ]. (i) The exact solution x(t) of (1.1) satisfies x(t ) = x0 (t ), t ∈ [t0 , ξ 1 ), where x0 (t) is the solution of the following equation:



x0 (t ) = f (t , x0 (t ), (t − τ )), t ∈ [t0 , ξ1 ], x0 (t0 ) = (t0 ).

(6.2)

(ii) x(t ) = x1 (t ), t ∈ [ξ 1 , ξ 2 ), where x1 (t) is the solution of the following equation:

⎧  (t , x1 (t ), (t − τ )), t ∈ [ξ1 , ξ2 ], ⎨x1 (t ) = f  z1 + Iw1 (z1 ), if ∃w1 ∈ S˜ such that ξ1 = τw1 , ⎩x1 (ξ1 ) = z1

(6.3)

otherwise,

where z1 = x0 (ξ1 ) which can be solved from (6.2). (iii) Similarly, x(t ) = xi (t ) for t ∈ [ξi , ξi+1 ), i = 2, 3, . . . , v, and xi (t) is the solution of the following equation:

⎧  ⎨xi (t ) = f(t , xi (t ), (t − τ )), t ∈ [ξi , ξi+1 ], z + Iwi (zi ), if ∃wi ∈ S˜ such that ξi = τwi , ⎩xi (ξi ) = i zi

where zi = xi−1 (ξi ).

otherwise,

(6.4)

18

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

(II) Consider the exact solution of (1.1) on the v + 1 subintervals of [t0 + τ , t0 + 2τ ]. For convenience, denote



X1,0 (t ) =

X1,i (t ) =

x0 (t )

xv+1 (t )

, t ∈ [t0 + τ , ξ1 + τ ],

xi (t ) , t ∈ [ξi + τ , ξi+1 + τ ], i = 1, 2, . . . , v, xv+1+i (t )



F1 (t, X1,0 (t )) =

F1 (t, X1,i (t )) =

f (t − τ , x0 (t ), (t − 2τ )) , t ∈ [t0 + τ , ξ1 + τ ], f (t, xv+1 (t ), x0 (t ))

f (t − τ , xi (t ), (t − 2τ )) , t ∈ [ξi + τ , ξi+1 + τ ], f (t, xv+1+i (t ), xi (t ))

i = 1, 2, . . . , v. (i) x(t ) = xv+1 for t ∈ [t0 + τ , ξ1 + τ ), where xv+1 (t ) is the last component of X1,0 (t ) = (x0 (t ), xv+1 (t ))T which is the solution of the following equation:



X1 ,0 (t ) = F1 (t , X1,0 (t )), t ∈ [t0 + τ , ξ1 + τ ], X1,0 (t0 + τ ) = Xˆ1,0 ,

where the initial data

Xˆ1,0 =

(6.5)

⎧ ⎨zv+1 + Iwv+1 (zv+1 ), ,xˆ1,0,1 = if ∃wv+1 ∈ S˜ such that ξv+1 = τwv+1 , ⎩z , otherwise,



(t0 ) xˆ1,0,1

v+1

and zv+1 = xv (ξv+1 ) = xv (t0 + τ ). (ii) The solution x(t) of IDDE (1.1) satisfies x(t ) = xv+1+i (t ) for t ∈ [ξi + τ , ξi+1 + τ ), where xv+1+i (t ) is the last component of X1,i (t ) = (xi (t ), xv+1+i (t ))T which is the solution of the following equation



X1 ,i (t ) = F1 (t , X1,i (t )), t ∈ [ξi + τ , ξi+1 + τ ], X1,i (ξi + τ ) = Xˆ1,i ,

(6.6)

where the initial data

Xˆ1,i =

xˆ1,i,0 xˆ1,i,1

⎧ ⎨z j(v+1)+i + Iw j(v+1)+i (z j(v+1)+i ), , xˆ1,i, j = if ∃w j (v+1)i ∈ S˜ such that ξ j (v+1)+i = τw j(v+1)+i , ⎩z , otherwise, j (v+1 )+i

where z j (v+1 )+i = x j (v+1 )+i−1 (ξi + jτ ) which can be obtained in the front process, j = 0, 1. (III) Consider the exact solution of (1.1) on the subintervals of [t0 + nτ , t0 + (n + 1 )τ ], n = 2, 3, . . . , nN . Similarly, denote

⎛ x (t ) 0 ⎜ xv+1 (t ) Xn,0 (t ) = ⎜ .. ⎝ ⎛

⎞

⎟ ⎟, t ∈ [t0 + nτ , ξ1 + nτ ], ⎠

. xn(v+1) (t ) xi (t )

⎞

⎜ xv+1+i (t ) ⎟ ⎟, t ∈ [ξi + nτ , ξi+1 + nτ ], i = 1, 2, . . . , v, .. ⎝ ⎠

Xn,i (t ) = ⎜

. xn(v+1)+i (t )

⎛ f (t − nτ , x (t ), (t − (n + 1 )τ ))⎞ 0 ⎜ f (t − (n − 1 )τ , xv+1 (t ), x0 (t )) ⎟ ⎟, t ∈ [t0 + nτ , ξ1 + nτ ], Fn (t, Xn,0 (t )) = ⎜ .. ⎝ ⎠ . f (t , xn(v+1) (t ), x(n−1)(v+1) (t ))

⎛ f (t − nτ , x (t ), (t − (n + 1 )τ ))⎞ i

⎜ f (t − (n − 1 )τ , xv+1+i (t ), xi (t )) ⎟ ⎟, Fn (t, Xn,i (t )) = ⎜ .. ⎝ ⎠ . f (t , xn(v+1)+i (t ), x(n−1)(v+1)+i (t ))

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

19

t ∈ [ξi + nτ , ξi+1 + nτ ], i = 1, 2, . . . , v. (i) The solution x(t) of (1.1) satisfies x(t ) = xn(v+1 ) (t ), t ∈ [t0 + nτ , ξ1 + nτ ), where xn(v+1 ) (t ) is the last component of Xn,0 (t ) = (x0 (t ), xv+1 (t ), . . . , xn(v+1 ) (t ))T which is the solution of the following equation£º



 (t ) = F (t , X (t )), t ∈ [t + nτ , ξ + nτ ], Xn, n n,0 0 1 0

Xn,0 (t0 + nτ ) = Xˆn,0 ,

(6.7)

where the initial data

⎛(t )⎞ 0 ⎧ ⎜ xˆn,0,1 ⎟ ⎨zi(v+1) + Iwi(v+1) (zi(v+1) ), ⎜ ⎟ Xˆn,0 = ⎜ . ⎟, xˆn,0,i = if ∃wi(v+1) ∈ S˜ such that t0 + iτ = τwi(v+1) , ⎩ ⎝ .. ⎠ zi(v+1) , otherwise,

xˆn,0,n

(ii) The solution x(t) of (1.1) satisfies x(t ) = xn(v+1 )+i (t ), t ∈ [ξi + nτ , ξi+1 + nτ ), i = 1, 2, . . . , v, where xn(v+1 )+i (t ) is the last component of Xn,i (t ) = (xi (t ), xv+1+i (t ), . . . , xn(v+1 )+i (t ))T which is the solution of the following equation£º



 (t ) = F (t , X (t )), t ∈ [ξ + nτ , ξ Xn,i n n,i i i+1 + nτ ],

Xn,i (ξi + nτ ) = Xˆn,i ,

(6.8)

where the initial data

⎛xˆ

⎞

⎧ ⎜xˆn,i,1 ⎟ ⎨z j(v+1)+i + Iw j(v+1)+i (z j(v+1)+i ), ⎜ ⎟ xˆn,i = ⎜ . ⎟, xˆn,i, j = if ∃w j (v+1)+i ∈ S˜ such that ξi + jτ = τw j(v+1)+i , ⎩ ⎝ .. ⎠ n,i,0

z j (v+1)+i , otherwise,

xˆn,i,n

where z j (v+1 )+i = x j (v+1 )+i−1 (ξi + jτ ) which can be obtained in the front process, j = 0, 1, . . . , n. 6.2. The numerical solutions’ relationship The relationship between the numerical solution obtained from the Runge–Kutta method (4.1) for (1.1) and the numerical solutions obtained from the same method for the above series of ODEs is studied as follows. (I) Obviously, for k ≤ v, the numerical solution obtained from the Runge–Kutta method (4.1) for (1.1) is the same as the numerical solutions obtained from the same method for the series of ODEs with another group initial data. (i)Obviously, for k = 0, the Runge–Kutta method (4.1) for (1.1) is the same as the same method for the following ODE:



x˜0 (t ) = f (t , x˜0 (t ), (t − τ )), t ∈ [t0 , ξ1 ], x˜0 (t0 ) = (t0 );

(6.9)

(ii) for k = 1, the same as the Runge–Kutta method for the following ODE:



x˜1 (t ) = f (t , x˜1 (t ), (t − τ )),

t ∈ [ ξ 1 , ξ2 ]

x˜1 (ξ1 ) = y1,0 ;

(6.10)

(iii) for k = i ≤ v, the same as the Runge–Kutta method for the following ODE:



x˜i (t ) = f (t , x˜i (t ), (t − τ )),

t ∈ [ξi , ξi+1 ],

x˜i (ξi ) = yi,0 .

(6.11)

(II) In fact, for v + 1 ≤ k < 2(v + 1 ), the numerical solution obtained from the Runge–Kutta method (4.1) for (1.1) equals the last components of the numerical solutions obtained from the same method for the series of ODEs as follows. (i) For k = v + 1, the numerical solution yv+1,l (l = 1, 2, . . . , m ) obtained from the Runge–Kutta method (4.1) for (1.1) is the last component of y˜1,0,l obtained from the following Runge–Kutta method for (6.13):

⎧ i s j j ⎨Y˜1,0,l = y˜1,0,l + hv+1 j=1 ai j F1 (t˜1,0,l , Y˜1,0,l ), i = 1, 2, . . . , s, s y˜ = y˜1,0,l + hv+1 i=1 bi F1 (t˜1i ,0,l , Y˜1i ,0,l ), l = 1, 2, . . . , m − 1, ⎩ 1,0,l+1 y˜1,0,0

(6.12)

= X¯ 1,0 ,

where = tv+1,l = t0 + τ + (l + c j )hv+1 . In fact, the relationship between the numerical solution obtained from (4.1) and the numerical solution obtained from (6.12) can be written as j t˜1,0,l

j



y˜1,0,l =

y0,l

yv+1,l

, l = 0, 1, . . . , m.

20

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

The following equation is the ODE (6.5) with another initial data:



X˜1 ,0 (t ) = F1 (t , X˜1,0 (t )), t ∈ [t0 + τ , ξ1 + τ ], X˜1,0 (t0 + τ ) = X¯ 1,0 ,

(6.13)

where the initial data



(t0 )

X¯ 1,0 =

.

yv+1,0

(ii) Similarly, for k = v + 1 + i, (i = 1, 2, . . . , v ), the numerical solution yv+1+i,l (l = 1, 2, . . . , m ) obtained from the Runge–Kutta method (4.1) for (1.1) is the last component of y˜1,i,l obtained from the following Runge–Kutta method for (6.15):

⎧ i s j j ⎨Y˜1,i,l = y˜1,i,l + hv+1+i j=1 ai j F1 (t˜1,i,l , Y˜1,i,l ), i = 1, 2, . . . , s, s y˜ = y˜1,i,l + hv+1+i i=1 bi F1 (t˜1i ,i,l , Y˜1i ,i,l ), l = 1, 2, . . . , m − 1, ⎩ 1,i,l+1

(6.14)

y˜1,i,0 = X¯ 1,i ,

 j j where t˜1,i,l = tv+1+i,l = t0 + τ + ik−1 m · hv+1+k + (l + c j )hv+1+i . In fact, the relationship between the numerical solu=0 tion obtained from (4.1) and the numerical solution obtained from (6.14) can be written as

y˜1,i,l =

yi,l

, i = 1, . . . , v, l = 0, 1, . . . , m.

yv+1+i,l

The following equation is the ODE (6.6) with another initial data:



X˜1 ,i (t ) = F1 (t , X˜1,i (t )), t ∈ [ξi + τ , ξi+1 + τ ], X˜1,i (ξi + τ ) = X¯ 1,i ,

where the initial data



X¯ 1,i =

yi,0

(6.15)

yv+1+i,0

.

(III) Similarly, for n(v + 1 ) ≤ k < (n + 1 )(v + 1 ), the numerical solution obtained from the Runge–Kutta method (4.1) for IDDE (1.1) equals the last components of the numerical solutions obtained from the same Runge–Kutta method for the series of ODEs as follows. (i) For k = n(v + 1 ), the numerical solution yn(v+1 ),l (l = 1, 2, . . . , m ) obtained from the Runge–Kutta method (4.1) for (1.1) is the last component of y˜n,0,l obtained from the following Runge–Kutta method for (6.17):

⎧ i s j j ⎨Y˜n,0,l = y˜n,0,l + hn(v+1) j=1 ai j Fn (t˜n,0,l , Y˜n,0,l ), i = 1, 2, . . . , s, s i y˜ = y˜n,0,l + hn(v+1) i=1 bi Fn (t˜n,0,l , Y˜n,i 0,l ), l = 1, 2, . . . , m − 1, ⎩ n,0,l+1 y˜n,0,0

(6.16)

= X¯ n,0 ,

where t˜n,0,l = tni (v+1 ),l = t0 + nτ + (l + c j )hn(v+1 ) . In fact, the relationship between the numerical solution obtained from (4.1) and the numerical solution obtained from (6.16) can be written as j

y˜n,0,l

⎛ y 0,l ⎜ yv+1,l ⎜ =⎜ ⎝ ...

⎞ ⎟ ⎟ ⎟, l = 0, 1, . . . , m. ⎠

yn(v+1),l The following equation is the ODE (6.7) with another initial data:



 (t ) = F (t , X˜ (t )), t ∈ [t + nτ , ξ + nτ ], X˜n, n n,0 0 1 0 X˜n,0 (t0 + nτ ) = X¯ n,0 ,

where the initial data

(6.17)

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

⎛ (t ) 0 ⎜ yv+1,0 ⎜ X¯ n,0 = ⎜ ⎝ ...

21

⎞ ⎟ ⎟ ⎟. ⎠

yn(v+1),0 (ii) Similarly, for k = n(v + 1 ) + i, (i = 1, 2, . . . , v ), the numerical solution yn(v+1 )+i,l (l = 1, 2, . . . , m ) obtained from the Runge–Kutta method (4.1) for IDDE (1.1) is the last component of y˜n,i,l obtained from the following Runge–Kutta method for ODE (6.19):

⎧ i s j j ⎨Y˜n,i,l = y˜n,i,l + hv+1+i j=1 ai j Fn (t˜n,i,l , Y˜n,i,l ), i = 1, 2, . . . , s, s i i y˜ = y˜n,i,l + hv+1+i i=1 bi Fn (t˜n,i,l , Y˜n,i,l ), l = 1, 2, . . . , m − 1, ⎩ n,i,l+1 y˜n,i,0

where

j t˜n,i,l

(6.18)

= X¯ n,i ,

= tn(v+1 )+i,l = t0 + nτ + j

i−1

k=0

m · hn(v+1 )+k + (l + c j )hn(v+1 )+i . In fact, the relationship between the numeri-

cal solution obtained from (4.1) and the numerical solution obtained from (6.18) can be written as

⎛

⎞

yi,l

⎜ yv+1+i,l ⎜ .. ⎝ .

⎟ ⎟ ⎟, i = 1, . . . , v, l = 0, 1, . . . , m. ⎠

y˜n,i,l = ⎜

yn(v+1)+i,l The following equation is the ODE (6.8) with another initial data:



 (t ) = F (t , X˜ (t )), t ∈ [ξ + nτ , ξ X˜n,i n n,i i i+1 + nτ ], X˜n,i (ξi + nτ ) = X¯ n,i ,

where the initial data

⎛

yi,0

(6.19)

⎞

⎜ y(v+1)+i,0 ⎟ ⎟ ⎟. .. ⎝ ⎠ .

⎜ X¯ n,i = ⎜

yn(v+1)+i,0 6.3. Global errors Consider the global errors between the numerical solution of (4.1) and the exact solution of (1.1) by mathematical induction. (I) Consider the global errors of Runge–Kutta method (4.1) for (1.1) as k = 0, 1, 2, . . . , v + 1. (i) When k = 0, there is a constant C0 such that the global errors satisfy:

x(t0,l ) − y0,l  = x0 (t0,l ) − y0,l  = x˜0 (t0,l ) − y0,l  ≤ C0 h0p ≤ C0 h p , for ∀l = 1, 2, · · · , m − 1, and

x(ξ1− ) − y0,m  = x˜0 (ξ1 ) − y0,m  ≤ C0 h p , this can be interpreted as the convergence of the Runge–Kutta method for ODE (6.9). Obviously, if ξ1 = τ1 , we can obtain that

x(ξ1 ) − y1,0  = x(ξ1− ) + I1 (x(ξ1− )) − y0,m − I1 (y0,m ) = x˜0 (ξ1 ) + I1 (x˜0 (ξ1 )) − y0,m − I1 (y0,m ) ≤ x˜0 (ξ1 ) − y0,m  + I1 (x˜0 (ξ1 )) − I1 (y0,m ) ≤ (1 + γ )x˜0 (ξ1 ) − y0,m  ≤ (1 + γ )C0 h p , otherwise ξ1 = τw for ∀w ∈ S˜ , we also can obtain

x(ξ1 ) − y1,0  = x(ξ1− ) − y1,0  = x˜0 (ξ1 ) − y1,0  = x˜0 (ξ1 ) − y0,m  ≤ C0 h p ≤ (1 + γ )C0 h p .

22

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

(ii) When k = 1, because f satisfies Lipschitz condition and the Runge–Kutta method is pth order, we can obtain that for ∀l = 1, 2, · · · , m − 1,

x(t1,l ) − y1,l  = x1 (t1,l ) − y1,l  ≤ x1 (t1,l ) − x˜1 (t1,l ) + x˜1 (t1,l ) − y1,l  ≤ x1 (t1,l ) − x˜1 (t1,l ) + C˜1 h1p ≤ eα (ξ2 −ξ1 ) x1 (ξ1 ) − x˜1 (ξ1 ) + C˜1 h1p ≤ eατ x1 (ξ1 ) − y1,0  + C˜1 h p ≤ C1 h p where C˜1 is a positive constant and C1 = eατ (1 + γ )C0 + C˜1 . Similarly, we can obtain

x(t1−,m ) − yk,m  = x˜1 (tk,m ) − yk,m  ≤ C1 h p ,

(6.20)

x(t2,0 ) − y2,0  ≤ (1 + γ )C1 h p .

(6.21)

and

(iii) Similarly, for ∀k = 2, 3, · · · , v and ∀l = 1, 2, · · · , m − 1, we have

x(tk,l ) − yk,l  ≤ Ck h p ,

(6.22)

− x(tk,m ) − yk,m  = xk (tk,m ) − yk,m  ≤ Ck h p ,

(6.23)

x(tk+1,0 ) − yk+1,0  ≤ (1 + γ )Ck h p .

(6.24)

and

(II) Consider the global errors of Runge–Kutta method (4.1) for (1.1) as k = v + 1, (v + 1 ) + 1, . . . , (v + 1 ) + v. Because f satisfies Lipschitz condition which implies that F1 also satisfies Lipschitz condition, there is a constant L1 such that

F1 (t, X1 ) − F1 (t, X2 ) ≤ L1 X1 − X2  for arbitrary t ∈ [t0 + τ , t0 + 2τ ] and X1 , X2 ∈ C2d . And because the Runge–Kutta method is pth order, by [10, pp. 109, Theorem 2.1.12], we can obtain that for ∀l = 1, 2, · · · , m − 1,

x(tv+1,l ) − y1,l  ≤ X1,0 (tv+1,l ) − y˜1,0,l  ≤ X1,0 (tv+1,l ) − X˜1,0 (tv+1,l ) + X˜1,0 (tv+1,l ) − y˜1,0,l  ≤ X0 (tv+1,l ) − X˜0 (tv+1,l ) + C˜v+1 hvp+1 ≤ eL1 (ξ1 −t0 ) Xˆ1,0 − X¯ 1,0  + C˜v+1 hvp+1 ≤ Cv+1 h p where C˜v+1 is a positive constant and Cv+1 = eL1 τ (1 + γ )Cv + C˜v+1 . Similarly, we can obtain

x(tv−+1,m ) − yv+1,m  = xv+1 (tv+1,m ) − yv+1,m  ≤ Cv+1 h p ,

(6.25)

x(tv+2,0 ) − yv+2,0  ≤ (1 + γ )Cv+1 h p .

(6.26)

and

Similarly, we can obtain that for ∀k = v + 2, v + 3, . . . , 2(v + 1 ) and ∀l = 1, 2, · · · , m − 1,

x(tk,l ) − yk,l  ≤ Ck h p ,

(6.27)

− x(tk,m ) − yk,m  = xk (tk,m ) − yk,m  ≤ Ck h p

(6.28)

x(tk+1,0 ) − yk+1,0  ≤ (1 + γ )Ck h p .

(6.29)

and

(III) Assume that inequalities (6.27), (6.28) and (6.29) hold for all k < n(v + 1 ). Because f satisfies Lipschitz condition which implies that Fn also satisfies Lipschitz condition, there is a constant Ln such that

Fn (t, X1 ) − Fn (t, X2 ) ≤ Ln X1 − X2  for arbitrary t ∈ [t0 + nτ , t0 + (n + 1 )τ ] and X1 , X2 ∈ C(n+1 )d . And because the Runge–Kutta method is pth order, by [10, pp. 109, Theorem 2.1.12], we can obtain that for ∀l = 1, 2, · · · , m − 1,

x(tn(v+1),l ) − yn(v+1),l  ≤ Xn,0 (tv+1,l ) − y˜n,0,l  ≤ Xn,0 (tn(v+1),l ) − X˜n,0 (tn(v+1),l ) + X˜n,0 (tn(v+1),l ) − yn,0,l  ≤ Xn,0 (tn(v+1),l ) − X˜n,0 (tn(v+1),l ) + C˜n(v+1) hnp(v+1) ≤ eLn (ξ1 −t0 ) Xˆn,0 − X¯ n,0  + C˜n(v+1) hnp(v+1)

G.-L. Zhang / Applied Mathematics and Computation 313 (2017) 12–23

23

≤ eLn τ Xˆn,0 − X¯ n,0  + C˜n(v+1) h p where C˜n(v+1 ) is a positive constant. By the inductive assumption, we have x(t0 + iτ ) − yi(n+1 ),0  ≤ (1 + γ )Ci(v+1)−1 h p , i = 1, 2, . . . , n. So there is a constant C˘n such that

Xˆn,0 − X¯n,0  ≤ C˘n h p . Hence for ∀l = 1, 2, · · · , m − 1,

x(tn(v+1),l ) − yn(v+1),l  ≤ Cn(v+1) h p where Cn(v+1 ) = eLn τ C˘n + C˜n(v+1 ) . And we have

x(tn−(v+1),m ) − yn(v+1),m  = xn(v+1) (tn(v+1),m ) − yn(v+1),m  ≤ Cn(v+1) h p and

x(tn(v+1)+1,0 ) − yn(v+1)+1,0  ≤ (1 + γ )Cn(v+1) h p . Similarly, for i = 1, 2, . . . , v, l = 1, 2, . . . , m − 1, we can obtain that there is a positive constant Cn(v+1 )+i such that

x(tn(v+1)+i,l ) − yn(v+1)+i,l  ≤ Cn(v+1)+i h p and

x(tn−(v+1)+i,m ) − yn(v+1)+i,m  ≤ Cn(v+1)+i h p . Consequently, by mathematical induction, for ∀k = 1, 2, . . . , (v + 1 )nN and ∀l = 1, 2, · · · , m − 1, we can obtain that

x(tk,l ) − yk,l  ≤ Ch p and − x(tk,m ) − yk,m  = xk (tk,m ) − yk,m  ≤ Ch p ,

where C = max{C1 , C2 , . . . , C(v+1 )nN }. Hence the Runge–Kutta method (4.1) for (1.1) is convergent of order p. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]

A. Anokhin, L. Berezansky, E. Braverman, Exponential stability of linear delay impulsive differential equations, J. Math. Anal. Appl. 193 (1995) 923–941. L. Berezansky, E. Braverman, Exponential boundedness of solutions for impulsive delay differential equations, Appl. Math. Lett. 9 (1996) 91–95. G. Ballinger, X.Z. Liu, Existence and uniqueness results for impulsive delay differential equations, Dyn. Contin. Discrete Impuls. Syst. 5 (1999) 579–591. X.H. Ding, K.N. Wu, M.Z. Liu, The Euler scheme and its convergence for impulsive delay differential equations, Appl. Math. Comput. 216 (2010) 1566–1570. K. Gopalsamy, B.G. Zhang, On delay differential equations with impulses, J. Math. Anal. Appl. 139 (1989) 110–122. E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems, Springer-Verlag, New York, 1993. V. Lakshmikantham, D.D. Bainov, P.S. Simeonov, Theory of Impulsive Differential Equations, World Scientific, Singapore, 1989. X.Z. Liu, G. Ballinger, Uniform asymptotic stability of impulsive delay differential equations, Comput. Math. Appl. 41 (2001) 903–915. B. Liu, X.Z. Liu, K. Teo, Q. Wang, Razumikhin-type theorem on exponential stability of impulsive delay systems, IMA J. Appl. Math. 71 (2006) 47–61. A.M. Stuart, A.R. Humphries, Dynamical System and Numerical Analysis, Cambridge University Press, Britain, 1998. J.H. Shen, J.R. Yan, Razumikhin-type stability theorems for impulsive function differential equations, Nonlinear Anal. 33 (1998) 519–531. Q. Wang, X.Z. Liu, Exponential stability for impulsive delay differential equations by Razumikhin method, J. Math. Anal. Appl. 309 (2005) 462–473. G.L. Zhang, M.H. Song, M.Z. Liu, Asymptotic stability of a class of impulsive delay differential equations, J. Appl. Math. (2012), doi:10.1155/2012/723893. G.L. Zhang, M.H. Song, M.Z. Liu, Exponential stability of the exact solutions and the numerical solutions for a class of linear impulsive delay differential equations, J. Comput. Appl. Math. 285 (2015) 32–44.