Nonlinear Stability and D-Convergence of Additive Runge-Kutta Methods for Multidelay-Integro-Differential Equations

Abstract

This paper is devoted to the stability and convergence analysis of the Additive Runge-Kutta methods with the Lagrangian interpolation (ARKLMs) for the numerical solution of multidelay-integro-differential equations (MDIDEs). GDN-stability and D-convergence are introduced and proved. It is shown that strongly algebraically stability gives D-convergence, DA- DAS- and ASI-stability give GDN-stability. A numerical example is given to illustrate the theoretical results.

1. Introduction

Delay differential equations arise in a variety of fields as biology, economy, control theory, electrodynamics (see, e.g., [1–5]). When considering the applicability of numerical methods for the solution of DDEs, it is necessary to analyze the stability of the numerical methods. In the last three decades, many works had dealt with these problems (see, e.g., [6]). For the case of nonlinear delay differential equations, this kind of methodology had been first introduced by Torelli [7] and then developed by [8–12].

In this paper, we consider the following nonlinear multidelay-integro-differential equations (MDIDEs) with 𝑚 delays:𝑦^(𝑡)=𝑓^[1]𝑡,𝑦(𝑡),𝑦𝑡−𝜏¹,^{𝑡𝑡−𝜏}₁+𝑓_[2]𝑔^[1](𝑡,𝑠,𝑦(𝑠))𝑑𝑠𝑡,𝑦(𝑡),𝑦𝑡−𝜏²,^{𝑡𝑡−𝜏}₂+⋯+𝑓_[𝑚]𝑔^[2](𝑡,𝑠,𝑦(𝑠))𝑑𝑠𝑡,𝑦(𝑡),𝑦𝑡−𝜏^𝑚,^{𝑡𝑡−𝜏}₂𝑦(𝑡)=𝜑(𝑡),𝑡∈𝑡₀𝑔^[𝑚](𝑡,𝑠,𝑦(𝑠))𝑑𝑠,𝑡∈𝑡₀,𝑇,−𝜏,𝑡₀,(1.1) where 𝜏₁≤𝜏₂≤⋯≤𝜏_𝑚=𝜏, 𝑓^[𝑣]∶[𝑡₀,𝑇]×𝐶^𝑁×𝐶^𝑁×𝐶^𝑁→𝐶^𝑁,𝑔^[𝑣]∶[𝑡₀,𝑇]×𝐶^𝑁×𝐶^𝑁→𝐶^𝑁𝑣=1,2,…,𝑚, and 𝜑∶[𝑡₀−𝜏,𝑡₀]→𝐶^𝑁 are continuous functions such that (1.1) has a unique solution. Moreover, we assume that there exist some inner product ⟨⋅,⋅⟩ and the induced norm ||⋅|| such that Re≤𝑓𝛼^[𝑣]_𝑣𝑡,𝑦₁,𝑢₁,𝑤₁−𝑓^[𝑣]𝑡,𝑦₂,𝑢₂,𝑤₂,𝑦₁−𝑦₂‖‖𝑓‖‖𝑦₁−𝑦₂‖‖²+𝛽_𝑣‖‖𝑢₁−𝑢₂‖‖²+𝜎_𝑣‖‖𝑤₁−𝑤₂‖‖²,𝑣=1,2,…,𝑚,𝑡≥𝑡₀,(1.2)_[𝑣]‖‖𝑔𝑡,𝑦,𝑢₁,𝑤−𝑓^[𝑣]𝑡,𝑦,𝑢₂,𝑤‖‖≤𝑟_𝑣‖‖𝑢₁−𝑢₂‖‖,(1.3)_[𝑣]𝑡,𝑠,𝑤₁−𝑔^[𝑣]𝑡,𝑠,𝑤₂‖‖≤𝑟_𝑣̃‖‖𝑤₁−𝑤₂‖‖,(𝑡,𝑠)∈𝐷(1.4)

forall𝑡∈[𝑡₀,𝑇],forall𝑦,𝑦₁,𝑦₂,𝑢,𝑢₁,𝑢₂,𝑤,𝑤₁,𝑤₂∈𝐶^𝑁, (−𝛼_𝑣),𝛽_𝑣,𝜎_𝑣,𝑟_𝑣,𝑟_𝑣̃ are all nonnegative constants. Throughout this paper, we assume that the problem (1.1) has unique exact solution 𝑦(𝑡). Space discretization of some time-dependent delay partial differential equations give rises to such delay differential equations containing additive terms with different stiffness properties. In these situations, additive Runge-Kutta (ARK) methods are used. Some recent works about ARK can refer to [13, 14]. For the additive MDIDEs (1.1), similar to the proof of Theorem 2.1 in [7], it is straightforward to prove that under the conditions (1.2)~(1.4), the analytic solutions satisfy‖𝑦(𝑡)−𝑧(𝑡)‖≤max_{𝑡0−𝜏≤𝑡≤𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖,(1.5) where 𝑧(𝑡) is the solution of the perturbed problem to (1.1).

To demand the discrete numerical solutions to preserve the convergence properties of the analytic solutions, Torelli [7] introduced a concept of RN-, GRN-stability for numerical methods applied to dissipative nonlinear systems of DDEs such as (1.1) when 𝑔^[𝑣](𝑡,𝑠,𝑦(𝑠))=0,𝑣=1,2,…,𝑚, which is the straightforward generalization of the well-known concept of BN-stability of numerical methods with respect to dissipative systems of ODEs (see also [9]). More recently, one has noticed a growing interesting the analysis of delay integro-differential equations (DIDEs). This type of equations have been investigated in various fields, such as mathematical biology and control theory (see [15–17]). The theory of computational methods for delay integro-differential equations (DIDEs) has been studied by many authors, and a great deal of interesting results have been obtained (see [18–22]). Koto [23] dealt with the linear stability of Runge-Kutta (RK) methods for systems of DIDEs; Huang and Vandewalle [24] gave sufficient and necessary stability conditions for exact and discrete solutions of linear Scalar DIDEs. However, little attention has been paid to nonlinear multidelay-integro-differential equations (MDIDEs).

So, the aim of this paper is the study of stability and convergence properties for ARK methods when they are applied to nonlinear multidelay-integro-differential equations (MDIDEs) with 𝑚 delays.

2. The GDN-Stability of the Additive Runge-Kutta Methods

An additive Runge-Kutta method with the Lagrangian interpolation (ARKLM) of 𝑠 stages and 𝑚 levels can be organized in the Butcher tableau:𝐶𝐴^[1]𝐴^[2]⋯𝐴^[𝑚]𝑏^[1]𝑇𝑏^[2]𝑇⋯𝑏^[𝑚]𝑇=𝑐¹𝑎^[1]11𝑎^[1]12⋯𝑎^[1]1𝑠𝑎^[𝑚]11𝑎^[𝑚1]12⋯𝑎^[𝑚]1𝑠𝑐₂𝑎^[1]21𝑎^[1]22⋯𝑎^[1]2𝑠𝑎^[𝑚]21𝑎^[𝑚]22⋯𝑎^[𝑚]2𝑠⋮𝑐⋮⋮⋱⋮⋯⋮⋮⋱⋮_𝑠𝑎^[1]𝑠1𝑎^[1]𝑠2⋯𝑎^[1]𝑠𝑠𝑎^[𝑚]𝑠1𝑎^[𝑚]𝑠2⋯𝑎^{[𝑚]𝑠𝑠}𝑏^1[1]𝑏^2[1]⋯𝑏^𝑠[1]⋯𝑏^1[𝑚]𝑏^2[𝑚]⋯𝑏^𝑠[𝑚],(2.1) where 𝐶=[𝑐₁,𝑐₂,…,𝑐_𝑠]^𝑇, 𝑏^[𝑣]=[𝑏^1[𝑣],𝑏^2[𝑣],…,𝑏^𝑠[𝑣]], and 𝐴^[𝑣]=(𝑎^{[𝑣]𝑖𝑗})^{𝑠𝑖,𝑗=1}.

The adoption of the method (2.1) for solving the problem (1.1) leads to𝑦_𝑛+1=𝑦_𝑛+ℎ^𝑚_𝑠𝑣=1_𝑗=1𝑦_𝑖(𝑛)𝑏^𝑗[𝑣]𝑓^[𝑣]𝑡^𝑗(𝑛),𝑦^𝑗(𝑛),𝑦̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)},=𝑦_𝑛+ℎ^𝑚_𝑠𝑣=1_𝑗=1𝑎^{[𝑣]𝑖𝑗}𝑓^[𝑣]𝑡^𝑗(𝑛),𝑦^𝑗(𝑛),𝑦̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)},(2.2) where 𝑡_𝑛=𝑡₀+𝑛ℎ,𝑡^𝑗(𝑛)=𝑡_𝑛+𝑐_𝑗ℎ, 𝑦_𝑛, and 𝑦^𝑗(𝑛), 𝑦̃^{𝑗[𝑣](𝑛)} are approximations to the analytic solution 𝑦(𝑡_𝑛), 𝑦(𝑡_𝑛+𝑐_𝑗ℎ), 𝑦(𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣) of (1.1), respectively, and the argument 𝑦̃^{𝑗[𝑣](𝑛)} is determined by𝑦̃^{𝑗[𝑣](𝑛)}=⎧⎪⎨⎪⎩𝜑𝑡^𝑛+𝑐_𝑗ℎ−𝜏_𝑣𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣≤0_𝑟𝑃_𝑣=−𝑑𝐿_𝑃𝑣𝛿_𝑣𝑦^{(𝑛−𝑚𝑣+𝑃𝑣)𝑗}𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣>0,(2.3) with 𝜏_𝑣=(𝑚_𝑣−𝛿_𝑣)ℎ, 𝛿_𝑣∈[0,1), integer 𝑚_𝑣≥𝑟+1, 𝑟,𝑑≥0, and𝐿_𝑃𝑣𝛿_𝑣=^𝑟_𝑘=_{−𝑑𝑘≠𝑃𝑣}𝛿^𝑣−𝑘𝑃_𝑣−𝑘,𝑃^𝑣=−𝑑,−𝑑+1,…,𝑟.(2.4) We assume 𝑚_𝑣≥𝑟+1 is to guarantee that no (unknown) values 𝑦^𝑗(𝑖) with 𝑖≥𝑛 are used in the interpolation procedure𝑤^{𝑗[𝑣](𝑛)}isanapproximationto𝑤𝑡^𝑗(𝑛)∶=^{𝑡𝑗(𝑛)}_{𝑡𝑣)𝑗(𝑛−𝑚}𝑔^[𝑣]𝑡^𝑗(𝑛),𝑠,𝑦(𝑠)𝑑𝑠,(2.5) which can be computed by a appropriate compound quadrature rule:𝑤^{𝑗[𝑣](𝑛)}=ℎ^𝑚𝑣_𝑞=0𝑑_𝑞𝑔^[𝑣]𝑡^𝑗(𝑛),𝑡^{𝑗(𝑛−𝑞)},𝑦^{𝑗(𝑛−𝑞)},𝑣=1,2,…,𝑚,𝑗=1,2,…,𝑠.(2.6) As for the quadrature rule (2.6), we usually adopt the compound trapezoidal rule, the compound Simpsons rule or the compound Newton-Cotes rule, and so forth according to the requirement of the convergence of the method (see [19]) and denote 𝑀=max_{1≤𝑣≤𝑚}{𝑚_𝑣} and 𝜂=max_{1≤𝑣≤𝑚}{𝜂_𝑣} with 𝜂_𝑣 satisfing ∑^𝑚𝑣_𝑞=0|𝑑_𝑞|<𝜂_𝑣, 𝑣=1,2,…,𝑚.

In addition, we always put𝑦^𝑗(𝑛)=𝜑(𝑡_𝑛+𝑐_𝑗ℎ), 𝑦_𝑛=𝜑(𝑡_𝑛) whenever 𝑛≤0.

In order to write (2.2), (2.3), (2.5), and (2.6) in a more compact way, we introduce some notations. The 𝑁×𝑁 identity matrix will be denoted by 𝐼_𝑁, 𝑒=(1,1,…,1)^𝑇∈𝑅^𝑆, 𝐺=𝐺⊗𝐼_𝑁 is the Kronecker product of matrix 𝐺 and 𝐼_𝑁. For 𝑢=(𝑢₁,𝑢₂,…,𝑢_𝑠)^𝑇, 𝑣=(𝑣₁,𝑣₂,…,𝑣_𝑠)^𝑇∈𝐶^𝑁𝑆, we define the inner product and the induced norm in 𝐶^𝑁𝑆 as follows: ⟨𝑢,𝑣⟩=^𝑠_𝑖=1⟨𝑢_𝑖,𝑣_𝑖⟩,‖𝑢‖=⎷^𝑠_𝑖=1‖‖𝑢_𝑖‖‖².(2.7)

Moreover, we also adopt that𝑦^(𝑛)=⎡⎢⎢⎢⎣𝑦𝑦^{1(𝑛)2(𝑛)}𝑦⋮_𝑠(𝑛)⎤⎥⎥⎥⎦,𝑦̃_[𝑣](𝑛)=⎡⎢⎢⎢⎣𝑦𝑦̃̃^{1[𝑣](𝑛)2[𝑣](𝑛)}𝑦̃⋮_{𝑠[𝑣](𝑛)}⎤⎥⎥⎥⎦,𝑤_[𝑣](𝑛)=⎡⎢⎢⎢⎣𝑤𝑤^{1[𝑣](𝑛)2[𝑣](𝑛)}𝑤⋮_{𝑠[𝑣](𝑛)}⎤⎥⎥⎥⎦,𝑇_(𝑛)=⎡⎢⎢⎢⎣𝑡𝑡^{1(𝑛)2(𝑛)}𝑡⋮_𝑠(𝑛)𝑓_[𝑣]⎤⎥⎥⎥⎦,𝑇^(𝑛),𝑦^(𝑛),𝑦̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)=⎡⎢⎢⎢⎢⎢⎣𝑓^[𝑣]𝑓_[𝑣]𝑡^1(𝑛),𝑦^1(𝑛),𝑦̃^{1[𝑣](𝑛)},𝑤^{1[𝑣](𝑛)}𝑓_[𝑣]𝑡^2(𝑛),𝑦^2(𝑛),⋮𝑦̃^{2[𝑣](𝑛)},𝑤^{2[𝑣](𝑛)}𝑡^𝑠(𝑛),𝑦^𝑠(𝑛),𝑦̃^{𝑠[𝑣](𝑛)},𝑤^{𝑠[𝑣](𝑛)}⎤⎥⎥⎥⎥⎥⎦.(2.8)

With the above notation, method (2.2),(2.3), (2.5), and (2.6) can be written as𝑦_𝑛+1=𝑦_𝑛+ℎ^𝑚_𝑣=1𝑦_(𝑛)𝑏̃^[𝑣]𝑇𝑓^[𝑣]𝑇^(𝑛),𝑦^(𝑛),𝑦̃^[𝑣](𝑛),𝑤^[𝑣](𝑛),=𝑒𝑦̃_𝑛+ℎ^𝑚_𝑣=1𝑦̃_{𝑗[𝑣](𝑛)}𝐴^[𝑣]𝑓^[𝑣]𝑇^(𝑛),𝑦^(𝑛),𝑦̃^[𝑣](𝑛),𝑤^[𝑣](𝑛),=⎧⎪⎨⎪⎩𝑒𝜑̃𝑡^𝑛+𝑐_𝑗ℎ−𝜏_𝑣,𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣≤𝑡₀,_𝑟𝑃_𝑣=−𝑑𝑤_{𝑗[𝑣](𝑛)}𝐿_𝑃𝑣𝛿_𝑣𝑦^{(𝑛−𝑚𝑣+𝑃𝑣)𝑗},𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣>𝑡₀,=ℎ^𝑚𝑣_𝑞=0𝑑_𝑞𝑔^[𝑣]𝑡^𝑛+𝑐_𝑗ℎ,𝑡_𝑛−𝑞+𝑐_𝑗ℎ,𝑦^{𝑗(𝑛−𝑞)}.(2.9) In 1997, Zhang and Zhou [25] introduced the extension of RN-stability to GDN-stability as follows.

Definition 2.1. An ARKLM (2.1) for DDEs is called GDN-stable if, numerical approximations 𝑦_𝑛 and 𝑧_𝑛 to the solution of (1.1) and its perturbed problem, respectively, satisfy ‖‖𝑦_𝑛−𝑧_𝑛‖‖≤𝐶max_{𝑡0−𝜏≤𝑡≺𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖,𝑛≥0,(2.10) where constant 𝐶>0 depends only on the method, the parameter 𝛼_𝑣,𝛽_𝑣,𝜎_𝑣,𝑟_𝑣,𝑟_𝑣̃, and the interval length 𝑇−𝑡₀, 𝜓(𝑡) is the initial function to the perturbed problem of (1.1).

Definition 2.2. An ARKLM (2.1) is called strongly algebraically stable if matrices 𝑀_𝛾𝜇 are nonnegative definite, where 𝑀_𝛾𝜇=𝐵^[𝛾]𝐴^[𝜇]+𝐴^𝑇[𝛾]𝐵^[𝜇]−𝑏^[𝛾]𝑏^𝑇[𝜇],𝐵^[𝛾]=diag𝑏^1[𝛾],𝑏^2[𝛾],…,𝑏^𝑠[𝛾],(2.11) for 𝜇,𝛾=1,2,…,𝑚.
Let 𝑦𝑧_𝑛𝑛,𝑦^𝑗(𝑛),𝑦̃^{𝑗[1](𝑛)},𝑦̃^{𝑗[2](𝑛)},…,𝑦̃^{𝑗[𝑚](𝑛)},𝑤^{𝑗[1](𝑛)},𝑤^{𝑗[2](𝑛)},…,𝑤^{𝑗[𝑚](𝑛)}^𝑠𝑗=1,,𝑧^𝑗(𝑛),𝑧̃^{𝑗[1](𝑛)},𝑧̃^{𝑗[2](𝑛)},…,𝑧̃^{𝑗[𝑚](𝑛)},𝑤^{𝑗[1](𝑛)},𝑤^{𝑗[2](𝑛)},…,𝑤^{𝑗[𝑚](𝑛)}^𝑠𝑗=1(2.12)be two sequences of approximations to problems (1.1) and its perturbed problem, respectively. From method (2.1) with the same step size ℎ, and write 𝑄𝑈^𝑖(𝑛)_{𝑖[𝑣](𝑛)}=𝑦^𝑖(𝑛)−𝑧^𝑖(𝑛),𝑈^{𝑖[𝑣](𝑛)}=𝑦̃^{𝑖[𝑣](𝑛)}−𝑧̃^{𝑖[𝑣](𝑛)},𝑈^0(𝑛)=𝑦_𝑛−𝑧_𝑛,=ℎ𝑓^[𝑣]𝑡^𝑖(𝑛),𝑦^𝑖(𝑛),𝑦̃^{𝑖[𝑣](𝑛)},𝑤^{𝑖[𝑣](𝑛)}𝑖=−1𝑓,2,^[𝑣]…𝑡,𝑠^𝑖(𝑛),,𝑧^𝑖(𝑛)𝑣,𝑧=̃1,^{𝑖[𝑣](𝑛)}2,,𝑤…,𝑚^{𝑖[𝑣](𝑛)}.,(2.13) Then (2.2) and (2.3) read 𝑈^0(𝑛+1)=𝑈^0(𝑛)+^𝑚_𝑠𝑣=1_𝑗=1𝑈𝑏^𝑗[𝑣]𝑄^{𝑗[𝑣](𝑛)},(2.14)_𝑖(𝑛)=𝑈^0(𝑛)+^𝑚_𝑠𝑣=1_𝑗=1𝑎^{[𝑣]𝑖𝑗}𝑄^{𝑗[𝑣](𝑛)}.(2.15)

Our main results about GDN-stability are contained in the following theorem.

Theorem 2.3. Assume ARK method (2.2) is strongly algebraically stable, and then the corresponding ARKLM (2.1) is GDN-stable, and satisfies ‖‖𝑦_𝑛−𝑧_𝑛‖‖≤𝐶max_{𝑡0−𝜏≤𝑡≤𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖²,𝑛≥0,(2.16) where𝐶=exp6𝑇−𝑡₀𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1max_{𝑡0−𝜏≤𝑡≤𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖²,𝐿_𝑣=max_{−𝑑≤𝑝𝑣≤𝑟}𝐿_𝑝𝜈,(2.17)

Proof. From (2.14) and (2.15) we get ‖‖𝑈^0(𝑛+1)‖‖²=𝑈^0(𝑛)+^𝑚_𝑠𝑣=1_𝑖=1𝑏^𝑖[𝑣]𝑄^{𝑖[𝑣](𝑛)},𝑈^0(𝑛)+^𝑚_𝑠𝑣=1_𝑖=1𝑏^𝑖[𝑣]𝑄^{𝑖[𝑣](𝑛)}=‖‖𝑈_0(𝑛)‖‖²+2^𝑚_𝑠𝑣=1_𝑖=1+_𝑚𝑏^𝑖[𝑣]Re𝑄^{𝑖[𝑣](𝑛)},𝑈^0(𝑛)_𝑠_𝑢,𝑣=1_𝑖,𝑗=1=‖‖𝑈_0(𝑛)𝑏^𝑖[𝑢]𝑏^𝑗[𝑣]𝑄^{𝑖[𝑢](𝑛)},𝑄^{𝑗[𝑣](𝑛)}‖‖²+2^𝑚_𝑠𝑣=1_𝑖=1𝑏^𝑖[𝑣]Re𝑄^{𝑖[𝑣](𝑛)},𝑈^𝑖(𝑛)−^𝑚_𝑠𝑣=1_𝑗=1𝑎^{[𝑣]𝑖𝑗}𝑄^{𝑗[𝑣](𝑛)}+_𝑚_𝑠_𝑢,𝑣=1_𝑖,𝑗=1=‖‖𝑈_0(𝑛)𝑏^𝑖[𝑢]𝑏^𝑗[𝑣]𝑄^{𝑖[𝑣](𝑛)},𝑄^{𝑗[𝑢](𝑛)}‖‖²+2^𝑚_𝑠𝑣=1_𝑖=1−_𝑚𝑏^𝑖[𝑣]Re𝑄^{𝑖[𝑣](𝑛)},𝑈^𝑖(𝑛)_𝑠_𝑢,𝑣=1_𝑖,𝑗=1𝑏^𝑖[𝑢]𝑎^{[𝑣]𝑖𝑗}+𝑏^𝑗[𝑣]𝑎^{[𝑢]𝑖𝑗}−𝑏^𝑖[𝑢]𝑏^𝑗[𝑣]𝑄^{𝑖[𝑣](𝑛)},𝑄^{𝑗[𝑢](𝑛)}.(2.18) If the matrices 𝑀_𝛾𝜇 are nonnegative definite, then ‖‖𝑈^0(𝑛+1)‖‖²≤‖‖𝑈^0(𝑛)‖‖²+2^𝑚_𝑠𝑣=1_𝑗=1𝑏^𝑗[𝑣]Re𝑄^{𝑗[𝑣](𝑛)},𝑈^𝑗(𝑛).(2.19) Furthermore, by conditions (1.2)~(1.4) and Schwartz inequality we have Re𝑄^{𝑗[𝑣](𝑛)},𝑈^𝑗(𝑛)=ℎRe−𝑓𝑓^[𝑣]_[𝑣]𝑡^𝑗(𝑛),𝑦^𝑗(𝑛),𝑦̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}≤ℎ𝛼_𝑣𝑡^𝑗(𝑛),𝑧^𝑗(𝑛),𝑧̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)},𝑈^𝑗(𝑛)‖‖𝑈^𝑗(𝑛)‖‖²+ℎ𝛽_𝑣‖‖𝑈^{𝑗[𝑣](𝑛)}‖‖+ℎ𝜎_𝑣‖‖𝑤^{𝑗[𝑣](𝑛)}−𝑤^{𝑗[𝑣](𝑛)}‖‖²≤ℎ𝛼_𝑣‖‖𝑈^𝑗(𝑛)‖‖²+ℎ𝛽_𝑣‖‖𝑈^{𝑗[𝑣](𝑛)}‖‖+2ℎ³𝑟^2𝑣̃𝜂^2𝑣𝜎_𝑣𝑚^𝑣_𝑞=0‖‖𝑈^{𝑗(𝑛−𝑞)}‖‖²=ℎ𝛼_𝑣(2.20)‖‖𝑈^𝑗(𝑛)‖‖²+2ℎ³𝑟^2𝑣̃𝜂^2𝑣𝜎_𝑣𝑚^𝑣_𝑞=0‖‖𝑈^{𝑗(𝑛−𝑞)}‖‖²,𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣≤𝑡₀=ℎ𝛼_𝑣‖‖𝑈^𝑗(𝑛)‖‖²+ℎ𝛽_𝑣‖‖‖‖^𝑟_𝑝𝑣_=−𝑑𝐿_𝑝𝑣𝛿_𝑣𝑈^{(𝑛−𝑚𝑣+𝑝𝑣)𝑗}‖‖‖‖²+2ℎ³𝑟^2𝑣̃𝜂^2𝑣𝜎_𝑣𝑚^𝑣_𝑞=0‖‖𝑈^{𝑗(𝑛−𝑞)}‖‖²,𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣>𝑡₀≤ℎ𝛼_𝑣(2.21)‖‖𝑈^𝑗(𝑛)‖‖²+2ℎ³𝑟^2𝑣̃𝜂^2𝑣𝜎_𝑣𝑚^𝑣_𝑞=0≤ℎ𝛼_𝑣‖‖𝑈^{𝑗(𝑛−𝑞)}‖‖²,𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣≤𝑡₀(2.22)‖‖𝑈^𝑗(𝑛)‖‖²+2ℎ𝛽_𝑣𝐿^2𝑣‖‖‖‖^𝑟_𝑝𝑣_=−𝑑𝑈^{(𝑛−𝑚𝑣+𝑝𝑣)𝑗}‖‖‖‖²+2ℎ³𝑟^2𝑣̃𝜂^2𝑣𝜎_𝑣𝑚^𝑣_𝑞=0‖‖𝑈^{𝑗(𝑛−𝑞)}‖‖²,𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣>𝑡₀.(2.23) For (2.23), we have (2.23)≤2ℎ𝛽_𝑣𝐿^2𝑣‖‖‖‖^𝑟_𝑝𝑣_=−𝑑𝑈^{(𝑛−𝑚𝑣+𝑝𝑣)𝑗}‖‖‖‖²+2ℎ³𝑟^2𝑣̃𝜂^2𝑣𝜎_𝑣𝑚^𝑣_𝑞=0‖‖𝑈^{𝑗(𝑛−𝑞)}‖‖²≤3ℎ𝛽_𝑣𝐿^2𝑣‖‖‖‖^𝑚𝑣_𝑝𝑣_=−𝑑𝑈^{(𝑛−𝑚𝑣+𝑝𝑣)𝑗}‖‖‖‖².(2.24) By the same way, we can also get (2.22)≤3ℎ𝛽_𝑣𝐿^2𝑣‖‖‖‖^𝑚𝑣_𝑃_𝑣=−𝑑𝑈^{(𝑛−𝑚𝑣+𝑝𝑣)𝑗}‖‖‖‖².(2.25) Substituting (2.25) and (2.24) into (2.19), yields ‖‖𝑈^0(𝑛+1)‖‖²≤‖‖𝑈^0(𝑛)‖‖²+2ℎ^𝑚_𝑠𝑣=1_𝑗=13𝛽_𝑣𝐿^2𝑣𝑏^𝑗[𝑣]‖‖‖‖^𝑚𝑣_𝑝𝑣_=−𝑑𝑈^{(𝑛−𝑚𝑣+𝑝𝑣)𝑗}‖‖‖‖²≤‖‖𝑈(2.26)_0(𝑛)‖‖²+6ℎ^𝑚_𝑠𝑣=1_𝑗=1𝛽_𝑣𝐿^2𝑣𝑏^𝑗[𝑣]𝑚_𝑣+𝑑+1max_{−𝑑≤𝑝𝑣≤𝑚𝑣}‖‖𝑈^{(𝑛−𝑚𝑣+𝑝𝑣)𝑗}‖‖²≤‖‖𝑈^0(𝑛)‖‖²+6ℎ𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1max_{(𝑗,𝑝𝑣)∈𝐸𝑣}‖‖𝑈^{(𝑛−𝑚𝑣+𝑝𝑣)𝑗}‖‖²≤1+6ℎ𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1max‖‖𝑈^0(𝑛)‖‖²,max_{(𝑗,𝑝𝑣)∈𝐸𝑣}‖‖𝑈^{(𝑛−𝑚𝑣+𝑝𝑣)𝑗}‖‖²,(2.27) where 𝐸_𝑣={(𝑗,𝑃_𝑣)1≤𝑗≤𝑠,−𝑑≤𝑃_𝑣≤𝑟}.
Similar to (2.27), the inequalities: ‖‖𝑈^𝑖(𝑛)‖‖²≤1+6ℎ𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1max‖‖𝑈^0(𝑛)‖‖²,max_{(𝑗,𝑃𝑣)∈𝐸}‖‖𝑈^{(𝑛−𝑚𝑣+𝑃𝑣)𝑗}‖‖²(2.28) follows for 𝑖=1,2,…,𝑠.
In the following, with the help of inequalities (2.27), (2.28), and induction we shall prove the inequalities: ‖‖𝑈^𝑖(𝑛)‖‖²≤1+6ℎ𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1^(𝑛+1)max_{𝑡0−𝜏≤𝑡≤𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖²,(2.29) for 𝑛≥0, 𝑖=1,2,…,𝑠.
In fact, it is clear from (2.27), (2.28), and 𝑚_𝑣≥𝑟+1 such that ‖‖𝑈^𝑖(0)‖‖²≤1+6ℎ𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1max_{𝑡0−𝜏≤𝑡≤𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖²,𝑖=0,1,2,…,𝑠.(2.30) Suppose for 𝑛≤𝑘(𝑘≥0)that ‖‖𝑈^𝑖(𝑛)‖‖²≤1+6ℎ𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1^(𝑛+1)max_{𝑡0−𝜏≤𝑡≤𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖²,𝑖=0,1,2,…,𝑠.(2.31) Then from (2.27) and (2.28), 𝑚_𝑣≥𝑟+1 and 1+6ℎ𝑚𝑠∑^𝑚𝑣=1𝛽_𝑣𝐿^2𝑣(𝑚_𝑣+𝑑+1)>1, we conclude that ‖‖𝑈^𝑖(𝑘+1)‖‖²≤1+6ℎ𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1^(𝑘+2)max_{𝑡0−𝜏≤𝑡≤𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖²,𝑖=0,1,2,…,𝑠.(2.32) This completes the proof of inequalities (2.29). In view of (2.29), we get for 𝑛≥0 that ‖‖𝑈^0(𝑛)‖‖²≤1+6ℎ𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1^(𝑛+1)max_{𝑡0−𝜏≤𝑡≤𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖²≤exp(𝑛+1)6ℎ𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1max_{𝑡0−𝜏≤𝑡≤𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖²≤exp6𝑇−𝑡₀𝑚𝑠^𝑚_𝑣=1𝛽_𝑣𝐿^2𝑣𝑚_𝑣+𝑑+1max_{𝑡0−𝜏≤𝑡≤𝑡0}‖𝜑(𝑡)−𝜓(𝑡)‖².(2.33)As a result, we know that method (2.1) is GDN-stable.

3. D-Convergence

In order to study the convergence of numerical methods for MDIDEs, we have to mention the concept of the convergence for stiff ODEs.

In 1981, Frank et al. [26] introduced the important concept of B-convergence for numerical methods applied to nonlinear stiff initial value problems of ordinary differential equations. Later, there have been rapid developments in the study of B-convergence, and a significant number of important results have already been found for Runge-Kutta methods. In fact, B-convergence result is nothing but a realistic global error estimate based on one-sided Lipschitz constant [27]. In this section, we start discussing the convergence of ARKLM (2.1) for MDIDEs (1.1) with conditions (1.2)–(1.4). The approach to the derivation of these estimates is similar to that used in [25]. We assume the analytic solution 𝑦(𝑡) of (1.1) is smooth enough, and its derivatives used later are bounded by‖‖𝐷^(𝑖)𝑦(𝑡)‖‖≤𝑀_𝑖,𝑡∈𝑡₀−𝜏,𝑇,(3.1)

where 𝐷^(𝑖)𝑦(𝑡)=𝑦^(𝑖)𝑦(𝑡),𝑡∈𝑡₀+(𝑗−1)ℎ,𝑡₀+𝑗ℎ_(𝑖),𝑡₀+𝑗ℎ−0,𝑡=𝑡₀+𝑗ℎ.(3.2) If we introduce some notations𝑌^(𝑛)=⎡⎢⎢⎢⎣𝑦𝑡^𝑛𝑦𝑡_𝑛+𝑐₁ℎ𝑦𝑡_𝑛⋮+𝑐₂ℎ+𝑐_𝑠ℎ⎤⎥⎥⎥⎦,𝑌^[𝑣](𝑛)=⎡⎢⎢⎢⎣𝑦𝑡^𝑛𝑦𝑡_𝑛+𝑐₁ℎ−𝜏_𝑣𝑦𝑡_𝑛+𝑐⋮₂ℎ−𝜏_𝑣+𝑐_𝑠ℎ−𝜏_𝑣⎤⎥⎥⎥⎦,𝑤^[𝑣](𝑛)=⎡⎢⎢⎢⎣𝑤𝑡^𝑛𝑤𝑡_𝑛+𝑐₁ℎ−𝜏_𝑣𝑤𝑡_𝑛+⋮𝑐₂ℎ−𝜏_𝑣+𝑐_𝑠ℎ−𝜏_𝑣⎤⎥⎥⎥⎦.(3.3) With the above notations, the local errors in (2.9) can be defined as𝑦𝑡_𝑛+1=𝑦𝑡_𝑛+ℎ^𝑚_𝑣=1𝑏̃^[𝑣]𝑇𝑓^[𝑣]𝑌𝑇^(𝑛),𝑌^(𝑛),𝑌^[𝑣](𝑛),𝑤^[𝑣](𝑛)+𝑄^𝑛,(3.4)_(𝑛)=𝑒𝑦̃𝑡_𝑛+ℎ^𝑚_𝑣=1𝐴^[𝑣]𝑓^[𝑣]𝑌𝑇^(𝑛),𝑌^(𝑛),𝑌^[𝑣](𝑛),𝑤^[𝑣](𝑛)+𝑟^𝑛,(3.5)_[𝑣](𝑛)=𝑌^{1[𝑣](𝑛)},𝑌^{2[𝑣](𝑛)},…,𝑌^{𝑠[𝑣](𝑛)}^𝑇,(3.6) with𝑌^{𝑗[𝑣](𝑛)}=⎧⎪⎨⎪⎩𝜑𝑡^𝑛+𝑐_𝑗ℎ−𝜏_𝑣,𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣≤𝑡₀,_𝑟𝑃_𝑣=−𝑑𝑤𝐿_𝑃𝑣𝛿_𝑣𝑦^{(𝑛−𝑚𝑣+𝑃𝑣)𝑗}+𝜌^{𝑗[𝑣](𝑛)},𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣>𝑡₀,(3.7)_{𝑗[𝑣](𝑛)}=ℎ^𝑚𝑣_𝑞=0𝑑_𝑞𝑔^[𝑣]𝑡^𝑛+𝑐_𝑗ℎ,𝑡_𝑛−𝑞+𝑐_𝑗ℎ,𝑦^{𝑗(𝑛−𝑞)}+𝑅^{𝑗[𝑣](𝑛)}.(3.8) If we take 𝑦̆_𝑛=𝑦(𝑡_𝑛), 𝑦̆^(𝑛)=𝑌^(𝑛), 𝑦̆^[𝑣](𝑛)=𝑌^[𝑣](𝑛), and 𝑤̆^[𝑣](𝑛)=𝑤^[𝑣](𝑛)

Then we can get the perturbed scheme of (2.9),𝑦̆_𝑛+1=𝑦̆_𝑛+ℎ^𝑚_𝑣=1𝑏̃^[𝑣]𝑇𝑓^[𝑣]𝑦̆_(𝑛)𝑇^(𝑛),𝑦̆^(𝑛),𝑦̆̃^[𝑣](𝑛),𝑤̆^[𝑣](𝑛)+𝑄^𝑛,(3.9)=𝑒𝑦̃̆_𝑛+ℎ^𝑚_𝑣=1𝐴^[𝑣]𝑓^[𝑣]𝑦̆̃_{𝑗[𝑣](𝑛)}𝑇^(𝑛),𝑦̆^(𝑛),𝑦̆̃^[𝑣](𝑛),𝑤̆^[𝑣](𝑛)+𝑟^𝑛,(3.10)=⎧⎪⎨⎪⎩𝑒𝜑̃𝑡^𝑛+𝑐_𝑗ℎ−𝜏_𝑣,𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣≤0,_𝑟𝑃_𝑣=−𝑑𝑤𝐿_𝑃𝑣𝛿_𝑣𝑦̆^{(𝑛−𝑚𝑣+𝑃𝑣)𝑗}+𝜌_{𝑗[𝑣](𝑛)},𝑡_𝑛+𝑐_𝑗ℎ−𝜏_𝑣>0,(3.11)_{𝑗[𝑣](𝑛)}=ℎ^𝑚𝑣_𝑞=0𝑑_𝑞𝑔^[𝑣]𝑡^𝑛+𝑐_𝑗ℎ,𝑡_𝑛−𝑞+𝑐_𝑗ℎ,𝑦^{𝑗(𝑛−𝑞)}+𝑅^{𝑗[𝑣](𝑛)}.(3.12) With perturbations, 𝑄_𝑛∈𝐶^𝑁, 𝑟_𝑛=(𝑟^(𝑛)𝑇₁,𝑟^(𝑛)𝑇₂,…,𝑟^(𝑛)𝑇_s)^𝑇,𝑅^[𝑣](𝑛)=(𝑅^{1[𝑣](𝑛)},𝑅^{2[𝑣](𝑛)},…,𝑅^{𝑠[𝑣](𝑛)})^𝑇, 𝜌^(𝑛)=(𝜌^(𝑛)𝑇₁,𝜌^(𝑛)𝑇₂,…,𝜌^(𝑛)𝑇_s)∈𝐶^𝑁𝑆, according to Taylor formula and the formula in [28, pages 205–212], 𝑄_𝑛,𝑟_𝑛 and 𝜌_𝑛 can be determined respectively, as follows:𝑄_𝑛=^𝑃_𝑙=1ℎ^𝑙(𝑙−1)!1𝑙−^𝑚_𝑠𝑣=1_𝑗=1𝑏^𝑗[𝑣]𝑐^𝑗𝑙−1𝐷^(𝑙)𝑟𝑦𝑡_𝑛+𝑅^0(𝑛),(3.13)_𝑖(𝑛)=^𝑃_𝑙=1ℎ^𝑙(𝑙−1)!1𝑙𝑐^𝑙𝑖−^𝑚_𝑠𝑣=1_𝑗=1𝑎^{[𝑣]𝑖𝑗}𝑐^𝑗𝑙−1𝐷^(𝑙)𝜌𝑦𝑡_𝑛+𝑅^𝑖(𝑛),(3.14)_𝑖(𝑛)=ℎ^𝑞+1(𝑞+1)!^𝑚_𝑟𝑣=1_{𝑃𝑣=−𝑑}𝛿_𝑣−𝑃_𝑣𝐷^(𝑞+1)𝑦𝜉^𝑖(𝑛),𝜉^𝑖(𝑛)∈𝑡_{𝑛−𝑚𝑣−𝑑}+𝑐_𝑖ℎ,𝑡_{𝑛−𝑚𝑣+𝑟}+𝑐_𝑖ℎ,(3.15) where 𝑞=𝑑+𝑟,𝑅^𝑖(𝑛), and 𝜉^𝑖(𝑛) satisfy ||𝑅^𝑖(𝑛)||≤𝑀_𝑖ℎ^𝑖+1,𝑖=0,1,2,…,𝑠, ℎ∈(0,ℎ₀], ℎ₀ depends only on the method, and 𝑀_𝑖(𝑖=0,1,2,…,𝑠) depends only on the method and some 𝑀_𝑖 in (3.2).

Combining (2.2), (2.3), (2.5), and (2.6) with (3.9), (3.10), (3.11), and (3.12) yields the following recursion scheme for the 𝜀^0(𝑛+1)=𝑦̆_𝑛+1−𝑦_𝑛+1:𝜀^0(𝑛+1)=𝜀^0(𝑛)+ℎ^𝑚_𝑣=1𝑏̃^[𝑣]T+𝑔𝑓^[𝑣]_𝑛[𝑣]𝑇^(𝑛),𝑦̆^(𝑛),𝑦̆̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)−𝑓^[𝑣]𝑇^(𝑛),𝑦̆^(𝑛),𝑦̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)𝜀𝜀_𝑛+𝐻^[𝑣](𝑛)𝑤̆^[𝑣](𝑛)−𝑤^[𝑣](𝑛)+𝑄𝑛,_𝑛=𝑒𝜀̃^0(𝑛)+ℎ^𝑚_𝑣=1𝐴^[𝑣]+𝑔𝑓^[𝑣]_𝑛[𝑣]𝑇^(𝑛),𝑦̆^(𝑛),𝑦̆̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)−𝑓^[𝑣]𝑇^(𝑛),𝑦̆^(𝑛),𝑦̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)𝜀_𝑛+𝐻^[𝑣](𝑛)𝑤̆^[𝑣](𝑛)−𝑤^[𝑣](𝑛)+𝑟^𝑛,(3.16) where 𝜀^0(𝑛+1)=𝑦̆_𝑛+1−𝑦_𝑛+1, 𝜀_𝑛=(𝜀^1𝑇(𝑛),𝜀^2𝑇(𝑛),…,𝜀^{𝑠𝑇(𝑛)})^𝑇=𝑦̆^(𝑛)−𝑦^(𝑛),𝑔^{𝑖[𝑣](𝑛)}𝐻=¹⁰𝑓₂𝑡^𝑛+𝑐_𝑖ℎ,𝑦^𝑖(𝑛)+𝜃𝑦̆^𝑖(𝑛)−𝑦^𝑖(𝑛),𝑦̆̃^{𝑖[𝑣](𝑛)},𝑤̆^{𝑖[𝑣](𝑛)}𝑑𝜃,𝑖=1,2,…,𝑠,_{𝑖[𝑣](𝑛)}=¹⁰𝑓₄𝑡^𝑛+𝑐_𝑖ℎ,𝑦^𝑖(𝑛),𝑦̆̃^{𝑖[𝑣](𝑛)},𝑤̆^{𝑖[𝑣](𝑛)}+𝜃𝑤̆^{𝑖[𝑣](𝑛)}−𝑤^{𝑖[𝑣](𝑛)}𝑑𝜃,(3.17) here, 𝑓_𝑖(𝑥₁,𝑥₂,𝑥₃,𝑥₄) is the Jacobian matrix (𝜕𝑓(𝑥₁,𝑥₂,𝑥₃,𝑥₄)/𝜕𝑥_𝑖)𝑖=1,2,3,4.𝐻^[𝑣](𝑛)𝑤̆^[𝑣](𝑛)−𝑤^[𝑣](𝑛)=ℎ𝐻^{𝑚[𝑣](𝑛)𝑣}_𝑞=1+ℎ𝐻_[𝑣](𝑛)𝑑_𝑞𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛−𝑞,𝑦̆^{(𝑛−𝑞)}−𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛−𝑞,𝑦^{(𝑛−𝑞)}𝑅^[𝑣](𝑛)=+ℎ𝐻ℎ𝐻^[𝑣](𝑛)𝑑_{0𝑚[𝑣](𝑛)𝑣}𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛,𝑦̆^(𝑛)−𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛,𝑦^(𝑛)_𝑞=1+ℎ𝐻_[𝑣](𝑛)𝑑_𝑞𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛−𝑞,𝑦̆^{(𝑛−𝑞)}−𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛−𝑞,𝑦^{(𝑛−𝑞)}𝑅^[𝑣](𝑛)+ℎ𝐻^[𝑣](𝑛)𝑑₀¹⁰𝑔^3[𝑣]𝑡_𝑛,𝑡_𝑛,𝑦^(𝑛)+𝜃𝑦̆^(𝑛)−𝑦^(𝑛)𝑑𝜃⋅𝜀_𝑛.(3.18) Assume that (𝐼_𝑠−ℎ∑^𝑚𝑣=1𝐴^[𝑣](𝑔^[𝑣](𝑛)+ℎ𝐻^[𝑣](𝑛)𝑑₀𝑔^{3[𝑣](𝑛)})) is regular, from (3.16) and (3.17), (3.18), we can get𝜀^0(𝑛+1)=⎧⎪⎨⎪⎩𝐼_𝑁+ℎ^𝑚_𝑣=1𝑏̃^[𝑣]𝑇𝐼^𝑠−ℎ^𝑚_𝑣=1𝐴^[𝑣]𝑔^[𝑣](𝑛)+ℎ𝐻^[𝑣](𝑛)𝑑₀𝑔^{3[𝑣](𝑛)}⁻¹×𝑒̃𝑔^[𝑣](𝑛)+ℎ𝐻^[𝑣](𝑛)𝑑₀𝑔^{3[𝑣](𝑛)}⎫⎪⎬⎪⎭𝜀^0(𝑛)+ℎ^𝑚_𝑣=1𝑏̃^[𝑣]𝑇𝑔^[𝑣](𝑛)+ℎ𝐻^[𝑣](𝑛)𝑑₀𝑔^{3[𝑣](𝑛)}𝐼^𝑠−ℎ^𝑚_𝑣=1𝐴^[𝑣]𝑔^[𝑣](𝑛)+ℎ𝐻^[𝑣](𝑛)𝑑₀𝑔^{3[𝑣](𝑛)}⁻¹×𝑟_𝑛+ℎ^2𝑚_𝑣=1𝐴^[𝑣]𝐻^[𝑣](𝑛)𝑅^[𝑣](𝑛)+ℎ_𝑚_𝑣=1𝑏̃^[𝑣]𝑇𝐼^𝑠−ℎ^𝑚_𝑣=1𝐴^[𝑣]𝑔^[𝑣](𝑛)+ℎ𝐻^[𝑣](𝑛)𝑑₀𝑔^{3[𝑣](𝑛)}⁻¹×ℎ^𝑚_𝑣=1𝐴^[𝑣]⋅𝑓_[𝑣]𝑔^[𝑣](𝑛)+ℎ𝐻^[𝑣](𝑛)𝑑₀𝑔^{3[𝑣](𝑛)}+ℎ𝐻𝑇^(𝑛),𝑦̆^(𝑛)_{𝑚[𝑣](𝑛)𝑣},𝑦̆̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)−𝑓^[𝑣]𝑇^(𝑛),𝑦̆^(𝑛),𝑦̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)_𝑞=1+𝑄_𝑛𝑑_𝑞𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛−𝑞,𝑦̆^{(𝑛−𝑞)}−𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛−𝑞,𝑦^{(𝑛−𝑞)}+ℎ^2𝑚_𝑣=1𝑏̃^[𝑣]𝑇𝐻^[𝑣](𝑛)𝑅^[𝑣](𝑛).(3.19) Now, we introduce the concept of D-convergence from [25].

Definition 3.1. An ARKLM (2.1) with𝑦_𝑛=𝑦(𝑡_𝑛)(𝑛≤0), 𝑦^𝑖(𝑛)=𝑦(𝑡_𝑛+𝑐_𝑖ℎ)(𝑛<0) and 𝑦̃^{𝑖[𝑣](𝑛)}=𝑦(𝑡_𝑛+𝑐_𝑖ℎ−𝜏_𝑣)(𝑛<0) is called D-convergence of order 𝑝 if this method, when applied to any given DIDEs (1.1) subject to (1.2)–(1.4); produce an approximation sequence 𝑦_𝑛 and the global error satisfies a bound of the form: ‖‖𝑦𝑡_𝑛−𝑦_𝑛‖‖≤𝐶𝑡_𝑛ℎ^𝑃,ℎ∈0,ℎ₀,(3.20) where the maximum stepsize ℎ₀ depends on characteristic parameter 𝛼_𝑣,𝛽_𝑣,𝜎_𝑣,𝑟_𝑣,𝑟_𝑣 and the method, the function 𝐶(𝑡) depends only on some 𝑀_𝑖 in (3.2), delay 𝜏_𝑣, characteristic parameters 𝛼_𝑣,𝛽_𝑣,𝜎_𝑣,𝑟_𝑣,𝑟_𝑣̃, 𝑣=1,2,…,𝑚, and the method.

Definition 3.2. The ARKLM (2.2), (2.3), (2.5), and (2.6) is said to be DA-stable if the matrix(𝐼_𝑠−∑^𝑚𝜈=1𝐴^[𝜈]𝜉) is regular for 𝜉∈𝐶⁻∶={𝜉∈𝐶∣Re𝜉≤0}, and |𝑅_𝑖(𝜉)|≤1forall𝜉∈𝐶⁻,𝑖=0,1,…,𝑠.
Where 𝑅_𝑖𝜀₁=1+^𝑚_𝑣=1𝐴^𝑖[𝑣]𝜀₁𝐼^𝑠−^𝑚_𝜈=1𝐴𝐴^[𝜈]𝜉⁻¹𝑒,_0[𝑣]=𝑏^[𝑣],𝐴^𝑖[𝑣]=𝑎^[𝑣]𝑖1,𝑎^[𝑣]𝑖2,…,𝑎^{[𝑣]𝑖𝑠}^𝑇,𝑖=0,1,…,𝑠.(3.21)

Definition 3.3. The ARKLM (2.2), (2.3), (2.5), and (2.6) is said to be ASI-stable if the matrix (𝐼_𝑠−∑^𝑀𝜈=1𝐴^[𝜈]𝜉) is regular for𝜉∈𝐶⁻, and (𝐼_𝑠−∑^𝑀𝜈=1𝐴^[𝜈]𝜉)⁻¹ is uniformly bounded for 𝜉∈𝐶⁻.

Definition 3.4. The ARKLM (2.2), (2.3), (2.5), and (2.6) is said to be DAS-stable if the matrix(𝐼_𝑠−∑^𝑀𝜈=1𝐴^[𝜈]𝜉) is regular for 𝜉∈𝐶⁻, and ∑^𝑚𝜈=1𝐴^[𝜈]𝑇_𝑖𝜉(𝐼_𝑠−∑^𝑀𝜈=1𝐴^[𝜈]𝜉)⁻¹(𝑖=0,1,…,s) is uniformly bounded for 𝜉∈𝐶⁻.

Lemma 3.5. Suppose the ARKLM (2.2), (2.3), (2.5), and (2.6) is DA- DAS- and ASI-stable, then there exist positive constants ℎ₀,𝛾₁,𝛾₂,𝛾₃, which depend only on the method and the parameter 𝛼_𝑣,𝛽_𝑣,𝜎_𝑣,𝑟_𝑣,𝑟_𝑣 such that ‖‖‖‖𝐼_𝑠−^𝑀_𝜈=1‖‖‖‖𝐼_𝑁𝐴^[𝜈]𝜉‖‖‖‖≤𝛾₁,+^𝑚_𝑣=1𝐴^[𝑣]𝑇_𝑖𝜉𝐼_𝑠−^𝑚_𝑣=1‖‖‖‖𝐴^[𝑣]𝜉⁻¹𝑒‖‖‖‖̃≤1+𝛾₂ℎ,_𝑚_𝑣=1𝐴^[𝑣]𝑇_𝑖𝜉𝐼_𝑠−^𝑚_𝑣=1𝐴ℎ∈^[𝑣]𝜉0,⁻¹ℎ𝑣₀‖‖‖‖≤𝛾₃‖𝑣‖,𝑣∈𝐶^𝑁𝑆,,𝑖=0,1,2…,𝑠.(3.22)

Proof. This Lemma can be proved in the similar way as that of in [29, Lemmas 3.5–3.7].

Theorem 3.6. Suppose the ARKLM (2.2), (2.3), (2.5), and (2.6) is DA- DAS- and ASI-stable, then there exist positive constants ℎ₀,𝛾₃,𝛾₄,𝛾₅, which depend only on the method and the parameters𝛼_𝑣,𝛽_𝑣,𝜎_𝑣,𝑟_𝑣,𝑟_𝑣̃, such that for ℎ∈(𝑜,ℎ₀], ‖‖𝜀^𝑖(𝑛)‖‖≤⎧⎪⎪⎪⎨⎪⎪⎪⎩1+ℎ𝛾⁴+‖‖𝑄_𝑛−1max‖‖𝜀^0(𝑛+1)‖‖,max_{(𝑖,𝑝𝑣)∈𝐸}‖‖𝜀^{(𝑛−𝑚𝑣+𝑝𝑣)𝑖}‖‖,max_{(𝑖,𝑞)∈𝐸𝑞}‖‖𝜀^{𝑖(𝑛−𝑞)}‖‖+ℎ𝛾₅max_{1≤𝑖≤𝑠}‖‖𝜌^{𝑖(𝑛−1)}‖‖1+ℎ𝛾₄‖‖+𝛾₃‖‖𝛾̃_𝑛−1‖‖,𝑖=0,+‖‖𝑄_𝑛max‖‖𝜀^0(𝑛+1)‖‖,max_{(𝑖,𝑝𝑣)∈𝐸}‖‖𝜀^{(𝑛−𝑚𝑣+𝑝𝑣)𝑖}‖‖,max_{(𝑖,𝑞)∈𝐸𝑞}‖‖𝜀^{𝑖(𝑛−𝑞)}‖‖+ℎ𝛾₅max_{1≤𝑖≤𝑠}‖‖𝜌^𝑖(𝑛)‖‖‖‖+𝛾₃‖‖𝛾̃_𝑛‖‖,𝑖=1,2,…,𝑠,(3.23) where 𝜀^0(𝑛)=𝑦̆_𝑛−𝑦_𝑛,𝜀^𝑖(𝑛)=𝑦̆^𝑖(𝑛)−𝑦^𝑖(𝑛),𝐸={(𝑖,𝑝_𝑣)∣1≤𝑖≤𝑠,−𝑑≤𝑝_𝑣≤𝛾}, 𝐸_𝑞={(𝑖,𝑞)∣1≤𝑖≤𝑠,1≤𝑞≤𝑚},𝑄_𝑛=𝑄_𝑛+ℎ²∑^𝑚𝑣=1𝑏̃^[𝑣]𝑇𝐻^[𝑣](𝑛)𝑅^[𝑣](𝑛), 𝑟_𝑛̃=𝑟_𝑛+ℎ²∑^𝑚𝑣=1𝐴^[𝑣]𝐻^[𝑣](𝑛)𝑅^[𝑣](𝑛).

Proof. Using (3.19) and Lemma 3.5, for ℎ∈(0,ℎ₀], we obtain that 𝜀^0(𝑛+1)≤+1ℎ𝛾+₃𝛾₂ℎ‖‖𝜀^0(𝑛)‖‖+𝛾₃‖‖𝑟_𝑛̃‖‖+‖‖𝑄_𝑛‖‖‖‖‖‖^𝑚_𝑣=1𝐴^[𝑣]𝑓^[𝑣]𝑇^(𝑛)+ℎ𝐻_{𝑚[𝑣](𝑛)𝑣},𝑦̆^(𝑛),𝑦̆̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)−𝑓^[𝑣]𝑇^(𝑛),𝑦̆^(𝑛),𝑦̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)_𝑞=1+ℎ‖‖‖‖_𝑚𝑑_𝑞𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛−𝑞,𝑦̆^{(𝑛−𝑞)}−𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛−𝑞,𝑦^{(𝑛−𝑞)}‖‖‖‖_𝑣=1𝑏̃^[𝑣]𝑇𝑓^[𝑣]𝑇^(𝑛)+ℎ𝐻_[𝑣](𝑛),𝑦̆^(𝑛),𝑦̆̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)−𝑓^[𝑣]𝑇^(𝑛),𝑦̆^(𝑛),𝑦̃^[𝑣](𝑛),𝑤^[𝑣](𝑛)⋅^𝑚𝑣_𝑞=1≤1+𝛾₂𝑑_𝑞𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛−𝑞,𝑦̆^{(𝑛−𝑞)}−𝑔^[𝑣]𝑡_𝑛,𝑡_𝑛−𝑞,𝑦^{(𝑛−𝑞)}‖‖‖‖(3.24)+ℎ𝛾_3𝑚ℎ‖‖𝜀^0(𝑛)‖‖+𝛾₃‖‖𝑟_𝑛̃‖‖+‖‖𝑄_𝑛‖‖(3.25a)_𝑣=1^𝑠_𝑖=1‖‖‖‖^𝑠_𝑗=1+ℎ𝐻𝑎^{[𝑣]𝑖𝑗}_{𝑗𝑚[𝑣](𝑛)𝑉}𝑓^[𝑣]𝑡^𝑗(𝑛),𝑦^𝑗(𝑛),𝑦̆̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}−𝑓^[𝑣]𝑡^𝑗(𝑛),𝑦^𝑗(𝑛),𝑦̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}_𝑞=1𝑑_𝑞𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦̆^{𝑗(𝑛−𝑞)}−𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦^{𝑗(𝑛−𝑞)}‖‖‖‖²⎫⎪⎬⎪⎭^1/2+ℎ_𝑚(3.25b)_𝑣=1‖‖‖‖^𝑠_𝑗=1𝑏^𝑗[𝑣]𝑓^[𝑣]𝑡^𝑗(𝑛)+ℎ𝐻_{𝑗𝑚[𝑣](𝑛)𝑉},𝑦^𝑗(𝑛),𝑦̆̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}−𝑓^[𝑣]𝑡^𝑗(𝑛),𝑦^𝑗(𝑛),𝑦̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}_𝑞=1𝑑_𝑞𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦̆^{𝑗(𝑛−𝑞)}−𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦^{𝑗(𝑛−𝑞)}‖‖‖‖.(3.25c)For ℎ𝛾_3𝑚_𝑣=1^𝑠_𝑖=1‖‖‖‖^𝑠_𝑗=1+ℎ𝐻𝑎^{[𝑣]𝑖𝑗}_{𝑗𝑚[𝑣](𝑛)𝑉}𝑓^[𝑣]𝑡^𝑗(𝑛),𝑦^𝑗(𝑛),𝑦̆̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}−𝑓^[𝑣]𝑡^𝑗(𝑛),𝑦^𝑗(𝑛),𝑦̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}_𝑞=1𝑑_𝑞𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦̆^{𝑗(𝑛−𝑞)}−𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦^{𝑗(𝑛−𝑞)}‖‖‖‖²⎫⎪⎬⎪⎭^1/2=(3.25b).(3.26) Then (3.25b)≤ℎ𝛾_3𝑚_𝑣=1^𝑠_𝑖=12^𝑠_𝑗=1||𝑎^{[𝑣]𝑖𝑗}||²‖‖𝑓^[𝑣]𝑡^𝑗(𝑛)−𝑓_[𝑣],𝑦^𝑗(𝑛),𝑦̆̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}𝑡^𝑗(𝑛),𝑦^𝑗(𝑛),𝑦̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}‖‖²+‖‖‖‖ℎ𝐻^{𝑗𝑚[𝑣](𝑛)𝑉}_𝑞=1𝑑−^𝑞𝑔𝑔^[𝑣]_[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦̆^{𝑗(𝑛−𝑞)}𝑡^𝑛,𝑡_𝑛−𝑞,𝑦^{𝑗(𝑛−𝑞)}‖‖‖‖²^1/2≤ℎ𝛾_3𝑚_𝑣=1^𝑠_𝑖=12^𝑠_𝑗=1||𝑎^{[𝑣]𝑖𝑗}||²𝑟^2𝑣‖‖𝑦̆̃^{𝑗[𝑣](𝑛)}−𝑦̃^{𝑗[𝑣](𝑛)}‖‖²+2^𝑠_𝑗=1||𝑎^{[𝑣]𝑖𝑗}||²ℎ²𝐻^{𝑗2[𝑣](𝑛)}2^𝑚𝑣_𝑞=1𝑑_𝑞⋅‖‖𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦̆^{𝑗(𝑛−𝑞)}−𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦^{𝑗(𝑛−𝑞)}‖‖²^1/2≤ℎ𝛾_3𝑚_𝑣=1⎷_𝑠_𝑖=12^𝑠_𝑗=1||𝑎^{[𝑣]𝑖𝑗}||²𝑟^2𝑣‖‖𝑦̆̃^{𝑗[𝑣](𝑛)}−𝑦̃^{𝑗[𝑣](𝑛)}‖‖²+2ℎ²𝛾_3𝑚_𝑣=1⎷_𝑠_𝑠𝑖=1_𝑗=1||𝑎^{[𝑣]𝑖𝑗}||²𝐻^{𝑗2𝑚[𝑣](𝑛)𝑣}_𝑞=1𝑑^2𝑞⋅‖‖𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦̆^{𝑗(𝑛−𝑞)}−𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦^{𝑗(𝑛−𝑞)}‖‖²≤2ℎ𝛾_3𝑚_𝑠𝑣=1_𝑠𝑖=1_𝑗=1+2ℎ₂||𝑎^{[𝑣]𝑖𝑗}||𝑟_𝑣‖‖𝑦̆̃^{𝑗[𝑣](𝑛)}−𝑦̃^{𝑗[𝑣](𝑛)}‖‖𝛾_3𝑚_𝑠𝑣=1_𝑠𝑖=1_𝑚𝑗=1^𝑣_𝑞=1≤2ℎ𝛾_3𝑚𝑑_𝑞||𝑎^{[𝑣]𝑖𝑗}||||𝐻^{𝑗[𝑣](𝑛)}||⋅‖‖𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦̆^{𝑗(𝑛−𝑞)}−𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦^{𝑗(𝑛−𝑞)}‖‖_𝑠𝑣=1_𝑖,𝑗=1+2ℎ₂||𝑎^{[𝑣]𝑖𝑗}||𝑟_𝑣‖‖𝑦̆̃^{𝑗[𝑣](𝑛)}−𝑦̃^{𝑗[𝑣](𝑛)}‖‖𝛾_3𝑚^𝑣_𝑞=1𝑑_𝑞𝑚_𝑠𝑣=1_𝑖,𝑗=1||𝑎^{[𝑣]𝑖𝑗}||||𝐻^{𝑗[𝑣](𝑛)}||𝑟_𝑣̃‖‖𝑦̆^{𝑗(𝑛−𝑞)}−𝑦^{𝑗(𝑛−𝑞)}‖‖.(3.27) For ℎ^𝑚_𝑣=1‖‖‖‖^𝑠_𝑗=1𝑏^𝑗[𝑣]𝑓^[𝑣]𝑡^𝑗(𝑛)+ℎ𝐻,_{𝑗𝑚[𝑣](𝑛)𝑉}𝑦^𝑗(𝑛),𝑦̆̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}−𝑓^[𝑣]𝑡^𝑗(𝑛),𝑦^𝑗(𝑛),𝑦̃^{𝑗[𝑣](𝑛)},𝑤^{𝑗[𝑣](𝑛)}_𝑞=1𝑑_𝑞−𝑔𝑔^[𝑣]_[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦̆^{𝑗(𝑛−𝑞)}𝑡^𝑛,𝑡_𝑛−𝑞,𝑦^{𝑗(𝑛−𝑞)}‖‖‖‖=(3.25c),(3.28)then (3.25c)≤ℎ^𝑚_𝑠𝑣=1_𝑗=1+ℎ_2𝑚||𝑏^𝑗(𝜈)||𝛾_𝑣‖‖𝑦̆̃^{𝑗[𝑣](𝑛)}−𝑦̃^{𝑗[𝑣](𝑛)}‖‖_𝑠𝑣=1_𝑗=1||𝑏^𝑗[𝑣]||||𝐻^{𝑗[𝑣](𝑛)}||^𝑚𝑣_𝑞=1≤ℎ_𝑚𝑑_𝑞‖‖𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦̆^{𝑗(𝑛−𝑞)}−𝑔^[𝑣]𝑡^𝑛,𝑡_𝑛−𝑞,𝑦^{𝑗(𝑛−𝑞)}‖‖_𝑠𝑣=1_𝑗=1||𝑏^𝑗(𝜈)||𝛾_𝑣‖‖𝑦̆̃^{𝑗[𝑣](𝑛)}−𝑦̃^{𝑗[𝑣](𝑛)}‖‖+ℎ^2𝑚_𝑠𝑣=1_𝑗=1||𝑏^𝑗[𝑣]||||𝐻^{𝑗[𝑣](𝑛)}||^𝑚𝑣_𝑞=1𝑑_𝑞𝑟_𝑣̃‖‖𝑦̆^{𝑗(𝑛−𝑞)}−𝑦^{𝑗(𝑛−𝑞)}‖‖.(3.29) Combine (3.27), (3.29), and (3.25a), we have 𝜀^0(𝑛+1)≤+1ℎ+_𝑚𝛾₂ℎ‖‖𝜀^0(𝑛)‖‖+𝛾₃‖‖𝑟_𝑛̃‖‖+‖‖𝑄_𝑛‖‖_𝑣=1𝛾_𝑣2𝛾^3𝑠_𝑖,𝑗=1||𝑎^{[𝑣]𝑖𝑗}||+^𝑠_𝑗=1+ℎ_2𝑚||𝑏^𝑗[𝑣]||‖‖𝑦̆̃^{𝑗[𝑣](𝑛)}−𝑦̃^{𝑗[𝑣](𝑛)}‖‖_𝑣=1𝑟̃_𝑣𝑚^𝑣_𝑞=1𝑑_𝑞𝑠_𝑗=12𝛾^3𝑠_𝑖=1||𝑎^{[𝑣]𝑖𝑗}||⋅||𝐻^{𝑗[𝑣](𝑛)}||+||𝑏^𝑗[𝑣]||⋅||𝐻^{𝑗[𝑣](𝑛)}||‖‖𝜀^{𝑗(𝑛−𝑞)}‖‖.(3.30) Moreover, it follows from (2.5) and (3.11) that ‖‖𝑦̆̃^{𝑗[𝑣](𝑛)}−𝑦̃^{𝑗[𝑣](𝑛)}‖‖≤sup_{𝛿𝑣∈[𝑟0,1)𝑃}_𝑣=−𝑑||𝐿_𝑝𝑣𝛿_𝑣||max_{−𝑑≤𝑃𝑟≤𝑟}‖‖𝜀^{(𝑛−𝑚𝑣+𝑃𝑣)𝑗}‖‖+‖‖𝜌^{𝑗[𝑣](𝑛)}‖‖.(3.31) Substituting (3.31) in (3.30), we get ‖‖𝜀^0(𝑛+1)‖‖≤1+𝛾⁴⁽⁰⁾ℎ+𝛾⁶⁽⁰⁾ℎ²max‖‖𝜀^0(𝑛)+𝛾₃‖‖,max_{(𝑗,𝑝𝑣)∈𝐸}‖‖𝜀^{(𝑛−𝑚𝑣+𝑃𝑣)𝑗}‖‖,max_{(𝑗.𝑞)∈𝐸𝑞}‖‖𝜀^{𝑗(𝑛−𝑞)}‖‖+‖‖𝑄_𝑛‖‖‖‖𝑟_𝑛̃‖‖+ℎ𝛾⁵⁽⁰⁾max_{(𝑗,𝑣)∈𝐸𝑚}‖‖𝜌^{𝑗[𝑣](𝑛)}‖‖,ℎ∈0,ℎ₀,(3.32) where 𝛾⁴⁽⁰⁾=𝛾₂+𝛾⁵⁽⁰⁾sup_{𝛿𝑣∈[0,1)}^𝑟_𝑃_𝑣=−𝑑||𝐿_𝑝𝑣𝛿_𝑣||,𝛾⁵⁽⁰⁾=^𝑚_𝑣=1𝑟_𝑣2𝛾^3𝑠_𝑖,𝑗=1||𝑎^{[𝑣]𝑖𝑗}||+^𝑠_𝑗=1𝛾₆₍₀₎||𝑏^𝑗[𝑣]||,=^𝑚_𝑣=1𝛾̃_𝑣𝑚^𝑣_𝑠𝑞=1_𝑗=1𝑑_𝑞2𝛾^3𝑠_𝑖=1𝐸=𝑗,𝑃_𝑣||𝑎_{𝑖𝑗[𝑣]}⋅𝐻^{𝑗[𝑣](𝑛)}||+||𝑏^𝑗[𝑣]⋅𝐻^{𝑗[𝑣](𝑛)}||,∣1≤𝑗≤𝑠𝐸,_𝑚−𝑑≤𝑃_𝑣≤𝛾,𝐸_𝑞=(𝑗,𝑞)∣1≤𝑗≤𝑠,1≤𝑞≤𝑚_𝑣,={(𝑗,𝑣)∣1≤𝑗≤𝑠,1≤𝑣≤𝑚}.(3.33) By Lemma 3.5, similar to (3.32), we can obtain the inequalities: ‖‖𝜀^𝑖(𝑛)‖‖≤1+ℎ𝛾^4(𝑖)+ℎ²𝛾^6(𝑖)max‖‖𝜀^0(𝑛)+ℎ𝛾_5(𝑖)‖‖,max_{(𝑗,𝑃𝑣∈𝐸)}‖‖𝜀^{(𝑛−𝑚𝑣+𝑃𝑣)𝑗}‖‖,max_{(𝑗,𝑞)∈𝐸𝑞}‖‖𝜀^{𝑗(𝑛−𝑞)}‖‖max_{(𝑗,𝑣)∈𝐸𝑚}‖‖𝜌^{𝑗[𝑣](𝑛)}‖‖+‖‖𝑄_𝑛‖‖+𝛾₃‖‖𝑟_𝑛̃‖‖,𝑖=1,2,…,𝑠,ℎ∈0,ℎ₀,(3.34) where 𝛾^4(𝑖)=𝛾₂+𝛾^5(𝑖)sup_{𝛿𝑣∈[0,1)}^𝑟_𝑃_𝑣=−𝑑||𝐿_𝑝𝑣𝛿_𝑣||,𝛾^5(𝑖)=^𝑚_𝑣=1𝑟_𝑣2𝛾^3𝑠_𝑖,𝑗=1||𝑎^{[𝑣]𝑖𝑗}||+^𝑠_𝑗=1𝛾_6(𝑖)||𝑎^{[𝑣]𝑖𝑗}||,=^𝑚_𝑣=1𝛾̃_𝑣𝑚^𝑣_𝑠𝑞=1_𝑗=1𝑑_𝑞2𝛾^3𝑠_𝑖=1||𝑎_{𝑖𝑗[𝑣]}⋅𝐻^{𝑗[𝑣](𝑛)}||+||𝑎^{[𝑣]𝑖𝑗}⋅𝐻^{𝑗[𝑣](𝑛)}||.(3.35) Setting 𝛾₄=max_{0≤𝑖≤𝑠}𝛾^4(𝑖), 𝛾₅=max_{0≤𝑖≤𝑠}𝛾^5(𝑖), 𝛾₆=max_{0≤𝑖≤𝑠}𝛾^6(𝑖).
Combining (3.32) with (3.34), we immediately obtain the conclusion of this theorem.
Now, we turn to study the convergence of ARKLM (2.1) for (1.1). It is always assumed that the analytic solution 𝑦(𝑡) of (1.1) is smooth enough on each internal of the form (𝑡₀+(𝑗−1)ℎ,𝑡₀+𝑗ℎ) (𝑗 is a positive integer) as (3.2) defined.

Theorem 3.7. Assume ARKLM (2.1) with stage order 𝑝 is DA-, DAS- and ASI-stable, then the ARKLM (2.1) is D-convergent of order min{𝑝,𝑞+1,𝑠+1}, where 𝑞=𝑑+𝑟.

Proof. By Theorem 3.6, we have for ℎ∈(0,ℎ₀]‖‖𝜀^𝑖(𝑛)‖‖≤⎧⎪⎪⎪⎨⎪⎪⎪⎩1+ℎ𝛾⁴++ℎ²𝛾₆max‖‖𝜀^0(𝑛−1)‖‖,max_{(𝑖,𝑝𝑣)∈𝐸}‖‖𝜀^{(𝑛−1−𝑚𝑣+𝑝𝑣)𝑖}‖‖,max_{(𝑗,𝑞∈𝐸𝑞)}‖‖𝜀^{𝑗(𝑛−𝑞)}‖‖1+T₁ℎℎ𝛾^𝑝+1₄+T₂ℎ^𝑞+2+T₃ℎ^𝑠+2,𝑖=0++ℎ²𝛾₆max‖‖𝜀^0(𝑛)‖‖,max_{(𝑖,𝑝𝑣)∈𝐸}‖‖𝜀^{(𝑛−𝑚𝑣+𝑝𝑣)𝑖}‖‖,max_{(𝑗,𝑞)∈𝐸𝑞}‖‖𝜀^{𝑗(𝑛−𝑞)}‖‖T₁ℎ^𝑝+1+T₂ℎ^𝑞+2+T₃ℎ^𝑠+2,𝑖=1,2,…,𝑠,(3.36) where T₁=𝑀₀+𝛾₃⎷_𝑠_𝑖=1𝑀^2𝑖,T₂=𝛾^𝑠(𝑞+1)!^𝑟_{𝑝𝑣=−𝑑}||𝛿_𝑣−𝑃_𝑣||𝑀_𝑞+1,T₃=^𝑚_𝑣=1𝑏̃^[𝑣]T𝐻^[𝑣]𝑅+^𝑚_𝑣=1𝐴^[𝑣]𝐻^[𝑣]𝑅𝑔^[𝑣](𝑠)(𝜉).(3.37) It follows from an induction to (3.36) for 𝑛 that ‖‖𝜀^𝑖(𝑛)‖‖≤⎧⎪⎪⎨⎪⎪⎩^𝑛_𝑗=01+2ℎ𝛾₄^𝑗T₁ℎ^𝑝+1+T₂ℎ^𝑞+2+T₃ℎ^𝑠+2,𝑖=0,_𝑛+1_𝑗=01+2ℎ𝛾₄^𝑗T₁ℎ^𝑝+1+T₂ℎ^𝑞+2+T₃ℎ^𝑠+2,𝑖=1,2,…,𝑠.(3.38)Hence, forℎ∈(0,ℎ₀], we arrive at ‖‖𝑦𝑡_𝑛−𝑦_𝑛‖‖=‖‖𝜀^0(𝑛)‖‖≤^𝑛_𝑗=01+2ℎ𝛾₄^𝑗=1+2ℎ𝛾₄T₁ℎ^𝑝+1+T₂ℎ^𝑞+2+T₃ℎ^𝑠+2^𝑛+1−12ℎ𝛾₄≤exp2(𝑛+1)ℎ𝛾₄T₁ℎ^𝑝+1+T₂ℎ^𝑞+2+T₃ℎ^𝑠+2−12ℎ𝛾₄≤exp2T₁ℎ^𝑝+1+T₂ℎ^𝑞+2+T₃ℎ^𝑠+2T−𝑡₀𝛾₄exp2ℎ₀𝛾₄−12𝛾₄T₁ℎ^𝑝+T₂ℎ^𝑞+1+T₃ℎ^𝑠+1.(3.39) Therefore, the ARKLM (2.1) is D-Convergent of order min{𝑝,𝑞+1,𝑠+1},(𝑞=𝑟+𝑑).

4. Some Examples

Consider the following initial value problem of multidelay-integro-differential equations:𝑦^(𝑡)=−10⁴+99𝑖[1+𝑦(𝑡)]²1+[1+𝑦(𝑡)]²+500_𝑡𝑡−2+500[𝑦(𝑡−1)−𝑦(𝑡−2)]−500^𝑡𝑡−1𝑦(𝑠)𝑑𝑠𝑦(𝑠)𝑑𝑠−exp(−𝑡)+10⁴−99𝑖×1+exp(−𝑡)²1+1+exp(−𝑡)²𝑦(𝑡)=exp(−𝑡),−2≤𝑡≤0,,0≤𝑡≤3,(4.1) and its perturbed problem:𝑧^(𝑡)=−10⁴+99𝑖[1+𝑧(𝑡)]²1+[1+𝑧(𝑡)]²+500_𝑡𝑡−2+500[𝑧(𝑡−1)−𝑧(𝑡−2)]−500^𝑡𝑡−1𝑧(𝑠)𝑑𝑠𝑧(𝑠)𝑑𝑠−exp(−𝑡)+10⁴−99𝑖×1+exp(−𝑡)²1+1+exp(−𝑡)²𝑧(𝑡)=exp(−𝑡)+0.01,−2≤𝑡≤0.,0≤𝑡≤3,(4.2) It an be easily verified that 𝑎₁=𝑎₂=−9.5×10³, 𝛽₁=𝛽₂=250,𝜎₁=𝜎₂=250, and 𝑟₁̃=𝑟₂̃=1, with analytic solution 𝑦(𝑡)=exp(−𝑡), where the inner product is standard inner product. We apply the two-stages and two-order additive R-K method:111101022101122101(4.3)

(1) to the problem (4.1) and its perturbed problem (4.2). Since the order of the method is 2, we adopt the compound trapezoidal rule for computing the integer part. According to the result of Theorems 2.3 and 3.7, the corresponding method for DIDEs is GDN-stable and D-convergent. We denote the numerical solution of problem (4.1) and (4.2)𝑦_𝑛 and 𝑧_𝑛, where 𝑦_𝑛 and 𝑧_𝑛 are approximations to 𝑦(𝑡_𝑛) and 𝑧(𝑡_𝑛), respectively. The values 𝑦_𝑛 and 𝑧_𝑛 with ℎ=0.1 are listed inFigure 1, (where the abscissa and ordinate denote variable 𝑛 and 𝑦_𝑛, resp.) and Figure 2, (where the abscissa and ordinate denote variable 𝑛 and 𝑧_𝑛, resp.). It is shown in Table 1 that the numerical solutions are toward to the exact solutions as ℎ→0.

Table 1 

Comparison between the numerical solutions for different h and the exact solution.

Figure 1 

Values 𝑦𝑛 with ℎ=0.1.

Figure 2 

Values 𝑧𝑛 with ℎ=0.1.

It is obvious that the corresponding method for MDIDEs is GDN-stable and D-convergent.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (11101109) and the Natural Science Foundationof Hei-long-jiang Province of China (A201107).