# Stability Analysis of Analytical and Numerical Solutions to Nonlinear Delay Differential Equations with Variable Impulses

Authors: X. Liu; Y. M. Zeng
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

A stability theory of nonlinear impulsive delay differential equations (IDDEs) is established. Existing algorithm may not converge when the impulses are variable. A convergent numerical scheme is established for nonlinear delay differential equations with variable impulses. Some stability conditions of analytical and numerical solutions to IDDEs are given by the properties of delay differential equations without impulsive perturbations.


---

## Body

## 1. Introduction

In this paper, we study the stability of nonlinear impulsive differential equations:(1)x′t=ft,xt,xt-τ,t≥0,t≠kτ,Δx=rkxt,t=kτ,k=0,1,…,xt=φt,-τ≤t≤0,where f:[0,+∞)×Cq×Cq→Cq, rk≠-1, τ is a positive constant, φ(t) is a continuous function on [-τ,0), and limt→0-⁡φ(t) exists. Assume that Δx=x(t+0)-x(t).Definition 1 (see [1]).
x
(
t
) is said to be the solution of (1), if(1)
x(0)=1+r0limt→0-⁡φ(t),(2)
x(t) is differentiable and x′(t)=f(t,x(t),x(t-τ)) for t∈(0,+∞),t>0, t≠kτ,(3)
x(t) is left continuous in (0,+∞) and x(kτ+)=(1+rk)x(kτ), k=0,1,2,….Impulsive differential equations are widely used in actual modeling such as epidemic, optimal control and population dynamics (see [2, 3]). Up to now, extensive work has been done in the area of qualitative theories of IDDEs (see [4–6]). The numerical properties of impulsive differential equations attracted attentions of scholars since Ran et al.’s work in [7], which showed that the explicit Euler method is stable for impulsive differential equations, while the implicit Euler method is not. Some excellent works have been done in existing literature. In [8], the convergence of Euler method for linear IDDEs is studied. Zhang et al. studied the numerical stability of Runge-Kutta methods to IDDEs in [9, 10] and proved that backward Euler method and 2-stage Lobatto IIIC method can preserve the stability of nonlinear IDDEs with fixed impulses. Unfortunately, the numerical scheme may not converge when the impulses are variable. It is necessary to find a convergent stable method for IDDEs with variable pulses.The rest of this paper is organized as follows. In Section2, we establish a theorem that enables us to transform the stability of IDDEs to corresponding delay differential equations without impulsive perturbations. In Sections 3 and 4, stability of exact solutions is studied. In Section 5, a convergent numerical scheme is constructed. In Section 6, the stability of numerical solutions is researched. In Section 7, numerical experiments are given to confirm the conclusions.


## 2. Preliminaries

Recently, a new technique has been introduced in the stability analysis of exact and numerical solutions to impulsive differential equations by constructing equivalent equations (see [9–12]). But the equivalent equations may not work when the impulses are variable. In this section, we will propose an improved equivalent equation.Hypothesis 1.
Assume that the functionα(t) satisfies the following:(1)
α(t) is continuously differentiable on (kτ,(k+1)τ), k=0,1,….(2)
α(kτ-)=(1+rk)α(kτ+), k=0,1,….(3)
α(t)≠0, ∀t>0.Remark 2.
It is easy to check that(2)αt=1+rkt-k+1τ/τ,t∈kτ,k+1τ,(3)αt=∏i=1k11+rk,t∈kτ,k+1τsatisfy Hypothesis 1.Consider the equation(4)y′t=α′tαtyt+αtft,ytαt,yt-ταt-τ,t>0,yt=φt,-τ≤t≤0,where α(t) has been defined in Hypothesis 1. Here α′(t) is the left derivative of α(t).Corollary 3.
In the case thatrk=r, k=0,1,2,…, the transformation function (2) can be reduced to α(t)=(1/1+r){t/τ}-1. Then (4) is the same as the equivalent equations proposed in [9, 13].Theorem 4.
Assume that Hypothesis1 holds.(1)
Ifx(t) is the solution of (1), then y(t)=α(t)x(t) is the solution of (4).(2)
Ify(t) is the solution of (4), then x(t)=y(t)/α(t) is the solution of (1).Proof.
(1) It is obvious thaty(t) is continuous on each interval (kτ,(k+1)τ). On the other hand, for any k=0,1,2,…, (5)ykτ+=xkτ+αkτ+=1+rkxkτ-αkτ+=αkτ-xkτ-=ykτ-,which implies that y(t) is continuous in [0,+∞). It is easy to verify that y(t)=α(t)x(t) is the solution of (4).
(2) On each interval(kτ,(k+1)τ),(6)x′t=y′tαt-ytα′tα2t=α′t/αtyt+αtft,yt/αt,yt-τ/αt-ταt-ytα′tα2t=ft,ytαt,yt-ταt-τ=ft,xt,xt-τ.On the other hand, at each impulsive point t=kτ,(7)xkτ+=limt→kτ+ytαt=ykτ+αkτ+=1+rkykτ-αkτ-=1+rkxkτ-.It follows from Definition 1 that x(t) is the solution of (1).


## 3. Stability Analysis of Exact Solutions

Let〈·,·〉 be an inner product on Cq and · is the corresponding norm. We assume that there exist real constants α and β such that the function f in (1) satisfies (8)Rx1-x2,ft,x1,y-ft,x2,y≤αx1-x22,ft,x,y1-ft,x,y2≤βy1-y2.To study the stability of (1), we also consider the equation with differential initial value: (9)x-′t=ft,x-t,x-t-τ,t≥0,t≠kτ,Δx-=rkx-t,t=kτ,k=0,1,…,x-t=ψt,-τ≤t≤0.Based on the transformation(10)zt=αtx-t,we can get an equivalent equation to (9) as follows: (11)z′t=α′tαtzt+αtft,ztαt,zt-ταt-τ,t>0,zt=ψt,-τ≤t≤0.Definition 5.
The zero solution of (1) is said to be stable, if there exists a constant C such that(12)xt-x-t≤Csup-τ≤s≤0φs-ψs.Definition 6.
The zero solution of (1) is said to be asymptotically stable, if (13)limt→+∞xt-x-t=0.Lemma 7.
Assume that(14)Rα′tαt+α≤α0,βαtαt-τ≤β0hold. Then(1)
ifα0+β0≤0,(15)yt-zt≤sup-τ≤s≤0φs-ψs,(2)
ifα0+β0<0,(16)limt→+∞yt-zt=0.Proof.
Note that(17)12ddtyt-zt2=Ryt-zt,α′tαtyt-zt+αtft,ytαt,yt-ταt-τ-ft,ztαt,zt-ταt-τ=Rα′tαtyt-zt2+Ryt-zt,αtft,ytαt,yt-ταt-τ-ft,ztαt,yt-ταt-τ+Ryt-zt,αtft,ztαt,yt-ταt-τ-ft,ztαt,zt-ταt-τ≤Rα′tαt+αyt-zt2+βαtαt-τ·yt-zt·yt-τ-zt-τ≤α0yt-zt2+β0·yt-zt·yt-τ-zt-τ.Therefore,(18)ddtyt-zt≤α0yt-zt+β0yt-τ-zt-τ.It follows that(19)yt-zt≤eα0τ+β0α0eα0τ-1ksup-τ≤s≤0φs-ψs,t∈kτ,k+1τ.On the other hand, eα0τ+β0/α0(eα0τ-1)∈(0,1] when α0+β0≤0 and |eα0τ+β0/α0(eα0τ-1)|<1 when α0+β0<0. The proof is complete.By Theorem4 and Lemma 7, we can obtain the stability of (1).Theorem 8.
Assume that1/α(t) is bounded. Then(1)
ifα0+β0≤0, then zero solution of (1) is stable,(2)
ifα0+β0<0, then zero solution of (1) is asymptotically stable.To investigate the stability of (1) by Theorem 8, we need to establish a suitable function α(t). Different stability conditions may be given under different function α(t). In the next section, some specific stability conditions are given.


## 4. Some Specific Stability Conditions

In this section, some specific stability conditions are given based on the choice ofα(t) proposed in (2) and (3).Case 1 (takeα(t) as in (2)).
Takingα(t) as in (2), we can obtain the following conclusions by Theorem 8.Theorem 9.
Assume that(20)β1+rk1+rk-1t-k+1τ/τ≤β0,t∈kτ,k+1τ,ln⁡1+rkτ+α≤α0,k=0,1,2,….Then the zero solution of (1) is asymptotically stable if β0<-α0.Remark 10.
The zero solution of (1) is asymptotically stable if β<-(ln⁡(1+r)/τ+α), when rk=r, k=0,1,…. The zero solution of (1) is asymptotically stable if β<-α, when rk=0, k=0,1,…. Therefore, in the case that impulses are fixed or absent, the conclusions in Theorem 8 can be reduced to the same as in [9, 14].Corollary 11.
Assume that(21)1≥1+rk≥1+rk-1,k=0,1,2,….Then the zero solution of (1) is asymptotically stable if β<-α.Case 2 (takeα(t) as in (3)).
Takingα(t) as in (3), we can obtain the following conclusion by Theorem 8.Corollary 12.
If there exists a constantM, the following hold:(22)∏i=1k1+ri≤M,t∈kτ,k+1τ,(23)β1+rk≤β0,k=0,1,2,….The zero solution of (1) is stable when β0≤-α0 and is asymptotically stable when β0<-α0.Note that it is difficult to find sequencesrk satisfying conditions (22) and (23) simultaneously. It is necessary to find other stability conditions in this choice of α(t).Theorem 13.
Assume that the following holds:(24)1+ri≤1,i=1,2,…,1αt≤∏i=1k1+ri,t∈kτ,k+1τ.Then the zero solution of (1) is stable when β≤-α and is asymptotically stable when β<-α.Proof.
Note that while we takeα(t) as in (3), inequality (17) becomes(25)ddtyt-zt≤αyt-zt+β1+rkyt-τ-zt-τ,t∈kτ,k+1τ.Hence,(26)yt-zt≤∏i=1keατ+β1+riαeατ-1sup-τ≤s≤0φs-ψs.
By Theorem4 and inequality (17), for t∈(kτ,(k+1)τ], we have(27)xt-x-t=yt-ztαt≤1αt·∏i=1keατ+β1+riαeατ-1sup-τ≤s≤0⁡φs-ψs≤eατ+βαeατ-1ksup-τ≤s≤0⁡φs-ψs.Note that 0<eατ+β/α(eατ-1)≤1 if β≤-α and 0<eατ+β/α(eατ-1)<1 if β<-α. The proof is complete.


## 5. Numerical Process

In this section, we establish a convergent numerical process ofθ methods for (1).


### 5.1.θ Methods for Ordinary Differential Equations

The application of one-stageθ method in case of an ordinary differential equation, (28)y′t=ft,yt,yields(29)Yn=yn+hθftn+θh,Yn,yn+1=yn+hftn+θh,Yn,where h denotes the step size and Yn is an approximation of y(tn+θh).


### 5.2.θ Methods for IDDEs

Leth=τ/m be a given step size with integer m≥1. The application of one-stage θ methods, in case of (4) on interval (kτ,(k+1)τ], yields(30)yn+1=yn+hα′tn+θhαtn+θhYn+αtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ,Yn=yn+hθα′tn+θhαtn+θhYn+αtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ,where yn denotes an approximation to y(tn) for n>0. Based on Theorem 4, the solution of (1) can be approximated by (30) and(31)xn+1=yn+1αtn+1.In the same way, we can define the numerical solutions of (3).(32)zn+1=zn+hα′tn+θhαtn+θhZn+αtn+θhftn+θh,Znαtn+θh,Zn-mαtn+θh-τ,Zn=zn+hθα′tn+θhαtn+θhZn+αtn+θhftn+θh,Znαtn+θh,Zn-mαtn+θh-τ,x-n+1=zn+1αtn+1.Remark 14.
If we takeα(t) as in (2), numerical processes (30) and (31) can be reduced to(33)yn+1=yn+hln⁡1+rkτYn+αtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ,Yn=yn+hθln⁡1+rkτYn+αtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ,when tn∈(kτ,(k+1)τ].Remark 15.
If we takeα(t) as in (3), numerical processes (30) and (31) can be reduced to(34)yn+1=yn+hαtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ,Yn=yn+hθαtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ.The convergence of numerical processes (30) and (31) is described in the following theorem.Theorem 16.
Assume thatα(t) is continuously differentiable on (kτ,(k+1)τ), k=0,1,2,…, and 1/α(t) is bounded. Then numerical processes (30) and (31) are of order 1 when θ≠1/2 and of order 2 when θ=1/2.Proof.
On each interval(kτ,(k+1)τ], (4) reduces to a DDE. Therefore yn-y(tn)=O(h) for θ≠1/2 and yn-y(tn)=O(h2) for θ=1/2. On the other hand, |xn-x(tn)|=|1/α(tn)(yn-y(tn))|. The conclusion follows because α(t) is bounded.


## 5.1.θ Methods for Ordinary Differential Equations

The application of one-stageθ method in case of an ordinary differential equation, (28)y′t=ft,yt,yields(29)Yn=yn+hθftn+θh,Yn,yn+1=yn+hftn+θh,Yn,where h denotes the step size and Yn is an approximation of y(tn+θh).


## 5.2.θ Methods for IDDEs

Leth=τ/m be a given step size with integer m≥1. The application of one-stage θ methods, in case of (4) on interval (kτ,(k+1)τ], yields(30)yn+1=yn+hα′tn+θhαtn+θhYn+αtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ,Yn=yn+hθα′tn+θhαtn+θhYn+αtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ,where yn denotes an approximation to y(tn) for n>0. Based on Theorem 4, the solution of (1) can be approximated by (30) and(31)xn+1=yn+1αtn+1.In the same way, we can define the numerical solutions of (3).(32)zn+1=zn+hα′tn+θhαtn+θhZn+αtn+θhftn+θh,Znαtn+θh,Zn-mαtn+θh-τ,Zn=zn+hθα′tn+θhαtn+θhZn+αtn+θhftn+θh,Znαtn+θh,Zn-mαtn+θh-τ,x-n+1=zn+1αtn+1.Remark 14.
If we takeα(t) as in (2), numerical processes (30) and (31) can be reduced to(33)yn+1=yn+hln⁡1+rkτYn+αtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ,Yn=yn+hθln⁡1+rkτYn+αtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ,when tn∈(kτ,(k+1)τ].Remark 15.
If we takeα(t) as in (3), numerical processes (30) and (31) can be reduced to(34)yn+1=yn+hαtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ,Yn=yn+hθαtn+θhftn+θh,Ynαtn+θh,Yn-mαtn+θh-τ.The convergence of numerical processes (30) and (31) is described in the following theorem.Theorem 16.
Assume thatα(t) is continuously differentiable on (kτ,(k+1)τ), k=0,1,2,…, and 1/α(t) is bounded. Then numerical processes (30) and (31) are of order 1 when θ≠1/2 and of order 2 when θ=1/2.Proof.
On each interval(kτ,(k+1)τ], (4) reduces to a DDE. Therefore yn-y(tn)=O(h) for θ≠1/2 and yn-y(tn)=O(h2) for θ=1/2. On the other hand, |xn-x(tn)|=|1/α(tn)(yn-y(tn))|. The conclusion follows because α(t) is bounded.


## 6. Stability Analysis of Numerical Solutions

In this section, we study the stability property of numerical solutions. As usual, we expect that the numerical solution can reproduce the property of the exact solutions. The stability of the numerical solutions can be defined according to Theorem4.Definition 17.
A numerical method is said to be stable for (1), if there exists a constant C, such that the numerical solutions xn and x-n of (1) and (4) satisfy(35)xn-x-n≤C·max-τ≤s≤0φs-ψsfor every step size h under the constraint hm=τ, where m is a positive integer.Definition 18.
A numerical method is said to be asymptotically stable for (1), if the numerical solutions xn and x-n of (1) and (4) satisfy (35) and(36)limn→+∞xn-x-n=0for every step size h under the constraint hm=τ, where m is a positive integer.Denote(37)f^t,y,z=α′tαty+αtft,yαt,zαt-τ,δn=yn-zn,ωn=Yn-Zn,Kn=f^tn+θh,Yn,Yn-m-f^tn+θh,Zn,Yn-m,Hn=f^tn+θh,Zn,Yn-m-f^tn+θh,Zn,Zn-m,Δn=f^tn+θh,Yn,Yn-m-f^tn+θh,Zn,Yn-m.It follows from (30) that(38)ωn=δn+hθΔn,δn+1=δn+hΔn.Lemma 19.
Assume that (8) holds. Then(39)δn+12≤δn2+h2α0+β0ωn2+hβ0ωn-m2holds if θ≥1/2. Furthermore, if α0+β0≤0, inequality(40)δn+12≤δn2+hα0ωn2+hβ0ωn-m2holds.Proof.
By difference equation (38), it is obvious that(41)δn+12=δn+hΔn2=δn2+2hRωn,Δn-2θ-1h2RΔn,Δn≤δn2+2hRωn,Kn+Hn=δn2+2hRωn,Kn+2hRωn,Hn≤δn2+2hRα′tn+θhαtn+θh+αωn2+2hβαtn+θhαtn+θh-τ·ωn·ωn-m≤δn2+2hRα′tn+θhαtn+θh+α+hβαtn+θhαtn+θh-τωn2+hβαtn+θhαtn+θh-τ·ωn-m2≤δn2+h2α0+β0ωn2+hβ0ωn-m2.The proof is complete.Using a similar technique to the one in [15, 16], we can obtain the numerical stability of θ methods to (1).Theorem 20.
Assume that1/α(t) is bounded and θ∈[1/2,1]. Then the numerical methods (30) and (31) are asymptotically stable for (1) with(42)β0<22θ-122θ-12+1-α0.Proof.
It follows from (38) that(43)wn=θδn+1+1-θδn.Therefore,(44)wn2=wn,wn=θδn+1+1-θδn,θδn+1+1-θδn=θ2δn+12+1-θ2δn2+2θ1-θδn+1,δn.Note that(45)-δn+12+δn2≤2δn+1,δn≤δn+12+δn2.Hence,(46)wn2≤θδn+12+1-θδn2,wn2≥θ2θ-1δn+12+1-θ1-2θδn2.In view of inequality (39), we have(47)δn+12≤δn2+2α0+β0hwn2+β0hwn-m2≤δn2+β0hθδn+1-m2+1-θδn-m2+2θ-12α0+β0hθδn+12-1-θδn2.Hence,(48)δn+12≤1+2α0+β0h1-θ1-2θ1-2α0+β0hθ2θ-1δn2+β0h1-2α0+β0hθ2θ-1θδn+1-m2+1-θδn-m2≤1+2α0+β0h1-θ1-2θ1-2α0+β0hθ2θ-1δn2+β0h1-2α0+β0hθ2θ-1max⁡δn+1-m2,δn-m2.By induction, we have(49)δn+12≤1+2α0+β0h1-θ1-2θ+β0h1-2α0+β0hθ2θ-1ksup-τ≤s≤0φs-ψs2,tn∈kτ,k+1τ.The conclusion follows because (42) implies that(50)1+2α0+β0h1-θ1-2θ+β0h1-2α0+β0hθ2θ-1∈0,1.On the other hand,(51)xn+1-x-n+1=yn+1-zn+1αtn+1.Therefore, the conclusion holds because 1/α(t) is bounded.Corollary 21.
The implicit Euler method can preserve the asymptotic stability of (1) with β0<-α0.Corollary 22.
Assume that1/α(t) is bounded and θ∈[1/2,1]. Then the numerical solution (33) is asymptotically stable for (1) satisfying (42). Furthermore, if(52)1≥1+rk≥1+rk-1,k=0,1,2,…,the numerical solution (33) is asymptotically stable for (1) satisfying (53)β<22θ-122θ-12+1-α.It is necessary to find other stability criteria of numerical process (34), since the criteria in Theorem 20 are too tough for numerical scheme (34).Theorem 23.
Assume that1+rk≤1 and(54)α+β1+rk≤0,k=0,1,2,….Then the numerical solution xn calculated by numerical process (34) is asymptotically stable for (1) satisfying (55)β≤-1-2θ2α.Proof.
Note that (54) implies that α0+β0≤0. Hence, by Lemma 19, (56)δn+12≤δn2+hα0ωn2+hβ0ωn-m2,where(57)α0=α,β0=β1+rk,tn∈kτ,k+1τ.By inequality (46), we have(58)δn+12≤δn2+αhwn2+β1+rkhwn-m2≤δn2+β1+rkhθδn+1-m2+1-θδn-m2+2θ-1αhθδn+12-1-θδn2.
Therefore,(59)δn+12≤1+αh1-θ1-2θ1-αhθ2θ-1δn2+β/1+rkh1-αhθ2θ-1θδn+1-m2+1-θδn-m2≤1+αh1-θ1-2θ1-αhθ2θ-1δn2+β/1+rkh1-αhθ2θ-1max⁡δn+1-m2,δn-m2,tn∈kτ,k+1τ.Then, by induction, we have(60)δn+12≤∏i=0k1+αh1-θ1-2θ+β/1+rih1-αhθ2θ-1sup-τ≤s≤0φs-ψs2.Therefore,(61)xn+1-x-n+12=yn+1-zn+1αtn+12≤∏i=0k1+ri∏i=0k1+αh1-θ1-2θ+β/1+rih1-αhθ2θ-1sup-τ≤s≤0φs-ψs2≤1+αh1-θ1-2θ+βh1-αhθ2θ-1ksup-τ≤s≤0φs-ψs2.The conclusion follows, since (55) implies that(62)1+αh1-θ1-2θ+βh1-αhθ2θ-1∈0,1.Corollary 24.
The implicit Euler method can preserve the asymptotic stability of (1) with (54).Corollary 25.
If(63)rk=0,k=0,1,2,…,then (54) holds when β<-α. Then the conclusion in this paper coincides with the numerical stability of delay differential equations in [14].


## 7. Numerical Experiments

In this section, we will do some numerical experiments to illustrate the conclusion.Consider the following equations:(64)x′t=-xt-x3t+131+x2t-1,t≥0,t≠k,Δx=1k+1-1xk-,t=k,k=1,…,xt=t+1,-1≤t≤0,(65)x-′t=-x-t-x-3t+131+x-2t-1,t≥0,t≠k,Δx-=1k+1-1x-k-,t=k,k=1,…,x-t=2,-1≤t≤0.We can choose α(t)=(k+1)!, t∈[k,k+1). Then (64) and (65) can be transformed as(66)y′t=-yt-y3tα2t+αt31+y2t-1α2t-1,t≥0,t≠k,yt=t+1,-1≤t≤0,(67)z′t=-zt-z3tα2t+αt31+z2t-1α2t-1,t≥0,t≠k,zt=2,-1≤t≤0.We chose h=1/10 and calculated the numerical solutions by numerical process (34) in case of implicit Euler method. Figure 1 shows the relationship of the solutions between (64) and (66). It can be seen from Figure 2 that the difference of solutions to (64) and (65) tends to zero as n→∞, which coincides with the conclusion in the paper.Figure 1
Numerical solutions of (64) and (66) on [0,3].Figure 2
The difference of numerical solutions to (64) and (65) on [0,30].



---

*Source: 6723491-2017-11-13.xml*