# Classical Theory of Linear Multistep Methods for Volterra Functional Differential Equations

Authors: Yunfei Li; Shoufu Li
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

Based on the linear multistep methods for ordinary differential equations (ODEs) and the canonical interpolation theory that was presented by Shoufu Li who is exactly the second author of this paper, we propose the linear multistep methods for general Volterra functional differential equations (VFDEs) and build the classical stability, consistency, and convergence theories of the methods. The methods and theories presented in this paper are applicable to nonneutral, nonstiff, and nonlinear initial value problems in ODEs, Volterra delay differential equations (VDDEs), Volterra integro-differential equations (VIDEs), Volterra delay integro-differential equations (VDIDEs), etc. At last, some numerical experiments verify the correctness of our theories.


---

## Body

## 1. Introduction

VFDEs contain many subtypes, such as VDDEs, VIDEs, VDIDEs, etc. Certainly ODEs are also a subtype of them. VFDEs are widely applied in many fields of science and technology (see [1–5] and their references), for which there have been remarkable theoretical and numerical analysis research results. As for VDDEs, please refer to [6–19], as for VIDEs, please refer to [20–23], and as for VDIDEs, please refer to [24–29]. In recent decades, Li [30–35] has carried on systematic research for stiff general VFDEs and the numerical methods for them. In 2014, Li [36] established the classical stability and convergence theories of Runge–Kutta methods for nonstiff, nonlinear general VFDEs. It is well known that in solving ODEs, linear multistep methods have significant advantages in the aspects of format simplification and computation cost, and experts have presented many famous linear multistep methods such as Backward Differentiation Formula (BDF) methods and Adams methods. Based on the linear multistep methods for ODEs and the canonical interpolation theory presented by Li [34], we propose the linear multistep methods for general VFDEs and build the classical stability, consistency, and convergence theories of the methods. The methods and theories presented in this paper are applicable to nonneutral, nonstiff, and nonlinear initial value problems in ODEs, VDDEs, VIDEs, VDIDEs, etc., and are interesting companions to the methods and theories in [36]. Furthermore, the paper may have some value for the study of numerical methods for fractional order VFDEs because many numerical methods for fractional-order differential equations are presented by extending the methods for integer-order differential equations.This paper is organized as follows. In Section2, the linear multistep methods for general VFDEs are introduced. The classical stability of the linear multistep methods is discussed in Section 3. The classical consistency and convergence analyses are carried out in Section 4. Some numerical experiments are carried out in Section 5, which verifies the correctness of the theories presented in this paper.


## 2. Derivation of the Numerical Methods

LetRm be the m-dimensional Euclidean space with standard inner product ⋅,⋅ and the corresponding norm ⋅. For any given closed interval I⊂R, let CmI denote the Banach space consisting of all continuous mappings x:I⟶Rm with the maximum norm x=maxt∈Ixt.Consider VFDEs of the form [30–32, 34, 36](1)y′t=ft,yt,y⋅,t∈0,T,yt=φt,t∈−τ,0,where T∈0,+∞ and τ∈0,+∞ are constants, φt∈Cm−τ,0 is an initial function, and f:0,T×Rm×Cm−τ,T⟶Rm is a given continuous mapping which satisfies the classical Lipschitz conditions [36]:(2)ft,u1,ψ⋅−ft,u2,ψ⋅≤L1u1−u2,∀t∈0,T,u1,u2∈Rm,ψ∈Cm−τ,T,ft,u,ψ1⋅−ft,u,ψ2⋅≤L2maxξ∈−τ,tψ1ξ−ψ2ξ,∀t∈0,T,u∈Rm,ψ1,ψ2∈Cm−τ,T,where L1 and L2 are classical Lipschitz constants. Equations (1) has a unique solution yt [37], and we further make the same assumption as [36] that the system y′t=ft,yt,y⋅t∈0,T in (1) is stable for the initial function φt∈Cm−τ,0, which is to say there exists a function Ct>0t∈0,T which satisfies maxt∈0,TCt has only moderate size such that(3)yt−zt≤Ctmaxξ∈−τ,0φξ−χξ,∀t∈0,T,where zt is the solution of the perturbed problem:(4)z′t=ft,zt,z⋅,t∈0,T,zt=χt,t∈−τ,0.We useDL1,L2 to denote the problem class consisting of all the equations of form (1) which satisfy condition (2) and all the other assumption conditions made above. The problem class DL1,L2 contains nonneutral initial value problems in ODEs, DDEs, IDEs, DIDEs, etc. [36].Combine the lineark-step method:(5)∑i=0kαiyn+i=h∑i=0kβiftn+i,yn+i,for ODEs with piecewise Lagrangian interpolation operators Πnih which are constructed based on canonical interpolation theory [34, 36], and we get the linear k-step method:(6)yniht=Πniht;ψ,y1,y2,…,yn+k,−τ≤t≤tn+i,i=0,1,…,k,∑i=0kαiyn+i=h∑i=0kβiftn+i,yn+i,ynih⋅,for the VFDEs (1), where k≥1, all αi and βi are real constants, αk≠0, α0+β0>0, ρξ=∑i=0kαiξi is called the characteristic polynomial of the method (6), h=T/N is the stepsize, N is a given appropriate positive integer, tn=nhn=0,1,…,N are grid points, ψ∈Cm−τ,0 is an approximation to the initial function φt, y0:=ψ0, yn∈Rm is an approximation to the true value ytn, and the piecewise interpolation operators Πnih: R×Cm−τ,0×Rmn+k⟶Cm−τ,tn+i based on p-degree piecewise Lagrangian interpolation polynomials are defined as follows: if p=0,(7)Πniht;ψ,y1,y2,…,yn+k=ψt,t∈−τ,0,yj,t∈tj−1,tj,j=1,2,…,n+i,and if p≥1,(8)Πniht;ψ,y1,y2,…,yn+k=ψt,t∈−τ,0,Lj,pt,t∈tj−1,tj,j=1,2,…,n+i,with(9)Lj,pt=∑v=0plq+vtyq+v,q=j−j¯,j=1,2,…,n+i,lq+vt=∏s=0,s≠vpt−tq+stq+v−tq+s,q=j−j¯,v=0,1,…,p,j=1,2,…,n+i,where the integer j¯ can be freely determined under the conditions(10)1≤j¯≤j,j+p−n−k≤j¯≤p.According to [34, 36], we know a piecewise interpolation operator Πnih defined by (7) or by (8) and (9) which satisfies the canonical condition(11)max−τ≤t≤tn+iΠniht;ψ,y1,y2,…,yn+k−Πniht;χ,z1,z2,…,zn+k≤cπmaxmax1≤j≤n+kyj−zj,max−τ≤t≤0ψt−χt,∀ψ,χ∈Cm−τ,0,yj,zj∈Rm,j=1,2,…,n+k,where the constant cπ≥1 only depends on the polynomial order p and does not depend on n,h,yj,zj,ψ and χ. In this paper, we always assume that cπ is of moderate size. In fact, Li [34, 36] has given the computational formula of cπ and found cπ, respectively, equals 1,1,1.250, and 1.631 when p, respectively, equals 0,1,2, and 3. We call method (5) as the mother method of method (6).Ifβk=0 in method (6), in order to reduce the amount of calculation, we define(12)yniht=Πniht;ψ,y1,y2,…,yn+k−1,−τ≤t≤tn+i,i=0,1,…,k−1,and now the method is an explicit linear k-step method which is a special case of (6).In particular, we point out that when the lineark-step method (6) for VFDEs starts, on the one hand, like method (5) for ODEs, it needs k starting values, and on the other hand, it may need several grid points as interpolation nodes to construct interpolation operator Πnih. Based on the above statements, in this paper, we always assume the k⌣−1k⌣≥k values y1,y2,…,yk⌣−1, which are called additional starting values, should be provided in advance by other ways, and method (6) starts with n=k⌣−k≥0. By the Banach contraction mapping principle, we conclude that, for any given starting data ψ∈Cm−τ,0,y1,y2,…,yk⌣−1∈Rm, method (6) can uniquely determine the sequence yn,n=k⌣,k⌣+1,…,N when h≤hφ, where hφ can be any number under the condition(13)0<hφ<αkβkL1+cπL2.


## 3. Stability Analysis

Firstly, we list some concepts and elementary facts about scalar linear difference equations as follows [35, 38, 39].Consider thekth order scalar linear difference equation with constant coefficients:(14)∑j=0kajy¯m+j=bm,m=q0,q0+1,…,here and later, q0≥0 is an integer constant, and where k≥1 is an integer constant, every aj is a complex constant, ak≠0, any bm is a complex number depending on m. If a complex number sequence y¯q0,y¯q0+1,…, satisfies equation (14), it is denoted by y¯m, in this paper, and is called a solution of (14). It is easily known that equation (14) has a unique solution when given k initial values y¯q0,y¯q0+1,…,y¯q0+k−1. If every bm is zero, equation (14) is a so-called homogeneous equation:(15)∑j=0kajy¯m+j=0,m=q0,q0+1,….A system ofk linearly independent solutions y¯m1,y¯m2,…,y¯mk of equation (15) is called a fundamental system of (15), and the general solution y¯m of equation (15) can be written by(16)y¯m=∑i=1kdiy¯mi,m=q0,q0+1,…,where each di is an arbitrary complex constant.Call the algebraic equation(17)∑j=0kajξj=0,which corresponds to (14) or (15) the characteristic equation.Proposition 1.
(see [35, 38]). Suppose ξ1,ξ2,…,ξk0k0≤k are all different roots of the characteristic equation (17), and let rj denote the multiple number of ξj. Then, when a0≠0, the k=∑j=1k0rj sequences with the elements y¯m which can be equal to(18)ξjm−q0,m−q0ξjm−q0,…,m−q0m−q0−1…m−q0−rj+2ξjm−q0,j=1,2,…,k0,form a fundamental system of equation (15), and when a0=0, which means equation (17) has zero root and we further suppose ξ1=0 is a r1-tuple root of the characteristic equation (17), the k=∑j=1k0rj sequences with the elements y¯m which can be equal to(19)δm−q0q,q=0,1,…,r1−1,(20)ξjm−q0,m−q0ξjm−q0,…,m−q0m−q0−1…m−q0−rj+2ξjm−q0,j=2,3,…,k0,form a fundamental system of equation (15), where δij is the Kronecker function:(21)δij=0,wheni≠j,1,wheni=j.Remark 1.
Some explanations about the proof of Proposition1 are as follows:(i)
Whena0≠0 and q0=0, P. Henrici has given the proof in the famous book [38].(ii)
Whena0≠0 and q0≠0, from the conclusion in the case that a0≠0 and q0=0, the proof is easily completed.(iii)
Whena0=0 and q0=0, the corresponding part of Proposition 1 is based on the proof [38] in the case that a0≠0 and q0=0 mentioned in the above (i) and some conclusions in the case that a0=0 and q0=0 in [35], and we offer the explanations about the proof as follows. Firstly, ξ1=0 is a r1-tuple root of the characteristic equation (17) which means that a0=a1=⋯=ar1−1=0, and it is easily further known as the r1 sequences with the elements y¯m=δm−q0qq=0,1,…,r1−1 are the solutions of equation (15). Secondly, based on the proof presented by P. Henrici mentioned in the above (i), we can know the sequences with the elements shown in (20) are solutions of equation (15) and that the solutions with the elements shown in (19) and (20) are linearly independent.(iv)
Whena0=0 and q0≠0, from the conclusion in the case that a0=0 and q0=0, the proof is easily completed.Proposition 2.
(see [35]). When the initial values y¯q0=y¯q0+1=⋯=y¯q0+k−1=0, the solution y¯m of equation (14) is given by(22)y¯m=∑s=q0m−kbsym−s−1∗,m=q0,q0+1,…,where each ym∗ satisfies the homogeneous equation:(23)∑j=0kajym+j∗=0,m=0,1,…,y0∗=y1∗=⋯=yk−2∗=0,yk−1∗=1ak.Remark 2.
In this paper, whenq<p, we define ∑i=pqQi≡0 and maxp≤i≤qQi≡0, where each Qi is any a given expression.Remark 3.
By Theorem 5.2 in [38], Proposition 2 is easily shown.Proposition 3.
(see [35]). Suppose y¯mii=1,2,…,k is a fundamental system of equation (15) and ym∗ is the solution of the homogeneous equation (23). Then, the general solution y¯m of equation (14) can be expressed by(24)y¯m=∑i=1kciy¯mi+∑s=q0m−kbsym−s−1∗,m=q0,q0+1,…,where each ci is an arbitrary constant.Proposition 4.
(discrete Bellman inequality). Letμ1,μ2≥0,h>0,ηq0, ηq0+1,…,ηN be a sequence of nonnegative real numbers which satisfy(25)ηm≤μ2+μ1h∑j=q0m−1ηj,m=q0,q0+1,…,N.
Then,(26)ηm≤μ2eμ1m−q0h,m=q0,q0+1,…,N.Definition 1.
We say method (6) is zero-stable if the sequence yk⌣,yk⌣+1,…,yN determined by method (6) applied to equation (1) ∈DL1,L2 with starting data ψ∈Cm−τ,0,y0:=ψ0,y1,y2,…,yk⌣−1∈Rm (the method (6) starts with n=k⌣−k) and the sequence zk⌣,zk⌣+1,…,zN determined by the perturbation equations(27)zniht=∏niht;χ,z1,z2,…,zn+k+γnit,−τ≤t≤tn+i,i=0,1,…,k,∑i=0kαizn+i=h∑i=0kβiftn+i,zn+i,znih⋅+εn+k,n=k⌣−k,k⌣−k+1,…,N−k,with starting data χ∈Cm−τ,0,z0:=χ0,z1,z2,…,zk⌣−1∈Rm, satisfy(28)maxk⌣≤n≤Nzn−yn≤Cmax1≤j≤k⌣−1zj−yj+maxt∈−τ,0χt−ψt+maxk⌣−k≤s≤N−kmax0≤i≤kmaxt∈−τ,ts+iγsit+maxk⌣≤j≤Nεj,where γnit∈Cm−τ,tn+i and εn+k∈Rm are both perturbations, 0<h≤h0, the constant h0>0 depends on method (6), L1 and L2, and the constant C is independent of h.Definition 2.
We say method (6) satisfies the root condition if for all roots of the characteristic polynomial ρξ=∑i=0kαiξi, the moduli of them are all no larger than 1, and among them any root whose modulus equals 1 is simple root.Theorem 1.
Method (6) is zero-stable if and only if it satisfies the root condition.Proof.
Firstly, we prove method (6) satisfies the root condition if it is zero-stable. Consider the ODEs:(29)y′t=ft,yt,t∈0,T,yt=η,t=0,which are a special case of (1), where T∈0,+∞, η∈Rm, and the continuous mapping f:0,T×Rm⟶Rm satisfies the classical Lipschitz condition:(30)ft,u1−ft,u2≤L1u1−u2,∀t∈0,T,u1,u2∈Rm,where L1 is a classical Lipschitz constant. For the ODEs (29), method (6) degrades into its mother method (5). Method (6) is zero-stable implies that its mother method (5) is zero-stable, which means method (6) should satisfy the root condition if it is zero-stable [35, 40].
We assume that method (6) satisfies the root condition, and now we prove it is zero-stable.
Denoteδn=zn−yn,n=0,1,…,N. From (6) and (27), we obtain(31)∑i=0kαiδn+i=b¯n,n=k⌣−k,k⌣−k+1,…,N−k,where(32)b¯n=h∑i=0kβiftn+i,zn+i,znih⋅−ftn+i,yn+i,ynih⋅+εn+k.
Using Proposition1 and considering that all αi are real constants, we can find a real fundamental system [38]ωnii=1,2,…,k of the following scalar homogeneous equation:(33)∑j=0kαjωn+j=0,n=k⌣−k,k⌣−k+1,…,N−k,and by applying Proposition 3 to (31), we obtain(34)δn=∑i=1kciωni+∑s=k⌣−kn−kb¯syn−s−1∗,n=k⌣−k,k⌣−k+1,…,N,where each ci∈Rm is a constant vector and each yn∗ satisfies the scalar homogeneous equation:(35)∑j=0kαjyn+j∗=0,n=0,1,…,y0∗=y1∗=⋯=yk−2∗=0,yk−1∗=1αk.
From the assumption that method (6) satisfies the root condition, we can conclude for all ωni and yn−s−1∗ in (34), there exists a constant M which is only dependent on the mother method (5) such that(36)ωni≤M,yn−s−1∗≤M.
It follows from (34) and (36) that(37)δn≤M∑i=1kci+∑s=k⌣−kn−kb¯s,n=k⌣−k,k⌣−k+1,…,N.
Settingn=k⌣−k,k⌣−k+1,…,k⌣−1 in (34), we thus find(38)ci=∑j=1kW−1ijδk⌣−k−1+j,i=1,2,…,k,where the matrix(39)W=ωk⌣−k1ωk⌣−k2⋯ωk⌣−kkωk⌣−k+11ωk⌣−k+12⋯ωk⌣−k+1k⋯⋯⋯⋯ωk⌣−11ωk⌣−12⋯ωk⌣−1k.
By (38), we know(40)∑i=1kci≤W−1maxk⌣−k≤j≤k⌣−1δj,where W−1=∑i,j=1kW−1ij. Combining (37) with (40), we have(41)δn≤MW−1maxk⌣−k≤j≤k⌣−1δj+M∑s=k⌣−kn−kb¯s,n=k⌣−k,k⌣−k+1,…,N.
It is known from (2), (6), (11), (27), and (32) that(42)b¯s≤h∑i=0kβifts+i,zs+i,zsih⋅−fts+i,ys+i,ysih⋅+hεs+k≤h∑i=0kβifts+i,zs+i,zsih⋅−fts+i,ys+i,zsih⋅+fts+i,ys+i,zsih⋅−fts+i,ys+i,ysih⋅+hεs+k≤h∑i=0kβifts+i,zs+i,zsih⋅−fts+i,ys+i,zsih⋅+fts+i,ys+i,zsih⋅−fts+i,ys+i,ysih⋅+hεs+k≤h∑i=0kβiL1δs+i+L2maxt∈−τ,ts+izsiht−ysiht+hεs+k≤h∑i=0kβiL1δs+i+L2maxt∈−τ,ts+i∏siht;χ,z1,z2,…,zs+k−∏siht;ψ,y1,y2,…,ys+k+L2maxt∈−τ,ts+iγsit+hεs+k≤h∑i=0kβiL1δs+i+L2cπmaxmax1≤i¯≤s+kδi¯,maxt∈−τ,0χt−ψt+L2maxt∈−τ,ts+iγsit+hεs+k,s=k⌣−k,k⌣−k+1,…,N−k.
Then, (42) yields(43)∑s=k⌣−kn−kb¯s≤h∑i=0k∑s=k⌣−kn−kβiL1δs+i+L2cπmaxmax1≤i¯≤s+kδi¯,maxt∈−τ,0χt−ψt+L2maxt∈−τ,ts+iγsit+h∑s=k⌣−kn−kεs+k≤h∑i=0k∑j=k⌣−knβiL1δj+L2cπh∑i=0k∑s=k⌣−kn−kβimaxmax1≤i¯≤s+kδi¯,maxt∈−τ,0χt−ψt+hL2∑i=0k∑s=k⌣−kn−kβimaxt∈−τ,ts+iγsit+h∑s=k⌣−kn−kϵs+k,n=k⌣,k⌣+1,…,N.
Settingβ=max0≤i≤kβi and noticing (43), we have(44)∑s=k⌣−kn−kb¯s≤k+1hβL1∑j=k⌣−knδj+k+1hβL2cπ∑s=k⌣−kn−kmaxmax1≤i¯≤s+kδi¯,maxt∈−τ,0χt−ψt+k+1hβL2∑s=k⌣−kn−kmax0≤i≤kmaxt∈−τ,ts+iγsit+tnmaxk⌣≤j≤nεj,n=k⌣,k⌣+1,…,N.
Denote(45)Q0=δ0,Qn=maxδn,max0≤i≤n−1δi≤δn+Qn−1,n≥1.
Combining (44) with (45), we conclude(46)∑s=k⌣−kn−kb¯s≤k+1hβL1∑j=k⌣−knQj+k+1hβL2cπ∑j=k⌣nQj+k+1hβL2cπ∑s=k⌣−kn−kmaxt∈−τ,0χt−ψt+k+1hβL2∑s=k⌣−kn−kmax0≤i≤kmaxt∈−τ,ts+iγsit+tnmaxk⌣≤j≤nεj,n=k⌣,k⌣+1,…,N.
From inequalities (41) and (46), noticing Remark 2, we obtain(47)δn≤MW−1maxk⌣−k≤j≤k⌣−1δj+k+1hβML1+L2cπ∑j=k⌣−knQj+k+1βL2cπMtnmaxt∈−τ,0χt−ψt+k+1βL2Mtnmaxk⌣−k≤s≤n−kmax0≤i≤kmaxt∈−τ,ts+iγsit+Mtnmaxk⌣≤j≤nεj,n=k⌣−k,k⌣−k+1,…,N.
By (45), we have(48)∑j=k⌣−knQj=∑j=k⌣−kn−1Qj+Qn≤2∑j=k⌣−kn−1Qj+δn,n=k⌣,k⌣+1,…,N.
It can be concluded from (47) and (48) that(49)δn≤MW−1maxk⌣−k≤j≤k⌣−1δj+k+1βL2cπMtnmaxt∈−τ,0χt−ψt+k+1βL2Mtnmaxk⌣−k≤s≤n−kmax0≤i≤kmaxt∈−τ,ts+iγsit+Mtnmaxk⌣≤j≤nεj+k+1hβML1+L2cπδn+2k+1hβML1+L2cπ∑j=k⌣−kn−1Qj,n=k⌣,k⌣+1,…,N.
Let(50)α˜n=MW−1maxk⌣−k≤j≤k⌣−1δj+k+1βL2cπMtnmaxt∈−τ,0χt−ψt+k+1βL2Mtnmaxk⌣−k≤s≤n−kmax0≤i≤kmaxt∈−τ,ts+iγsit+Mtnmaxk⌣≤j≤nεj+maxt∈−τ,0χt−ψt+max1≤j≤k⌣−1δj,n=0,1,…,N.
From (49) and (50), we know that, for n=0,1,…,N,(51)δn≤α˜n+k+1hβML1+L2cπδn+2k+1hβML1+L2cπ∑j=k⌣−kn−1Qj,(52)1−k+1hβML1+L2cπδn≤α˜n+2k+1hβML1+L2cπ∑j=k⌣−kn−1Qj.
For any givenv1∈0,1, set(53)h0=min1−v1k+1βML1+L2cπ,hφ,which depends on method (6), L1 and L2, where hφ is defined in Section 2. Noticing k⌣−k≥0, when 0<h≤h0, from (52), we obtain(54)δn≤v1−1α˜n+2v1−1k+1hβML1+L2cπ∑j=k⌣−kn−1Qj≤v1−1α˜n+2v1−1k+1hβML1+L2cπ∑j=0n−1Qj,n=0,1,…,N.
Becauseα˜n does not decrease as n increases, when n≥1 and 0≤i≤n−1, from (54), we know(55)δi≤v1−1α˜n+2v1−1k+1hβML1+L2cπ∑j=0n−1Qj,0<h≤h0.
It is concluded from (54) and (55) that(56)Qn≤v1−1α˜n+2v1−1k+1hβML1+L2cπ∑j=0n−1Qj,0<h≤h0,n=0,1,…,N,which implies(57)Qi≤v1−1α˜n+2v1−1k+1hβML1+L2cπ∑j=0i−1Qj,0<h≤h0,i=0,1,…,n,n=0,1,…,N.
From (57) and Proposition 4, we have(58)Qn≤v1−1α˜ne2v1−1k+1hnβML1+L2cπ,0<h≤h0,n=0,1,…,N,so(59)δn≤Qn≤v1−1α˜ne2v1−1k+1tnβML1+L2cπ,0<h≤h0,n=0,1,…,N,(60)max0≤i≤Nδi≤v1−1α˜Ne2v1−1k+1TβML1+L2cπ≤v1−1e2v1−1k+1TβML1+L2cπMW−1+1max1≤j≤k⌣−1δj+k+1βL2cπMT+MW−1+1maxt∈−τ,0χt−ψt+k+1βL2MTmaxk⌣−k≤s≤N−kmax0≤i≤kmaxt∈−τ,ts+iγsit+MTmaxk⌣≤j≤Nεj≤Cmax1≤j≤k⌣−1δj+maxt∈−τ,0χt−ψt+maxk⌣−k≤s≤N−kmax0≤i≤kmaxt∈−τ,ts+iγsit+maxk⌣≤j≤Nεj,0<h≤h0,where(61)C=v1−1e2v1−1k+1TβML1+L2cπmaxMW−1+1,k+1βL2cπMT+MW−1+1,k+1βL2MT,MT,which is independent of h. This completes the proof.


## 4. Consistency and Convergence Analyses

In this section, for equation (1), we assume the mapping f possesses sufficiently high-order continuous partial derivatives, and its true solution yt possesses sufficiently high-order continuous derivatives on the interval0,T and the constants Mi used later which are defined as(62)Mi=sup0≤t≤Tdiytdtiare all of moderate size except for some special cases such as in the transient phase of a stiff problem.Definition 3.
It is said the piecewise Lagrangian interpolation operatorΠnih in (6) is consistent of order pp≥1 if, for any given function u:−τ,T⟶Rm which is sufficiently differentiable on the subinterval 0,T, there exists(63)max−τ≤t≤tn+iut−∏niht;ϕ,ut1,ut2,…,utn+k≤g^hp,where the function ϕ is a restriction of the function ut on the subinterval −τ,0 and the constant g^ only depends on some M˜j defined as M˜j=supt∈0,Tdjut/dtj.Remark 4.
By Definition3 and the Lagrangian interpolation theorem, we know the piecewise interpolation operators Πnih in (6) which based on pp≥0-degree Lagrangian interpolation polynomial are all consistent of order p+1.Definition 4.
It is said method (6) is consistent of order p if its mother method (5) and the piecewise interpolation operators Πnihn=k⌣−k,k⌣−k+1,…,N−k,i=0,1,…,k are all consistent of order p.Remark 5.
As for the consistency of the mother method (5), please see references [35, 38, 40].Definition 5.
Method (6) is said to be convergent of order p if, for the sequence yk⌣,yk⌣+1,…,yN determined by method (6) applied to any given equations (1) in DL1,L2 with the starting data ψ∈Cm−τ,0,y0:=ψ0,y1,y2,…,yk⌣−1∈Rm (the method (6) starts with n=k⌣−k), there exists a sufficiently small positive number h˜, such that(64)ytn−yn≤C0tnmaxmax−τ≤t≤0φt−ψt,max1≤j≤k⌣−1ytj−yj+C1tnhp,n=k⌣,k⌣+1,…,N,0<h≤h˜,where the continuous functions C0t and C1t depend on t, method (6), L1, L2, and some Mi .Theorem 2.
If method (6) is consistent of order p and satisfies the root condition, it is convergent of order p.Proof.
Assume method (6) is consistent of order p and satisfies the root condition. In this proving course, we always assume 0<h≤h0, where h0 is defined by (53). Using method (6) to solve any given equations (1) ∈DL1,L2, we easily know the true solution yt obviously satisfies the perturbation equations:(65)ynit=∏niht;φ,yt1,yt2,…,ytn+k+γnit,−τ≤t≤tn+i,i=0,1,…,k,∑i=0kαiytn+i=h∑i=0kβiftn+i,ytn+i,yni⋅+1hTn+k,n=k⌣−k,k⌣−k+1,…,N−k,where(66)γnit=yt−∏niht;φ,yt1,yt2,…,ytn+k,−τ≤t≤tn+i,i=0,1,…,k,(67)Tn+k=∑i=0kαiytn+i−hβiftn+i,ytn+i,yni⋅.
From (1) and (67), we conclude(68)Tn+k=∑i=0kαiytn+i−hβiy′tn+i,n=k⌣−k,k⌣−k+1,…,N−k.
Since the mother method (5) is consistent of order p, equation (68) implies [35, 38, 40](69)Tn+k≤C˜hp+1,n=k⌣−k,k⌣−k+1,…,N−k,where the constant C˜ only depends on the mother method (5) and Mp+1 defined above. Because all the interpolation operators Πnihn=k⌣−k,k⌣−k+1,…,N−k,i=0,1,…,k are consistent of order p, by Definition 3, we obtain(70)maxk⌣−k≤n≤N−kmax0≤i≤kmax−τ≤t≤tn+iγnit≤g^hp,where the constant g^ only depends on some Mi defined above. By Theorem 1 and assumption that method (6) satisfies the root condition, we know method (6) is zero-stable, which tells us that inequality (59) holds. From (6), (27), (50), (59), (65), (69), and (70), we further conclude that(71)ytn−yn≤v1−1e2v1−1k+1tnβML1+L2cπMW−1maxk⌣−k≤j≤k⌣−1ytj−yj+k+1βL2cπMtnmaxt∈−τ,0φt−ψt+k+1βL2Mtnmaxk⌣−k≤s≤n−kmax0≤i≤kmaxt∈−τ,ts+iγsit+Mtnhmaxk⌣≤j≤nTj+maxt∈−τ,0φt−ψt+max1≤j≤k⌣−1ytj−yj≤v1−1e2v1−1k+1tnβML1+L2cπmaxt∈−τ,0φt−ψt+max1≤j≤k⌣−1ytj−yj+MW−1maxk⌣−k≤j≤k⌣−1ytj−yj+k+1βL2cπMtnmaxt∈−τ,0φt−ψt+k+1βL2Mtng^hp+MtnC˜hp≤C^0tnmaxmaxt∈−τ,0φt−ψt,max1≤j≤k⌣−1ytj−yj+C^1tnhp,n=k⌣,k⌣+1,…,N,where(72)C^0t=v1−1MW−1+k+1βL2cπMt+2e2v1−1k+1tβML1+L2cπ,C^1t=v1−1k+1βL2Mg^+MC˜te2v1−1k+1tβML1+L2cπ,and v1,M,W−1,andβ defined in the course of proving Theorem 1 in Section 3. It can be seen that the continuous functions C^0t and C^1t depend on t, method (6), L1, L2, and some Mi. This completes the proof.Remark 6.
In this paper, the strict condition (2) can be generally weakened as the inequalities in (2) hold only in some neighborhood of the true solution yt of the equations [35, 36, 38, 41].


## 5. Numerical Examples

In this section, the linear multistep methods formed as (6) are used to solve VFDEs formed as (1) ∈DL1,L2, and for convenience, we give the starting data ψt,y0,y1,y2,…,yk⌣−1 according to true solution, i.e., let(73)ψt=φt,−τ≤t≤0,y0=φ0,yi=yti,i=1,2,…,k⌣−1.For a method formed as (6) whose stepsize is h, we denote the maximum global error:(74)Eh=max0≤n≤Nmax1≤i≤myitn−yni,where yitn and yni are the ith components of the true solution ytn and the approximation solution yn, respectively, and furthermore, according to [36, 41], the observing order p˜ of the numerical method is defined as(75)p˜=log2EhEh/2.In order to adequately verify the correctness of the theory presented in this paper, we use various different zero-stable mother methods including the Adams–Bashforth methods [38, 42], the Adams–Moulton method [38, 43], the Backward Differentiation Formula (BDF) method, and the Milne-Simpson method [44, 45] with corresponding piecewise Lagrangian interpolation operator to build the methods formed as (6), where each mother method and its corresponding interpolation operator are both consistent of the same order p2≤p≤4. According to Definition 4, these methods are consistent of order p, and by Theorem 2, we know they are convergent of order p. For convenience, these methods are named the same as their corresponding mother methods. For simplicity, let AB p denote the Adams–Bashforth method which is consistent of order p, let AM p denote the Adams–Moulton method which is consistent of order p, let B p denote the BDF method which is consistent of order p, and let MS4 denote the Milne-Simpson method which is consistent of order 4.Remark 7.
We list the mother methods used in this section as follows:(i)
The Adams–Bashforth method which is consistent of order 2 takes the form(76)−yn+1+yn+2=h−12ftn,yn+32ftn+1,yn+1.(ii)
The Adams–Bashforth method which is consistent of order 3 takes the form(77)−yn+2+yn+3=h512ftn,yn−1612ftn+1,yn+1+2312ftn+2,yn+2.(iii)
The Adams–Moulton method which is consistent of order 2 takes the form(78)−yn+yn+1=h12ftn,yn+12ftn+1,yn+1.(iv)
The BDF method which is consistent of order 2 takes the form(79)12yn−2yn+1+32yn+2=hftn+2,yn+2.(v)
The Milne-Simpson method takes the form(80)−yn+yn+2=h13ftn,yn+43ftn+1,yn+1+13ftn+2,yn+2.Example 1.
Consider the pantograph equations [36]:(81)y1′t=−2.8y1t+3.6y2t−0.4sin0.5ty10.5t+0.8sin0.5ty20.5t+1+10sint,0≤t≤1,y2′t=3.6y1t−8.2y2t+0.8sin0.5ty10.5t−1.6sin0.5ty20.5t−2−20sint,0≤t≤1,y10=2,y20=1,whose unique true solution is y1t=sint+2e−t,y2t=−2sint+e−t. Equation (81) is a special case ∈DL1,L2, where L1=10andL2=2sin0.5≈0.9588[36]. We use the methods AB2, AM2, and B2 with h=1/80×2ii=0,1,…,5 to solve (81), respectively. The maximum global errors Eh and the observing orders p˜ are listed in Table 1.
From Table1, we can see for the pantograph equation (81) the observing orders of the three methods AB2, AM2, and B2 are all very close to the theoretical convergent order 2.Table 1
The maximum global errorsEh and the observing orders p˜ produced by AB2, AM2, and B2 when applied to equation (81).
hAB2AM2B2Ehp˜Ehp˜Ehp˜1/805.179403×10−59.921322×10−64.066459×10−51/1601.297365×10−51.9972012.497546×10−61.9900211.021635×10−51.9928931/3203.246682×10−61.9985466.265474×10−71.9950162.560343×10−61.9964711/6408.120817×10−71.9992691.569068×10−71.9975166.408666×10−71.9982411/12802.030719×10−71.9996343.926050×10−81.9987571.603144×10−71.9991201/25605.077445×10−81.9998169.819351×10−91.9993794.009081×10−81.999561Example 2.
Consider the nonlinear VIDEs [46]:(82)y1′t=2y2t−13t4+cosy1t−1+∫0t2ssiny1s+sty2sds,0≤t≤1,y2′t=1−tsiny2t−12t2siny1t+∫0tst2cosy1s+tcosy2sds,0≤t≤1,y10=0,y20=0.
The unique true solution of equation (82) is y1t=t2andy2t=t. Let(83)yt=y1t,y2tT,ft,u,ψ⋅=2u2−13t4+cosu1−1+∫0t2ssinψ1s+stψ2sds,1−tsinu2−12t2sinu1+∫0tst2cosψ1s+tcosψ2sds,where 0≤t≤1, u=u1,u2T∈R2,andψ∈C20,1. Then, equation (82) can be written in form (1)∈DL1,L2, where L1 and L2 are both of moderate size. We use the method AB2 with h=1/80×2ii=0,1,…,6 to solve (82), respectively. The maximum global errors Eh and the observing orders p˜ are listed in Table 2.
From Table2, we can see, for the nonlinear VIDEs (82), the observing orders of the method AB2 are all very close to the theoretical convergent order 2.Table 2
The maximum global errorsEh and the observing orders p˜ produced by AB2 when applied to equation (82).
hEhp˜1/802.567537×10−61/1606.421382×10−71.9994291/3201.605505×10−71.9998571/6404.013861×10−81.9999651/12801.003472×10−81.9999901/25602.508682×10−91.9999991/51206.271714×10−101.999998Example 3.
Consider the nonlinear VDIDEs [36]:(84)y′t=y2t+2yt−2yt−π2+2yt2yt−π2+∫t−πt−π/2yzdz+cos2t−1,0≤t≤6,yt=sint,−π≤t≤0.
The unique true solution of equation (84) is yt=sint. Equation (84) ∈DL1,L2 in the sense of Remark 6, where L1 and L2 are both of moderate size [36]. We use MS4 with h=0.01/2ii=0,1,2,3 to solve (84), respectively, and list the maximum global errors Eh and the observing orders p˜ in Table 3.
From Table3, we can see, for the nonlinear VDIDEs (84), the observing orders of the method MS4 are all very close to the theoretical convergent order 4.Table 3
The maximum global errorsEh and the observing orders p˜ produced by MS4 when applied to equation (84).
hEhp˜0.016.757335×10−70.0054.041475×10−84.0635010.00252.483531×10−94.0244170.001251.524486×10−104.025998



---

*Source: 6633554-2021-03-12.xml*