# Nonlinear Stability and Convergence of Two-Step Runge-Kutta Methods for Neutral Delay Differential Equations

Authors: Haiyan Yuan; Cheng Song; Peichen Wang
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper is devoted to the stability and convergence analysis of the two-step Runge-Kutta (TSRK) methods with the Lagrange interpolation of the numerical solution for nonlinear neutral delay differential equations. Nonlinear stability and D-convergence are introduced and proved. We discuss theGR(l)-stability,  GAR(l)-stability, and the weak  GAR(l)-stability on the basis of (k,l)-algebraically stable of the TSRK methods; we also discuss the D-convergence properties of TSRK methods with a restricted type of interpolation procedure.


---

## Body

## 1. Introduction

Neutral delay differential equations (NDDEs) arise in a variety of fields as biology, economy, control theory, and electrodynamics (see, e.g., [1–5]). The stability and convergence properties of numerical methods for linear NDDEs have been widely researched by many authors (see, e.g., [6–11]). For the case of nonlinear delay differential equations, this kind of methodology had been first introduced by Bellen and Zennaro [12] and Torelli [13] and then developed by Torelli [14], Bellen [15], and Zennaro [16, 17]. In 1997, Koto proved the asymptotic stability of natural Runge-Kutta method for a class of nonlinear delay differential equations in [18]. Bellen et al. [19] gave a discussion of the stability of continuous numerical methods for a special class of nonlinear neutral delay differential equations. In particular, Jackiewicz [20–22] systematically investigated the convergence of various numerical methods for more general neutral functional differential equations (NFDEs). In 2009, Yang et al. gave a novel robust stability criteria for stochastic Hopfield neural networks with time delays in [23]. Yang et al. [24] studied the exponential stability on stochastic neural networks with discrete interval and distributed delays in 2010. In 2011, Liu [25] gave the robust stability for Neutral time-varying delay systems with non-linear peturbations. On the stability, Tanikawa studied the values of random zero-sum games in [26], and in [27] Basin and Calderon-Alyarez gave the delay-dependent stability studies for vector nonlinear stochastic systems with multiple delays.However, these important convergence results are based on the classical Lipschitz conditions. The studies focusing on the stability and convergence of the numerical method for nonlinear NDDEs based on a one-sided Lipschitz condition have not yet been seen in literature until now. By means of a one-sided Lipschitz condition, in the present paper we discuss the stability and convergence of two-step Runge-Kutta (TSRK) methods for nonlinear NDDEs. Thanks to the one-sided nature of the Lipschitz condition, the error bounds obtained in the present paper are sharper than those given in the references mentioned.


## 2. Two-Step Runge-Kutta Methods for NDDEs

It is the purpose of this paper to investigate the nonlinear stability and convergence properties of the following NDDEs:(1)y′(t)=f(t,y(t),y(t-τ),y′(t-τ)),t∈[0,T],y(t)=φ(t),t∈[-τ,0],
where f:[0,T]×CN×CN×CN→CN is a given mapping, τ is a positive delay term, and φ:[-τ,0]→CN is a continuous function. Moreover, we assume that there exist some inner product 〈·,·〉 and the induced norm ∥·∥ in CN, such that
(2)Re〈f(t,y1,u,v)-f(t,y2,u,v),y1-y2〉≤α∥y1-y2∥2,∥f(t,y,u1,v)-f(t,y,u2,v)∥≤β∥u1-u2∥,∥f(t,y,u,v1)-f(t,y,u,v2)∥≤γ∥v1-v2∥,∥K(t,y,u1,v,w)-K(t,y,u2,v,w)∥≤q∥u1-u2∥,
where K(t,y,u,v,w)=f(t,y,u,f(t-τ,u,v,w)), t∈[0,T], ∀y,y1,y2,u,u1,u2,v,v1,v2,w∈CN, and α, β, γ, q are constants.In order to make the error analysis feasible, we always assume that problem (1) has a unique solution y(t) which is sufficiently differentiable and satisfies
(3)∥diy(t)dti∥≤Mi,i=1,2,…,t∈[-τ,T],
and denotes the problem class R(α,β,γ,q) that consists of all NDDEs with (2).Many numerical methods have been proposed for the numerical solution of problem (1).In this paper, we are concerned with two-step Runge-Kutta (TSRK) method of the form(4a)Yi(n)=yn+h×∑j=1saijf(tn+cjh,Yj(n)),i=1,2,…,s,(4b)Yi(n-1)=yn-1+h×∑j=1saijf(tn-1+cjh,Yj(n-1)),i=1,2,…,s,(4c)yn+1=(1-θ)yn+θyn-1+h∑i=1sbif(tn+cih,Yi(n))+h∑i=1sb~if(tn-1+cih,Yi(n-1)),where ∑i=1sbi+∑i=1sb~i=1+θ, ci=∑j=1saij, yn is the numerical approximation at tn=nh to the analytic solution y(tn), h>0 is a step size, and 0≤θ≤1. The above methods are studied in [11]. Now we consider the adaptation of the two-step Runge-Kutta method to (1):(5a)Yi(n)=yn+h×∑j=1saijf(tn+cjh,Yj(n),Y-j(n),Y~j(n)),i=1,2,…,s,(5b)Yi(n-1)=yn-1+h×∑j=1saijf(Y~j(n-1)tn-1+cjh,Yj(n-1),Y-j(n-1),Y~j(n-1)),i=1,2,…,s,(5c)yn+1=(1-θ)yn+θyn-1+h∑i=1sbif(tn+cih,Yi(n),Y-i(n),Y~i(n))+h∑i=1sb~if(tn-1+cih,Yi(n-1),Y-i(n-1),Y~i(n-1)),where yn is the numerical approximation to the analytic solution y(tn) with tn=nh.In particular,y0=φ(0). The argument Y-j(n) denotes an approximation to y(tn+cjh-τ) and the argument Y~j(n) denotes an approximation to y′(tn+cjh-τ) which are obtained by a specific interpolation procedure at the point t=tn+cjh-τ. Using values Yj(k)=φ(tk+cjh) with k<1, tk+cjh≤0.Letτ=(m-δ)h with integer m and δ∈[0,1), v,μ≥0 be integers. Define
(6)Y-j(n)={φ(tn+cjh-τ)tn+cjh-τ≤0,∑i=-μvL~i(δ)Yj(n-m+i)tn+cjh-τ>0,v+1≤m,
where
(7)L~i(δ)=∏k=-μk≠iv(δ-ki-k),δ∈[0,1),Y~j(n)=f(tn+cjh-τ,Y-j(n),Y-j(n+1),Y~j(n+1)).
We assume m≥v+1 is to guarantee that no (unknown) values Yj(i) with i≥n are used in the interpolation procedure.It should be pointed out that the adopted interpolation procedures (6) is only a class of  interpolation procedure for Y-j(n); there also exist some other types of interpolation procedures, such as numerical schemes which use Hermite interpolation between grid points (see [28–30]). It is the aim of our future research to investigate the future adaptation of two-step Runge-Kutta methods to NDDEs by means of other interpolation procedures.


## 3. The Nonlinear Stability Analysis

In this section, we will investigate the stability of the two-step Runge-Kutta methods for NDDEs.In order to consider the stability property, we also need to consider the perturbed problem of (1):
(8)z′(t)=f(t,z(t),z(t-τ),z′(t-τ)),t∈[0,T],z(t)=ψ(t),t∈[-τ,0],
where ψ:[-τ,0]→CN is a continuous function. The unique exact solution of the problem (8) is denoted as z(t).Applying the two-step Runge-Kutta method (4a), (4b), and (4c) to (8) leads to(9a)Zi(n)=zn+h×∑j=1saijf(tn+cjh,Zj(n),Z-j(n),Z~j(n)),i=1,2,…,s,(9b)Zi(n-1)=zn-1+h×∑j=1saijf(tn-1+cjh,Zj(n-1),Z-j(n-1),Z~j(n-1)),i=1,2,…,s,(9c)zn+1=(1-θ)zn+θzn-1+h∑i=1sbif(tn+cih,Zi(n),Z-i(n),Z~i(n))+h∑i=1sb~if(tn-1+cih,Zi(n-1),Z-i(n-1),Z~i(n-1)).


### 3.1. Some Concepts

Definition 1.
Letl be a real constant, a two-step Runge-Kutta method with an interpolation procedure is said to be R(l)-stable if there exists a constant C dependent only on the method and l such that (α+β+γq)h≤l and
(10)∥yn-zn∥2≤[2+C(β+γq)]maxt≤0∥φ(t)-ψ(t)∥2,n≥1,
with step size h satisfying hm=τ, where m is a positive integer.
GR
(
l
)-stability is defined by dropping the restriction hm=τ.Definition 2.
Letl be a real constant, a two-step Runge-Kutta method with an interpolation procedure is said to be AR(l)-stable if
(11)limn→∞∥yn-zn∥=0
with step size h satisfying (α+β+γq)h≤l and hm=τ, where m is a positive integer.
GAR
(
l
)-stability is defined by dropping the restriction hm=τ.Definition 3.
Letl be a real constant; a two-step Runge-Kutta method with an interpolation procedure is said to be weak AR(l)-stable if, under the conditions of Definition 2, (11) holds when f further satisfies
(12)∥f(t,u1,v,w)-f(t,u2,v,w)∥≤L∥u1-u2∥p0
with p0 being a positive real number and L being a nonnegative real number.
WeakGAR(l)-stability is defined by dropping the restriction hm=τ.


### 3.2. The Stability of TSRK Methods

Letwn=yn-zn, Wi(n)=Yi(n)-Zi(n), W-i(n)=Y-i(n)-Z-i(n), W~i(n)=Y~i(n)-Z~i(n),  Qi(n)=h[f(tn+cih,Yi(n),Y-i(n),Y~i(n))-f(tn+cih,Zi(n),Z-i(n),Z~i(n))], and i=1,2,…,s.It follows from (5a), (5b), (5c), (9a), (9b), and (9c) that(13a)Wi(n)=wn+∑j=1saijQj(n),i=1,2,…,s,(13b)Wi(n-1)=wn-1+∑j=1saijQj(n-1),i=1,2,⋯,s,(13c)wn+1=(1-θ)wn+θwn-1+∑i=1sbiQi(n)+∑i=1sb~iQi(n-1).Now we will write the s-stage TSRK methods (4a), (4b), and (4c) as a general linear method.LetVi(n)=(Wi(n),Wi(n-1))T be the internal stages, μn+1=(wn+1,wn)T the external vectors, and Pi(n)=(Qi(n),Qi(n-1))T. Then we have a 2(s+1)-stage partitioned general linear method:(14a)Vi(n)=∑j=1sCij11Pi(n)+∑j=1sCij12μn,i=1,2,…,s,(14b)μn+1=∑j=1sCij21Pi(n)+∑j=1sCij22μn,i=1,2,…,s,where
(15)C11=(Cij11)=(A00A),C12=(Cij12)=(e00e),C21=(Cij21)=(bb~00),C22=(Cij22)=(1-θθ00),A=(aij),e=(1,1,…1),b=(b1,b2,…,bs),b~=(b~1,b~2,…,b~s).Definition 4.
Letk,l be real constants; a TSRK method is said to be (k,l)-algebraically stable if there exists a diagonal nonnegative matrix G and D=diag(d1,d2,…,d2s) such that M=(mij) is nonnegative, where(16)
M
(
k
,
l
)
=
(
k
G
-
C
22
T
G
C
22
-
2
l
C
12
T
D
C
12
C
12
T
D
-
C
22
T
G
C
21
-
2
l
C
12
T
D
C
11
D
C
12
-
C
21
T
G
C
22
-
2
l
C
11
T
D
C
12
C
11
T
D
+
D
C
11
-
C
21
T
G
C
21
-
2
l
C
11
T
D
C
11
)
.In this paper, we use the linear interpolation procedure. Letτ=(m-δ)h with integer m≥1 and δ∈[0,1).Define(17a)Y-j(n)=δYj(n-m+1)+(1-δ)Yj(n-m),(17b)Z-j(n)=δZj(n-m+1)+(1-δ)Zj(n-m),j=1,2,…,s,where Yj(i)=φ(ti+cjh) and Zj(i)=ψ(ti+cjh) for i<0. When the step size h satisfies τ=mh, we have
(18)Y-j(n)=Yj(n-m),Z-j(n)=Zj(n-m),j=1,2,…,s.Theorem 5.
Assume that a TSRK method is(k,l)-algebraically stable. Then
(19)∥μn+1∥2≤k∥μn∥2+∑j=12sdj[((2α+β+γq)h-2l)∥Vj(n)∥2+h(β+γq)×(δ∥Vj(n-m)∥2+(1-δ)∥Vj(n-m-1)∥2)].Proof.
It is well known that(20)∥μn+1∥2-k∥μn∥2-2∑j=12sdjRe〈Vj(n),Pj(n)-lVj(n)〉=-∑i=12s+2∑j=12s+2Mij〈ri,rj〉,
where
(21)r1=wn+1,r2=wn,rj=Qj-2(n),j=3,4,…,s+2,rj=Qj-s-2(n-1),j=s+3,s+4,…,2s+2.
By means of (k,l)-algebraical stability of the method, we have
(22)∥μn+1∥2≤k∥μn∥2+2∑j=12sdjRe〈Vj(n),Pj(n)-lVj(n)〉.
It follows from (2) and (6) that
(23)2Re〈Wj(n),Qj(n)〉=2hRe〈Wj(n),f(tn+cjh,Yj(n),Y-j(n),Y~j(n))-f(tn+cjh,Zj(n),Y-j(n),Y~j(n))〉+2hRe〈Wj(n),f(tn+cjh,Zj(n),Y-j(n),Y~j(n))-f(tn+cjh,Zj(n),Z-j(n),Y~j(n))〉+2hRe〈Wj(n),f(tn+cjh,Zj(n),Z-j(n),Y~j(n))-f(tn+cjh,Zj(n),Z-j(n),Z~j(n))〉≤2hα∥Wj(n)∥2+2hβ∥Wj(n)∥·∥W-j(n)∥+2hγ∥Wj(n)∥·∥W~j(n)∥≤2hα∥Wj(n)∥2+2hβ∥Wj(n)∥·∥W-j(n)∥+2hγq∥Wj(n)∥·∥W-j(n)∥≤2hα∥Wj(n)∥2+2h(β+γq)∥Wj(n)∥·∥W-j(n)∥≤2hα∥Wj(n)∥2+2h(β+γq)∥Wj(n)∥·(δ∥Wj(n-m)∥+(1-δ)∥Wj(n-1-m)∥)≤h(2α+β+γq)∥Wj(n)∥2+h(β+γq)(δ∥Wj(n-m)∥2+(1-δ)∥Wj(n-1-m)∥2).
Substitution into (20) gives (19).Theorem 6.
Assume that a TSRK method is(k,l)-algebraically stable and k≤1. Then the method with linear interpolation procedure is GR(l)-stable.Proof.
The inequality(α+β+γq)h≤l and Theorem 5 lead to
(24)∥μn+1∥2≤∥μn∥2+∑j=12sdj[-(β+γq)h∥Vj(n)∥2+h(β+γq)ai×(δ∥Vj(n-m)∥2+(1-δ)∥Vj(n-m-1)∥2)].
By induction, we have
(25)∥μn+1∥2≤∥μ0∥2+h∑j=12sdj(β+γq)(δ∑i=-m-1∥Vj(i)∥2+(1-δ)∑i=-m-1-1∥Vj(i)∥2)≤∥μ0∥2+h∑j=12sdj(β+γq)×[δ(m-1)+(1-δ)m]max-m-1≤i≤-1∥Vj(i)∥2=∥μ0∥2+(β+γq)τ∑j=12sdjmax-m-1≤i≤-1∥Vj(i)∥2≤∥μ0∥2+(β+γq)τ∑j=12sdjmaxt≤0∥φ(t)-ψ(t)∥2.
Because μn+1=(wn+1,wn)T, so we have
(26)∥wn∥2≤∥w0∥2+∥w-1∥2+(β+γq)τ×∑j=12sdjmaxt≤0∥φ(t)-ψ(t)∥2≤2maxt≤0∥φ(t)-ψ(t)∥2+(β+γq)τ×∑j=12sdjmaxt≤0∥φ(t)-ψ(t)∥2≤[2+C(β+γq)τ]maxt≤0∥φ(t)-ψ(t)∥2.
Therefore, it isGR(l)-stable, where C=∑j=12sdj.Theorem 7.
Assume that a TSRK method is(k,l)-algebraically stable and k<1. Then the method with linear interpolation procedure is GAR(l)-stable.Proof.
Letμ=(2α+β+γq)h-2l and k-=max{k,(β+γq)h/(-μ)1/(m+1)}.
Then, when(α+β+γq)h<l, we have μ<-(β+γq)h and 0<k-<1.
The application of Theorem5 yields
(27)∥μn+1∥2≤k-∥μn∥2+∑j=12sdj[μ∥Vj(n)∥2+h(β+γq)×(δ∥Vj(n-m)∥2+(1-δ)∥Vj(n-m-1)∥2)].
By induction, we have
(28)∥μn+1∥2≤k-n+1∥μ0∥2+∑i=0nk-n-i∑j=12sdj[μ∥Vj(i)∥2+h(β+γq)i×(δ∥Vj(i-m)∥2+(1-δ)i×∥Vj(i-m-1)∥2)]≤k-n+1∥μ0∥2+∑j=12sdj[∑i=0n-m(k-n-iμ+h(β+γq)δk-n-m-iii+h(β+γq)(1-δ)k-n-m-i-1)ii×∥Vj(i)∥2+h(β+γq)ii×(δ∑i=-m-1k-n-m-i∥Vj(i)∥2ii+(1-δ)∑i=-m-1-1k-n-m-i-1∥Vj(i)∥2)]≤k-n+1∥μ0∥2+∑j=12sdj[∑i=0n-mk-n-m-i-1(k-m+1μ+h(β+γq)δk-ai+h(β+γq)(1-δ)k-m+1)∥Vj(i)∥2ii+k-n-mh(β+γq)∑i=-m-1-1∥Vj(i)∥2].
On the other hand,
(29)μk-m+1+hk-(β+γq)δ+h(β+γq)(1-δ)≤h(β+γq)+hk-(β+γq)δ+h(β+γq)(1-δ)≤0.
Considering dj≥0, j=1,2,…,s and 0<k-<1, we have
(30)limn→∞∥μn∥=0.
Because μn+1=(wn+1,wn)T, so we have
(31)limn→∞∥wn∥=0;
which shows that the method isGAR(l)-stable.Theorem 8.
Assume that a TSRK method is(k,l)-algebraically stable, k<1, and dj>0, j=1,2…,s. Then the method with linear interpolation procedure is weak GAR(l)-stable.Proof.
It follows from Theorem5 that
(32)∥μn+1∥2≤∥μn∥2+∑j=12sdj×[-σ∥Vj(n)∥2+h(β+γq)i×(δ∥Vj(n-m)∥2+(1-δ)∥Vj(n-m-1)∥2)],
where  σ=2l-(2α+β+γq)h. When (α+β+γq)h<l, we have σ>0. Analogous to Theorem 6, we can easily obtain
(33)∥μn+1∥2≤∥μ0∥2+∑j=12sdj[∑i=1n∥Vj(i)∥2(β+γq)τmax-m-1≤i≤-1∥Vj(i)∥2-σ∑i=1n∥Vj(i)∥2],
which shows
(34)limn→∞∑j=12sdj∥Vj(n)∥2=0,limn→∞∑j=12sdj∥Wj(n)∥2=0,(35)limn→∞∑j=12sdj∥V-j(n)∥2=0,limn→∞∑j=12sdj∥W-j(n)∥2=0.
On the other hand,
(36)∥Qi(n)∥≤h∥f(tn+cih,Yi(n),Y-i(n),Y~i(n))-f(tn+cih,Zi(n),Y-i(n),Y~i(n))∥+h∥f(tn+cih,Zi(n),Y-i(n),Y~i(n))-f(tn+cih,Zi(n),Z-i(n),Z~i(n))∥+h∥f(tn+cih,Zi(n),Z-i(n),Y~i(n))-f(tn+cih,Zi(n),Z-i(n),Z~i(n))∥≤h∥f(tn+cih,Yi(n),Y-i(n),Y~i(n))-f(tn+cih,Zi(n),Y-i(n),Y~i(n))∥+hβ∥W-j(n)∥+hγ∥W~j(n)∥≤h∥f(tn+cih,Yi(n),Y-i(n),Y~i(n))-f(tn+cih,Zi(n),Y-i(n),Y~i(n))∥+hβ∥W-j(n)∥+hγq∥W-j(n)∥.
In view of (12) and (35), we have
(37)limn→∞∑j=12sdj∥Qj(n)∥2=0
Considering (9a), (9b), (9c), (22), and (37) with dj>0, we have
(38)limn→∞∥μn∥=0.
Because μn+1=(wn+1,wn)T, so we have
(39)limn→∞∥wn∥=0,
which shows that the method is weakGAR(l)-stable.


## 3.1. Some Concepts

Definition 1.
Letl be a real constant, a two-step Runge-Kutta method with an interpolation procedure is said to be R(l)-stable if there exists a constant C dependent only on the method and l such that (α+β+γq)h≤l and
(10)∥yn-zn∥2≤[2+C(β+γq)]maxt≤0∥φ(t)-ψ(t)∥2,n≥1,
with step size h satisfying hm=τ, where m is a positive integer.
GR
(
l
)-stability is defined by dropping the restriction hm=τ.Definition 2.
Letl be a real constant, a two-step Runge-Kutta method with an interpolation procedure is said to be AR(l)-stable if
(11)limn→∞∥yn-zn∥=0
with step size h satisfying (α+β+γq)h≤l and hm=τ, where m is a positive integer.
GAR
(
l
)-stability is defined by dropping the restriction hm=τ.Definition 3.
Letl be a real constant; a two-step Runge-Kutta method with an interpolation procedure is said to be weak AR(l)-stable if, under the conditions of Definition 2, (11) holds when f further satisfies
(12)∥f(t,u1,v,w)-f(t,u2,v,w)∥≤L∥u1-u2∥p0
with p0 being a positive real number and L being a nonnegative real number.
WeakGAR(l)-stability is defined by dropping the restriction hm=τ.


## 3.2. The Stability of TSRK Methods

Letwn=yn-zn, Wi(n)=Yi(n)-Zi(n), W-i(n)=Y-i(n)-Z-i(n), W~i(n)=Y~i(n)-Z~i(n),  Qi(n)=h[f(tn+cih,Yi(n),Y-i(n),Y~i(n))-f(tn+cih,Zi(n),Z-i(n),Z~i(n))], and i=1,2,…,s.It follows from (5a), (5b), (5c), (9a), (9b), and (9c) that(13a)Wi(n)=wn+∑j=1saijQj(n),i=1,2,…,s,(13b)Wi(n-1)=wn-1+∑j=1saijQj(n-1),i=1,2,⋯,s,(13c)wn+1=(1-θ)wn+θwn-1+∑i=1sbiQi(n)+∑i=1sb~iQi(n-1).Now we will write the s-stage TSRK methods (4a), (4b), and (4c) as a general linear method.LetVi(n)=(Wi(n),Wi(n-1))T be the internal stages, μn+1=(wn+1,wn)T the external vectors, and Pi(n)=(Qi(n),Qi(n-1))T. Then we have a 2(s+1)-stage partitioned general linear method:(14a)Vi(n)=∑j=1sCij11Pi(n)+∑j=1sCij12μn,i=1,2,…,s,(14b)μn+1=∑j=1sCij21Pi(n)+∑j=1sCij22μn,i=1,2,…,s,where
(15)C11=(Cij11)=(A00A),C12=(Cij12)=(e00e),C21=(Cij21)=(bb~00),C22=(Cij22)=(1-θθ00),A=(aij),e=(1,1,…1),b=(b1,b2,…,bs),b~=(b~1,b~2,…,b~s).Definition 4.
Letk,l be real constants; a TSRK method is said to be (k,l)-algebraically stable if there exists a diagonal nonnegative matrix G and D=diag(d1,d2,…,d2s) such that M=(mij) is nonnegative, where(16)
M
(
k
,
l
)
=
(
k
G
-
C
22
T
G
C
22
-
2
l
C
12
T
D
C
12
C
12
T
D
-
C
22
T
G
C
21
-
2
l
C
12
T
D
C
11
D
C
12
-
C
21
T
G
C
22
-
2
l
C
11
T
D
C
12
C
11
T
D
+
D
C
11
-
C
21
T
G
C
21
-
2
l
C
11
T
D
C
11
)
.In this paper, we use the linear interpolation procedure. Letτ=(m-δ)h with integer m≥1 and δ∈[0,1).Define(17a)Y-j(n)=δYj(n-m+1)+(1-δ)Yj(n-m),(17b)Z-j(n)=δZj(n-m+1)+(1-δ)Zj(n-m),j=1,2,…,s,where Yj(i)=φ(ti+cjh) and Zj(i)=ψ(ti+cjh) for i<0. When the step size h satisfies τ=mh, we have
(18)Y-j(n)=Yj(n-m),Z-j(n)=Zj(n-m),j=1,2,…,s.Theorem 5.
Assume that a TSRK method is(k,l)-algebraically stable. Then
(19)∥μn+1∥2≤k∥μn∥2+∑j=12sdj[((2α+β+γq)h-2l)∥Vj(n)∥2+h(β+γq)×(δ∥Vj(n-m)∥2+(1-δ)∥Vj(n-m-1)∥2)].Proof.
It is well known that(20)∥μn+1∥2-k∥μn∥2-2∑j=12sdjRe〈Vj(n),Pj(n)-lVj(n)〉=-∑i=12s+2∑j=12s+2Mij〈ri,rj〉,
where
(21)r1=wn+1,r2=wn,rj=Qj-2(n),j=3,4,…,s+2,rj=Qj-s-2(n-1),j=s+3,s+4,…,2s+2.
By means of (k,l)-algebraical stability of the method, we have
(22)∥μn+1∥2≤k∥μn∥2+2∑j=12sdjRe〈Vj(n),Pj(n)-lVj(n)〉.
It follows from (2) and (6) that
(23)2Re〈Wj(n),Qj(n)〉=2hRe〈Wj(n),f(tn+cjh,Yj(n),Y-j(n),Y~j(n))-f(tn+cjh,Zj(n),Y-j(n),Y~j(n))〉+2hRe〈Wj(n),f(tn+cjh,Zj(n),Y-j(n),Y~j(n))-f(tn+cjh,Zj(n),Z-j(n),Y~j(n))〉+2hRe〈Wj(n),f(tn+cjh,Zj(n),Z-j(n),Y~j(n))-f(tn+cjh,Zj(n),Z-j(n),Z~j(n))〉≤2hα∥Wj(n)∥2+2hβ∥Wj(n)∥·∥W-j(n)∥+2hγ∥Wj(n)∥·∥W~j(n)∥≤2hα∥Wj(n)∥2+2hβ∥Wj(n)∥·∥W-j(n)∥+2hγq∥Wj(n)∥·∥W-j(n)∥≤2hα∥Wj(n)∥2+2h(β+γq)∥Wj(n)∥·∥W-j(n)∥≤2hα∥Wj(n)∥2+2h(β+γq)∥Wj(n)∥·(δ∥Wj(n-m)∥+(1-δ)∥Wj(n-1-m)∥)≤h(2α+β+γq)∥Wj(n)∥2+h(β+γq)(δ∥Wj(n-m)∥2+(1-δ)∥Wj(n-1-m)∥2).
Substitution into (20) gives (19).Theorem 6.
Assume that a TSRK method is(k,l)-algebraically stable and k≤1. Then the method with linear interpolation procedure is GR(l)-stable.Proof.
The inequality(α+β+γq)h≤l and Theorem 5 lead to
(24)∥μn+1∥2≤∥μn∥2+∑j=12sdj[-(β+γq)h∥Vj(n)∥2+h(β+γq)ai×(δ∥Vj(n-m)∥2+(1-δ)∥Vj(n-m-1)∥2)].
By induction, we have
(25)∥μn+1∥2≤∥μ0∥2+h∑j=12sdj(β+γq)(δ∑i=-m-1∥Vj(i)∥2+(1-δ)∑i=-m-1-1∥Vj(i)∥2)≤∥μ0∥2+h∑j=12sdj(β+γq)×[δ(m-1)+(1-δ)m]max-m-1≤i≤-1∥Vj(i)∥2=∥μ0∥2+(β+γq)τ∑j=12sdjmax-m-1≤i≤-1∥Vj(i)∥2≤∥μ0∥2+(β+γq)τ∑j=12sdjmaxt≤0∥φ(t)-ψ(t)∥2.
Because μn+1=(wn+1,wn)T, so we have
(26)∥wn∥2≤∥w0∥2+∥w-1∥2+(β+γq)τ×∑j=12sdjmaxt≤0∥φ(t)-ψ(t)∥2≤2maxt≤0∥φ(t)-ψ(t)∥2+(β+γq)τ×∑j=12sdjmaxt≤0∥φ(t)-ψ(t)∥2≤[2+C(β+γq)τ]maxt≤0∥φ(t)-ψ(t)∥2.
Therefore, it isGR(l)-stable, where C=∑j=12sdj.Theorem 7.
Assume that a TSRK method is(k,l)-algebraically stable and k<1. Then the method with linear interpolation procedure is GAR(l)-stable.Proof.
Letμ=(2α+β+γq)h-2l and k-=max{k,(β+γq)h/(-μ)1/(m+1)}.
Then, when(α+β+γq)h<l, we have μ<-(β+γq)h and 0<k-<1.
The application of Theorem5 yields
(27)∥μn+1∥2≤k-∥μn∥2+∑j=12sdj[μ∥Vj(n)∥2+h(β+γq)×(δ∥Vj(n-m)∥2+(1-δ)∥Vj(n-m-1)∥2)].
By induction, we have
(28)∥μn+1∥2≤k-n+1∥μ0∥2+∑i=0nk-n-i∑j=12sdj[μ∥Vj(i)∥2+h(β+γq)i×(δ∥Vj(i-m)∥2+(1-δ)i×∥Vj(i-m-1)∥2)]≤k-n+1∥μ0∥2+∑j=12sdj[∑i=0n-m(k-n-iμ+h(β+γq)δk-n-m-iii+h(β+γq)(1-δ)k-n-m-i-1)ii×∥Vj(i)∥2+h(β+γq)ii×(δ∑i=-m-1k-n-m-i∥Vj(i)∥2ii+(1-δ)∑i=-m-1-1k-n-m-i-1∥Vj(i)∥2)]≤k-n+1∥μ0∥2+∑j=12sdj[∑i=0n-mk-n-m-i-1(k-m+1μ+h(β+γq)δk-ai+h(β+γq)(1-δ)k-m+1)∥Vj(i)∥2ii+k-n-mh(β+γq)∑i=-m-1-1∥Vj(i)∥2].
On the other hand,
(29)μk-m+1+hk-(β+γq)δ+h(β+γq)(1-δ)≤h(β+γq)+hk-(β+γq)δ+h(β+γq)(1-δ)≤0.
Considering dj≥0, j=1,2,…,s and 0<k-<1, we have
(30)limn→∞∥μn∥=0.
Because μn+1=(wn+1,wn)T, so we have
(31)limn→∞∥wn∥=0;
which shows that the method isGAR(l)-stable.Theorem 8.
Assume that a TSRK method is(k,l)-algebraically stable, k<1, and dj>0, j=1,2…,s. Then the method with linear interpolation procedure is weak GAR(l)-stable.Proof.
It follows from Theorem5 that
(32)∥μn+1∥2≤∥μn∥2+∑j=12sdj×[-σ∥Vj(n)∥2+h(β+γq)i×(δ∥Vj(n-m)∥2+(1-δ)∥Vj(n-m-1)∥2)],
where  σ=2l-(2α+β+γq)h. When (α+β+γq)h<l, we have σ>0. Analogous to Theorem 6, we can easily obtain
(33)∥μn+1∥2≤∥μ0∥2+∑j=12sdj[∑i=1n∥Vj(i)∥2(β+γq)τmax-m-1≤i≤-1∥Vj(i)∥2-σ∑i=1n∥Vj(i)∥2],
which shows
(34)limn→∞∑j=12sdj∥Vj(n)∥2=0,limn→∞∑j=12sdj∥Wj(n)∥2=0,(35)limn→∞∑j=12sdj∥V-j(n)∥2=0,limn→∞∑j=12sdj∥W-j(n)∥2=0.
On the other hand,
(36)∥Qi(n)∥≤h∥f(tn+cih,Yi(n),Y-i(n),Y~i(n))-f(tn+cih,Zi(n),Y-i(n),Y~i(n))∥+h∥f(tn+cih,Zi(n),Y-i(n),Y~i(n))-f(tn+cih,Zi(n),Z-i(n),Z~i(n))∥+h∥f(tn+cih,Zi(n),Z-i(n),Y~i(n))-f(tn+cih,Zi(n),Z-i(n),Z~i(n))∥≤h∥f(tn+cih,Yi(n),Y-i(n),Y~i(n))-f(tn+cih,Zi(n),Y-i(n),Y~i(n))∥+hβ∥W-j(n)∥+hγ∥W~j(n)∥≤h∥f(tn+cih,Yi(n),Y-i(n),Y~i(n))-f(tn+cih,Zi(n),Y-i(n),Y~i(n))∥+hβ∥W-j(n)∥+hγq∥W-j(n)∥.
In view of (12) and (35), we have
(37)limn→∞∑j=12sdj∥Qj(n)∥2=0
Considering (9a), (9b), (9c), (22), and (37) with dj>0, we have
(38)limn→∞∥μn∥=0.
Because μn+1=(wn+1,wn)T, so we have
(39)limn→∞∥wn∥=0,
which shows that the method is weakGAR(l)-stable.


## 4. The Convergence of TSRK Method for NDDEs

### 4.1. Some Concepts

In order to study the convergence of the method, we define(40)Y(n)=(Y1(n)!T,…,Ys(n)!T,Y1(n-1)!T,…,Ys(n-1)!T)T∈C2sN,Y-(n)=(Y-1(n)!T,…,Y-s(n)!T,Y-1(n-1)!T,…,Y-s(n-1)!T)T∈C2sN,Y~(n)=(Y~1(n)!T,…,Y~s(n)!T,Y~1(n-1)!T,…,Y~s(n-1)!T)T∈C2sN,ζ(n+1)=(y(n+1)!T,y(n)!T)T∈C2N,F(tn,Y(n),Y-(n),Y~(n))=(f(tn+c1h,Y1(n),Y-1(n),Y~1(n))T,…,f(tn+csh,Ys(n),Y-s(n),Y~s(n))T,f(tn-1+c1h,Y1(n-1),Y-1(n-1),Y~1(n-1))T,…,f(tn-1+csh,Ys(n-1),Y-s(n-1),Y~s(n-1))T)T∈C2sNThus, process (5a), (5b), and (5c) can be written in the more compact form(41a)Y(n)=hC11F(tn,Y(n),Y-(n),Y~(n))+C12ζ(n),(41b)ζ(n+1)=hC21F(tn,Y(n),Y-(n),Y~(n))+C22ζ(n).Definition 9.
Method (4a), (4b), and (4c) with an interpolation procedure is said to be D-convergent of order p if the global error satisfies a bound of the form
(42)∥H(tn)-y(n)∥≤ρ1(tn)(∥H(t0)-y(0)∥+maxk≤0∥Y(tk)-y(k)∥+hp),n≥1,h∈(0,h0],
where H(t) is defined by
(43)H(t)=(y(t+h),y(t))∈C2rN,Y(t)=(y(t+c1h),…,y(t+csh),iiy(t-h+c1h),…,y(t-h+csh))∈C2sN.ρ1(t) and h0 depend on Mi, α, β, γ, and τ.Definition 10.
TSRK method (4a), (4b), and (4c) is said to be algebraically stable if there exist a real symmetric, positive definite 2r×2r matrix G and a nonnegative diagonal 2s×2s matrix D such that the matrix
(44)M=[G-C22TGC22C12TD-C22TGC21DC12-C21TGC22DC11+C11TD-C21TGC21]
is nonnegative definite.Definition 11.
TSRK method (4a), (4b), and (4c) is said to be diagonally stable if there exist an 2s×2s diagonal matrix Q>0 such that the matrix QC11+C11TQ is positive definite.Remark 12.
The concepts of algebraic stability and diagonal stability of TSRK method are the generalizations of corresponding concepts of Runge-Kutta methods. Although it is difficult to examine these conditions, many results have been found; in particular, there exist algebraically stable and diagonally stable multistep formulas of arbitrarily high order (cf. [31]).Definition 13.
TSRK method (4a), (4b), and (4c) is said to have generalized stage order P if P is the largest integer which possesses the following properties.
For any given problem (1) and h∈[0,h-0], there exists an abstract function Hh(t),
(45)Hh(t)=(H1h(t),H2h(t))∈C2N,
such that
(46)∥H(t)-Hh(t)∥≤p1hp,∥Δh(t)∥≤p1hp+1,∥δh(t)∥≤p1hp+1,
where the maximum step size h-0>0 and the constant p1 depend only on the method and the bounds Mi, Δh(t), and δh(t); they are defined by;(47a)Y(t)=hC11Y′(t)+C12Hh(t)+Δh(t),(47b)Hh(t+h)=hC21Y′(t)+C22Hh(t)+δh(t).
The functionY′(t) is defined by(48a)Y(t)=(y′(t+c1h),…,y′(t+csh),iiy′(t-h+c1h),…,y′(t-h+csh))∈C2sN.
In particular whenH(t)=Hh(t), generalized stage order is called stage order.


### 4.2. D-Convergence and Proofs

In this section, we focus on the error analysis of TSRK method for (1). For the sake of simplicity, we always assume that all constants hi, γi, di, and Li are dependent only on the method, some of the bounds Mi, and the parameters α, β, γ, and τ.First, we give a preliminary result which will later be used several times. For anyY, Y-, Y~, Z, Z-, Z~∈C2sN for Y(n), Y-(n), Y~(n), Z(n), Z-(n), Z~(n), and ζ, ϖ∈C2N for ζ(n+1), ϖ(n+1) where ϖ(n+1)=(z(n+1),z(n)). Define Δ~ and δ~ by(49a)Δ~=Y-Z-hC11(F(t,Y,Y-,Y~)-F(t,Z,Z-,Z~)),(49b)δ~=ζ-ϖ-hC21(F(t,Y,Y-,Y~)-F(t,Z,Z-,Z~)).Theorem 14.
Suppose that method (4a), (4b), and (4c) is diagonally stable. Then there exist constants h1, p2, and p3 such that(50a)∥Y-Z∥≤p2(∥Δ~∥+h∥Y--Z-∥+h∥Y~-Z~∥),h∈(0,h1],(50b)∥ζ-ϖ∥≤p3(∥Δ~∥+∥δ~∥+h∥Y--Z-∥+h∥Y~-Z~∥),h∈(0,h1].Proof.
Since the method (4a), (4b), and (4c) is diagonally stable, there exists a positive definite diagonal matrix Q such that the matrix E=QC11+C11TQ is positive definite. Therefore, the matrix C11 is obviously nonsingular, and there exists an l>0 which depends only on the method such that the matrix
(51)El=C11-TEC11-2lQ
is also positive definite.
Define(52a)W=Y-Z,W-=Y--Z-,W~=Y~-Z~,(52b)K1=h(F(t,Y,Y-,Y~)-F(t,Z,Y-,Y~)),(52c)K2=h(F(t,Z,Y-,Y~)-F(t,Z,Z-,Y~)),(52d)K3=h(F(t,Z,Z-,Y~)-F(t,Z,Z-,Z~)).
Then(53a)Δ~=W-C11(K1+K2+K3),(53b)δ~=ζ-ϖ-C21(K1+K2+K3).
Using (2), (6), (53a), and (53b), we have, for hα≤l,
(54)0≤〈W,2lQW〉+2Re〈K1,-QW〉=-〈W,ElW〉+2Re〈W,Q(C11-1-K1)〉≤-λl∥W∥2+2Re〈W,Q(C11-1Δ~+K2+K3)〉≤-λl∥W∥2+2∥QC11-1∥∥W∥∥Δ~∥+2hβ∥Q∥∥W∥∥W-∥+2hγ∥Q∥∥W∥∥W~∥,
where λl is the minimum eigenvalue of El. Therefore,
(55)∥Y-Z∥≤p1(∥Δ~∥+h∥Y--Z-∥+h∥Y~-Z~∥),h∈(0,h1],
where(56a)p1=2λlmax(∥QC11-1∥,β∥Q∥,γ∥Q∥),(56b)h1={1α≤0,min(1,lα)α>0.From (50a), (56a), and (56b), it follows that
(57)∥ζ-ϖ∥=∥δ~+C21C11-1(W-Δ~)∥≤∥δ~∥+p1h∥C21C11-1∥(∥W-∥+∥W~∥)+(1+p1)∥C21C11-1∥·∥Δ~∥≤p3(∥Δ~∥+∥δ~∥+h∥Y--Z-∥+h∥Y~-Z~∥),h∈(0,h1],
where p3=max{1,(1+p2)∥C21C11-1∥}, which completes the proof of Theorem 14.Consider the compact form of (9a), (9b), and (9c):(58a)Z(n)=hC11F(tn,Z(n),Z-(n),Z~(n))+C12ϖ(n),(58b)ϖ(n+1)=hC21F(tn,Z(n),Z-(n),Z~(n))+C22ϖ(n),where(59)Z(n)=(Z1(n),…,Zs(n),Z1(n-1),…,Zs(n-1))∈C2sN,Z-(n)=(Z-1(n),…,Z-s(n),Z-1(n-1),…,Z-s(n-1))∈C2sN,Z~(n)=(Z~1(n),…,Z~s(n),Z~1(n-1),…,Z~s(n-1))∈C2sN,ϖ(n+1)=(z(n+1),z(n))∈C2N.Theorem 15.
Suppose that the method (4a), (4b), and (4c) is algebraically stable for the matrices G and D. Then for (41a), (41b), (58a), and (58b) we have
(60)∥ζ(n+1)-ϖ(n+1)∥G2≤∥ζ(n)-ϖ(n)∥G2+p4h(∥Y(n)-Z(n)∥2+∥Y(n-m)-Z(n-m)∥2+∥Y~(n)-Z~(n)∥2),
where p4=∥D∥·max{(2α+β+γ),β,γ}, ∥·∥ is a norm on C2sNdefined by
(61)∥U∥G=〈U,GU〉1/2=(∑i,j=12sgij〈ui,uj〉)1/2,U=(u1,u2,…,u2s)∈C2sN,ui∈CN.Proof.
Defineu(n)=ζ(n)-ϖ(n). We get from (52a), (52b), (52c), and (52d) that
(62)K(n)=h(F(tn,Y(n),Y-(n),Y~(n))-F(tn,Z(n),Z-(n),Z~(n))),W(n)=C11K(n)+C12u(n),u(n+1)=C21K(n)+C22u(n).
With algebraic stability, the matrix
(63)M=[G-C22TGC22C12TD-C22TGC21DC12-C21TGC22DC11+C11TD-C21TGC21]
is nonnegative definite. As in [32], we have
(64)〈u(n+1),Gu(n+1)〉-〈u(n),Gu(n)〉-2Re〈W(n),DK(n)〉=-〈〈u(n),K(n)〉,M〈u(n),K(n)〉〉≤0.
Using (2), we further obtain
(65)〈u(n+1),Gu(n+1)〉≤〈u(n),Gu(n)〉+2Re〈W(n),DK(n)〉≤〈u(n),Gu(n)〉+2Re〈W(n),DK1(n)〉+2Re〈W(n),DK2(n)+DK3(n)〉≤〈u(n),Gu(n)〉+2hα〈W(n),DW(n)〉+2hβ∥D1/2W(n)∥·∥D1/2W-(n)∥+2hγ∥D1/2W(n)∥·∥D1/2W~(n)∥≤〈u(n),Gu(n)〉+2hα〈W(n),DW(n)〉+hβ(∥D∥·∥W(n)∥2+∥D∥·∥W-(n)∥2)+hγ(∥D∥·∥W(n)∥2+∥D∥·∥W~(n)∥2)≤〈u(n),Gu(n)〉+(2hα+hβ+hγ)∥D∥·∥W(n)∥2+∥D∥(hβ∥W(n-m)∥2+hγ∥W~(n)∥2)
which gives (60). The proof is completed.In the following, we assume that the method (4a), (4b), and (4c) has generalized stage order P; that is, there exists a function Hh(t) such that (46) holds. For any n>0, we define Y^(n) and y^(n+1) by(66a)Y^(n)hC11F(tn,Y^(n),Y(tn-τ),Y~^(n))+C12Hh(tn),(66b)y^(n+1)hC21F(tn,Y^(n),Y(tn-τ),Y~^(n))+C22Hh(tn),where(67)Y~^(n)=(f(tn+c1h-τ,Y-1(n),Y-1(n+1),H1h(tn+c1h-τ)),iif(tn+c2h-τ,Y-2(n),Y-2(n+1),H1h(tn+c2h-τ)),…,iif(tn+csh-τ,Y-s(n),Y-s(n+1),H1h(tn+csh-τ)),iif(tn+c1h-h-τ,Y-1(n-1),Y-1(n),iiiH1h(tn+c1h-h-τ)Y-1(n-1)),…,iif(tn+csh-h-τ,Y-s(n-1),Y-s(n),H1h(tn+csh-h-τ))).Theorem 16.
Suppose that the method (4a), (4b), and (4c) is diagonally stable and its generalized stage order is P. Then there exist constants p5 and h2 such that
(68)∑k=1n∥Y-(k)-Y(tk-τ)∥2≤p4(∑k=1n∥ζ(k)-Hh(tk)∥2+nh(p+1)+nh(v+μ+1)+∑k=-m+10∥Y(k)-Y(tk)∥2),h∈(0,h2].Proof.
It follows from (6) that
(69)∥Y-j(k)-y(tk+cjh-τ)∥≤∥∑i=-μvL~i(δ)(Yj(k-m+i)-y(tk-m+i+cjh))∥+∥∑i=-μvL~i(δ)y(tk-m+i+cjh)-y(tk-m+cjh+δh)∥.
From the remainder estimation of Lagrange interpolation formula, we have
(70)∥∑i=-μvL~i(δ)y(tk-m+i+cjh)-y(tk-m+cjh+δh)∥≤Mμ+v+1hμ+v+1.
Using Cauchy inequality, we further obtain
(71)∥Y-(k)-Y(tk-τ)∥2≤2s(μ+v+2)(L~02∑i=-μv∥Y(k-m+i)-Y(tk-m+i)∥2+Mμ+v+1h(v+μ+1)∑∑i=-μv∥Y(k-m+i)-Y(tk-m+i)∥2),
where L~02=max-μ≤i≤vSupx∈[0,1]|L~i(x)|. Hence, there exists a constant d1 such that
(72)∑k=1n∥Y-(k)-Y(tk-τ)∥2≤d1(∑k=-μ-m+1n∥Y(k)-Y(tk)∥2+nh2(μ+v+1)).
On the other hand, a combination of (41a) and (47a) leads to
(73)Y(k)-Y(tk)=hC11(F(tk,Y(k),Y-(k),Y~(k))-F(tk,Y(tk),Y(tk-τ),Y~(tk-τ))(tk,Y(k),Y-(k),Y~(k)))+C12(y(k-1)-Hh(tk-h))-Δh(tk).
It follows from Theorem 14 that
(74)∥Y(k)-Y(tk)∥=p2{∥C12(y(k-1)-Hh(tk-h)-Δh(tk))∥ii+h∥Y-(k)-Y(tk-τ)∥+h∥Y~(k)-Y~(tk-τ)∥(y(k-1))-Hh(tk-h)-Δh(tk))},
which on substitution into (72) gives
(75)∑k=1n∥Y-(k)-Y(tk-τ)∥2≤d1(p12h2∑k=1n∥Y-(k)-Y(tk-τ)∥2+p12h2∑k=1n∥Y~(k)-Y~(tk-τ)∥2+p12∥C12∥2∑k=1n∥y(k-1)-Hh(tk-1)∥2+p12∥C12∥2∥Δh(tk)∥2+nh2(μ+v+1)+∑k=-μ-m+10∥Y(k)-Y(tk)∥2).
Therefore, there existp4 and h2 such that (68) holds. The proof of Theorem 16 is completed.Theorem 17.
Suppose that method (4a), (4b), and (4c) is algebraically stable and diagonally stable and its generalized stage order is p. Then the method with interpolation procedure (6) is D-convergent of order at least min{p,μ+v+1}.Proof.
In view of (41a), (41b), (66a), and (66b), it follows from Theorem 15 that
(76)∥y(n)-y^(n)∥G2≤p3h(∥Y(n)-Y^(n)∥2+∥Y-(n)-Y(tn-τ)∥2+∥Y~(n)-Y~(tn-τ)∥2)+∥y(n-1)-Hh(tn-1)∥G2.
Using Theorem 14, we have
(77)∥Y(n)-Y^(n)∥≤p1{∥C12∥∥y(n-1)-Hh(tn-1)∥i+h∥Y-(n)-Y(tn-τ)∥+h∥Y~(n)-Y~(tn-τ)∥},h∈(0,h1],
which on substitution into (76) gives
(78)∥y(n)-y^(n)∥G2≤p3h(1+2p12h2)(1+2sr2q2)×∥Y-(n)-Y(tn-τ)∥2+(1+2hp3p12∥C12∥2λ2)×∥y(n-1)-Hh(tn-1)∥G2,h∈(0,h1],
where λ2 is the minimum characteristic value of G. On the other hand,
(79)∥y(n)-Hh(tn)∥G2≤∥y(n)-y^(n)∥G2+∥y^(n)-Hh(tn)∥G2+2∥y(n)-y^(n)∥G∥y^(n)-Hh(tn)∥G≤(1+h)∥y(n)-y^(n)∥G2+(1+1h)∥y^(n)-Hh(tn)∥G2.
In view of (47a), (47b), (48a), (66a), and (66b), the application of Theorem 14 leads to
(80)∥y^(n)-Hh(tn)∥≤p2(∥Δh(tn)∥+∥δh(tn)∥),h∈(0,h-0],
which gives
(81)∥y^(n)-Hh(tn)∥G2≤λ1p12(∥Δh(tn)∥+∥δh(tn)∥)2,h∈(0,h-0],
where λ1 denotes the maximum eigenvalue of the matrix G. A combination of (46) and (76)–(81) leads to
(82)∥y(n)-Hh(tn)∥G2≤(1+hγ1)∥y(n-1)-Hh(tn-1)∥G2+γ2h∥Y-(n)-Y(tn-τ)∥2+γ3h2p+1,h∈(0,h3],
where
(83)h3=min{h-0,h1}≤1,γ1=1+4p3p12∥C12∥2λ2,γ2=2p3(1+2p12)(1+2sr2q2),(84)γ3=8λ1p02p22.
Therefore,
(85)∥y(n)-Hh(tn)∥G2≤∥y(0)-Hh(t0)∥G2+h∑i=1n(γ1∥y(i-1)-Hh(ti-1)∥G2+γ2∥Y-(i)-Y(ti-τ)∥2)+γ3tnh2p,h∈(0,h3].
Considering Theorem 16, we further obtain
(86)∥y(n)-Hh(tn)∥G2≤∥y(0)-Hh(t0)∥G2+h(γ1+γ2p4)×∑i=0n-1∥y(i)-Hh(ti)∥G2+γ3tnh2p+γ2p4tn(h2(p+1)+h2(μ+v+1))+p4γ2h×∑k=-μ-m+10∥Y(k)-Y(tk)∥2,h∈(0,h0],
where h0=min{h2,h3}≤1. Using discrete Bellman inequality, we have
(87)∥y(n)-Hh(tn)∥G2≤(∥y(0)-Hh(t0)∥G2+p4γ2hoo×∑k=-μ-m+10∥Y(k)-Y(tk)∥2oo+(γ3+γ2p4)tnh2p+γ2p4tnh2(μ+v+1)∥y(0)-Hh(t0)∥G2)×exp((γ1+γ2p4)tn),h∈(0,h0].
Considering ∥H(t)-Hh(tn)∥≤p0hp, we obtain
(88)∥y(n)-H(tn)∥≤∥y(n)-Hh(tn)∥+∥Hh(tn)-H(tn)∥≤∥y(n)-Hh(tn)∥+p0hp.
Considering (87) and (88), we can easily conclude that method (4a), (4b), and (4c) with interpolation procedure (6) is D-convergent of order at least min{p,μ+v+1}. The proof is completed.


## 4.1. Some Concepts

In order to study the convergence of the method, we define(40)Y(n)=(Y1(n)!T,…,Ys(n)!T,Y1(n-1)!T,…,Ys(n-1)!T)T∈C2sN,Y-(n)=(Y-1(n)!T,…,Y-s(n)!T,Y-1(n-1)!T,…,Y-s(n-1)!T)T∈C2sN,Y~(n)=(Y~1(n)!T,…,Y~s(n)!T,Y~1(n-1)!T,…,Y~s(n-1)!T)T∈C2sN,ζ(n+1)=(y(n+1)!T,y(n)!T)T∈C2N,F(tn,Y(n),Y-(n),Y~(n))=(f(tn+c1h,Y1(n),Y-1(n),Y~1(n))T,…,f(tn+csh,Ys(n),Y-s(n),Y~s(n))T,f(tn-1+c1h,Y1(n-1),Y-1(n-1),Y~1(n-1))T,…,f(tn-1+csh,Ys(n-1),Y-s(n-1),Y~s(n-1))T)T∈C2sNThus, process (5a), (5b), and (5c) can be written in the more compact form(41a)Y(n)=hC11F(tn,Y(n),Y-(n),Y~(n))+C12ζ(n),(41b)ζ(n+1)=hC21F(tn,Y(n),Y-(n),Y~(n))+C22ζ(n).Definition 9.
Method (4a), (4b), and (4c) with an interpolation procedure is said to be D-convergent of order p if the global error satisfies a bound of the form
(42)∥H(tn)-y(n)∥≤ρ1(tn)(∥H(t0)-y(0)∥+maxk≤0∥Y(tk)-y(k)∥+hp),n≥1,h∈(0,h0],
where H(t) is defined by
(43)H(t)=(y(t+h),y(t))∈C2rN,Y(t)=(y(t+c1h),…,y(t+csh),iiy(t-h+c1h),…,y(t-h+csh))∈C2sN.ρ1(t) and h0 depend on Mi, α, β, γ, and τ.Definition 10.
TSRK method (4a), (4b), and (4c) is said to be algebraically stable if there exist a real symmetric, positive definite 2r×2r matrix G and a nonnegative diagonal 2s×2s matrix D such that the matrix
(44)M=[G-C22TGC22C12TD-C22TGC21DC12-C21TGC22DC11+C11TD-C21TGC21]
is nonnegative definite.Definition 11.
TSRK method (4a), (4b), and (4c) is said to be diagonally stable if there exist an 2s×2s diagonal matrix Q>0 such that the matrix QC11+C11TQ is positive definite.Remark 12.
The concepts of algebraic stability and diagonal stability of TSRK method are the generalizations of corresponding concepts of Runge-Kutta methods. Although it is difficult to examine these conditions, many results have been found; in particular, there exist algebraically stable and diagonally stable multistep formulas of arbitrarily high order (cf. [31]).Definition 13.
TSRK method (4a), (4b), and (4c) is said to have generalized stage order P if P is the largest integer which possesses the following properties.
For any given problem (1) and h∈[0,h-0], there exists an abstract function Hh(t),
(45)Hh(t)=(H1h(t),H2h(t))∈C2N,
such that
(46)∥H(t)-Hh(t)∥≤p1hp,∥Δh(t)∥≤p1hp+1,∥δh(t)∥≤p1hp+1,
where the maximum step size h-0>0 and the constant p1 depend only on the method and the bounds Mi, Δh(t), and δh(t); they are defined by;(47a)Y(t)=hC11Y′(t)+C12Hh(t)+Δh(t),(47b)Hh(t+h)=hC21Y′(t)+C22Hh(t)+δh(t).
The functionY′(t) is defined by(48a)Y(t)=(y′(t+c1h),…,y′(t+csh),iiy′(t-h+c1h),…,y′(t-h+csh))∈C2sN.
In particular whenH(t)=Hh(t), generalized stage order is called stage order.


## 4.2. D-Convergence and Proofs

In this section, we focus on the error analysis of TSRK method for (1). For the sake of simplicity, we always assume that all constants hi, γi, di, and Li are dependent only on the method, some of the bounds Mi, and the parameters α, β, γ, and τ.First, we give a preliminary result which will later be used several times. For anyY, Y-, Y~, Z, Z-, Z~∈C2sN for Y(n), Y-(n), Y~(n), Z(n), Z-(n), Z~(n), and ζ, ϖ∈C2N for ζ(n+1), ϖ(n+1) where ϖ(n+1)=(z(n+1),z(n)). Define Δ~ and δ~ by(49a)Δ~=Y-Z-hC11(F(t,Y,Y-,Y~)-F(t,Z,Z-,Z~)),(49b)δ~=ζ-ϖ-hC21(F(t,Y,Y-,Y~)-F(t,Z,Z-,Z~)).Theorem 14.
Suppose that method (4a), (4b), and (4c) is diagonally stable. Then there exist constants h1, p2, and p3 such that(50a)∥Y-Z∥≤p2(∥Δ~∥+h∥Y--Z-∥+h∥Y~-Z~∥),h∈(0,h1],(50b)∥ζ-ϖ∥≤p3(∥Δ~∥+∥δ~∥+h∥Y--Z-∥+h∥Y~-Z~∥),h∈(0,h1].Proof.
Since the method (4a), (4b), and (4c) is diagonally stable, there exists a positive definite diagonal matrix Q such that the matrix E=QC11+C11TQ is positive definite. Therefore, the matrix C11 is obviously nonsingular, and there exists an l>0 which depends only on the method such that the matrix
(51)El=C11-TEC11-2lQ
is also positive definite.
Define(52a)W=Y-Z,W-=Y--Z-,W~=Y~-Z~,(52b)K1=h(F(t,Y,Y-,Y~)-F(t,Z,Y-,Y~)),(52c)K2=h(F(t,Z,Y-,Y~)-F(t,Z,Z-,Y~)),(52d)K3=h(F(t,Z,Z-,Y~)-F(t,Z,Z-,Z~)).
Then(53a)Δ~=W-C11(K1+K2+K3),(53b)δ~=ζ-ϖ-C21(K1+K2+K3).
Using (2), (6), (53a), and (53b), we have, for hα≤l,
(54)0≤〈W,2lQW〉+2Re〈K1,-QW〉=-〈W,ElW〉+2Re〈W,Q(C11-1-K1)〉≤-λl∥W∥2+2Re〈W,Q(C11-1Δ~+K2+K3)〉≤-λl∥W∥2+2∥QC11-1∥∥W∥∥Δ~∥+2hβ∥Q∥∥W∥∥W-∥+2hγ∥Q∥∥W∥∥W~∥,
where λl is the minimum eigenvalue of El. Therefore,
(55)∥Y-Z∥≤p1(∥Δ~∥+h∥Y--Z-∥+h∥Y~-Z~∥),h∈(0,h1],
where(56a)p1=2λlmax(∥QC11-1∥,β∥Q∥,γ∥Q∥),(56b)h1={1α≤0,min(1,lα)α>0.From (50a), (56a), and (56b), it follows that
(57)∥ζ-ϖ∥=∥δ~+C21C11-1(W-Δ~)∥≤∥δ~∥+p1h∥C21C11-1∥(∥W-∥+∥W~∥)+(1+p1)∥C21C11-1∥·∥Δ~∥≤p3(∥Δ~∥+∥δ~∥+h∥Y--Z-∥+h∥Y~-Z~∥),h∈(0,h1],
where p3=max{1,(1+p2)∥C21C11-1∥}, which completes the proof of Theorem 14.Consider the compact form of (9a), (9b), and (9c):(58a)Z(n)=hC11F(tn,Z(n),Z-(n),Z~(n))+C12ϖ(n),(58b)ϖ(n+1)=hC21F(tn,Z(n),Z-(n),Z~(n))+C22ϖ(n),where(59)Z(n)=(Z1(n),…,Zs(n),Z1(n-1),…,Zs(n-1))∈C2sN,Z-(n)=(Z-1(n),…,Z-s(n),Z-1(n-1),…,Z-s(n-1))∈C2sN,Z~(n)=(Z~1(n),…,Z~s(n),Z~1(n-1),…,Z~s(n-1))∈C2sN,ϖ(n+1)=(z(n+1),z(n))∈C2N.Theorem 15.
Suppose that the method (4a), (4b), and (4c) is algebraically stable for the matrices G and D. Then for (41a), (41b), (58a), and (58b) we have
(60)∥ζ(n+1)-ϖ(n+1)∥G2≤∥ζ(n)-ϖ(n)∥G2+p4h(∥Y(n)-Z(n)∥2+∥Y(n-m)-Z(n-m)∥2+∥Y~(n)-Z~(n)∥2),
where p4=∥D∥·max{(2α+β+γ),β,γ}, ∥·∥ is a norm on C2sNdefined by
(61)∥U∥G=〈U,GU〉1/2=(∑i,j=12sgij〈ui,uj〉)1/2,U=(u1,u2,…,u2s)∈C2sN,ui∈CN.Proof.
Defineu(n)=ζ(n)-ϖ(n). We get from (52a), (52b), (52c), and (52d) that
(62)K(n)=h(F(tn,Y(n),Y-(n),Y~(n))-F(tn,Z(n),Z-(n),Z~(n))),W(n)=C11K(n)+C12u(n),u(n+1)=C21K(n)+C22u(n).
With algebraic stability, the matrix
(63)M=[G-C22TGC22C12TD-C22TGC21DC12-C21TGC22DC11+C11TD-C21TGC21]
is nonnegative definite. As in [32], we have
(64)〈u(n+1),Gu(n+1)〉-〈u(n),Gu(n)〉-2Re〈W(n),DK(n)〉=-〈〈u(n),K(n)〉,M〈u(n),K(n)〉〉≤0.
Using (2), we further obtain
(65)〈u(n+1),Gu(n+1)〉≤〈u(n),Gu(n)〉+2Re〈W(n),DK(n)〉≤〈u(n),Gu(n)〉+2Re〈W(n),DK1(n)〉+2Re〈W(n),DK2(n)+DK3(n)〉≤〈u(n),Gu(n)〉+2hα〈W(n),DW(n)〉+2hβ∥D1/2W(n)∥·∥D1/2W-(n)∥+2hγ∥D1/2W(n)∥·∥D1/2W~(n)∥≤〈u(n),Gu(n)〉+2hα〈W(n),DW(n)〉+hβ(∥D∥·∥W(n)∥2+∥D∥·∥W-(n)∥2)+hγ(∥D∥·∥W(n)∥2+∥D∥·∥W~(n)∥2)≤〈u(n),Gu(n)〉+(2hα+hβ+hγ)∥D∥·∥W(n)∥2+∥D∥(hβ∥W(n-m)∥2+hγ∥W~(n)∥2)
which gives (60). The proof is completed.In the following, we assume that the method (4a), (4b), and (4c) has generalized stage order P; that is, there exists a function Hh(t) such that (46) holds. For any n>0, we define Y^(n) and y^(n+1) by(66a)Y^(n)hC11F(tn,Y^(n),Y(tn-τ),Y~^(n))+C12Hh(tn),(66b)y^(n+1)hC21F(tn,Y^(n),Y(tn-τ),Y~^(n))+C22Hh(tn),where(67)Y~^(n)=(f(tn+c1h-τ,Y-1(n),Y-1(n+1),H1h(tn+c1h-τ)),iif(tn+c2h-τ,Y-2(n),Y-2(n+1),H1h(tn+c2h-τ)),…,iif(tn+csh-τ,Y-s(n),Y-s(n+1),H1h(tn+csh-τ)),iif(tn+c1h-h-τ,Y-1(n-1),Y-1(n),iiiH1h(tn+c1h-h-τ)Y-1(n-1)),…,iif(tn+csh-h-τ,Y-s(n-1),Y-s(n),H1h(tn+csh-h-τ))).Theorem 16.
Suppose that the method (4a), (4b), and (4c) is diagonally stable and its generalized stage order is P. Then there exist constants p5 and h2 such that
(68)∑k=1n∥Y-(k)-Y(tk-τ)∥2≤p4(∑k=1n∥ζ(k)-Hh(tk)∥2+nh(p+1)+nh(v+μ+1)+∑k=-m+10∥Y(k)-Y(tk)∥2),h∈(0,h2].Proof.
It follows from (6) that
(69)∥Y-j(k)-y(tk+cjh-τ)∥≤∥∑i=-μvL~i(δ)(Yj(k-m+i)-y(tk-m+i+cjh))∥+∥∑i=-μvL~i(δ)y(tk-m+i+cjh)-y(tk-m+cjh+δh)∥.
From the remainder estimation of Lagrange interpolation formula, we have
(70)∥∑i=-μvL~i(δ)y(tk-m+i+cjh)-y(tk-m+cjh+δh)∥≤Mμ+v+1hμ+v+1.
Using Cauchy inequality, we further obtain
(71)∥Y-(k)-Y(tk-τ)∥2≤2s(μ+v+2)(L~02∑i=-μv∥Y(k-m+i)-Y(tk-m+i)∥2+Mμ+v+1h(v+μ+1)∑∑i=-μv∥Y(k-m+i)-Y(tk-m+i)∥2),
where L~02=max-μ≤i≤vSupx∈[0,1]|L~i(x)|. Hence, there exists a constant d1 such that
(72)∑k=1n∥Y-(k)-Y(tk-τ)∥2≤d1(∑k=-μ-m+1n∥Y(k)-Y(tk)∥2+nh2(μ+v+1)).
On the other hand, a combination of (41a) and (47a) leads to
(73)Y(k)-Y(tk)=hC11(F(tk,Y(k),Y-(k),Y~(k))-F(tk,Y(tk),Y(tk-τ),Y~(tk-τ))(tk,Y(k),Y-(k),Y~(k)))+C12(y(k-1)-Hh(tk-h))-Δh(tk).
It follows from Theorem 14 that
(74)∥Y(k)-Y(tk)∥=p2{∥C12(y(k-1)-Hh(tk-h)-Δh(tk))∥ii+h∥Y-(k)-Y(tk-τ)∥+h∥Y~(k)-Y~(tk-τ)∥(y(k-1))-Hh(tk-h)-Δh(tk))},
which on substitution into (72) gives
(75)∑k=1n∥Y-(k)-Y(tk-τ)∥2≤d1(p12h2∑k=1n∥Y-(k)-Y(tk-τ)∥2+p12h2∑k=1n∥Y~(k)-Y~(tk-τ)∥2+p12∥C12∥2∑k=1n∥y(k-1)-Hh(tk-1)∥2+p12∥C12∥2∥Δh(tk)∥2+nh2(μ+v+1)+∑k=-μ-m+10∥Y(k)-Y(tk)∥2).
Therefore, there existp4 and h2 such that (68) holds. The proof of Theorem 16 is completed.Theorem 17.
Suppose that method (4a), (4b), and (4c) is algebraically stable and diagonally stable and its generalized stage order is p. Then the method with interpolation procedure (6) is D-convergent of order at least min{p,μ+v+1}.Proof.
In view of (41a), (41b), (66a), and (66b), it follows from Theorem 15 that
(76)∥y(n)-y^(n)∥G2≤p3h(∥Y(n)-Y^(n)∥2+∥Y-(n)-Y(tn-τ)∥2+∥Y~(n)-Y~(tn-τ)∥2)+∥y(n-1)-Hh(tn-1)∥G2.
Using Theorem 14, we have
(77)∥Y(n)-Y^(n)∥≤p1{∥C12∥∥y(n-1)-Hh(tn-1)∥i+h∥Y-(n)-Y(tn-τ)∥+h∥Y~(n)-Y~(tn-τ)∥},h∈(0,h1],
which on substitution into (76) gives
(78)∥y(n)-y^(n)∥G2≤p3h(1+2p12h2)(1+2sr2q2)×∥Y-(n)-Y(tn-τ)∥2+(1+2hp3p12∥C12∥2λ2)×∥y(n-1)-Hh(tn-1)∥G2,h∈(0,h1],
where λ2 is the minimum characteristic value of G. On the other hand,
(79)∥y(n)-Hh(tn)∥G2≤∥y(n)-y^(n)∥G2+∥y^(n)-Hh(tn)∥G2+2∥y(n)-y^(n)∥G∥y^(n)-Hh(tn)∥G≤(1+h)∥y(n)-y^(n)∥G2+(1+1h)∥y^(n)-Hh(tn)∥G2.
In view of (47a), (47b), (48a), (66a), and (66b), the application of Theorem 14 leads to
(80)∥y^(n)-Hh(tn)∥≤p2(∥Δh(tn)∥+∥δh(tn)∥),h∈(0,h-0],
which gives
(81)∥y^(n)-Hh(tn)∥G2≤λ1p12(∥Δh(tn)∥+∥δh(tn)∥)2,h∈(0,h-0],
where λ1 denotes the maximum eigenvalue of the matrix G. A combination of (46) and (76)–(81) leads to
(82)∥y(n)-Hh(tn)∥G2≤(1+hγ1)∥y(n-1)-Hh(tn-1)∥G2+γ2h∥Y-(n)-Y(tn-τ)∥2+γ3h2p+1,h∈(0,h3],
where
(83)h3=min{h-0,h1}≤1,γ1=1+4p3p12∥C12∥2λ2,γ2=2p3(1+2p12)(1+2sr2q2),(84)γ3=8λ1p02p22.
Therefore,
(85)∥y(n)-Hh(tn)∥G2≤∥y(0)-Hh(t0)∥G2+h∑i=1n(γ1∥y(i-1)-Hh(ti-1)∥G2+γ2∥Y-(i)-Y(ti-τ)∥2)+γ3tnh2p,h∈(0,h3].
Considering Theorem 16, we further obtain
(86)∥y(n)-Hh(tn)∥G2≤∥y(0)-Hh(t0)∥G2+h(γ1+γ2p4)×∑i=0n-1∥y(i)-Hh(ti)∥G2+γ3tnh2p+γ2p4tn(h2(p+1)+h2(μ+v+1))+p4γ2h×∑k=-μ-m+10∥Y(k)-Y(tk)∥2,h∈(0,h0],
where h0=min{h2,h3}≤1. Using discrete Bellman inequality, we have
(87)∥y(n)-Hh(tn)∥G2≤(∥y(0)-Hh(t0)∥G2+p4γ2hoo×∑k=-μ-m+10∥Y(k)-Y(tk)∥2oo+(γ3+γ2p4)tnh2p+γ2p4tnh2(μ+v+1)∥y(0)-Hh(t0)∥G2)×exp((γ1+γ2p4)tn),h∈(0,h0].
Considering ∥H(t)-Hh(tn)∥≤p0hp, we obtain
(88)∥y(n)-H(tn)∥≤∥y(n)-Hh(tn)∥+∥Hh(tn)-H(tn)∥≤∥y(n)-Hh(tn)∥+p0hp.
Considering (87) and (88), we can easily conclude that method (4a), (4b), and (4c) with interpolation procedure (6) is D-convergent of order at least min{p,μ+v+1}. The proof is completed.


## 5. Numerical Experiments

Consider the following nonlinear neutral delay differential equations:(89)y′(t)=-y(t)+0.9cos(y(t-1)+y′(t-1))+cos(t)e-t-0.9cos(cos(t-τ)e-(t-τ)),10≥t≥0y(t)=e-tsin(t),-1≤t≤0,
and its perturbed problem
(90)z′(t)=-z(t)+0.9cos(z(t-1)+z′(t-1))+cos(t)e-t-0.9cos(cos(t-τ)e-(t-τ)),10≥t≥0,z(t)=e-tsin(t)+0.2,-1≤t≤0.We can calculateα=-1, β=0.9, γ=0.9, and q=0. Equation (89) has a unique true solution:
(91)y(t)=e-tsin(t),t≥-1.Apply the two-step Runge-Kutta method induced by the GL method in [12](92)to (89) and (90), Y~j(n) is computed by the Lagrange interpolation procedure with n=5 (see Figures 1 and 2 and Table 1).Table 1
Convergence order of two-step Runge-Kutta methods for system (89).
m
5
10
20
40
Convergence order
4.0012
3.9880
3.9864
3.9860Figure 1
The numerical solution of (92) for (89) with h=0.1.Figure 2
The numerical solution of (92) for (90) with h=0.1.It is obvious that the corresponding method for NDDEs is stable and convergent, and the convergence order ismin{4,5}.


## 6. Conclusions

In this paper we gave the stability and convergence results of two-step Runge-Kutta methods with linear interpolation procedure for solving nonlinear NDDEs (1). First, we gave the definitions of (k,l)-algebraically stable, algebraically stable and diagonally stable. Then we proved that if a TSRK method is (k,l)-algebraically stable, k<1 and dj>0  j=1,2,…,s, then the method with linear interpolation procedure is weak GAR(l)-stable. We also proved that if a TSRK is algebraically stable and diagonally stable and its generalized stage order is p, then the method with interpolation procedure is D-convergent of order at least min{p,μ+v+1}.We believe that the results presented in this paper can be extended to other general NDDEs. However, it is difficult to extend these results to more general neutral functional differential equations. Results extending the results presented in this paper to more general neutral functional differential equations and other delay differential equations such as delay integral differential equations will be discussed elsewhere.



---

*Source: 683137-2013-05-23.xml*