# On the One-Leg Methods for Solving Nonlinear Neutral Differential Equations with Variable Delay

Authors: Wansheng Wang; Shoufu Li
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

Based onA-stable one-leg methods and linear interpolations, we introduce four algorithms for solving neutral differential equations with variable delay. A natural question is which algorithm is better. To answer this question, we analyse the error behavior of the four algorithms and obtain their error bounds under a one-sided Lipschitz condition and some classical Lipschitz conditions. After extensively numerically experimenting, we give a positive conclusion.


---

## Body

## 1. Introduction

In this paper, we focus on establishing convergence properties of one-leg methods for the numerical solution of the neutral differential equations with variable time delay (NDDEs). These are equations of the form(1.1)y′(t)=f(t,y(t),y(η(t)),y′(η(t))),t∈I0=[0,T],(1.2)y(t)=ϕ(t),y′(t)=ϕ′(t),t∈Iτ=[-τ,0],
where -τ=inft∈I0η(t), τ>0, T>0, f and ϕ are complex smooth vector functions, and η is scalar function with η(t)≤t for t∈I0. Variations of (1.1) include equations with multiple time delays as follows:
(1.3)y′(t)=f(t,y(t),y(η1(t)),…,y(ηm(t)),y′(η1(t)),…,y′(ηm(t))),t∈I0,
and equations with state-dependent time delay as follows:
(1.4)y′(t)=f(t,y(t),y(η(t,y(t))),y′(η(t,y(t)))),t∈I0.The initial-value problems for ordinary differential equations (ODEs) and delay differential equations (DDEs) are special cases of (1.1)-(1.2).Neutral delay differential equations has found applications in many areas of science (see, e.g., [1–4]). The stability and convergence of numerical methods for NDDEs have drawn much attention in the last decades. A number of authors have investigated the linear stability of numerical methods for NDDEs (see, e.g., [5–12]). In 2000, based on a one-sided Lipschitz condition and some classical Lipschitz conditions, Bellen et al. [13] discussed the contractivity and asymptotic stability of Runge-Kutta methods for nonlinear NDDEs with special form. Following this paper, the nonlinear stability of numerical methods for NDDEs of the “Hale’ form” [14, 15] and for NDDEs (1.1) [16–22] has been extensively examined.On the other hand, as far back as the 1970s, Castleton and Grimm [23] investigated the convergence of a first-order numerical method for state-dependent delay NDDEs. Jackiewicz [24–28] and Jackiewicz and Lo [29, 30] investigated the convergence of numerical methods for more general neutral functional differential equations (NFDEs) which contains problem (1.1) as a particular case. Enright and Hayashi in [31] investigated the convergence of continuous Runge-Kutta methods for state-dependent delay NDDEs. Baker [32] and Zennaro [33] discussed the convergence of numerical methods for NDDEs. Jackiewicz et al. [34] and Bartoszewski and Kwapisz [35] gave the convergence results of waveform relaxation methods for NFDEs. Observe that these convergence results were based on the assumption that the right-hand function f in (1.1) satisfies a classical Lipschitz condition with respect to the second variable and the error bounds were obtained by using directly the classical Lipschitz constant. Recently, Wang et al. [36] gave some new error bounds of a class of linear multi-step methods for the NDDEs (1.1) by means of the one-sided Lipschitz constant; Wang and Li [37] investigated the convergence of waveform relaxation methods for the NDDEs (1.1) in which the right-hand function f satisfies a one-sided Lipschitz condition with respect to the second variable.The main contribution of this paper is that we introduce four algorithms for solving NDDEs (1.1) based on a one-leg method and present their convergence results. To accomplish this, we in Section 2 consider applying a one-leg method to NDDEs (1.1) and introduce two NDDEs numerical algorithms based on direct estimation. The convergence of the two algorithms then is analysized in Section 3. Noting that this class of algorithms may create some implementation problems, we introduce another two algorithms based on interpolation and analyze their convergence in Section 4. The application of the four algorithms for NDDEs is illustrated by means of two examples in Section 5. These numerical results confirm our theoretical analysis. We end with some concluding remarks in Section 6.


## 2. One-Leg Methods Discretization

Leth>0 be a fixed step size, and E be the translation operator Eyn=yn+1. The one-leg k-step method gives(2.1)ρ(E)yn=hf(σ(E)tn,σ(E)yn),n=0,1,…,
for ODEs will be said to be determined by “the method (ρ,σ),” where ρ and σ are defined by
(2.2)ρ(x)≡∑j=0k‍αjxj,σ(x)≡∑j=0k‍βjxj;|α0|+|β0|>0,αk≠0,αi,βi∈R,
with the additional requirements that ρ(x) and σ(x) are coprime polynomials and
(2.3)ρ(1)=0,ρ′(1)=σ(1)=1.
The application of the method (ρ,σ) to (1.1) yields
(2.4)ρ(E)yn=hf(σ(E)tn,σ(E)yn,yh(η(σ(E)tn)),Yh(η(σ(E)tn))),n=0,1,…,
where yh(η(σ(E)tn)) and Yh(η(σ(E)tn)) are approximations to y(η(σ(E)tn)) and y′(η(σ(E)tn)), respectively. For -m≤n≤0, yn=ϕ(tn). When -τ≤η(σ(E)tn)≤0, y(η(σ(E)tn))=ϕ(η(σ(E)tn, and y′(η(σ(E)tn))=ϕ′(η(σ(E)tn)); when 0≤η(σ(E)tn)≤T, y(η(σ(E)tn)) is obtained by a specific interpolation, and y′(η(σ(E)tn)) is obtained by a specific approximation at the point t=η(σ(E)tn).Therefore, the NDDEs method (2.4) is determined completely by the method (ρ,σ), the interpolation procedure for yh(η(σ(E)tn)) and the approximation procedure for Yh(η(σ(E)tn)).


### 2.1. The First Class of Algorithms Derived from the Direct Evaluation: OLIDE

Now we consider the different interpolation procedures foryh(η(σ(E)tn)) and the different approximation procedures for Yh(η(σ(E)tn)). It is well known that any A-stable one-leg method for ODEs has order of at most 2. So we can use the linear interpolation procedure for yh(η(σ(E)tn)). Let us define
(2.5)η(0)(t)=t,η(1)(t)=η(t),…,η(i)(t)=η(η(i-1)(t)),…
and suppose that
(2.6)η(i)(σ(E)tn)=(mn(i)+δn(i))h,
where mn(i) is an integer such that tmn(i)≥-τ, δn(i)∈[0,1). Then define
(2.7)yh(η(i)(σ(E)tn))={δn(i)ymn(i)+1+(1-δn(i))ymn(i),tmn(i)+δn(i)h>0,ϕ(tmn(i)+δn(i)h),tmn(i)+δn(i)h≤0,
where yj=ϕ(jh) for j≤0.In approximating nonneutral DDEs with constant delay, another interpolation procedure is used to approximateyh(η(σ(E)tn)) (see, e.g., [38]). Slightly modifying this interpolation procedure, we have
(2.8)yh(η(i)(σ(E)tj))=δ-j(i)σ(E)yj-m-j(i)+1+(1-δ-j(i))σ(E)yj-m-j(i),0≤j≤n,
provided
(2.9)σ(E)tj-η(i)(σ(E)tj)=(m-j(i)-δ-j(i))h,
where m-j(i)≥0 is an integer, δ-j(i)∈[0,1).ForYh(η(σ(E)tn)), we first consider the direct evaluation:
(2.10)Yh(η(i)(σ(E)tn))=f(η(i)(σ(E)tn),yh(η(i)(σ(E)tn)),yh(η(i+1)(σ(E)tn)),Yh(η(i+1)(σ(E)tn))).
Then we can obtain two algorithms for solving NDDEs: one is (2.4)–(2.7)–(2.10), which is simply denoted by OLIDE(I), the other is (2.4)–(2.8)–(2.10), which is simply denoted by OLIDE(II). In [18], Wang et al. have shown that the algorithm OLIDE(II) have better stability properties than the algorithm OLIDE(I) for DDEs with constant delay. It is interesting to prove that the methods OLIDE(I) and OLIDE(II) have the same convergence properties.


## 2.1. The First Class of Algorithms Derived from the Direct Evaluation: OLIDE

Now we consider the different interpolation procedures foryh(η(σ(E)tn)) and the different approximation procedures for Yh(η(σ(E)tn)). It is well known that any A-stable one-leg method for ODEs has order of at most 2. So we can use the linear interpolation procedure for yh(η(σ(E)tn)). Let us define
(2.5)η(0)(t)=t,η(1)(t)=η(t),…,η(i)(t)=η(η(i-1)(t)),…
and suppose that
(2.6)η(i)(σ(E)tn)=(mn(i)+δn(i))h,
where mn(i) is an integer such that tmn(i)≥-τ, δn(i)∈[0,1). Then define
(2.7)yh(η(i)(σ(E)tn))={δn(i)ymn(i)+1+(1-δn(i))ymn(i),tmn(i)+δn(i)h>0,ϕ(tmn(i)+δn(i)h),tmn(i)+δn(i)h≤0,
where yj=ϕ(jh) for j≤0.In approximating nonneutral DDEs with constant delay, another interpolation procedure is used to approximateyh(η(σ(E)tn)) (see, e.g., [38]). Slightly modifying this interpolation procedure, we have
(2.8)yh(η(i)(σ(E)tj))=δ-j(i)σ(E)yj-m-j(i)+1+(1-δ-j(i))σ(E)yj-m-j(i),0≤j≤n,
provided
(2.9)σ(E)tj-η(i)(σ(E)tj)=(m-j(i)-δ-j(i))h,
where m-j(i)≥0 is an integer, δ-j(i)∈[0,1).ForYh(η(σ(E)tn)), we first consider the direct evaluation:
(2.10)Yh(η(i)(σ(E)tn))=f(η(i)(σ(E)tn),yh(η(i)(σ(E)tn)),yh(η(i+1)(σ(E)tn)),Yh(η(i+1)(σ(E)tn))).
Then we can obtain two algorithms for solving NDDEs: one is (2.4)–(2.7)–(2.10), which is simply denoted by OLIDE(I), the other is (2.4)–(2.8)–(2.10), which is simply denoted by OLIDE(II). In [18], Wang et al. have shown that the algorithm OLIDE(II) have better stability properties than the algorithm OLIDE(I) for DDEs with constant delay. It is interesting to prove that the methods OLIDE(I) and OLIDE(II) have the same convergence properties.


## 3. Convergence Analysis of Algorithms OLIDE

### 3.1. A General Assumption

In the current section, we shall consider the convergence of the algorithms OLIDE for NDDEs (1.1)-(1.2). To do this, we always make the following assumptions unless otherwise stated: (𝒜1) There exist constants, α, ϱ, β, and γ for which the following inequalities hold:
(3.1)ℜe〈y1-y2,f(t,y1,u,v)-f(t,y2,u,v)〉≤α∥y1-y2∥2,(3.2)∥f(t,y,u1,v1)-f(t,y,u2,v2)∥≤β∥u1-u2∥+γ∥v1-v2∥,(3.3)∥f(t,y,u1,f(η(t),u1,v,w))-f(t,y,u2,f(η(t),u2,v,w))∥≤ϱ∥u1-u2∥,
where 〈·,·〉 is the inner product in complex Euclid space with the corresponding norm ∥·∥.(𝒜2)  α0=max{α,0} and ϱ/(1-γ) are of moderate size.(𝒜3) We assume the existence, uniqueness, and stability of solutions to the mathematical problems under consideration. For example, it is required that γ<1 (see, e.g., [31]). This restriction ensures that when η(t)=t, the matrix
(3.4)[I-∂f∂y′]|t,y(t),y′(t)
is full rank so that y′(t) is uniquely determined as a function of t and y(t) at this point.(𝒜4) The unique true solution y(t) of problem (1.1)-(1.2) is also sufficiently differentiable on the interval [-τ,T], and all its derivatives used later are continuous and satisfy
(3.5)∥diy(t)dti∥≤Mi,-τ≤t≤T.Remarks
(1) In general, the regularity assumption(𝒜4) does not hold for NDDEs (1.1)-(1.2). Even if the functions f and ϕ are sufficiently smooth, the solution has low regularity at the discontinuity points (see, e.g., [11]). This problem has attracted some researchers’ attention. Some approaches to handling discontinuities have been proposed, for example, discontinuity tracking [39, 40], discontinuity detection [27, 30, 31, 41, 42], perturbing the initial function [43]. Since the main purpose of this paper is to present some convergence results but not to discuss the treatment of discontinuities, we will think that the regularity assumption (𝒜4) holds.
(2) The regularity assumption(𝒜4) does not hold for general NDDEs (1.1)-(1.2) but holds for NDDEs (1.1)-(1.2) with proportional delay, that is, η(t)=qt,q∈(0.1), also holds for NDDEs (1.1)-(1.2) with η(t)=t, that is, implicit ODEs.
(3) For the case that
𝒞
1 the function η(t) satisfies such conditions that we can divide the time interval I0 into some subintervals Ik=[ξk-1,ξk] for k≥1, where ξ0=0 and ξk+1 is the unique solution of η(ξ)=ξk, then the analysis of the error behaviour of the solutions can be done interval-by-interval since the regularity assumption (𝒜4) generally holds on each subinterval Ik. For example, the function η(t) satisfies (𝒞1) if it is continuous and increasing and there exists a positive constant τ0 such that τ(t)=t-η(t)≥τ0, for all t∈I0,The aim of this section is to derive theoretical estimates for the errors generated by the algorithms OLIDE(I) and OLIDE(II). To obtain the theoretical results, we need the following definition.Definition 3.1.
The one-leg method(ρ,σ) with two approximation procedures is said to be EB-convergent of order p if this method when applied to any given problem (1.1) satisfying the assumptions (𝒜1)–(𝒜4) with initial values y0,y1,…,yk-1, produces an approximation {yn}, and the global error satisfies a bound of the form
(3.6)∥y(tn)-yn∥≤C(tn)(hp+max0≤i≤k-1∥y(ti)-yi∥),n≥k,h∈(0,h0],
where the function C(t) and the maximum step-size h0 depend only on the method, some of the bounds Mi, the parameters α and ϱ/(1-γ).Proposition 3.2.
E
B-convergence of the method (2.4) implies B-convergence of the method (ρ,σ).Remark 3.3.
For nonneutral DDEs, Huang et al. [44] introduced the D-convergence of one-leg methods. Obviously, when an EB-convergent NDDEs method is applied to nonneutral DDEs, it is D-convergent.


### 3.2. Convergence Analysis of OLIDE(I)

First, we consider the convergence of the algorithm OLIDE(I). For this purpose, let us consider(3.7)ρ(E)y^n+αken=hf(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn))),n=0,1,…,
where
(3.8)y^n=y(tn)+c1h2y′′(tn),(3.9)y-h(η(i)(σ(E)tn))={y(η(i)(σ(E)tn)),ifη(i)(σ(E)tn)≤tn+k-1,δn(i)(y^n+k+en)+(1-δn(i))y(tn+k-1),iftn+k-1<η(i)(σ(E)tn)≤tn+k,(3.10)Y-h(η(i)(σ(E)tn))=f((η(i)(σ(E)tn))Y-h(η(i)(σ(E)tn))η(i)(σ(E)tn),y-h(η(i)(σ(E)tn)),y-h(η(i+1)(σ(E)tn)),Y-h(η(i)(σ(E)tn))),(3.11)c1=-12∑j=0k-1‍(βj-βkαkαj)j2-βkαk∑j=0k‍jβj+12(∑j=0k‍jβj)2.
It can be seen that for n≥0, when the step-size h satisfies certain conditions, en is defined completely by (3.7).For any givenk×k real symmetric positive definite matrix G=[gij], the norm ∥·∥G is defined as
(3.12)∥U∥G=(∑i,j=1k‍gij〈ui,uj〉)1/2,∀U=[u1,u2,…,us]∈CNk.Lemma 3.4.
If the method(ρ,σ) is A-stable, then the numerical solution produced by the algorithm OLIDE(I) satisfies
(3.13)∥εn+1∥G2≤(1+d1h)max0≤i≤n∥εi∥G2+d2h5+d3h-1∥en∥2,n=0,1,…,h∈(0,h1],
where εn=[(yn-y^n)T,(yn+1-y^n+1)T,…,(yn+k-1-y^n+k-1)T]T.Proof.
SinceA-stability is equivalent to G-stability (see [45]), it follows from the definition of G-stability that there exists a k×k real symmetric positive definite matrix G such that for any real sequence {ai}i=0k(3.14)A1TGA1-A0TGA0≤2σ(E)a0ρ(E)a0,
where Ai=(ai,ai+1,…,ai+k-1)T(i=0,1). Therefore, we can easily obtain (see [45, 46])
(3.15)∥εn+1∥G2-∥εn∥G2≤2ℜe〈σ(E)(yn-y^n),ρ(E)(yn-y^n)〉.
Write ε^n+1=[(yn+1-y^n+1)T,…,(yn+k-1-y^n+k-1)T,(yn+k-y^n+k-en)T]T, and note the difference between it and εn+1. Then use the argument technique above to obtain
(3.16)∥ε^n+1∥G2≤∥εn∥G2+2ℜe〈σ(E)(yn-y^n)-βken,ρ(E)(yn-y^n)-αken〉.
By means of condition (3.1), it follows from (3.16) that
(3.17)∥ε^n+1∥G2≤∥εn∥G2+2hℜe〈(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))σ(E)(yn-y^n)-βken,f(σ(E)tn,σ(E)yn,yh(η(σ(E)tn)),Yh(η(σ(E)tn)))-f(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))〉≤∥εn∥G2+2h[(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))α∥σ(E)(yn-y^n)-βken∥2+∥σ(E)(yn-y^n)-βken∥×∥(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))f(σ(E)tn,σ(E)y^n+βken,yh(η(σ(E)tn)),Yh(η(σ(E)tn)))-f(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))∥].
If there exists a positive integer l such that η(l)(σ(E)tn)≤0, then in view of yh(η(l)(σ(E)tn))-y-h(η(l)(σ(E)tn))=0, we have
(3.18)∥(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))f(σ(E)tn,σ(E)y^n+βken,yh(η(σ(E)tn)),Yh(η(σ(E)tn)))-f(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))∥≤ϱ∥yh(η(σ(E)tn))-y-h(η(σ(E)tn))∥+γ∥(η(σ(E)tn),y-h(η(σ(E)tn)),y-h(η(2)(σ(E)tn)),Y-h(η(2)(σ(E)tn)))f(η(σ(E)tn),y-h(η(σ(E)tn)),yh(η(2)(σ(E)tn)),Yh(η(2)(σ(E)tn)))-f(η(σ(E)tn),y-h(η(σ(E)tn)),y-h(η(2)(σ(E)tn)),Y-h(η(2)(σ(E)tn)))∥≤ϱ1-γmaxi≥1∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥.
Conversely, if for all l, η(l)(σ(E)tn)>0, we similarly have (3.18). In short, it follows from (3.16) that
(3.19)∥ε^n+1∥G2≤∥εn∥G2+2h[ϱ1-γmaxi≥1∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥α∥σ(E)(yn-y^n)-βken∥2+∥σ(E)(yn-y^n)-βken∥×ϱ1-γmaxi≥1∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥].
As an important step toward the proof of this lemma, we show that
(3.20)∥ε^n+1∥G2≤∥εn∥G2+c2h∥σ(E)(yn-y^n)-βken∥2+2ϱh1-γ(1+|c1|)2M22h4+2ϱh1-γmax{∥yn+k-y^n+k-en∥2,max1≤i≤n+k-1∥yi-y^i∥2},
where
(3.21)c2={0,2α+ϱ1-γ≤0,2α+ϱ1-γ,2α+ϱ1-γ>0.
For this purpose, we consider the following two cases successively.
Case1.  tn+k-1≤η(i)(σ(E)tn)≤tn+k. In this case, from (2.7) and (3.9) we have
(3.22)∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥≤(1-δn(i))∥yn+k-1-y-h(tn+k-1)∥+δn(i)∥yn+k-y^n+k-en∥≤(1-δn(i))[∥yn+k-1-y^n+k-1∥+|c1|M2h2]+δn(i)∥yn+k-y^n+k-en∥≤max{∥yn+k-1-y^n+k-1∥,∥yn+k-y^n+k-en∥}+|c1|M2h2.Case2.  η(i)(σ(E)tn)<tn+k-1. For this case, similarly, it follows from (2.7) and (3.9) that
(3.23)∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥≤∥(1-δn(i))yn-mn(i)+δn(i)yn-mn(i)+1-y(tn-mn(i)+δn(i)h)∥≤(1-δn(i))∥yn-mn(i)-y^n-mn(i)∥+δn(i)∥yn-mn(i)+1-y^n-mn(i)+1∥+∥δn(i)y^n-mn(i)+1+(1-δn(i))y^n-mn(i)-y(tn-mn(i)+δn(i)h)∥≤max{∥yn-mn(i)-y^n-mn(i)∥,∥yn-mn(i)+1-y^n-mn(i)+1∥}+(1+|c1|)M2h2.
Substituting (3.22) and (3.23) into (3.19), and using the Cauchy inequality yields (3.20). From (3.20) one further gets
(3.24)∥ε^n+1∥G2≤∥εn∥G2+c2h∥σ(E)(yn-y^n)-βken∥2+2ϱh1-γ(1+|c1|)2M22h4+2ϱhλminG(1-γ)∥ε^n+1∥G2+2ϱhλminG(1-γ)max0≤i≤n∥εi∥G2,
where λminG denotes the minimum eigenvalues of the matrix G. On the other hand, it is easy to verify that there exists a constant d4 such that
(3.25)∥σ(E)(yn-y^n)∥2=∥∑j=0k‍βj(yn+j-y^n+j)∥2≤d4∑j=0k‍∥yn+j-y^n+j∥2≤d4λminG[∥εn∥G2+∥εn+1∥G2].
Thus, we get
(3.26)∥σ(E)(yn-y^n)-βken∥2≤2∥σ(E)(yn-y^n)∥2+2βk2∥en∥2≤2d4λminG[∥εn∥G2+∥εn+1∥G2]+2βk2∥en∥2.
Substituting the above inequality into (3.24) yields
(3.27)∥ε^n+1∥G2≤(1+2c2d4hλminG)∥εn∥G2+2c2d4hλminG∥εn+1∥G2+2c2βk2h∥en∥2+2ϱh1-γ(1+|c1|)2M22h4+2ϱhλminG(1-γ)∥ε^n+1∥G2+2ϱhλminG(1-γ)max0≤i≤n∥εi∥G2.
For any given c0∈(0,1), we let c-2=(1-c0)λminG/2. Then for (ϱ/(1-γ))h≤c-2, the above inequality leads to
(3.28)∥ε^n+1∥G2≤11-(2ϱh/λminG(1-γ))[(1+2c2d4hλminG)∥εn∥G2+2c2d4hλminG∥εn+1∥G2+2c2βk2h∥en∥2+2ϱh51-γ(1+|c1|)2M22]+2ϱhλminG(1-γ)-2ρhmax0≤i≤n∥εi∥G2≤(1+d5h)max0≤i≤n∥εi∥G2+2c2d4hc0λminG∥εn+1∥G2+2c2βk2hc0∥en∥2+2ϱh5c0(1-γ)(1+|c1|)2M22,
where d5=2c2d4/c0λminG+4ϱ/c0λminG(1-γ). Noting that
(3.29)∥εn+1∥G2≤∥ε^n+1∥G2+λmaxG∥en∥2+2λmaxG∥en∥∥ε^n+1∥G≤(1+h)∥ε^n+1∥G2+(1+1h)λmaxG∥en∥2,
where λmaxG denotes the maximum eigenvalues of the matrix G, we have
(3.30)∥εn+1∥G2≤(1+h)[(1+d5h)max0≤i≤n∥εi∥G2+2c2d4hc0λminG∥εn+1∥G2+2c2βk2hc0∥en∥2+2ϱh5c0(1-γ)(1+|c1|)2M22]+(1+1h)λmaxG∥en∥2.
Now for any given c~0∈(0,1), when (4c2d4/c0λminG)h≤1-c~0, the above inequality leads to (3.13) with
(3.31)h1=min{1,(1-γ)c-2ϱ,(1-c~0)c0λminG4c2d4},d1=1+2d5c~0+4c2d4c~0c0λminG,d2=4ϱ(1+|c1|)2M22(1-γ)c~0c0,d3=4c2βk2c~0c0+2λmaxGc~0.
This completes the proof of this lemma.Compared to the nonneutral DDEs with a constant delay considered in [38, 44], the problems considered in this paper are more complex such that a series of new difficulties need to be overcome. Nevertheless, we note that some basic proof ideas are related to the ones used in [38, 44].Lemma 3.5.
If the method(ρ,σ) is A-stable, then there exist two constants d6 and h2, which depend only on the method, some of the bounds Mi, and the parameters α and ϱ/(1-γ), such that
(3.32)∥en∥≤d6hp+1,h∈(0,h2],n=0,1,…,
where p is the consistent order in the classical sense, p=1,2.Proof.
The idea is a generalization of that used in [38, 44]. The A-stability of the method (ρ,σ) implies that βk/αk>0, and it is consistent of order p=1,2 in the classical sense (see [45, 46]). Consider the following scheme:
(3.33)y(σ(E)tn)=∑j=0k-1‍(βj-βkαkαj)y^n+j+βkαkhy′(σ(E)tn)+R1(n),(3.34)ρ(E)y^n=hy′(σ(E)tn)+R2(n).
By Taylor expansion, we can assert that there exists constant c3 such that
(3.35)R1(n)≤c3h3,R2(n)≤c3hp+1.
From (3.7) and (3.33), it is easy to obtain
(3.36)y(σ(E)tn)-σ(E)y^n-βken=R1(n)+βkαkh[y′(σ(E)tn)-f(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))].
If η(i)(σ(E)tn)∈(tn+k-1,tn+k], by Taylor expansion, we have
(3.37)∥y(η(i)(σ(E)tn))-δn(i)(y^n+k+en)-(1-δn(i))y(tn+k-1)∥≤c4h2M2+∥en∥.
This inequality and (3.36) imply that
(3.38)∥y(σ(E)tn)-σ(E)y^n-βken∥2≤βkαkh∥y(σ(E)tn)-σ(E)y^n-βken∥×[α∥y(σ(E)tn)-σ(E)y^n-βken∥+ϱ1-γ(c4h2M2+∥en∥)]+R1(n)∥y(σ(E)tn)-σ(E)y^n-βken∥.
For any given c-0∈(0,1), we let c-3=(1-c-0)αk/βk. Then for αh≤c-3, the above inequality leads to
(3.39)∥y(σ(E)tn)-σ(E)y^n-βken∥≤1c-0[βkϱhαk(1-γ)(c4h2M2+∥en∥)+∥R1(n)∥]≤c5h3+βkϱhαkc-0(1-γ)∥en∥,
where
(3.40)c5=1c-0[βkϱc4M2αk(1-γ)+c3].
Substitute (3.39) into (3.36) to obtain
(3.41)∥hy′(σ(E)tn)-hf(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))∥≤αkβk∥y(σ(E)tn)-σ(E)y^n-βken∥+αkβk∥R1(n)∥≤αkβk(c3+c5)h3+ϱhc-0(1-γ)∥en∥.
On the other hand, since
(3.42)αken=hf(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))-hy′(σ(E)tn)-R2(n)
it follows from (3.41) that
(3.43)|αk|∥en∥≤αkβk(c3+c5)h3+ϱhc-0(1-γ)∥en∥+c3hp+1.
For any given c^0∈(0,1), let
(3.44)h2={min{h1,(1-γ)(1-c^0)c-0|αk|ϱ},ifα≤0,min{h1,(1-γ)(1-c^0)c-0|αk|ϱ,c-3α},ifα>0.
Then the required inequality (3.32) follows from (3.43), where
(3.45)d6=1c^0(c3+c5|βk|+c3|αk|).
This completes the proof of the lemma.Lemma3.4, together with Lemma 3.5, implies the following theorem.Theorem 3.6.
The algorithm OLIDE(I) extended by the method(ρ,σ) is EB-convergent of order p if and only if the method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2.Proof.
Firstly, we observe that theEB-convergence of order p of the method OLIDE(I) implies that the method (ρ,σ) is B-convergent of order p, and it is A-stable and consistent of order p in the classical sense for ODEs, p=1 or 2 (see, e.g., [38]).
On the other hand, it follows from Lemmas3.4 and 3.5 that
(3.46)∥εn+1∥G2≤(1+d1h)max0≤i≤n∥εi∥G2+d2h5+d3d6h2p+1,n=0,1,…,h∈(0,h2]
which implies that
(3.47)max0≤i≤n+1∥εn+1∥G2≤(1+d1h)max0≤i≤n∥εi∥G2+d2h5+d3d6h2p+1,n=0,1,…,h∈(0,h2].
By induction, we have
(3.48)∥εn∥G2≤max0≤i≤n∥εi∥G2≤(1+d1h)n∥ε0∥G2+(d2h5+d3d6h2p+1)[1+(1+d1h)+⋯+(1+d1h)n-1]≤(1+d1h)n[∥ε0∥G2+n(d2h5+d3d6h2p+1)]≤exp(d1tn)[∥ε0∥G2+tn(d2+d3d6)h2p]
and therefore
(3.49)∥yn-y^n∥2≤exp(d1tn)λminG[λmaxG∑j=0k-1‍∥yj-y^j∥2+tn(d2+d3d6)h2p].
Consequently, the following holds:
(3.50)∥yn-y(tn)∥=∥yn-y^n+y^n-y(tn)∥≤|c1|M2h2+[kλmaxGλminG(max0≤i≤k-1∥yi-y(ti)∥+|c1|M2h2)+tn(d2+d3d6)λminGhp]×exp(12d1tn),n=0,1,…,h∈(0,h2]
which implies the algorithm OLIDE(I) is EB-convergent of order p,p=1 or 2. The proof of Theorem 3.6 is completed.In [44], Huang et al. proved that A-stable one-leg methods with the interpolation (2.7) for constant-delay DDEs are D-convergent under a condition that η(i)(σ(E)tn)≤tn+k-1. As a special case, we have the following corollary.Corollary 3.7.
A one-leg method(ρ,σ) with the interpolation (2.7) for DDEs with any variable delay is D-convergent of order p if and only if the method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2.


### 3.3. Convergence Analysis of OLIDE(II)

In [18], Wang et al. have shown that the algorithm OLIDE(II) has better stability properties than the algorithm OLIDE(I) for DDEs with constant delay. Here we will prove that the algorithms OLIDE(II) and OLIDE(I) have the same convergence properties. To obtain the convergence result of the algorithm OLIDE(II), we also consider (3.7). But in this case, y-h(η(i)(σ(E)tn))  (i=1,2,…) are determined by
(3.51)y-h(η(i)(σ(E)tn))={βk(y^n+k+en)+∑j=0k-1‍βjy(tn+j),iftn+k-1<η(i)(σ(E)tn)≤tn+k,m-n(i)=0,δ-n(i)[βk(y^n+k+en)+∑j=0k-1‍βjy(tn+j)]+(1-δ-n(i))∑j=0k‍βjy(tn-1+j),iftn+k-1<η(i)(σ(E)tn)≤tn+k,m-n(i)=1,y(η(i)(σ(E)tn)),othercase.Whenη(i)(σ(E)tn)≤tn+k-1, we have
(3.52)∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥=∥(1-δ-n(i))σ(E)yn-m-n(i)+δ-n(i)σ(E)yn-m-n(i)+1-y(η(i)(σ(E)tn))∥≤(1-δ-n(i))∥σ(E)(yn-m-n(i)-y^n-m-n(i))∥+δ-n(i)∥σ(E)(yn-m-n(i)+1-y^n-mn(i)+1)∥+∥(1-δ-n(i))σ(E)y^n-m-n(i)+δ-n(i)σ(E)y^n-m-n(i)+1-y(η(i)(σ(E)tn))∥≤(1-δ-n(i))∑j=0k‍|βj|∥yn-m-n(i)+j-y^n-m-n(i)+j∥+δ-n(i)∑j=0k‍|βj|∥yn-m-n(i)+j+1-y^n-m-n(i)+j+1∥+∥(1-δ-n(i))∑j=0k‍βjY^n-m-n(i)+j+δ-n(i)∑j=0k‍βjY^n-m-n(i)+j+1-y′(η(i)(σ(E)tn))∥,
where
(3.53)Y^i=f(ti,y^i,y-h(η(ti)),Y-h(η(ti))),i=1,2,…,n+k-1.
By Taylor expansion, there exist constants c-1 and h3 such that
(3.54)∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥≤c-1[max0≤j≤k+1∥yn+j-m-n(i)-y^n+j-m-n(i)∥+h2],h∈(0,h3].
Noting that for tn+k-1≤η(σ(E)tn)≤tn+k we have similar inequality, we find that the proof of the following theorem is similar to that of Theorem 3.6.Theorem 3.8.
The algorithm OLIDE(II) extended by the method(ρ,σ) is EB-convergent of order p if and only if the method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2.It should be pointed out that Huang et al. [38] proved that A-stable one-leg methods with the interpolation (2.8) for constant-delay DDEs are D-convergent. As a special case, in this paper we obtain the convergence result of A-stable one-leg methods with the interpolation (2.8) for DDEs with any variable delay.Since the algorithm OLIDE(II) has better stability properties than the algorithm OLIDE(I) [18] and they have the same convergence properties, we believe that the algorithm OLIDE(II) is more effective than the algorithm OLIDE(I) for solving delay problem.Although we can show theoretically the convergence of the algorithms OLIDE for solving NDDEs with any variable delay, in practical implementation, these algorithms are generally expensive and will create a serious storage problem since they always require to trace back the recursion until the initial interval. So for general variable delay NDDEs, we need other algorithms, which are based on interpolation.


## 3.1. A General Assumption

In the current section, we shall consider the convergence of the algorithms OLIDE for NDDEs (1.1)-(1.2). To do this, we always make the following assumptions unless otherwise stated: (𝒜1) There exist constants, α, ϱ, β, and γ for which the following inequalities hold:
(3.1)ℜe〈y1-y2,f(t,y1,u,v)-f(t,y2,u,v)〉≤α∥y1-y2∥2,(3.2)∥f(t,y,u1,v1)-f(t,y,u2,v2)∥≤β∥u1-u2∥+γ∥v1-v2∥,(3.3)∥f(t,y,u1,f(η(t),u1,v,w))-f(t,y,u2,f(η(t),u2,v,w))∥≤ϱ∥u1-u2∥,
where 〈·,·〉 is the inner product in complex Euclid space with the corresponding norm ∥·∥.(𝒜2)  α0=max{α,0} and ϱ/(1-γ) are of moderate size.(𝒜3) We assume the existence, uniqueness, and stability of solutions to the mathematical problems under consideration. For example, it is required that γ<1 (see, e.g., [31]). This restriction ensures that when η(t)=t, the matrix
(3.4)[I-∂f∂y′]|t,y(t),y′(t)
is full rank so that y′(t) is uniquely determined as a function of t and y(t) at this point.(𝒜4) The unique true solution y(t) of problem (1.1)-(1.2) is also sufficiently differentiable on the interval [-τ,T], and all its derivatives used later are continuous and satisfy
(3.5)∥diy(t)dti∥≤Mi,-τ≤t≤T.Remarks
(1) In general, the regularity assumption(𝒜4) does not hold for NDDEs (1.1)-(1.2). Even if the functions f and ϕ are sufficiently smooth, the solution has low regularity at the discontinuity points (see, e.g., [11]). This problem has attracted some researchers’ attention. Some approaches to handling discontinuities have been proposed, for example, discontinuity tracking [39, 40], discontinuity detection [27, 30, 31, 41, 42], perturbing the initial function [43]. Since the main purpose of this paper is to present some convergence results but not to discuss the treatment of discontinuities, we will think that the regularity assumption (𝒜4) holds.
(2) The regularity assumption(𝒜4) does not hold for general NDDEs (1.1)-(1.2) but holds for NDDEs (1.1)-(1.2) with proportional delay, that is, η(t)=qt,q∈(0.1), also holds for NDDEs (1.1)-(1.2) with η(t)=t, that is, implicit ODEs.
(3) For the case that
𝒞
1 the function η(t) satisfies such conditions that we can divide the time interval I0 into some subintervals Ik=[ξk-1,ξk] for k≥1, where ξ0=0 and ξk+1 is the unique solution of η(ξ)=ξk, then the analysis of the error behaviour of the solutions can be done interval-by-interval since the regularity assumption (𝒜4) generally holds on each subinterval Ik. For example, the function η(t) satisfies (𝒞1) if it is continuous and increasing and there exists a positive constant τ0 such that τ(t)=t-η(t)≥τ0, for all t∈I0,The aim of this section is to derive theoretical estimates for the errors generated by the algorithms OLIDE(I) and OLIDE(II). To obtain the theoretical results, we need the following definition.Definition 3.1.
The one-leg method(ρ,σ) with two approximation procedures is said to be EB-convergent of order p if this method when applied to any given problem (1.1) satisfying the assumptions (𝒜1)–(𝒜4) with initial values y0,y1,…,yk-1, produces an approximation {yn}, and the global error satisfies a bound of the form
(3.6)∥y(tn)-yn∥≤C(tn)(hp+max0≤i≤k-1∥y(ti)-yi∥),n≥k,h∈(0,h0],
where the function C(t) and the maximum step-size h0 depend only on the method, some of the bounds Mi, the parameters α and ϱ/(1-γ).Proposition 3.2.
E
B-convergence of the method (2.4) implies B-convergence of the method (ρ,σ).Remark 3.3.
For nonneutral DDEs, Huang et al. [44] introduced the D-convergence of one-leg methods. Obviously, when an EB-convergent NDDEs method is applied to nonneutral DDEs, it is D-convergent.


## 3.2. Convergence Analysis of OLIDE(I)

First, we consider the convergence of the algorithm OLIDE(I). For this purpose, let us consider(3.7)ρ(E)y^n+αken=hf(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn))),n=0,1,…,
where
(3.8)y^n=y(tn)+c1h2y′′(tn),(3.9)y-h(η(i)(σ(E)tn))={y(η(i)(σ(E)tn)),ifη(i)(σ(E)tn)≤tn+k-1,δn(i)(y^n+k+en)+(1-δn(i))y(tn+k-1),iftn+k-1<η(i)(σ(E)tn)≤tn+k,(3.10)Y-h(η(i)(σ(E)tn))=f((η(i)(σ(E)tn))Y-h(η(i)(σ(E)tn))η(i)(σ(E)tn),y-h(η(i)(σ(E)tn)),y-h(η(i+1)(σ(E)tn)),Y-h(η(i)(σ(E)tn))),(3.11)c1=-12∑j=0k-1‍(βj-βkαkαj)j2-βkαk∑j=0k‍jβj+12(∑j=0k‍jβj)2.
It can be seen that for n≥0, when the step-size h satisfies certain conditions, en is defined completely by (3.7).For any givenk×k real symmetric positive definite matrix G=[gij], the norm ∥·∥G is defined as
(3.12)∥U∥G=(∑i,j=1k‍gij〈ui,uj〉)1/2,∀U=[u1,u2,…,us]∈CNk.Lemma 3.4.
If the method(ρ,σ) is A-stable, then the numerical solution produced by the algorithm OLIDE(I) satisfies
(3.13)∥εn+1∥G2≤(1+d1h)max0≤i≤n∥εi∥G2+d2h5+d3h-1∥en∥2,n=0,1,…,h∈(0,h1],
where εn=[(yn-y^n)T,(yn+1-y^n+1)T,…,(yn+k-1-y^n+k-1)T]T.Proof.
SinceA-stability is equivalent to G-stability (see [45]), it follows from the definition of G-stability that there exists a k×k real symmetric positive definite matrix G such that for any real sequence {ai}i=0k(3.14)A1TGA1-A0TGA0≤2σ(E)a0ρ(E)a0,
where Ai=(ai,ai+1,…,ai+k-1)T(i=0,1). Therefore, we can easily obtain (see [45, 46])
(3.15)∥εn+1∥G2-∥εn∥G2≤2ℜe〈σ(E)(yn-y^n),ρ(E)(yn-y^n)〉.
Write ε^n+1=[(yn+1-y^n+1)T,…,(yn+k-1-y^n+k-1)T,(yn+k-y^n+k-en)T]T, and note the difference between it and εn+1. Then use the argument technique above to obtain
(3.16)∥ε^n+1∥G2≤∥εn∥G2+2ℜe〈σ(E)(yn-y^n)-βken,ρ(E)(yn-y^n)-αken〉.
By means of condition (3.1), it follows from (3.16) that
(3.17)∥ε^n+1∥G2≤∥εn∥G2+2hℜe〈(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))σ(E)(yn-y^n)-βken,f(σ(E)tn,σ(E)yn,yh(η(σ(E)tn)),Yh(η(σ(E)tn)))-f(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))〉≤∥εn∥G2+2h[(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))α∥σ(E)(yn-y^n)-βken∥2+∥σ(E)(yn-y^n)-βken∥×∥(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))f(σ(E)tn,σ(E)y^n+βken,yh(η(σ(E)tn)),Yh(η(σ(E)tn)))-f(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))∥].
If there exists a positive integer l such that η(l)(σ(E)tn)≤0, then in view of yh(η(l)(σ(E)tn))-y-h(η(l)(σ(E)tn))=0, we have
(3.18)∥(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))f(σ(E)tn,σ(E)y^n+βken,yh(η(σ(E)tn)),Yh(η(σ(E)tn)))-f(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))∥≤ϱ∥yh(η(σ(E)tn))-y-h(η(σ(E)tn))∥+γ∥(η(σ(E)tn),y-h(η(σ(E)tn)),y-h(η(2)(σ(E)tn)),Y-h(η(2)(σ(E)tn)))f(η(σ(E)tn),y-h(η(σ(E)tn)),yh(η(2)(σ(E)tn)),Yh(η(2)(σ(E)tn)))-f(η(σ(E)tn),y-h(η(σ(E)tn)),y-h(η(2)(σ(E)tn)),Y-h(η(2)(σ(E)tn)))∥≤ϱ1-γmaxi≥1∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥.
Conversely, if for all l, η(l)(σ(E)tn)>0, we similarly have (3.18). In short, it follows from (3.16) that
(3.19)∥ε^n+1∥G2≤∥εn∥G2+2h[ϱ1-γmaxi≥1∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥α∥σ(E)(yn-y^n)-βken∥2+∥σ(E)(yn-y^n)-βken∥×ϱ1-γmaxi≥1∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥].
As an important step toward the proof of this lemma, we show that
(3.20)∥ε^n+1∥G2≤∥εn∥G2+c2h∥σ(E)(yn-y^n)-βken∥2+2ϱh1-γ(1+|c1|)2M22h4+2ϱh1-γmax{∥yn+k-y^n+k-en∥2,max1≤i≤n+k-1∥yi-y^i∥2},
where
(3.21)c2={0,2α+ϱ1-γ≤0,2α+ϱ1-γ,2α+ϱ1-γ>0.
For this purpose, we consider the following two cases successively.
Case1.  tn+k-1≤η(i)(σ(E)tn)≤tn+k. In this case, from (2.7) and (3.9) we have
(3.22)∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥≤(1-δn(i))∥yn+k-1-y-h(tn+k-1)∥+δn(i)∥yn+k-y^n+k-en∥≤(1-δn(i))[∥yn+k-1-y^n+k-1∥+|c1|M2h2]+δn(i)∥yn+k-y^n+k-en∥≤max{∥yn+k-1-y^n+k-1∥,∥yn+k-y^n+k-en∥}+|c1|M2h2.Case2.  η(i)(σ(E)tn)<tn+k-1. For this case, similarly, it follows from (2.7) and (3.9) that
(3.23)∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥≤∥(1-δn(i))yn-mn(i)+δn(i)yn-mn(i)+1-y(tn-mn(i)+δn(i)h)∥≤(1-δn(i))∥yn-mn(i)-y^n-mn(i)∥+δn(i)∥yn-mn(i)+1-y^n-mn(i)+1∥+∥δn(i)y^n-mn(i)+1+(1-δn(i))y^n-mn(i)-y(tn-mn(i)+δn(i)h)∥≤max{∥yn-mn(i)-y^n-mn(i)∥,∥yn-mn(i)+1-y^n-mn(i)+1∥}+(1+|c1|)M2h2.
Substituting (3.22) and (3.23) into (3.19), and using the Cauchy inequality yields (3.20). From (3.20) one further gets
(3.24)∥ε^n+1∥G2≤∥εn∥G2+c2h∥σ(E)(yn-y^n)-βken∥2+2ϱh1-γ(1+|c1|)2M22h4+2ϱhλminG(1-γ)∥ε^n+1∥G2+2ϱhλminG(1-γ)max0≤i≤n∥εi∥G2,
where λminG denotes the minimum eigenvalues of the matrix G. On the other hand, it is easy to verify that there exists a constant d4 such that
(3.25)∥σ(E)(yn-y^n)∥2=∥∑j=0k‍βj(yn+j-y^n+j)∥2≤d4∑j=0k‍∥yn+j-y^n+j∥2≤d4λminG[∥εn∥G2+∥εn+1∥G2].
Thus, we get
(3.26)∥σ(E)(yn-y^n)-βken∥2≤2∥σ(E)(yn-y^n)∥2+2βk2∥en∥2≤2d4λminG[∥εn∥G2+∥εn+1∥G2]+2βk2∥en∥2.
Substituting the above inequality into (3.24) yields
(3.27)∥ε^n+1∥G2≤(1+2c2d4hλminG)∥εn∥G2+2c2d4hλminG∥εn+1∥G2+2c2βk2h∥en∥2+2ϱh1-γ(1+|c1|)2M22h4+2ϱhλminG(1-γ)∥ε^n+1∥G2+2ϱhλminG(1-γ)max0≤i≤n∥εi∥G2.
For any given c0∈(0,1), we let c-2=(1-c0)λminG/2. Then for (ϱ/(1-γ))h≤c-2, the above inequality leads to
(3.28)∥ε^n+1∥G2≤11-(2ϱh/λminG(1-γ))[(1+2c2d4hλminG)∥εn∥G2+2c2d4hλminG∥εn+1∥G2+2c2βk2h∥en∥2+2ϱh51-γ(1+|c1|)2M22]+2ϱhλminG(1-γ)-2ρhmax0≤i≤n∥εi∥G2≤(1+d5h)max0≤i≤n∥εi∥G2+2c2d4hc0λminG∥εn+1∥G2+2c2βk2hc0∥en∥2+2ϱh5c0(1-γ)(1+|c1|)2M22,
where d5=2c2d4/c0λminG+4ϱ/c0λminG(1-γ). Noting that
(3.29)∥εn+1∥G2≤∥ε^n+1∥G2+λmaxG∥en∥2+2λmaxG∥en∥∥ε^n+1∥G≤(1+h)∥ε^n+1∥G2+(1+1h)λmaxG∥en∥2,
where λmaxG denotes the maximum eigenvalues of the matrix G, we have
(3.30)∥εn+1∥G2≤(1+h)[(1+d5h)max0≤i≤n∥εi∥G2+2c2d4hc0λminG∥εn+1∥G2+2c2βk2hc0∥en∥2+2ϱh5c0(1-γ)(1+|c1|)2M22]+(1+1h)λmaxG∥en∥2.
Now for any given c~0∈(0,1), when (4c2d4/c0λminG)h≤1-c~0, the above inequality leads to (3.13) with
(3.31)h1=min{1,(1-γ)c-2ϱ,(1-c~0)c0λminG4c2d4},d1=1+2d5c~0+4c2d4c~0c0λminG,d2=4ϱ(1+|c1|)2M22(1-γ)c~0c0,d3=4c2βk2c~0c0+2λmaxGc~0.
This completes the proof of this lemma.Compared to the nonneutral DDEs with a constant delay considered in [38, 44], the problems considered in this paper are more complex such that a series of new difficulties need to be overcome. Nevertheless, we note that some basic proof ideas are related to the ones used in [38, 44].Lemma 3.5.
If the method(ρ,σ) is A-stable, then there exist two constants d6 and h2, which depend only on the method, some of the bounds Mi, and the parameters α and ϱ/(1-γ), such that
(3.32)∥en∥≤d6hp+1,h∈(0,h2],n=0,1,…,
where p is the consistent order in the classical sense, p=1,2.Proof.
The idea is a generalization of that used in [38, 44]. The A-stability of the method (ρ,σ) implies that βk/αk>0, and it is consistent of order p=1,2 in the classical sense (see [45, 46]). Consider the following scheme:
(3.33)y(σ(E)tn)=∑j=0k-1‍(βj-βkαkαj)y^n+j+βkαkhy′(σ(E)tn)+R1(n),(3.34)ρ(E)y^n=hy′(σ(E)tn)+R2(n).
By Taylor expansion, we can assert that there exists constant c3 such that
(3.35)R1(n)≤c3h3,R2(n)≤c3hp+1.
From (3.7) and (3.33), it is easy to obtain
(3.36)y(σ(E)tn)-σ(E)y^n-βken=R1(n)+βkαkh[y′(σ(E)tn)-f(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))].
If η(i)(σ(E)tn)∈(tn+k-1,tn+k], by Taylor expansion, we have
(3.37)∥y(η(i)(σ(E)tn))-δn(i)(y^n+k+en)-(1-δn(i))y(tn+k-1)∥≤c4h2M2+∥en∥.
This inequality and (3.36) imply that
(3.38)∥y(σ(E)tn)-σ(E)y^n-βken∥2≤βkαkh∥y(σ(E)tn)-σ(E)y^n-βken∥×[α∥y(σ(E)tn)-σ(E)y^n-βken∥+ϱ1-γ(c4h2M2+∥en∥)]+R1(n)∥y(σ(E)tn)-σ(E)y^n-βken∥.
For any given c-0∈(0,1), we let c-3=(1-c-0)αk/βk. Then for αh≤c-3, the above inequality leads to
(3.39)∥y(σ(E)tn)-σ(E)y^n-βken∥≤1c-0[βkϱhαk(1-γ)(c4h2M2+∥en∥)+∥R1(n)∥]≤c5h3+βkϱhαkc-0(1-γ)∥en∥,
where
(3.40)c5=1c-0[βkϱc4M2αk(1-γ)+c3].
Substitute (3.39) into (3.36) to obtain
(3.41)∥hy′(σ(E)tn)-hf(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))∥≤αkβk∥y(σ(E)tn)-σ(E)y^n-βken∥+αkβk∥R1(n)∥≤αkβk(c3+c5)h3+ϱhc-0(1-γ)∥en∥.
On the other hand, since
(3.42)αken=hf(σ(E)tn,σ(E)y^n+βken,y-h(η(σ(E)tn)),Y-h(η(σ(E)tn)))-hy′(σ(E)tn)-R2(n)
it follows from (3.41) that
(3.43)|αk|∥en∥≤αkβk(c3+c5)h3+ϱhc-0(1-γ)∥en∥+c3hp+1.
For any given c^0∈(0,1), let
(3.44)h2={min{h1,(1-γ)(1-c^0)c-0|αk|ϱ},ifα≤0,min{h1,(1-γ)(1-c^0)c-0|αk|ϱ,c-3α},ifα>0.
Then the required inequality (3.32) follows from (3.43), where
(3.45)d6=1c^0(c3+c5|βk|+c3|αk|).
This completes the proof of the lemma.Lemma3.4, together with Lemma 3.5, implies the following theorem.Theorem 3.6.
The algorithm OLIDE(I) extended by the method(ρ,σ) is EB-convergent of order p if and only if the method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2.Proof.
Firstly, we observe that theEB-convergence of order p of the method OLIDE(I) implies that the method (ρ,σ) is B-convergent of order p, and it is A-stable and consistent of order p in the classical sense for ODEs, p=1 or 2 (see, e.g., [38]).
On the other hand, it follows from Lemmas3.4 and 3.5 that
(3.46)∥εn+1∥G2≤(1+d1h)max0≤i≤n∥εi∥G2+d2h5+d3d6h2p+1,n=0,1,…,h∈(0,h2]
which implies that
(3.47)max0≤i≤n+1∥εn+1∥G2≤(1+d1h)max0≤i≤n∥εi∥G2+d2h5+d3d6h2p+1,n=0,1,…,h∈(0,h2].
By induction, we have
(3.48)∥εn∥G2≤max0≤i≤n∥εi∥G2≤(1+d1h)n∥ε0∥G2+(d2h5+d3d6h2p+1)[1+(1+d1h)+⋯+(1+d1h)n-1]≤(1+d1h)n[∥ε0∥G2+n(d2h5+d3d6h2p+1)]≤exp(d1tn)[∥ε0∥G2+tn(d2+d3d6)h2p]
and therefore
(3.49)∥yn-y^n∥2≤exp(d1tn)λminG[λmaxG∑j=0k-1‍∥yj-y^j∥2+tn(d2+d3d6)h2p].
Consequently, the following holds:
(3.50)∥yn-y(tn)∥=∥yn-y^n+y^n-y(tn)∥≤|c1|M2h2+[kλmaxGλminG(max0≤i≤k-1∥yi-y(ti)∥+|c1|M2h2)+tn(d2+d3d6)λminGhp]×exp(12d1tn),n=0,1,…,h∈(0,h2]
which implies the algorithm OLIDE(I) is EB-convergent of order p,p=1 or 2. The proof of Theorem 3.6 is completed.In [44], Huang et al. proved that A-stable one-leg methods with the interpolation (2.7) for constant-delay DDEs are D-convergent under a condition that η(i)(σ(E)tn)≤tn+k-1. As a special case, we have the following corollary.Corollary 3.7.
A one-leg method(ρ,σ) with the interpolation (2.7) for DDEs with any variable delay is D-convergent of order p if and only if the method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2.


## 3.3. Convergence Analysis of OLIDE(II)

In [18], Wang et al. have shown that the algorithm OLIDE(II) has better stability properties than the algorithm OLIDE(I) for DDEs with constant delay. Here we will prove that the algorithms OLIDE(II) and OLIDE(I) have the same convergence properties. To obtain the convergence result of the algorithm OLIDE(II), we also consider (3.7). But in this case, y-h(η(i)(σ(E)tn))  (i=1,2,…) are determined by
(3.51)y-h(η(i)(σ(E)tn))={βk(y^n+k+en)+∑j=0k-1‍βjy(tn+j),iftn+k-1<η(i)(σ(E)tn)≤tn+k,m-n(i)=0,δ-n(i)[βk(y^n+k+en)+∑j=0k-1‍βjy(tn+j)]+(1-δ-n(i))∑j=0k‍βjy(tn-1+j),iftn+k-1<η(i)(σ(E)tn)≤tn+k,m-n(i)=1,y(η(i)(σ(E)tn)),othercase.Whenη(i)(σ(E)tn)≤tn+k-1, we have
(3.52)∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥=∥(1-δ-n(i))σ(E)yn-m-n(i)+δ-n(i)σ(E)yn-m-n(i)+1-y(η(i)(σ(E)tn))∥≤(1-δ-n(i))∥σ(E)(yn-m-n(i)-y^n-m-n(i))∥+δ-n(i)∥σ(E)(yn-m-n(i)+1-y^n-mn(i)+1)∥+∥(1-δ-n(i))σ(E)y^n-m-n(i)+δ-n(i)σ(E)y^n-m-n(i)+1-y(η(i)(σ(E)tn))∥≤(1-δ-n(i))∑j=0k‍|βj|∥yn-m-n(i)+j-y^n-m-n(i)+j∥+δ-n(i)∑j=0k‍|βj|∥yn-m-n(i)+j+1-y^n-m-n(i)+j+1∥+∥(1-δ-n(i))∑j=0k‍βjY^n-m-n(i)+j+δ-n(i)∑j=0k‍βjY^n-m-n(i)+j+1-y′(η(i)(σ(E)tn))∥,
where
(3.53)Y^i=f(ti,y^i,y-h(η(ti)),Y-h(η(ti))),i=1,2,…,n+k-1.
By Taylor expansion, there exist constants c-1 and h3 such that
(3.54)∥yh(η(i)(σ(E)tn))-y-h(η(i)(σ(E)tn))∥≤c-1[max0≤j≤k+1∥yn+j-m-n(i)-y^n+j-m-n(i)∥+h2],h∈(0,h3].
Noting that for tn+k-1≤η(σ(E)tn)≤tn+k we have similar inequality, we find that the proof of the following theorem is similar to that of Theorem 3.6.Theorem 3.8.
The algorithm OLIDE(II) extended by the method(ρ,σ) is EB-convergent of order p if and only if the method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2.It should be pointed out that Huang et al. [38] proved that A-stable one-leg methods with the interpolation (2.8) for constant-delay DDEs are D-convergent. As a special case, in this paper we obtain the convergence result of A-stable one-leg methods with the interpolation (2.8) for DDEs with any variable delay.Since the algorithm OLIDE(II) has better stability properties than the algorithm OLIDE(I) [18] and they have the same convergence properties, we believe that the algorithm OLIDE(II) is more effective than the algorithm OLIDE(I) for solving delay problem.Although we can show theoretically the convergence of the algorithms OLIDE for solving NDDEs with any variable delay, in practical implementation, these algorithms are generally expensive and will create a serious storage problem since they always require to trace back the recursion until the initial interval. So for general variable delay NDDEs, we need other algorithms, which are based on interpolation.


## 4. The Second Class of Algorithms Derived from the Interpolation: OLIIT

In this section, forYh(η(σ(E)tn)) we consider the consistent interpolation schemes with for yh(η(σ(E)tn)). Firstly, if (2.6) holds, the consistent interpolation with (2.7) is
(4.1)Yh(η(i)(σ(E)tn))={δn(i)Ymn(i)+1+(1-δn(i))Ymn(i),tmn(i)+δn(i)h>0,ϕ′(tmn(i)+δn(i)h),tmn(i)+δn(i)h≤0,
where Yj is computed by the following formula:
(4.2)Yj=f(tj,yj,yh(η(tj)),Yh(η(tj))).
Of course, when j≤0, Yj is determined by Yj=ϕ′(jh). The algorithm (2.4)–(2.7)–(4.1) is simply denoted by OLIIT(I). Obviously, when the methods OLIDE(I) and OLIIT(I) are applied to non-neutral DDEs, they are identical and have been used in [11].Given (2.9), the corresponding interpolation scheme with (2.8) is
(4.3)Yh(η(i)(σ(E)tj))=δ-j(i)σ(E)Yj-m-j(i)+1+(1-δ-j(i))σ(E)Yj-m-j(i).
Similarly, Yj is produced by (4.2) and Yj=ϕ′(jh) for j≤0. The algorithm (2.4)–(2.8)–(4.3) is denoted by OLIIT(II).Now we discuss the convergence of OLIIT(I) and OLIIT(II). In order to do this, the condition (3.3) will be replaced by the following:
(4.4)∥f(t,y1,u,v)-f(t,y2,u,v)∥≤L∥y1-y2∥.
If we allow the error bound to depend on (β+γL)/(1-γ), we can also obtain the convergence results.


### 4.1. Convergence Analysis of OLIIT(I)

Below, we establish the convergence result of the algorithm OLIIT(I). To accomplish this, we need to consider (3.7) with Y-h(η(i)(σ(E)tn))   defined by the following formula:
(4.5)Y-h(η(i)(σ(E)tn))={y′(η(i)(σ(E)tn)),ifη(i)(σ(E)tn)≤tn+k-1,δn(i)Y^n+k+(1-δn(i))y′(tn+k-1),iftn+k-1<η(i)(σ(E)tn)≤tn+k,
where
(4.6)Y^n+k=f(tn+k,y^n+k+en,y-h(η(tn+k)),Y-h(η(tn+k))).
Then we have the following lemma.Lemma 4.1.
If the method(ρ,σ) is A-stable, then the numerical solution produced by the algorithm OLIIT(I) satisfies
(4.7)∥εn+1∥G2≤(1+d1h)max0≤i≤n∥εi∥G2+d2h5+d3h-1∥en∥2,n=0,1,…,h∈(0,h1].
Here the constants d1, d2, d3, and h1 depend on the method (ρ,σ), α, (β+γL)/(1-γ), and bounds Mi for certain derivatives of the true solution y(t).Proof.
Similar to the proof of Lemma3.4, we have
(4.8)∥ε^n+1∥G2≤∥εn∥G2+2h[∥Yh(η(σ(E)tn))-Y-h(η(σ(E)tn))∥α∥σ(E)(yn-y^n)-βken∥2+∥σ(E)(yn-y^n)-βken∥×(∥Yh(η(σ(E)tn))-Y-h(η(σ(E)tn))∥β∥yh(η(σ(E)tn))-y-h(η(σ(E)tn))∥+γ∥Yh(η(σ(E)tn))-Y-h(η(σ(E)tn))∥)].
If η(σ(E)tn)≤tn+k-1, it follows from (2.7) and (4.5) that
(4.9)β∥yh(η(σ(E)tn))-y-h(η(σ(E)tn))∥+γ∥Yh(η(σ(E)tn))-Y-h(η(σ(E)tn))∥=β∥(1-δn(1))ymn(1)+δn(1)ymn(1)+1-y(η(σ(E)tn))∥+γ∥(1-δn(1))Ymn(1)+δn(1)Ymn(1)+1-y′(η(σ(E)tn))∥≤β(1-δn(1))∥ymn(1)-y^mn(1)∥+βδn(1)∥ymn(1)+1-y^mn(1)+1∥+β(1+|c1|)M2h2+γ(1-δn(1))∥Ymn(1)-Y^mn(1)∥+γδn(1)∥Ymn(1)+1-Y^mn(1)+1∥+γL(1+|c1|)M2h2≤(β+γL)[(1-δn(1))∥ymn(1)-y^mn(1)∥+δn(1)∥ymn(1)+1-y^mn(1)+1∥+(1+|c1|)M2h2]+γ(β∥yh(η(tmn(1)))-y-h(η(tmn(1)))∥+γ∥Yh(η(tmn(1)))-Y-h(η(tmn(1)))∥)+γ(β∥yh(η(tmn(1)+1))-y-h(η(tmn(1)+1))∥+γ∥Yh(η(tmn(1)+1))-Y-h(η(tmn(1)+1))∥)≤β+γL1-γmax1≤i≤n+k-1∥yi-y^i∥+β+γL1-γ(1+|c1|)M2h2.
If tn+k-1≤η(σ(E)tn)≤tn+k, then from (2.7) and (4.5) we have
(4.10)β∥yh(η(σ(E)tn))-y-h(η(σ(E)tn))∥+γ∥Yh(η(σ(E)tn))-Y-h(η(σ(E)tn))∥=β(1-δn(1))∥yn+k-1-y-h(tn+k-1)∥+βδn(1)∥yn+k-y^n+k-en∥+γ(1-δn(1))∥Yn+k-1-Y-h(tn+k-1)∥+γδn(1)∥Yn+k-Y^n+k∥≤β(1-δn(1))∥yn+k-1-y^n+k-1∥+β(1+|c1|)M2h2+βδn(1)∥yn+k-y^n+k-en∥+γ(1-δn(1))∥Yn+k-1-Y^n+k-1∥+γL(1+|c1|)M2h2+γδn(1)∥Yn+k-Y^n+k∥≤β+γL1-γmax{max1≤i≤n+k-1∥yi-y^i∥,∥yn+k-y^n+k-en∥}+β+γL1-γ(1+|c1|)M2h2.
Consequently, we have a similar inequality to (3.20). The remaining part of this proof is analogous with that of Lemma 3.4, and we omit it here. This completes the proof.Similarly, we can give the same lemma as Lemma3.5 and obtain the following theorem.Theorem 4.2.
The algorithm OLIIT(I) extended by the method(ρ,σ) is convergent of order p if and only if the method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2.


### 4.2. Convergence Analysis of OLIIT(II)

In this subsection, we establish the convergence result of the algorithm OLIIT(II).Theorem 4.3.
The algorithm OLIIT(II) extended by the method(ρ,σ) is convergent of order p if and only if the method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2.There is no essential difference between the proofs of the previous three theorems and the proof of this theorem, and therefore we omit it here.We conclude this section with a few remarks about our results.Our first remark is that if we allow the error bound to depend on(β+γL)/(1-γ), the algorithms OLIDE(I) and OLIDE(II) are also convergent. Moreover, since ϱ≤β+γL, we have reason to believe that the algorithms OLIDE(I) and OLIDE(II) whose error bounds depend on α0 and ϱ/(1-γ) are more efficient than the algorithms OLIIT(I) and OLIIT(II) whose error bounds depend on α0 and (β+γL)/(1-γ).Second, a key property of NDDEs is that its solution does not possess the smoothing property that the solution becomes smoother with increasing values oft (see, [1]), that is, if
(4.11)ϕ′(t0)=f(t0,ϕ(t0),ϕ(η(t0)),ϕ′(η(t0)))
does not hold, the solution to NFDEs will in general have a discontinuous derivative at a finite set of points of discontinuity of the first kind. Our results cannot be applied to this case, but they would provide better insight to numerical researchers in the field of numerical methods for solving NDDEs. On the other hand, for the problem which satisfies the condition (𝒞1), the convergence analysis of the four algorithms can be done interval-by-interval if its true solution is sufficiently differentiable on each subinterval Ik.


## 4.1. Convergence Analysis of OLIIT(I)

Below, we establish the convergence result of the algorithm OLIIT(I). To accomplish this, we need to consider (3.7) with Y-h(η(i)(σ(E)tn))   defined by the following formula:
(4.5)Y-h(η(i)(σ(E)tn))={y′(η(i)(σ(E)tn)),ifη(i)(σ(E)tn)≤tn+k-1,δn(i)Y^n+k+(1-δn(i))y′(tn+k-1),iftn+k-1<η(i)(σ(E)tn)≤tn+k,
where
(4.6)Y^n+k=f(tn+k,y^n+k+en,y-h(η(tn+k)),Y-h(η(tn+k))).
Then we have the following lemma.Lemma 4.1.
If the method(ρ,σ) is A-stable, then the numerical solution produced by the algorithm OLIIT(I) satisfies
(4.7)∥εn+1∥G2≤(1+d1h)max0≤i≤n∥εi∥G2+d2h5+d3h-1∥en∥2,n=0,1,…,h∈(0,h1].
Here the constants d1, d2, d3, and h1 depend on the method (ρ,σ), α, (β+γL)/(1-γ), and bounds Mi for certain derivatives of the true solution y(t).Proof.
Similar to the proof of Lemma3.4, we have
(4.8)∥ε^n+1∥G2≤∥εn∥G2+2h[∥Yh(η(σ(E)tn))-Y-h(η(σ(E)tn))∥α∥σ(E)(yn-y^n)-βken∥2+∥σ(E)(yn-y^n)-βken∥×(∥Yh(η(σ(E)tn))-Y-h(η(σ(E)tn))∥β∥yh(η(σ(E)tn))-y-h(η(σ(E)tn))∥+γ∥Yh(η(σ(E)tn))-Y-h(η(σ(E)tn))∥)].
If η(σ(E)tn)≤tn+k-1, it follows from (2.7) and (4.5) that
(4.9)β∥yh(η(σ(E)tn))-y-h(η(σ(E)tn))∥+γ∥Yh(η(σ(E)tn))-Y-h(η(σ(E)tn))∥=β∥(1-δn(1))ymn(1)+δn(1)ymn(1)+1-y(η(σ(E)tn))∥+γ∥(1-δn(1))Ymn(1)+δn(1)Ymn(1)+1-y′(η(σ(E)tn))∥≤β(1-δn(1))∥ymn(1)-y^mn(1)∥+βδn(1)∥ymn(1)+1-y^mn(1)+1∥+β(1+|c1|)M2h2+γ(1-δn(1))∥Ymn(1)-Y^mn(1)∥+γδn(1)∥Ymn(1)+1-Y^mn(1)+1∥+γL(1+|c1|)M2h2≤(β+γL)[(1-δn(1))∥ymn(1)-y^mn(1)∥+δn(1)∥ymn(1)+1-y^mn(1)+1∥+(1+|c1|)M2h2]+γ(β∥yh(η(tmn(1)))-y-h(η(tmn(1)))∥+γ∥Yh(η(tmn(1)))-Y-h(η(tmn(1)))∥)+γ(β∥yh(η(tmn(1)+1))-y-h(η(tmn(1)+1))∥+γ∥Yh(η(tmn(1)+1))-Y-h(η(tmn(1)+1))∥)≤β+γL1-γmax1≤i≤n+k-1∥yi-y^i∥+β+γL1-γ(1+|c1|)M2h2.
If tn+k-1≤η(σ(E)tn)≤tn+k, then from (2.7) and (4.5) we have
(4.10)β∥yh(η(σ(E)tn))-y-h(η(σ(E)tn))∥+γ∥Yh(η(σ(E)tn))-Y-h(η(σ(E)tn))∥=β(1-δn(1))∥yn+k-1-y-h(tn+k-1)∥+βδn(1)∥yn+k-y^n+k-en∥+γ(1-δn(1))∥Yn+k-1-Y-h(tn+k-1)∥+γδn(1)∥Yn+k-Y^n+k∥≤β(1-δn(1))∥yn+k-1-y^n+k-1∥+β(1+|c1|)M2h2+βδn(1)∥yn+k-y^n+k-en∥+γ(1-δn(1))∥Yn+k-1-Y^n+k-1∥+γL(1+|c1|)M2h2+γδn(1)∥Yn+k-Y^n+k∥≤β+γL1-γmax{max1≤i≤n+k-1∥yi-y^i∥,∥yn+k-y^n+k-en∥}+β+γL1-γ(1+|c1|)M2h2.
Consequently, we have a similar inequality to (3.20). The remaining part of this proof is analogous with that of Lemma 3.4, and we omit it here. This completes the proof.Similarly, we can give the same lemma as Lemma3.5 and obtain the following theorem.Theorem 4.2.
The algorithm OLIIT(I) extended by the method(ρ,σ) is convergent of order p if and only if the method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2.


## 4.2. Convergence Analysis of OLIIT(II)

In this subsection, we establish the convergence result of the algorithm OLIIT(II).Theorem 4.3.
The algorithm OLIIT(II) extended by the method(ρ,σ) is convergent of order p if and only if the method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2.There is no essential difference between the proofs of the previous three theorems and the proof of this theorem, and therefore we omit it here.We conclude this section with a few remarks about our results.Our first remark is that if we allow the error bound to depend on(β+γL)/(1-γ), the algorithms OLIDE(I) and OLIDE(II) are also convergent. Moreover, since ϱ≤β+γL, we have reason to believe that the algorithms OLIDE(I) and OLIDE(II) whose error bounds depend on α0 and ϱ/(1-γ) are more efficient than the algorithms OLIIT(I) and OLIIT(II) whose error bounds depend on α0 and (β+γL)/(1-γ).Second, a key property of NDDEs is that its solution does not possess the smoothing property that the solution becomes smoother with increasing values oft (see, [1]), that is, if
(4.11)ϕ′(t0)=f(t0,ϕ(t0),ϕ(η(t0)),ϕ′(η(t0)))
does not hold, the solution to NFDEs will in general have a discontinuous derivative at a finite set of points of discontinuity of the first kind. Our results cannot be applied to this case, but they would provide better insight to numerical researchers in the field of numerical methods for solving NDDEs. On the other hand, for the problem which satisfies the condition (𝒞1), the convergence analysis of the four algorithms can be done interval-by-interval if its true solution is sufficiently differentiable on each subinterval Ik.


## 5. Numerical Experiments

In this section in order to illustrate the one-leg methods (2.4) for solving the NDDEs (1.1)-(1.2) we will consider two examples.


### 5.1. Example 1: A Test Equation

First of all, we consider a simple test equation(5.1)y′(t)=a[y(t)-g(t)]+g′(t)+b[y(t-τ)-g(t-τ)]+c[y′(t-τ)-g′(t-τ)],t∈[0,10],y(t)=g(t),t∈[-τ,0]
which is a modification of a test equation for stiff ODEs (see, e.g., [47, 48]), where τ=1, g:R→R is a given smooth function and a,b, and c are given real constants. Its exact solution is y(t)=g(t). Now let us take a=-108,b=0.9·108,c=0.9,andg(t)=sin(t). Then α=-108, β=0.9·108, γ=0.9, and ϱ=0. Observe that when the constant coefficient linear NDDE (5.1) is discretized on a constrained mesh, that is, h=0.1,h=0.01,andh=0.001, the four algorithms, OLIDE(I), OLIIT(I), OLIDE(II), and OLIIT(II), extended by the midpoint rule are identical. The same conclusion holds for the algorithms extended by second-order backward differentiation formula (BDF2). So we consider only applying the algorithms OLIDE(I) extended by the midpoint rule (OLIDE(I)-MP)
(5.2)yn+1=yn+hf(tn+12h,yn+yn+12),n=0,1,…
and by second-order backward differentiation formula (OLIDE(I)-BDF2)
(5.3)yn+2=43yn+1-13yn+23hf(tn+2,yn+2),n=0,1,…
to the problem (5.1). Table 1 shows the numerical errors at t=10.Table 1
The numerical errors att=10 when the algorithm OLIDE(I) is applied to the problem (5.1).
h
OLIDE(I)-MP
OLIDE(I)-BDF2
0.1
6.807355
e
-
004
2.922906
e
-
011
0.01
6.800335
e
-
006
2.811085
e
-
013
0.001
6.800264
e
-
008
2.775558
e
-
015Note that for this problem, the two algorithms both are convergent of order2 and the algorithm OLIDE(I)-BDF2 is more effective than the algorithm OLIDE(I)-MP.


### 5.2. Example 2: A Partial Functional Differential Equation

The next problem is defined fort∈[0,10] and x∈[0,1], as follows:
(5.4)∂∂tu(t,x)=a(t)∂2∂x2u(t,x)+b(t)u(η(t),x)+c(t)∂∂tu(η(t),x)+g(t,x).
The function g(t,x) and the initial and boundary values are selected in such a way that the exact solution becomes
(5.5)u(t,x)=(x-x2)exp(-t).
After application of the numerical method of lines, we obtain the following neutral-delay differential equations of the form
(5.6)ui′(t)=a(t)ui-1(t)-2ui(t)+ui+1(t)Δx2+b(t)ui(η(t))+c(t)ui′(η(t))+gi(t),t∈[0,10],u0(t)=uNx(t)=0,t∈[0,10],ui(t)=iΔx(1-iΔx)exp(-t),t≤0,i=1,2,…,Nx-1,
where Δx is the spatial step, Nx is a natural number such that ΔxNx=1, xi=iΔx,i=0,1,2,…,Nx, gi(t)=g(t,xi), and ui(t) is meant to approximate the solution of (5.4) at the point (t,xi).We takeΔx=0.1 or Δx=0.01 for the numerical method of lines and use the midpoint rule connecting the different interpolation approximations for the numerical integration of the problem (5.6).For the purpose of comparison, we consider three groups of coefficients as follows:Coefficient I:a(t)=sin2t, b(t)=0, c(t)=-0.0001;Coefficient II:a(t)=sin2t, b(t)=0, c(t)=-0.9;Coefficient III:a(t)=100sin2t, b(t)=0, c(t)=-0.9. It is easy to verify that α=-π2mint∈[0,10]a(t), L=4Nx2maxt∈[0,10]a(t), β=maxt∈[0,10]b(t), γ=maxt∈[0,10]|c(t)|, and ϱ=4Nx2maxt∈[0,10]{a(t)|c(t)|}. Then only for coefficients I, ϱ/(1-γ) or (β+γL)/(1-γ) is of moderate size.Let(5.7)ϵ=max1≤i≤Nx-1|ui(10)-u(10,iΔx)|
denote the error of an algorithm when applied to problem (5.6).(a) Constant Delay Problem
First of all, we consider the problem (5.6) with a constant delay τ(t)=t-η(t)=1. When we choose the step-sizes h=0.1, h=0.01, and h=0.001, the two algorithms OLIDE(I) and OLIDE(II) extended by the midpoint rule are identical (OLIDE-MP), and the two algorithms OLIIT(I) and OLIIT(II) extended by the midpoint rule are identical (OLIIT-MP). But since (5.6) is a variable-coefficient system, the algorithm OLIDE-MP differs from the algorithm OLIIT-MP. The errors ϵ at t=10 are listed in Tables 2, 3, and 4 when the methods are applied to the problem (5.6) with three groups of coefficients, respectively.
Observe that for the constant delay problem with the coefficients I, both algorithms, OLIDE-MP and OLIIT-MP, are convergent of order 2. But for the problem with the coefficients II, the numerical results of OLIIT-MP is not ideal whenNx=100 and h=0.1. This situation has further deteriorated when the coefficients become III. On the one hand, this implies that OLIDE-MP are stronger than OLIIT-MP, which confirm our theoretical analysis. On the other hand, this implies that the coefficients, the spatial step-size Δx and the time step size h affect the efficiency of the algorithm OLIIT-MP. It is well-known that the midpoint rule is A-stable and is convergent for stiff ODEs. But now the numerical results become worse when the time step-size is larger and the grid is finer, which further confirms our theoretical results.Table 2
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group I, where τ(t)=1.
h
OLIDE-MP
OLIIT-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
2.804687
e
-
008
2.795105
e
-
008
2.808233
e
-
008
2.798658
e
-
008
0.01
2.802893
e
-
010
2.793375
e
-
010
2.806421
e
-
010
2.796910
e
-
010
0.001
2.802876
e
-
012
2.793359
e
-
012
2.806405
e
-
012
2.796895
e
-
012Table 3
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group II, where τ(t)=1.
h
OLIDE-MP
OLIIT-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
3.500122
e
-
007
3.433799
e
-
007
3.345386
e
-
006
5.307093
e
-
002
0.01
2.658399
e
-
009
2.575017
e
-
009
3.215464
e
-
008
3.198151
e
-
008
0.001
2.649853
e
-
011
2.566321
e
-
011
3.214185
e
-
010
3.196886
e
-
010Table 4
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group III, where τ(t)=1.
h
OLIDE-MP
OLIIT-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
1.932937
e
-
006
2.085638
e
-
006
5.303484
e
+
000
5.282569
e
+
001
0.01
1.432539
e
-
010
1.432852
e
-
010
2.926272
e
-
010
9.189176
e
-
004
0.001
1.429940
e
-
012
1.429875
e
-
012
2.938189
e
-
012
2.941581
e
-
012(b) Bounded Variable Delay Problem
Next we consider the problem (5.6) with a bounded variable delay τ(t)=t-η(t)=sint+1. Observe that the function η(t)=t-sin(t)-1 satisfies the condition (𝒞1) such that the convergence analysis and numerically solving this equation can be done interval-by-interval. Because we have known that the true solution is sufficiently differentiable on the whole interval, the step sizes h=0.1, h=0.01, and h=0.001 are chosen.In this case, we do not consider OLIDE(I) and OLIDE(II) since the two algorithms will produce implementation and computational complexity issues. We explore only the algorithm OLIIT(I) extended by the midpoint rule (OLIIT(I)-MP) and the algorithm OLIIT(II) extended by the midpoint rule (OLIIT(II)-MP). The errorsϵ at t=10 are listed in Tables 5, 6, and 7 when the algorithms are applied to the problem (5.6) with three groups of coefficients, respectively.Table 5
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group I, where τ(t)=sint+1.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
2.809033
e
-
008
2.799451
e
-
008
2.808259
e
-
008
2.798684
e
-
008
0.01
2.807181
e
-
010
2.797665
e
-
010
2.806578
e
-
010
2.797068
e
-
010
0.001
2.807207
e
-
012
2.797692
e
-
012
2.806491
e
-
012
2.796982
e
-
012Table 6
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients  group II, where τ(t)=sint+1.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
1.126902
e
-
006
9.145590
e
-
004
4.938563
e
-
006
1.656698
e
-
002
0.01
1.106792
e
-
008
1.102020
e
-
008
6.263681
e
-
009
6.274124
e
-
009
0.001
1.142587
e
-
010
1.136783
e
-
010
6.251279
e
-
011
6.255985
e
-
011Table 7
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group III, where τ(t)=sint+1.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
9.066809
e
-
002
1.199248
e
+
000
1.651731
e
+
000
3.112600
e
+
001
0.01
2.328686
e
-
009
6.123243
e
-
005
2.292876
e
-
009
2.828239
e
-
004
0.001
2.261214
e
-
011
4.619492
e
-
009
2.239975
e
-
011
2.251171
e
-
011From these numerical data, we also see that both algorithms, OLIIT(I)-MP and OLIIT(II)-MP, are convergent of order 2 for this equation with the coefficients I, but such a good situation is destroyed when the coefficients become II or III.(c) Proportional Delay Problem
Finally, we consider the problem (5.6) with a proportional variable delay τ(t)=0.5t. We still choose the step size h=0.1, h=0.01, and h=0.001. Similar to the case of the bounded variable delay, we use only the algorithm OLIIT(I)-MP and the algorithm OLIIT(II)-MP to solve the problem (5.6). Tables 8, 9, and 10 show the errors ϵ at t=10.Table 8
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group I, where τ(t)=0.5t.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
3.007486
e
-
008
2.998332
e
-
008
2.985691
e
-
008
2.976663
e
-
008
0.01
3.005176
e
-
010
2.996115
e
-
010
2.983553
e
-
010
2.974621
e
-
010
0.001
3.005155
e
-
012
2.996096
e
-
012
2.983534
e
-
012
2.974603
e
-
012Table 9
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group II, where τ(t)=0.5t.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
3.794966
e
-
005
2.539557
e
-
001
1.627824
e
-
005
2.019670
e
-
004
0.01
4.905939
e
-
007
4.870288
e
-
007
1.643376
e
-
007
1.649717
e
-
007
0.001
5.026856
e
-
009
4.977278
e
-
009
1.770753
e
-
009
1.770849
e
-
009Table 10
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group III, where τ(t)=0.5t.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
4.279370
e
+
001
3.962658
e
+
001
2.324247
e
-
002
6.345608
e
-
002
0.01
5.989558
e
-
008
3.347081
e
-
004
6.024260
e
-
008
6.023323
e
-
008
0.001
5.951239
e
-
010
5.951353
e
-
010
5.842018
e
-
011
5.842442
e
-
010From these numerical data, we still observe that both algorithms, OLIIT(I)-MP and OLIIT(II)-MP, are convergent of order 2 for this equation with the coefficients I, but for the same equation with the coefficients II or III, the numerical results become worse when(β+γL)/(1-γ) becomes larger.


## 5.1. Example 1: A Test Equation

First of all, we consider a simple test equation(5.1)y′(t)=a[y(t)-g(t)]+g′(t)+b[y(t-τ)-g(t-τ)]+c[y′(t-τ)-g′(t-τ)],t∈[0,10],y(t)=g(t),t∈[-τ,0]
which is a modification of a test equation for stiff ODEs (see, e.g., [47, 48]), where τ=1, g:R→R is a given smooth function and a,b, and c are given real constants. Its exact solution is y(t)=g(t). Now let us take a=-108,b=0.9·108,c=0.9,andg(t)=sin(t). Then α=-108, β=0.9·108, γ=0.9, and ϱ=0. Observe that when the constant coefficient linear NDDE (5.1) is discretized on a constrained mesh, that is, h=0.1,h=0.01,andh=0.001, the four algorithms, OLIDE(I), OLIIT(I), OLIDE(II), and OLIIT(II), extended by the midpoint rule are identical. The same conclusion holds for the algorithms extended by second-order backward differentiation formula (BDF2). So we consider only applying the algorithms OLIDE(I) extended by the midpoint rule (OLIDE(I)-MP)
(5.2)yn+1=yn+hf(tn+12h,yn+yn+12),n=0,1,…
and by second-order backward differentiation formula (OLIDE(I)-BDF2)
(5.3)yn+2=43yn+1-13yn+23hf(tn+2,yn+2),n=0,1,…
to the problem (5.1). Table 1 shows the numerical errors at t=10.Table 1
The numerical errors att=10 when the algorithm OLIDE(I) is applied to the problem (5.1).
h
OLIDE(I)-MP
OLIDE(I)-BDF2
0.1
6.807355
e
-
004
2.922906
e
-
011
0.01
6.800335
e
-
006
2.811085
e
-
013
0.001
6.800264
e
-
008
2.775558
e
-
015Note that for this problem, the two algorithms both are convergent of order2 and the algorithm OLIDE(I)-BDF2 is more effective than the algorithm OLIDE(I)-MP.


## 5.2. Example 2: A Partial Functional Differential Equation

The next problem is defined fort∈[0,10] and x∈[0,1], as follows:
(5.4)∂∂tu(t,x)=a(t)∂2∂x2u(t,x)+b(t)u(η(t),x)+c(t)∂∂tu(η(t),x)+g(t,x).
The function g(t,x) and the initial and boundary values are selected in such a way that the exact solution becomes
(5.5)u(t,x)=(x-x2)exp(-t).
After application of the numerical method of lines, we obtain the following neutral-delay differential equations of the form
(5.6)ui′(t)=a(t)ui-1(t)-2ui(t)+ui+1(t)Δx2+b(t)ui(η(t))+c(t)ui′(η(t))+gi(t),t∈[0,10],u0(t)=uNx(t)=0,t∈[0,10],ui(t)=iΔx(1-iΔx)exp(-t),t≤0,i=1,2,…,Nx-1,
where Δx is the spatial step, Nx is a natural number such that ΔxNx=1, xi=iΔx,i=0,1,2,…,Nx, gi(t)=g(t,xi), and ui(t) is meant to approximate the solution of (5.4) at the point (t,xi).We takeΔx=0.1 or Δx=0.01 for the numerical method of lines and use the midpoint rule connecting the different interpolation approximations for the numerical integration of the problem (5.6).For the purpose of comparison, we consider three groups of coefficients as follows:Coefficient I:a(t)=sin2t, b(t)=0, c(t)=-0.0001;Coefficient II:a(t)=sin2t, b(t)=0, c(t)=-0.9;Coefficient III:a(t)=100sin2t, b(t)=0, c(t)=-0.9. It is easy to verify that α=-π2mint∈[0,10]a(t), L=4Nx2maxt∈[0,10]a(t), β=maxt∈[0,10]b(t), γ=maxt∈[0,10]|c(t)|, and ϱ=4Nx2maxt∈[0,10]{a(t)|c(t)|}. Then only for coefficients I, ϱ/(1-γ) or (β+γL)/(1-γ) is of moderate size.Let(5.7)ϵ=max1≤i≤Nx-1|ui(10)-u(10,iΔx)|
denote the error of an algorithm when applied to problem (5.6).(a) Constant Delay Problem
First of all, we consider the problem (5.6) with a constant delay τ(t)=t-η(t)=1. When we choose the step-sizes h=0.1, h=0.01, and h=0.001, the two algorithms OLIDE(I) and OLIDE(II) extended by the midpoint rule are identical (OLIDE-MP), and the two algorithms OLIIT(I) and OLIIT(II) extended by the midpoint rule are identical (OLIIT-MP). But since (5.6) is a variable-coefficient system, the algorithm OLIDE-MP differs from the algorithm OLIIT-MP. The errors ϵ at t=10 are listed in Tables 2, 3, and 4 when the methods are applied to the problem (5.6) with three groups of coefficients, respectively.
Observe that for the constant delay problem with the coefficients I, both algorithms, OLIDE-MP and OLIIT-MP, are convergent of order 2. But for the problem with the coefficients II, the numerical results of OLIIT-MP is not ideal whenNx=100 and h=0.1. This situation has further deteriorated when the coefficients become III. On the one hand, this implies that OLIDE-MP are stronger than OLIIT-MP, which confirm our theoretical analysis. On the other hand, this implies that the coefficients, the spatial step-size Δx and the time step size h affect the efficiency of the algorithm OLIIT-MP. It is well-known that the midpoint rule is A-stable and is convergent for stiff ODEs. But now the numerical results become worse when the time step-size is larger and the grid is finer, which further confirms our theoretical results.Table 2
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group I, where τ(t)=1.
h
OLIDE-MP
OLIIT-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
2.804687
e
-
008
2.795105
e
-
008
2.808233
e
-
008
2.798658
e
-
008
0.01
2.802893
e
-
010
2.793375
e
-
010
2.806421
e
-
010
2.796910
e
-
010
0.001
2.802876
e
-
012
2.793359
e
-
012
2.806405
e
-
012
2.796895
e
-
012Table 3
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group II, where τ(t)=1.
h
OLIDE-MP
OLIIT-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
3.500122
e
-
007
3.433799
e
-
007
3.345386
e
-
006
5.307093
e
-
002
0.01
2.658399
e
-
009
2.575017
e
-
009
3.215464
e
-
008
3.198151
e
-
008
0.001
2.649853
e
-
011
2.566321
e
-
011
3.214185
e
-
010
3.196886
e
-
010Table 4
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group III, where τ(t)=1.
h
OLIDE-MP
OLIIT-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
1.932937
e
-
006
2.085638
e
-
006
5.303484
e
+
000
5.282569
e
+
001
0.01
1.432539
e
-
010
1.432852
e
-
010
2.926272
e
-
010
9.189176
e
-
004
0.001
1.429940
e
-
012
1.429875
e
-
012
2.938189
e
-
012
2.941581
e
-
012(b) Bounded Variable Delay Problem
Next we consider the problem (5.6) with a bounded variable delay τ(t)=t-η(t)=sint+1. Observe that the function η(t)=t-sin(t)-1 satisfies the condition (𝒞1) such that the convergence analysis and numerically solving this equation can be done interval-by-interval. Because we have known that the true solution is sufficiently differentiable on the whole interval, the step sizes h=0.1, h=0.01, and h=0.001 are chosen.In this case, we do not consider OLIDE(I) and OLIDE(II) since the two algorithms will produce implementation and computational complexity issues. We explore only the algorithm OLIIT(I) extended by the midpoint rule (OLIIT(I)-MP) and the algorithm OLIIT(II) extended by the midpoint rule (OLIIT(II)-MP). The errorsϵ at t=10 are listed in Tables 5, 6, and 7 when the algorithms are applied to the problem (5.6) with three groups of coefficients, respectively.Table 5
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group I, where τ(t)=sint+1.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
2.809033
e
-
008
2.799451
e
-
008
2.808259
e
-
008
2.798684
e
-
008
0.01
2.807181
e
-
010
2.797665
e
-
010
2.806578
e
-
010
2.797068
e
-
010
0.001
2.807207
e
-
012
2.797692
e
-
012
2.806491
e
-
012
2.796982
e
-
012Table 6
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients  group II, where τ(t)=sint+1.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
1.126902
e
-
006
9.145590
e
-
004
4.938563
e
-
006
1.656698
e
-
002
0.01
1.106792
e
-
008
1.102020
e
-
008
6.263681
e
-
009
6.274124
e
-
009
0.001
1.142587
e
-
010
1.136783
e
-
010
6.251279
e
-
011
6.255985
e
-
011Table 7
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group III, where τ(t)=sint+1.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
9.066809
e
-
002
1.199248
e
+
000
1.651731
e
+
000
3.112600
e
+
001
0.01
2.328686
e
-
009
6.123243
e
-
005
2.292876
e
-
009
2.828239
e
-
004
0.001
2.261214
e
-
011
4.619492
e
-
009
2.239975
e
-
011
2.251171
e
-
011From these numerical data, we also see that both algorithms, OLIIT(I)-MP and OLIIT(II)-MP, are convergent of order 2 for this equation with the coefficients I, but such a good situation is destroyed when the coefficients become II or III.(c) Proportional Delay Problem
Finally, we consider the problem (5.6) with a proportional variable delay τ(t)=0.5t. We still choose the step size h=0.1, h=0.01, and h=0.001. Similar to the case of the bounded variable delay, we use only the algorithm OLIIT(I)-MP and the algorithm OLIIT(II)-MP to solve the problem (5.6). Tables 8, 9, and 10 show the errors ϵ at t=10.Table 8
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group I, where τ(t)=0.5t.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
3.007486
e
-
008
2.998332
e
-
008
2.985691
e
-
008
2.976663
e
-
008
0.01
3.005176
e
-
010
2.996115
e
-
010
2.983553
e
-
010
2.974621
e
-
010
0.001
3.005155
e
-
012
2.996096
e
-
012
2.983534
e
-
012
2.974603
e
-
012Table 9
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group II, where τ(t)=0.5t.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
3.794966
e
-
005
2.539557
e
-
001
1.627824
e
-
005
2.019670
e
-
004
0.01
4.905939
e
-
007
4.870288
e
-
007
1.643376
e
-
007
1.649717
e
-
007
0.001
5.026856
e
-
009
4.977278
e
-
009
1.770753
e
-
009
1.770849
e
-
009Table 10
The errorsϵ at t=10 when the algorithms are applied to the problem (5.6) with the coefficients group III, where τ(t)=0.5t.
h
OLIIT(I)-MP
OLIIT(II)-MP
N
x
=
10
N
x
=
100
N
x
=
10
N
x
=
100
0.1
4.279370
e
+
001
3.962658
e
+
001
2.324247
e
-
002
6.345608
e
-
002
0.01
5.989558
e
-
008
3.347081
e
-
004
6.024260
e
-
008
6.023323
e
-
008
0.001
5.951239
e
-
010
5.951353
e
-
010
5.842018
e
-
011
5.842442
e
-
010From these numerical data, we still observe that both algorithms, OLIIT(I)-MP and OLIIT(II)-MP, are convergent of order 2 for this equation with the coefficients I, but for the same equation with the coefficients II or III, the numerical results become worse when(β+γL)/(1-γ) becomes larger.


## 6. Concluding Remarks

In this paper we introduced four algorithms connecting with an ODEs method to solve nonlinear NDDEs with general variable delay, established their convergence properties and compared their numerical results by means of extensive numerical data. This paper can be regarded as an extension from the nonneutral DDEs with a constant delay [38, 44] to neutral DDEs with general variable delay. Although some basic proof ideas are related to the ones used in [38, 44], the problems considered in this paper are more complex such that some new proof techniques were introduced to overcome a series of new difficulties encountered in theoretical analysis.From theoretical analysis given in Sections3 and 4 and numerical results shown in Section 5, we come to the following remarks:(1)
Ifα0 and ϱ/(1-γ) are of moderate size, the algorithms OLIDE(I) and OLIDE(II) based on A-stable one-leg methods (ρ,σ) are convergent of order p, where p=1,2 is consistent of order in the classical sense. When α0 and (β+γL)/(1-γ) are of moderate size, the four algorithms introduced in this paper are convergent of order p if the one-leg method (ρ,σ) is A-stable and consistent of order p in the classical sense, where p=1,2. But if (β+γL)/(1-γ) is very large, the algorithms OLIIT(I) and OLIIT(II) may produce bad numerical results when the time step size is large even if the ODEs method is A-stable. This revels the difference between numerically solving ODEs and NDDEs.(2)
If using the direct estimation (2.10) does not create implementation or computational complexity problem, we prefer the algorithms OLIDE to the algorithms OLIIT. Furthermore, considering the algorithm OLIDE(II) has better stability properties than the algorithm OLIDE(I) (see [18] and the numerical Example 2 in Section 5.2), it is our belief that the algorithm OLIDE(II) could become an effective numerical method for solving this class of problem, for example, NDDEs with constant delay, if the algorithm is easy to implement. Of course, for general NDDEs with variable delay, we have to use the algorithms OLIIT(I) or OLIIT(II).(3)
The results we established in this paper can be extended in a straightforward way to the case of NDDEs with multiple time delay (1.3)-(1.2).



---

*Source: 989867-2012-09-18.xml*