Nonlinear Stability and Convergence of Two-Step Runge-Kutta Methods for Volterra Delay Integro-Differential Equations

Abstract

This paper introduces the stability and convergence of two-step Runge-Kutta methods with compound quadrature formula for solving nonlinear Volterra delay integro-differential equations. First, the definitions of (𝑘,𝑙)-algebraically stable and asymptotically stable are introduced; then the asymptotical stability of a (𝑘,𝑙)-algebraically stable two-step Runge-Kutta method with 0<𝑘<1 is proved. For the convergence, the concepts of 𝐷-convergence, diagonally stable, and generalized stage order are firstly introduced; then it is proved by some theorems that if a two-step Runge-Kutta method is algebraically stable and diagonally stable and its generalized stage order is 𝑝, then the method with compound quadrature formula is 𝐷-convergent of order at least min{𝑝,𝜈}, where 𝜈 depends on the compound quadrature formula.

1. Introduction

Volterra delay integro-differential equations (VDIDEs) arise widely in the mathematical modeling of physical and biological phenomena. Significant advances in the theoretical analysis and in the numerical analysis for these problems have been made in the last few decades (see, e.g., [1, 2]). For the case of linear stability and convergence for these equations, numerical time-integration techniques of one-step collocation and Runge-Kutta type were investigated in [3–8]. Linear multistep-based methods were studied in [9–12]. De la Sen and Luo studied the uniform exponential stability of a wide class of linear time-delay systems in [13]; De la Sen considered the stability of impulsive time-varying systems in [14].

For the case of nonlinear stability and convergence, stability results were obtained in [15, 16], where the authors investigated the nonlinear stability of continuous Runge-Kutta methods, discrete Runge-Kutta methods, and backward differentiation (BDF) methods, respectively. However, most of these important results are based on the classical Lipschitz conditions, while the classical Lipschitz conditions are so strong that there are few equations satisfying them. Most of the Volterra delay integro-differential equations satisfy the one-sided Lipschitz condition, while the studies focusing on the stability and convergence of the numerical method for nonlinear VDIDEs based on a one-sided Lipschitz condition have not yet been seen in the literature. By means of a one-sided Lipschitz condition, we will discuss the stability and convergence of two-step Runge-Kutta (TSRK) methods for nonlinear VDIDEs in the present paper.

The paper is organized as follows. In Section 2, a fairly general class of VDIDEs is defined. We present a stability criterion for such problems, which generalizes the criteria in the above references. A class of two-step Runge-Kutta methods is also derived for solving VDIDEs. They are obtained by compound quadrature rules. In Sections 3 and 4, nonlinear stability and convergence of TSRK method for NDDEs are derived and proved. In Section 5 we present some numerical examples in order to illustrate the nonlinear stability and convergence of a two-step Runge-Kutta method. These numerical results show that the new methods are quite effective.

2. A Class of VDIDES and the Two-Step Runge-Kutta Methods

It is the purpose of this paper to investigate the nonlinear stability and convergence properties of the following initial-value problem VDIDEs: 𝑦^′(𝑡)=𝑓𝑡,𝑦(𝑡),𝑦(𝑡−𝜏),^𝑡_𝑡−𝜏𝑔(𝑣),[],𝑦(𝑡)(𝑡)[],𝑡,𝑣,𝑦𝑑𝑣𝑡∈0,Τ=𝜑,𝑡∈−𝜏,0 (1) where 𝑓, 𝑔, and 𝜑 are smooth enough such that (1) has a unique solution 𝑦(𝑡) and 𝜏 is a positive delay term. We assume that there exist some inner product ⟨⋅,⋅⟩ and the induced norm ||⋅|| in 𝐶^𝑁, such that 𝜑∶[−𝜏,0]→𝐶^𝑁 and 𝑓∶[0,Τ]×𝐶^𝑁×𝐶^𝑁×𝐶^𝑁→𝐶^𝑁 satisfy the following conditions: Re𝑓𝑡,𝑢₁,𝑣,𝑤−𝑓𝑡,𝑢₂,𝑣,𝑤,𝑢₁−𝑢₂≤𝛼𝑢₁−𝑢₂²,𝑓𝑡,𝑢,𝑣₁,𝑤₁−𝑓𝑡,𝑢,𝑣₂,𝑤₂≤𝛽𝑣₁−𝑣₂+𝜎𝑤₁−𝑤₂,𝑔𝑡,𝑣,𝑢₁−𝑔𝑡,𝑣,𝑢₂≤𝛾𝑢₁−𝑢₂,(𝑡,𝑣)∈𝐷, (2)(3)(4) for 𝑡∈[0,Τ], 𝐷={(𝑡,𝑣)∶𝑡∈[0,Τ),𝑣∈[𝑡−𝜏,𝑡]}, and for all 𝑢₁,𝑢₂,𝑣₁,𝑣₂,𝑤₁,𝑤₂∈𝐶^𝑁(−𝛼),𝛽,𝜎, and 𝛾 are nonnegative constants.

Throughout this paper, we assume that (1) has unique solution 𝑦(𝑡) and denote the problem class 𝑅(𝛼,𝛽,𝜎,𝛾) that consist problem of type (1) with (2)–(4). In order to make the error analysis feasible, we always assume that (1) has a unique solution 𝑦(𝑡) which is sufficiently differentiable and satisfies 𝑑^𝑖𝑦(𝑡)𝑑𝑡^𝑖≤𝑀_𝑖,𝑖=1,2,…,𝑡∈[−𝜏,𝑇]. (5) For function 𝑔(𝑡,𝜃,𝑦(𝜃)), we make the following assumptions, unless otherwise stated; all its partial derivatives used later exist and satisfy the following: 𝜕^𝑖𝑔𝑡,𝜃,𝑦(𝜃)𝜕𝜃^𝑖≤𝑁_𝑖,𝑡∈[0,𝑇],𝜃∈[−𝜏,𝑇]. (6) Many numerical methods have been proposed for the numerical solution of (1).

In this paper, we are concerned with two-step Runge-Kutta (TSRK) methods of the form𝑌_𝑖^(𝑛)=𝑦_𝑛+ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑌_𝑗^(𝑛)𝑌,𝑖=1,2,…,𝑠,_{𝑖(𝑛−1)}=𝑦_𝑛−1+ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑓𝑡_𝑛−1+𝑐_𝑗ℎ,𝑌_𝑗^(𝑛−1),𝑦𝑖=1,2,…,𝑠,_𝑛+1=(1−𝜃)𝑦_𝑛+𝜃𝑦_𝑛−1+ℎ^𝑠_𝑖=1𝑏_𝑖𝑓𝑡_𝑛+𝑐_𝑖ℎ,𝑌_𝑖^(𝑛)+ℎ^𝑠_𝑖=1𝑏̃_𝑖𝑓𝑡_𝑛−1+𝑐_𝑖ℎ,𝑌_𝑖^(𝑛−1), (7a)(7b)(7c)where ∑^𝑠_𝑖=1𝑏_𝑖+∑^𝑠_𝑖=1𝑏̃_𝑖=1+𝜃, 𝑐_𝑖=∑^𝑠_𝑗=1𝑎_𝑖𝑗, ℎ=𝜏/𝑚 is a stepsize, 𝑚 is an arbitrarily given positive integer, and 0≤𝜃≤1. The above methods are studied in [17]. Now we consider the adaptation of the two-step Runge-Kutta method to (1):𝑌_𝑖^(𝑛)=𝑦_𝑛+ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑌_𝑗^(𝑛),𝑌_𝑗^(𝑛),𝑌_𝑗^(𝑛),𝑌𝑖=1,2,…,𝑠,_{𝑖(𝑛−1)}=𝑦_𝑛−1+ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑓𝑡_𝑛−1+𝑐_𝑗ℎ,𝑌_𝑗^(𝑛−1),𝑌_𝑗^(𝑛−1),𝑌_𝑗^(𝑛−1),𝑦𝑖=1,2,…,𝑠,_𝑛+1=(1−𝜃)𝑦_𝑛+𝜃𝑦_𝑛−1+ℎ^𝑠_𝑖=1𝑏_𝑖𝑓𝑡_𝑛+𝑐_𝑖ℎ,𝑌_𝑖^(𝑛),𝑌_𝑖^(𝑛),𝑌_𝑖^(𝑛)+ℎ^𝑠_𝑖=1𝑏̃_𝑖𝑓𝑡_𝑛−1+𝑐_𝑖ℎ,𝑌_𝑖^(𝑛−1),𝑌_𝑖^(𝑛−1),𝑌_𝑖^(𝑛−1). (8a)(8b)(8c)In particular, 𝑦₀=𝜑(0), 𝑦_𝑛 is the numerical approximation at 𝑡_𝑛=𝑛ℎ to the analytic solution 𝑦(𝑡_𝑛), the argument 𝑌_𝑗^(𝑛)=𝑌_𝑗^{(𝑛−𝑚)} denotes an approximation to 𝑦(𝑡_𝑛+𝑐_𝑗ℎ−𝜏), and the argument 𝑌_𝑗^(𝑛) denotes an approximation to ∫^{𝑡𝑛+𝑐𝑗ℎ}_{𝑡𝑛+𝑐𝑗ℎ−𝜏}𝑔(𝑡_𝑛+𝑐_𝑗ℎ,𝜉,𝑦(𝜉))𝑑𝜉 which are obtained by a convergent compound quadrature (CQ) formula: 𝑌_𝑗^(𝑛)=ℎ^𝑚_𝑞=0𝑝_𝑞𝑔𝑡_𝑛+𝑐_𝑗ℎ,𝑡_𝑛−𝑞+𝑐_𝑗ℎ,𝑌_𝑗^{(𝑛−𝑞)},𝑗=1,2,…,𝑠, (9) using values 𝑌_𝑗^{(𝑛−𝑞)}=𝜑(𝑡_𝑛−𝑞+𝑐_𝑗ℎ) with 𝑛≤𝑞, 𝑡_𝑛−𝑞+𝑐_𝑗ℎ≤0.

The class of Runge-Kutta methods with CQ formula has been applied to delay-integro-differential equations by many authors (c.f. [18, 19]). For the CQ formula (9), we usually adopt the repeated trapezoidal rule, the repeated Simpson’s rule, or the repeated Newton-cotes rule, and so forth, denote 𝜂=max{𝑝₀,𝑝₁,…,𝑝_𝑚}. It should be pointed out that the adopted quadrature formula (9) is only a class of quadrature formula for 𝑌_𝑖^(𝑛), there also exist some other types of quadrature formula, such as Pouzet quadrature (PQ) formula and the quadrature formula based on Laguerre-Radau interpolations [20, 21]. It is the aim of our future research to investigate the adaptation of PQ formula and the quadrature formula based on Laguerre-Radau interpolations to VDIDEs.

3. The Nonlinear Stability Analysis

In this section, we will investigate the stability of the two-step Runge-Kutta methods for VDIDEs. In order to consider the stability property, we also need to consider the perturbed problem of (1): 𝑧^′(𝑡)=𝑓𝑡,𝑧(𝑡),𝑧(𝑡−𝜏),^𝑡_𝑡−𝜏𝑔((𝜉,[],𝑧(𝑡)(𝑡)[],𝑡,𝜉,𝑧))𝑑𝜉𝑡∈0,Τ=𝜓,𝑡∈−𝜏,0 (10) where 𝜓∶[−𝜏,0]→𝐶^𝑁 is a given continuous function. The unique exact solution of the problem (10) is denoted as 𝑧(𝑡).

Applying the two-step Runge-Kutta method (7a)–(7c) to (10) leads to𝑍_𝑖^(𝑛)=𝑧_𝑛+ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑍_𝑗^(𝑛),𝑍_𝑗^(𝑛),𝑍_𝑗^(𝑛),𝑍𝑖=1,2,…,𝑠,_{𝑖(𝑛−1)}=𝑧_𝑛−1+ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑓𝑡_𝑛−1+𝑐_𝑗ℎ,𝑍_𝑗^(𝑛−1),𝑍_𝑗^(𝑛−1),𝑍_𝑗^(𝑛−1),𝑧𝑖=1,2,…,𝑠,_𝑛+1=(1−𝜃)𝑧_𝑛+𝜃𝑧_𝑛−1+ℎ^𝑠_𝑖=1𝑏_𝑖𝑓𝑡_𝑛+𝑐_𝑖ℎ,𝑍_𝑖^(𝑛),𝑍_𝑖^(𝑛),𝑍_𝑖^(𝑛)+ℎ^𝑠_𝑖=1𝑏̃_𝑖𝑓𝑡_𝑛−1+𝑐_𝑖ℎ,𝑍_𝑖^(𝑛−1),𝑍_𝑖^(𝑛−1),𝑍_𝑖^(𝑛−1), (11a)(11b)(11c)where 𝑧_𝑛 and 𝑍_𝑖^(𝑛) denote approximations to 𝑧(𝑡_𝑛) and 𝑧(𝑡_𝑛+𝑐_𝑖ℎ), respectively; the argument 𝑍_𝑗^(𝑛) denotes an approximation to ∫^{𝑡𝑛+𝑐𝑗ℎ}_{𝑡𝑛+𝑐𝑗ℎ−𝜏}𝑔(𝑡_𝑛+𝑐_𝑗ℎ,𝜉,𝑧(𝜉))𝑑𝜉 which are obtained by a convergent compound quadrature (CQ) formula: 𝑍_𝑖^(𝑛)=ℎ^𝑚_𝑞=0𝑝_𝑞𝑔𝑡_𝑛+𝑐_𝑖ℎ,𝑡_𝑛−𝑞+𝑐_𝑖ℎ,𝑍_𝑖^{(𝑛−𝑞)},𝑖=1,2,…,𝑠, (12) and 𝑍_𝑗^(𝑛)=𝑍_𝑗^{(𝑛−𝑚)} denotes an approximation to 𝑧(𝑡_𝑛+𝑐_𝑗ℎ−𝜏), and 𝑧_𝑛=𝜓(𝑡_𝑛) for 𝑛=0, 𝑍_𝑖^{(𝑛−𝑞)}=𝜓𝑡_𝑛−𝑞+𝑐_𝑖ℎ,for𝑡_𝑛−𝑞+𝑐_𝑖ℎ≤0. (13)

3.1. Some Concepts

For the stability analysis, we need the compound quadrature formula (9) that satisfies the following condition: ℎ²(𝑚+1)^2𝑚_𝑞=0𝑝_𝑞²<𝑝², (14) where 𝑝 is a positive constant.

Let 𝑤_𝑛=𝑦_𝑛−𝑧_𝑛,𝑊_𝑗^(𝑛)=𝑌_𝑗^(𝑛)−𝑍_𝑗^(𝑛),𝑊_𝑗^(𝑛)=𝑌_𝑗^(𝑛)−𝑍_𝑗^(𝑛),𝑄_𝑗^(𝑛)=ℎ𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑌_𝑗^(𝑛),𝑌_𝑗^(𝑛),𝑌_𝑗^(𝑛)−𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑍_𝑗^(𝑛),𝑍_𝑗^(𝑛),𝑍_𝑗^(𝑛),𝑗=1,2,…,𝑠. (15) It follows from (8a)–(8c) and (11a)–(11c) that𝑊_𝑖^(𝑛)=𝑤_𝑛+^𝑠_𝑗=1𝑎_𝑖𝑗𝑄_𝑗^(𝑛)𝑊,𝑖=1,2,…,𝑠,_{𝑖(𝑛−1)}=𝑤_𝑛−1+^𝑠_𝑗=1𝑎_𝑖𝑗𝑄_𝑗^(𝑛−1)𝑤,𝑖=1,2,…,𝑠,_𝑛+1=(1−𝜃)𝑤_𝑛+𝜃𝑤_𝑛−1+^𝑠_𝑖=1𝑏_𝑖𝑄_𝑖^(𝑛)+^𝑠_𝑖=1𝑏̃_𝑖𝑄_𝑖^(𝑛−1). (16a)(16b)(16c)Now we will write the 𝑠-stage TSRK methods (7a)–(7c) as a general linear method.

Let 𝑉_𝑖^(𝑛)=(𝑊_𝑖^Τ(𝑛),𝑊_𝑖^Τ(𝑛−1))^Τ be the internal stages and 𝜇_𝑛=(𝑤^Τ_𝑛+1,𝑤^Τ_𝑛)^Τ the external vectors and 𝑃_𝑖^(𝑛)=(𝑄_𝑖^Τ(𝑛),𝑄_𝑖^Τ(𝑛−1))^Τ. Then we have a 2(𝑠+1)-stage partitioned general linear method:𝑉_𝑖^(𝑛)=^𝑠_𝑗=1𝐶¹¹_𝑖𝑗𝑃_𝑖^(𝑛)+^𝑠_𝑗=1𝐶¹²_𝑖𝑗𝜇_𝑛𝜇,𝑖=1,2,…,𝑠,_𝑛+1=^𝑠_𝑗=1𝐶²¹_𝑖𝑗𝑃_𝑖^(𝑛)+^𝑠_𝑗=1𝐶²²_𝑖𝑗𝜇_𝑛,𝑖=1,2,…,𝑠, (17a)(17b)where 𝐶¹¹_𝑖𝑗=𝑎_𝑖𝑗,𝑖,𝑗=1,2,⋯,𝑠, 𝐶¹¹_𝑖𝑗=0,𝑖=1,2,⋯,𝑠,𝑗=𝑠+1,𝑠+2,⋯,2𝑠, 𝐶¹¹_𝑖𝑗=0, 𝑖=𝑠+1,𝑠+2,⋯,2𝑠, 𝑗=1,2,⋯,𝑠, 𝐶¹¹_𝑖𝑗=𝑎_𝑖𝑗, 𝑖,𝑗=𝑠+1,𝑠+2,⋯,2𝑠, 𝐶¹²_𝑖𝑗=1, 𝑖,𝑗=1,2,⋯,𝑠, 𝐶¹²_𝑖𝑗=0, 𝑖=1,2,⋯,𝑠, 𝑗=𝑠+1,𝑠+2,⋯,2𝑠, 𝐶¹²_𝑖𝑗=0, 𝑖=𝑠+1,𝑠+2,⋯,2𝑠, 𝑗=1,2,⋯,𝑠, 𝐶¹²_𝑖𝑗=1, 𝑖,𝑗=𝑠+1,𝑠+2,⋯,2𝑠, 𝐶²¹_𝑖𝑗=𝑏_𝑗, 𝑖,𝑗=1,2,⋯,𝑠, 𝐶²¹_𝑖𝑗=𝑏̃_𝑗, 𝑖=1,2,⋯,𝑠, 𝑗=𝑠+1,𝑠+2,⋯,2𝑠, 𝐶²¹_𝑖𝑗=0, 𝑖=𝑠+1,𝑠+2,⋯,2𝑠, 𝑗=1,2,⋯,𝑠, 𝐶²¹_𝑖𝑗=0, 𝑖,𝑗=𝑠+1,𝑠+2,⋯,2𝑠, 𝐶²²_𝑖𝑗=1−𝜃, 𝑖,𝑗=1,2,⋯,𝑠, 𝐶²²_𝑖𝑗=𝜃, 𝑖=1,2,⋯,𝑠, 𝑗=𝑠+1,𝑠+2,⋯,2𝑠, 𝐶²²_𝑖𝑗=1, 𝑖=𝑠+1,𝑠+2,⋯,2𝑠, 𝑗=1,2,⋯,𝑠, 𝐶²²_𝑖𝑗=0, 𝑖,𝑗=𝑠+1,𝑠+2,⋯,2𝑠, and 𝐴=(𝑎_𝑖𝑗), 𝑒=(1,1,…1), 𝑏=(𝑏₁,𝑏₂,…,𝑏_𝑠), 𝑏̃=(𝑏̃₁,𝑏̃₂,…,𝑏̃_𝑠).

We can identify the coefficient matrices 𝐶₁₁=𝐶¹¹_𝑖𝑗=𝐴00𝐴,𝐶₁₂=𝐶¹²_𝑖𝑗=,𝐶𝑒00𝑒₂₁=𝐶²¹_𝑖𝑗=𝑏𝑏̃00,𝐶₂₂=𝐶²²_𝑖𝑗=.1−𝜃𝜃10 (18)

Definition 1 (see [22]). Let 𝑘,𝑙 be real constants, a TSRK method is said to be (𝑘,𝑙) algebraically stable if there exists a diagonal matrix 𝐷=diag(𝑑₁,𝑑₂,…,𝑑_2𝑠) and a diagonal nonnegative matrix 𝐺 such that 𝑀=(𝑚_𝑖𝑗) is nonnegative, where𝑀(𝑘,𝑙)=𝑘𝐺−𝐶^Τ₂₂𝐺𝐶₂₂−2𝑙𝐶^Τ₁₂𝐷𝐶₁₂𝐶^Τ₁₂𝐷−𝐶^Τ₂₂𝐺𝐶₂₁−2𝑙𝐶^Τ₁₂𝐷𝐶₁₁𝐷𝐶₁₂−𝐶^Τ₂₁𝐺𝐶₂₂−2𝑙𝐶^Τ₁₁𝐷𝐶₁₂𝐶^Τ₁₁𝐷+𝐷𝐶₁₁−𝐶^Τ₂₁𝐺𝐶₂₁−2𝑙𝐶^Τ₁₁𝐷𝐶₁₁. (19)In particular, a (1,0)-algebraically stable method is called algebraically stable.

Definition 2. The TSRK method (7a)–(7c) is called asymptotically stable for (1) with (2)–(4) if lim_𝑛→∞‖𝑤_𝑛‖=0.

3.2. Numerical Stability of the Methods

Numerical stability is an important feature of an effective numerical method. An unstable numerical method may be consistent of high order, yet arbitrarily small perturbations will eventually cause large deviations from the true solution. In this section, we will focus on the asymptotic stability of the TSRK method.

Theorem 3. Assume that the TSRK method (7a)–(7c) is (𝑘,𝑙)-algebraically stable with 0<𝑘≤1, suppose that the quadrature formula (9) satisfies the condition (14) and the conditions (2)–(4) hold, then, method (17a) and (17b) satisfies the following: 𝑦_𝑛−𝑧_𝑛≤𝐶max_{𝑡∈[−𝜏,0]}𝜑(𝑡)−𝜓(𝑡),𝑛=1,2,…, (20) when the following condition holds: ℎ𝛼+𝛽+𝜎+𝜎𝛾²𝑝²≤𝑙, (21) where 𝐶 depends only on the method, the parameters 𝛼, 𝛽, 𝜎, 𝛾, 𝜂, and 𝜏.

Proof. It follows from a fairly straightforward (but tedious) computation and (𝑘,𝑙) algebraically stability that (compare also [22, 23]) 𝜇_𝑛+1²−𝑘𝜇_𝑛²−2^2𝑠_𝑗=1𝑑_𝑗Re𝑉_𝑗^(𝑛),𝑃_𝑗^(𝑛)−𝑙𝑉_𝑗^(𝑛)=−^2𝑠+2_𝑖=1^2𝑠+2_𝑗=1𝑀_𝑖𝑗𝑟_𝑖,𝑟_𝑗, (22) where 𝑟₁=𝑤_𝑛+1,𝑟₂=𝑤_𝑛,𝑟_𝑗=𝑄^(𝑛)_𝑗−2𝑟,𝑗=3,4,…,𝑠+2,_𝑗=𝑄^(𝑛−1)_{𝑗−𝑠−2},𝑗=𝑠+3,𝑠+4,…,2𝑠+2. (23) By means of (𝑘,𝑙) algebraical stability of the method, we have the following: 𝜇_𝑛+1²≤𝑘𝜇_𝑛²+2^2𝑠_𝑗=1𝑑_𝑗Re𝑉_𝑗^(𝑛),𝑃_𝑗^(𝑛)−𝑙𝑉_𝑗^(𝑛). (24) It follows from (2)–(4) that 2Re𝑊_𝑗^(𝑛),𝑄_𝑗^(𝑛)=2ℎRe𝑊_𝑗^(𝑛),𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑌_𝑗^(𝑛),𝑌_𝑗^(𝑛),𝑌_𝑗^(𝑛)−𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑍_𝑗^(𝑛),𝑍_𝑗^(𝑛),𝑍_𝑗^(𝑛)=2ℎRe𝑊_𝑗^(𝑛),𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑌_𝑗^(𝑛),𝑌_𝑗^(𝑛),𝑌_𝑗^(𝑛)−𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑍_𝑗^(𝑛),𝑌_𝑗^(𝑛),𝑌_𝑗^(𝑛)+2ℎRe𝑊_𝑗^(𝑛),𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑍_𝑗^(𝑛),𝑌_𝑗^(𝑛),𝑌_𝑗^(𝑛)−𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑍_𝑗^(𝑛),𝑍_𝑗^(𝑛),𝑍_𝑗^(𝑛)≤2ℎ𝛼𝑊_𝑗^(𝑛)2+2ℎ𝛽𝑊_𝑗^(𝑛)⋅𝑊_𝑗^(𝑛)+2ℎ𝜎𝑊_𝑗^(𝑛)⋅𝑊_𝑗^(𝑛)≤2ℎ𝛼𝑊_𝑗^(𝑛)2+ℎ𝛽𝑊_𝑗^(𝑛)2+𝑊_𝑗^(𝑛)2+ℎ𝜎𝑊_𝑗^(𝑛)2+𝑊_𝑗^(𝑛)2=2ℎ𝛼+ℎ𝛽+ℎ𝜎𝑊_𝑗^(𝑛)2+ℎ𝛽𝑊_𝑗^(𝑛)2+ℎ𝜎𝑊_𝑗^(𝑛)2=2ℎ𝛼+ℎ𝛽+ℎ𝜎𝑊_𝑗^(𝑛)2+ℎ𝛽𝑊_𝑗^(𝑛)2+ℎ𝜎ℎ^𝑚_𝑞=0𝑝_𝑞𝑔𝑡_𝑛+𝑐_𝑗ℎ,𝑡_𝑛−𝑞+𝑐_𝑗ℎ,𝑌_𝑗^{(𝑛−𝑞)}−ℎ^𝑚_𝑞=0𝑝_𝑞𝑔𝑡_𝑛+𝑐_𝑗ℎ,𝑡_𝑛−𝑞+𝑐_𝑗ℎ,𝑍_𝑗^{(𝑛−𝑞)2}≤2ℎ𝛼+ℎ𝛽+ℎ𝜎𝑊_𝑗^(𝑛)2+ℎ𝛽𝑊_𝑗^{(𝑛−𝑚)2}+ℎ𝜎(𝑚+1)ℎ²𝑝²_𝑞×^𝑚_𝑞=0𝑔𝑡_𝑛+𝑐_𝑗ℎ,𝑡_𝑛−𝑞+𝑐_𝑗ℎ,𝑌_𝑗^{(𝑛−𝑞)}−𝑔𝑡_𝑛+𝑐_𝑗ℎ,𝑡_𝑛−𝑞+𝑐_𝑗ℎ,𝑍_𝑗^{(𝑛−𝑞)2}≤2ℎ𝛼+ℎ𝛽+ℎ𝜎𝑊_𝑗^(𝑛)2+ℎ𝛽𝑊_𝑗^{(𝑛−𝑚)2}+ℎ𝜎(𝑚+1)ℎ²𝛾²×^𝑚_𝑞=0𝑝²_𝑞⋅𝑊_𝑗^{(𝑛−𝑞)2}. (25) Substituting (25) into (24), using (14) we get the following: 𝜇_𝑛+1²≤𝜇_𝑛²+^2𝑠_𝑗=1𝑑_𝑗2ℎ𝛼+ℎ𝛽+ℎ𝜎−2𝑙𝑊_𝑗^(𝑛)2+ℎ𝛽𝑊_𝑗^{(𝑛−𝑚)}+ℎ𝜎𝛾²𝑝^2𝑚_𝑞=0𝑊_𝑗^{(𝑛−𝑞)2}. (26) By induction, we get the following: 𝜇_𝑛+1²≤𝜇₋₁²+^𝑛_𝑖=−1^2𝑠_𝑗=1𝑑_𝑗2ℎ𝛼+2ℎ𝛽+ℎ𝜎−2𝑙+ℎ𝜎𝛾²𝑝²𝑊_𝑗^(𝑖)2+⁻²_𝑖=−𝑚^2𝑠_𝑗=1𝑑_𝑗ℎ𝛽+ℎ𝜎𝛾²𝑝²𝑊_𝑗^(𝑖)2≤𝜇₋₁²+2^𝑛_𝑖=−1^2𝑠_𝑗=1𝑑_𝑗ℎ𝛼+𝛽+𝜎+𝜎𝛾²𝑝²−𝑙𝑊_𝑗^(𝑖)2+⁻²_{𝑖=−𝑚−1}^2𝑠_𝑗=1𝑑_𝑗ℎ𝛽+ℎ𝜎𝛾²𝑝²𝑊_𝑗^(𝑖)2≤𝜇₋₁²+𝑚ℎ^2𝑠_𝑗=1𝑑_𝑗𝛽+𝜎𝛾²𝑝²max_{−𝑚−1≤𝑖≤−2}𝑊_𝑗^(𝑖)2≤𝑤₀²+𝑤₋₁²+𝑚ℎ^2𝑠_𝑗=1𝑑_𝑗𝛽+𝜎𝛾²𝑝²max_{−𝑚−1≤𝑖≤−2}𝑊_𝑗^(𝑖)2≤2+𝜏𝛽+𝜎𝛾²𝑝^22𝑠_𝑗=1𝑑_𝑗max_{[]𝑡∈−𝜏,0}𝜑(𝑡)−𝜓(𝑡)². (27) Hence, ‖𝜇_𝑛+1‖≤𝐶max_{𝑡∈[−𝜏,0]}‖𝜑(𝑡)−𝜓(𝑡)‖, where 𝐶=2+𝜏(𝛽+𝜎𝛾²𝑝²)∑^2𝑠_𝑗=1𝑑_𝑗.
The proof of Theorem 3 is completed.

Theorem 4. Assume that a TSRK method (7a)–(7c) is (𝑘,𝑙)-algebraically stable with 0<𝑘<1, then the TSRK method (7a)–(7c) with (9), (12) and (14) is asymptotically stable for the problem (1) with (2)–(4), when the following condition holds: ℎ1𝛼+𝛽+21𝜎+2𝜎𝛾²𝑝²<𝑙. (28)

Proof. Let 𝑢=2𝛼+𝛽+𝜎ℎ−2𝑙,𝑘=max𝑘,𝛽+𝜎𝛾²𝑝²ℎ(−𝑢)^1/𝑚. (29) Then when ℎ(𝛼+𝛽+(1/2)𝜎+(1/2)𝜎𝛾²𝑝²)<𝑙, we have 𝑢<−(𝛽+𝜎𝛾²𝑝²)ℎ and 0<𝑘<1.
The application of Theorem 3 yields the following: 𝜇_𝑛+1²≤𝑘𝜇_𝑛²+^2𝑠_𝑗=1𝑑_𝑗2ℎ𝛼+ℎ𝛽+ℎ𝜎−2𝑙𝑊_𝑗^(𝑛)2+ℎ𝛽𝑊_𝑗^{(𝑛−𝑚)}+ℎ𝜎𝛾²𝑝^2𝑚_𝑞=0𝑊_𝑗^{(𝑛−𝑞)2}≤𝑘𝜇_𝑛²+^2𝑠_𝑗=1𝑑_𝑗𝑢𝑊_𝑗^(𝑛)2+ℎ𝛽𝑊_𝑗^{(𝑛−𝑚)}+ℎ𝜎𝛾²𝑝^2𝑚_𝑞=0𝑊_𝑗^{(𝑛−𝑞)2}. (30) By induction, we get the following: 𝜇_𝑛+1²≤𝑘^𝑛+2𝜇₋₁²+^𝑛_𝑖=−1𝑘^{𝑛−𝑚−𝑖}×^2𝑠_𝑗=1𝑑_𝑗𝑢𝑘^𝑚+ℎ𝛽+ℎ𝜎𝛾²𝑝²⋅𝑊_𝑗^(𝑖)2+⁻²_{𝑖=−𝑚−1}^2𝑠_𝑗=1𝑑_𝑗𝑘^{𝑛−𝑚−𝑖}ℎ𝛽+ℎ𝜎𝛾²𝑝²⋅𝑊_𝑗^(𝑖)2≤𝑘^𝑛+2𝜇₋₁²+𝑘^{𝑛−𝑚+2}⋅⋅𝑚ℎ_2𝑠𝑗=1𝑑_𝑗𝛽+𝜎𝛾²𝑝²⋅max_{−𝑚−1≤𝑖≤−2}𝑊_𝑗^(𝑖)2≤𝑘^𝑛+2𝑤₋₁²+𝑤₀²+𝑘^{𝑛−𝑚+2}⋅𝜏⋅^2𝑠_𝑗=1𝑑_𝑗𝛽+𝜎𝛾²𝑝²⋅max_{−𝑚−1≤𝑖≤−2}𝑊_𝑗^(𝑖)2≤2𝑘^𝑛+2+𝑘^{𝑛−𝑚+2}⋅𝜏⋅^2𝑠_𝑗=1𝑑_𝑗𝛽+𝜎𝛾²𝑝²⋅max_{[]𝑡∈−𝜏,0}𝜑(𝑡)−𝜓(𝑡)². (31)(32) The inequality together with the knowledge 0<𝑘<1 leads to the following: lim_𝑛→∞𝜇_𝑛+1=0. (33) Because 𝜇_𝑛=(𝑤_𝑛^Τ,𝑤^Τ_𝑛−1)^Τ, we can get from (33) that lim_𝑛→∞𝑤_𝑛=0. (34) The proof of the theorem is completed.

4. The Convergence of TSRK Method for NDDEs

4.1. Some Concepts

In order to study the convergence of the method, we define the following: 𝑌^(𝑛)=𝑌₁^Τ(𝑛),…,𝑌_𝑠^Τ(𝑛),𝑌₁^Τ(𝑛−1),…,𝑌_𝑠^{Τ(𝑛−1)Τ}∈𝐶^2𝑠𝑁,𝑌^(𝑛)=𝑌₁^Τ(𝑛),…,𝑌_𝑠^Τ(𝑛),𝑌₁^Τ(𝑛−1),…,𝑌_𝑠^{Τ(𝑛−1)Τ}∈𝐶^2𝑠𝑁,𝑌^(𝑛)=𝑌₁^Τ(𝑛),…,𝑌_𝑠^Τ(𝑛),𝑌₁^Τ(𝑛−1),…,𝑌_𝑠^{Τ(𝑛−1)Τ}∈𝐶^2𝑠𝑁,𝜁^(𝑛+1)=𝑦^Τ(𝑛+1),𝑦^Τ(𝑛)Τ∈𝐶^2𝑁,𝐹𝑡_𝑛,𝑌^(𝑛),𝑌^(𝑛),𝑌^(𝑛)=𝑓𝑡_𝑛+𝑐₁ℎ,𝑌₁^(𝑛),𝑌₁^(𝑛),𝑌₁^(𝑛)Τ𝑓𝑡,…,_𝑛+𝑐_𝑠ℎ,𝑌_𝑠^(𝑛),𝑌_𝑠^(𝑛),𝑌_𝑠^(𝑛)Τ,𝑓𝑡_𝑛−1+𝑐₁ℎ,𝑌₁^(𝑛−1),𝑌₁^(𝑛−1),𝑌₁^(𝑛−1)Τ𝑓𝑡,…,_𝑛−1+𝑐_𝑠ℎ,𝑌_𝑠^(𝑛−1),𝑌_𝑠^(𝑛−1),𝑌_𝑠^{(𝑛−1)ΤΤ}∈𝐶^2𝑠𝑁. (35) Thus, process (8a)–(8c) can be written in the more compact form:𝑌^(𝑛)=ℎ𝐶₁₁𝐹𝑡_𝑛,𝑌^(𝑛),𝑌^(𝑛),𝑌^(𝑛)+𝐶₁₂𝜁^(𝑛),𝜁^(𝑛+1)=ℎ𝐶₂₁𝐹𝑡_𝑛,𝑌^(𝑛),𝑌^(𝑛),𝑌^(𝑛)+𝐶₂₂𝜁^(𝑛). (36a)(36b)

Definition 5. Method (7a)–(7c) with an approximation procedure (9) is said to be 𝐷 convergent of order 𝑝 if the global error satisfies a bound of the form 𝐻𝑡_𝑛−𝜁^(𝑛)≤𝜌₁𝑡_𝑛𝐻𝑡₀−𝜁⁽⁰⁾+max_𝑘≤0𝑌𝑡_𝑘−𝑌^(𝑘)+ℎ^𝑝,𝑛≥1,ℎ∈0,ℎ₀, (37) where 𝐻(𝑡) and 𝑌(𝑡) are defined by 𝐻(𝑡)=𝑦(𝑡+ℎ),𝑦(𝑡)∈𝐶^2𝑁,𝑌(𝑡)=𝑦𝑡+𝑐₁ℎ,…,𝑦𝑡+𝑐_𝑠ℎ,𝑦𝑡−ℎ+𝑐₁ℎ𝑦,…,𝑡−ℎ+𝑐_𝑠ℎ∈𝐶^2𝑠𝑁, (38) where 𝜌₁(𝑡) and ℎ₀ depend on 𝑀_𝑖, 𝛼,𝛽,𝛾,𝜎, and 𝜏.

Definition 6. TSRK Method (7a)–(7c) is said to be diagonally stable if there exist a 2𝑠×2𝑠 diagonal matrix 𝑄>0 such that the matrix 𝑄𝐶₁₁+𝐶^Τ₁₁𝑄 is positive definite.

Remark 7. The concepts of algebraic stability and diagonal stability of TSRK method are the generalizations of corresponding concepts of Runge-Kutta methods. Although it is difficult to examine these conditions, many results have been found, especially, there exist algebraically stable and diagonally stable multistep formulas of arbitrarily high order (cf. [24]).

Definition 8. TSRK Method (7a)–(7c) is said to have generalized stage order 𝑝 if 𝑝 is the largest integer which possesses the following properties.
For any given problem (1) and ℎ∈[0,ℎ₀], there exists an abstract function 𝐻^ℎ(𝑡), 𝐻^ℎ(𝑡)=𝐻^ℎ₁(𝑡),𝐻^ℎ₂(𝑡)∈𝐶^2𝑁 (39) such that 𝐻(𝑡)=𝐻^ℎ(𝑡)≤𝑝₁ℎ^𝑝,Δ^ℎ(𝑡)≤𝑝₁ℎ^𝑝+1,𝛿^ℎ(𝑡)≤𝑝₁ℎ^𝑝+1, (40) where the maximum stepsize ℎ₀>0 and the constant 𝑝₁ depend only on the method and the bounds 𝑀_𝑖, Δ^ℎ(𝑡) and 𝛿^ℎ(𝑡), they are defined by the following equations:𝑌(𝑡)=ℎ𝐶₁₁𝑌^′(𝑡)+𝐶₁₂𝐻^ℎ(𝑡)+Δ^ℎ(𝑡),𝐻^ℎ(𝑡+ℎ)=ℎ𝐶₂₁𝑌^′(𝑡)+𝐶₂₂𝐻^ℎ(𝑡)+𝛿^ℎ(𝑡). (41a)(41b) The function 𝑌^′(𝑡) is defined by the following: 𝑌(𝑡)=𝑦^′𝑡+𝑐₁ℎ,…,𝑦^′𝑡+𝑐_𝑠ℎ,𝑦^′𝑡−ℎ+𝑐₁ℎ𝑦,…,_′𝑡−ℎ+𝑐_𝑠ℎ∈𝐶^2𝑠𝑁. (41c)Particularly, when 𝐻(𝑡)=𝐻^ℎ(𝑡), generalized stage order is called stage order.

4.2. 𝐷-Convergence and Proofs

In this section, we focus on the error analysis of TSRK method for (1). For the sake of simplicity, we always assume that all constants ℎ_𝑖, 𝛾_𝑖, 𝑑_𝑖, and 𝐿_𝑖 are dependent on the method, the bounds 𝑀_𝑖, 𝑁_𝑖, the parameters 𝛼,𝛽,𝛾,𝜎, and 𝜏.

First, we give a preliminary result which will later be used several times. To simplify, we denote 𝑌,𝑌,𝑌,𝑍,𝑍,𝑍∈𝐶^2𝑠𝑁 for 𝑌^(𝑛),𝑌^(𝑛),𝑌^(𝑛),𝑍^(𝑛),𝑍^(𝑛)  𝑍^(𝑛), and 𝜁,𝜛∈𝐶^2𝑁 for 𝜁^(𝑛+1), 𝜛^(𝑛+1), where 𝜛^(𝑛+1)=(𝑧^(𝑛+1),𝑧^(𝑛)). Define Δ and 𝛿 by the following:Δ=𝑌−𝑍−ℎ𝐶₁₁𝐹𝑡,𝑌,𝑌,𝑌−𝐹𝑡,𝑍,𝑍,𝑍,𝛿=𝜁−𝜛−ℎ𝐶₂₁𝐹𝑡,𝑌,𝑌,𝑌−𝐹𝑡,𝑍,𝑍,𝑍. (42a)(42b)

Theorem 9. Suppose method (7a)–(7c) is diagonally stable, then there exist constants ℎ₁, 𝑝₂, and 𝑝₃ such that‖𝑌−𝑍‖≤𝑝₂Δ+ℎ𝑌−𝑍+ℎ𝑌−𝑍,ℎ∈0,ℎ₁,𝜁−𝜛≤𝑝₃Δ+𝛿+ℎ𝑌−𝑍+ℎ𝑌−𝑍,ℎ∈0,ℎ₁. (43a)(43b)

Proof. Since the method (7a)–(7c) is diagonally stable, there exists a positive definite diagonal matrix 𝑄 such that the matrix 𝐸=𝑄𝐶₁₁+𝐶^Τ₁₁𝑄 is positive definite. Therefore, the matrix 𝐶₁₁ is obviously nonsingular and there exists a 𝑙>0 depends only on the method such that the matrix 𝐸_𝑙=𝐶^−Τ₁₁𝐸𝐶₁₁−2𝑙𝑄 (44) is also positive definite.
Define𝑊=𝑌−𝑍,𝑊=𝑌−𝐾𝑍,𝑊=𝑌−𝑍,₁=ℎ𝐹𝑡,𝑌,𝑌,𝑌−𝐹𝑡,𝑍,𝑌𝐾𝑌,,₂=ℎ𝐹𝑡,𝑍,𝑌,𝑌−𝐹𝑡,𝑍,𝑌𝐾𝑍,,₃=ℎ𝐹𝑡,𝑍,𝑍,𝑌−𝐹𝑡,𝑍,𝑍,𝑍; (45a)(45b)(45c)(45d)thenΔ=𝑊−𝐶₁₁𝐾₁+𝐾₂+𝐾₃,𝛿=𝜁−𝜛−𝐶₂₁𝐾₁+𝐾₂+𝐾₃. (46a)(46b)Using (2)–(4), (44), (46a), and (46b), we have for ℎ𝛼≤𝑙, 0≤⟨𝑊,2𝑙𝑄𝑊⟩+2Re𝐾₁,−𝑄𝑊=−𝑊,𝐸_𝑙𝑊+2Re𝑊,𝑄𝐶⁻¹₁₁𝑊−𝐾₁≤−𝜆₁‖𝑊‖²+2Re𝑊,𝑄𝐶⁻¹₁₁Δ+𝐾₂+𝐾₃≤−𝜆_𝑙‖𝑊‖²+2𝑄𝐶⁻¹₁₁𝑊‖⋅Δ‖𝑄‖⋅‖𝑊‖⋅‖+2ℎ𝛽𝑊+2ℎ𝜎‖𝑄‖⋅‖𝑊‖⋅𝑊, (47) where 𝜆_l is the minimum eigenvalue of 𝐸_𝑙. Therefore, ‖𝑌−𝑍‖≤𝑝₂Δ+ℎ𝑌−𝑍+ℎ𝑌−𝑍,ℎ∈0,ℎ₁, (48) where𝑝₂=2𝜆₁max𝑄𝐶⁻¹₁₁‖𝑄‖‖𝑄,ℎ,𝛽,𝜎‖₁=𝑙1𝛼≤0,min1,𝛼𝛼>0. (49a)(49b)From (43a), (49a), and (49b), it follows that 𝜁−𝜛=𝛿+𝐶₂₁𝐶₁₁⁻¹Δ≤𝛿𝑊−+𝑝₁ℎ𝐶₂₁𝐶⁻¹₁₁𝑊+𝑊+1+𝑝₁𝐶₂₁𝐶⁻¹₁₁⋅Δ≤𝑝₃Δ+𝛿+ℎ𝑌−𝑍+ℎ𝑌−𝑍,ℎ∈0,ℎ₁, (50) where 𝑝₃=max{1,(1+𝑝₂)‖𝐶₂₁𝐶⁻¹₁₁‖}, which completes the proof of Theorem 9.

Consider the compact form of (11a)–(11c):𝑍^(𝑛)=ℎ𝐶₁₁𝐹𝑡_𝑛,𝑍^(𝑛),𝑍^(𝑛),𝑍^(𝑛)+𝐶₁₂𝜛^(𝑛),𝜛^(𝑛+1)=ℎ𝐶₂₁𝐹𝑡_𝑛,𝑍^(𝑛),𝑍^(𝑛),𝑍^(𝑛)+𝐶₂₂𝜛^(𝑛), (51a)(51b)where 𝑍^(𝑛)=𝑍₁^(𝑛),…,𝑍_𝑠^(𝑛),𝑍₁^(𝑛−1),…,𝑍_𝑠^(𝑛−1)∈𝐶^2𝑠𝑁,𝑍^(𝑛)=𝑍₁^(𝑛),…,𝑍_𝑠^(𝑛),𝑍₁^(𝑛−1),…,𝑍_𝑠^(𝑛−1)∈𝐶^2𝑠𝑁,𝑍^(𝑛)=𝑍₁^(𝑛),…,𝑍_𝑠^(𝑛),𝑍₁^(𝑛−1),…,𝑍_𝑠^(𝑛−1)∈𝐶^2𝑠𝑁,𝜛^(𝑛+1)=𝑧^(𝑛+1),𝑧^(𝑛)∈𝐶^2𝑁. (52)

Theorem 10. Suppose the method (7a)–(7c) is algebraically stable for the matrices 𝐺 and 𝐷, then for (36a), (36b), (51a), and (51b), one has the following: 𝜁^(𝑛+1)−𝜛^(𝑛+1)2_𝐺≤𝜁^(𝑛)−𝜛^(𝑛)2_𝐺+𝑝₄ℎ𝑌^(𝑛)−𝑍^(𝑛)2+𝑌^{(𝑛−𝑚)}−𝑍^{(𝑛−𝑚)2}+𝑌^(𝑛)−𝑍^(𝑛)2, (53) where 𝑝₄=‖𝐷‖⋅max{(2𝛼+𝛽+𝜎),𝛽,𝜎}, ‖⋅‖ is a norm on 𝐶^2𝑠𝑁 defined by the following: ‖𝑈‖_𝐺=⟨𝑈,𝐺𝑈⟩^1/2=^2𝑠_𝑖,𝑗=1𝑔_𝑖𝑗𝑢_𝑖,𝑢_𝑗^1/2,𝑈=𝑢₁,𝑢₂,…,𝑢_2𝑠∈𝐶^2𝑠𝑁,𝑢_𝑖∈𝐶^𝑁. (54)

Proof. Define 𝑢^(𝑛)=𝜁^(𝑛)−𝜛^(𝑛), we get from (45a)–(45d) that 𝐾^(𝑛)=ℎ𝐹𝑡_𝑛,𝑌^(𝑛),𝑌^(𝑛),𝑌^(𝑛)−𝐹𝑡_𝑛,𝑍^(𝑛),𝑍^(𝑛),𝑍^(𝑛),𝑊^(𝑛)=𝐶₁₁𝐾^(𝑛)+𝐶₁₂𝑢^(𝑛),𝑢^(𝑛+1)=𝐶₂₁𝐾^(𝑛)+𝐶₂₂𝑢^(𝑛). (55) With algebraic stability, the matrix 𝑀=𝐺−𝐶^Τ₂₂𝐺𝐶₂₂𝐶^Τ₁₂𝐷−𝐶^Τ₂₂𝐺𝐶₂₁𝐷𝐶₁₂−𝐶^Τ₂₁𝐺𝐶₂₂𝐷𝐶₁₁+𝐶^Τ₁₁𝐷−𝐶^Τ₂₁𝐺𝐶₂₁ (56) is nonnegative definite. As in [25], we have 𝑢^(𝑛+1),𝐺𝑢^(𝑛+1)−𝑢^(𝑛),𝐺𝑢^(𝑛)−2Re𝑊^(𝑛),𝐷𝐾^(𝑛)=−𝑢^(𝑛),𝐾^(𝑛),𝑀𝑢^(𝑛),𝐾^(𝑛)≤0. (57) Using (2)–(4), we further obtain the following: 𝑢^(𝑛+1),𝐺𝑢^(𝑛+1)≤𝑢^(𝑛),𝐺𝑢^(𝑛)+2Re𝑊^(𝑛),𝐷𝐾^(𝑛)≤𝑢^(𝑛),𝐺𝑢^(𝑛)+2Re𝑊^(𝑛),𝐷𝐾₁^(𝑛)+2Re𝑊^(𝑛),𝐷𝐾₂^(𝑛)+𝐷𝐾₃^(𝑛)≤𝑢^(𝑛),𝐺𝑢^(𝑛)+2ℎ𝛼𝑊^(𝑛),𝐷𝑊^(𝑛)+2ℎ𝛽𝐷^1/2𝑊^(𝑛)⋅𝐷^1/2𝑊^(𝑛)+2ℎ𝜎𝐷^1/2𝑊^(𝑛)⋅𝐷^1/2𝑊^(𝑛)≤𝑢^(𝑛),𝐺𝑢^(𝑛)+2ℎ𝛼𝑊^(𝑛),𝐷𝑊^(𝑛)+ℎ𝛽‖𝐷‖⋅𝑊^(𝑛)2+‖𝐷‖⋅𝑊^(𝑛)2+ℎ𝜎‖𝐷‖⋅𝑊^(𝑛)2+‖𝐷‖⋅𝑊^(𝑛)2≤𝑢^(𝑛),𝐺𝑢^(𝑛)+2ℎ𝛼+ℎ𝛽+ℎ𝜎‖𝐷‖⋅𝑊^(𝑛)2+‖𝐷‖ℎ𝛽𝑊^{(𝑛−𝑚)2}+ℎ𝜎𝑊^(𝑛)2, (58) which gives (53). The proof is completed.

In the following, we assume that the method (7a)–(7c) has generalized stage order 𝑝; that is, there exists a function 𝐻^ℎ(𝑡) such that (40) holds. For any >0, we define 𝑌^(𝑛) and 𝑦^(𝑛+1) by the following:𝑌^(𝑛)=ℎ𝐶₁₁𝐹𝑡_𝑛,𝑌^(𝑛),𝑌𝑡_𝑛−𝜏,𝑌^(𝑛)+𝐶₁₂𝐻^ℎ𝑡_𝑛𝑦_(𝑛+1),=ℎ𝐶₂₁𝐹𝑡_𝑛,𝑌^(𝑛),𝑌𝑡_𝑛−𝜏,𝑌^(𝑛)+𝐶₂₂𝐻^ℎ𝑡_𝑛, (59a)(59b)where 𝑌^(𝑛)=^{𝑡𝑛+𝑐1ℎ}_{𝑡𝑛−𝜏+𝑐1ℎ}𝑔𝑡,𝜉,𝐻^ℎ₁(𝜉)𝑑𝜉,…,_{𝑡𝑛+𝑐𝑠ℎ𝑡𝑛−𝜏+𝑐𝑠ℎ}𝑔𝑡,𝜉,𝐻^ℎ₁(𝜉)𝑑𝜉,…,_{𝑡𝑛+𝑐1ℎ𝑡𝑛−𝜏+𝑐1ℎ}𝑔𝑡,𝜉,𝐻^ℎ₂(𝜉)𝑑𝜉,…,_{𝑡𝑛+𝑐𝑠ℎ𝑡𝑛−𝜏+𝑐𝑠ℎ}𝑔𝑡,𝜉,𝐻^ℎ₂(𝜉)=𝑑𝜉_{𝑡𝑛+𝑐1ℎ𝑡𝑛−𝜏+𝑐1ℎ}𝑔𝑡,𝜉,𝐻^ℎ₁(𝜉)𝑑𝜉,…,_{𝑡𝑛+𝑐𝑠ℎ𝑡𝑛−𝜏+𝑐𝑠ℎ}𝑔𝑡,𝜉,𝐻^ℎ₁(𝜉)𝑑𝜉,…,_{𝑡𝑛−ℎ+𝑐1ℎ𝑡𝑛−𝜏−ℎ+𝑐1ℎ}𝑔𝑡,𝜉,𝐻^ℎ₁(𝜉)𝑑𝜉,…,_{𝑡𝑛−ℎ+𝑐𝑠ℎ𝑡𝑛−ℎ−𝜏+𝑐𝑠ℎ}𝑔𝑡,𝜉,𝐻^ℎ₁(𝜉)𝑑𝜉. (60)

Theorem 11. Suppose the method (7a)–(7c) is diagonally stable and its generalized stage order is 𝑝, then there exist constants 𝑝₅ and ℎ₂ such that ^𝑛_𝑘=1𝑌^(𝑘)−𝑌(𝑡_𝑘−𝜏)²=^𝑛_𝑘=1𝑌^{(𝑘−𝑚)}−𝑌(𝑡_𝑘−𝜏)²≤𝑝₅^𝑛_𝑘=1𝜁^(𝑘)−𝐻^ℎ(𝑡_𝑘)²+𝐶₁ℎ^(𝑝+1)+𝐶₂ℎ^2(𝜈+1)+𝐶₃⁰_{𝑘=−𝑚+1}𝑌^(𝑘)−𝑌(𝑡_𝑘)²,ℎ∈0,ℎ₂, (61) where 𝑝₅,𝐶₁,𝐶₂, and 𝐶₃ depend on 𝑝₂,𝑠,𝑑₁,ℎ₁, and ‖𝐶₁₂‖.

Proof. A combination of (36a) and (41a) leads to the following: 𝑌^(𝑘)−𝑌𝑡_𝑘=ℎ𝐶₁₁𝐹𝑡_𝑘,𝑌^(𝑘),𝑌^(𝑘),𝑌^(𝑘)−𝐹𝑡_𝑘,𝑌𝑡_𝑘,𝑌𝑡_𝑘,𝑌𝑡_𝑘+𝐶₁₂𝜁^(𝑘)−𝐻^ℎ𝑡_𝑘−Δ^ℎ(𝑡). (62) It follows from Theorem 9 that 𝑌^(𝑘)−𝑌𝑡_𝑘≤𝑝₂ℎ𝑌^{(𝑘−𝑚)}−𝑌𝑡_𝑘−𝜏+ℎ𝑌^(𝑘)−𝑌𝑡_𝑘+𝐶₁₂𝜁^(𝑘)−𝐻^ℎ𝑡_𝑘−Δ^ℎ(𝑡). (63) A combination of (4) and (9) gives the following: 𝑌^(𝑘)−𝑌𝑡_𝑘≤ℎ^𝑚_𝑞=0𝑝_𝑞𝑔𝑡_𝑘+𝑐₁ℎ,𝑡_𝑘−𝑞+𝑐₁ℎ,𝑌₁^{(𝑘−𝑞)}ℎ,…,_𝑚𝑞=0𝑝_𝑞𝑔𝑡_𝑘+𝑐_𝑠ℎ,𝑡_𝑘−𝑞+𝑐_𝑠ℎ,𝑌_𝑠^{(𝑘−𝑞)},ℎ^𝑚_𝑞=0𝑝_𝑞𝑔𝑡_𝑘−1+𝑐₁ℎ,𝑡_{𝑘−𝑞−1}+𝑐₁ℎ,𝑌₁^{(𝑘−𝑞−1)}ℎ,…,_𝑚𝑞=0𝑝_𝑞𝑔𝑡_𝑘−1+𝑐_𝑠ℎ,𝑡_{𝑘−𝑞−1}+𝑐_𝑠ℎ,𝑌_𝑠^{(𝑘−𝑞−1)}−^{𝑡𝑛+𝑐1ℎ}_{𝑡𝑛−𝜏+𝑐1ℎ}𝑔𝑡,𝜉,𝑦(𝜉)𝑑𝜉,…,_{𝑡𝑛+𝑐𝑠ℎ𝑡𝑛−𝜏+𝑐𝑠ℎ}𝑔𝑡,𝜉,𝑦(𝜉)𝑑𝜉,…,_{𝑡𝑛+𝑐1ℎ𝑡𝑛−𝜏+𝑐1ℎ}𝑔𝑡,𝜉,𝑦(𝜉−ℎ)𝑑𝜉,…,_{𝑡𝑛+𝑐𝑠ℎ𝑡𝑛−𝜏+𝑐𝑠ℎ}𝑔()𝑑𝑡,𝜉,𝑦𝜉−ℎ𝑑𝜉≤2𝑠₁ℎ^𝜈, (64) where 𝑑₁ and 𝜈 depend only on the compound quadrature (CQ) formula (9), which on substitution into (63) gives the following: ^𝑛_𝑘=1𝑌^(𝑘)−𝑌(𝑡_𝑘)²≤3𝑝²₁ℎ^2𝑛_𝑘=1𝑌^{(𝑘−𝑚)}−𝑌(𝑡_𝑘−𝜏)²+𝑝²₂𝐶₁₂^2𝑛_𝑘=1𝜁^(𝑘)−𝐻^ℎ(𝑡_𝑘)²+𝑝²₁⋅Δ^ℎ𝑡_𝑘²+𝑝²₂4𝑠²𝑑²₁ℎ^2𝜈+2,_𝑛𝑘=1𝑌^{(𝑘−𝑚)}−𝑌(𝑡_𝑘−𝜏)²=^𝑛_𝑘=1𝑌^{(𝑘−𝑚)}−𝑌(𝑡_𝑘−𝑚ℎ)²=^𝑛−𝑚_{𝑖=1−𝑚}𝑌^(𝑖)−𝑌(𝑡_𝑖)²≤^𝑛_𝑘=1𝑌^(𝑘)−𝑌(𝑡_𝑘)²+⁰_{𝑘=1−𝑚}𝑌^(𝑘)−𝑌(𝑡_𝑘)². (65) Therefore, there exist 𝑝₅ and ℎ₂ such that (61) holds. The proof of Theorem 11 is completed.

Theorem 12. Suppose method (7a)–(7c) is algebraically stable and diagonally stable and its generalized stage order is 𝑝. Then the method with quadrature formula (9) is 𝐷-convergent of order at least min{𝑝,𝜈}.

Proof. A combination of (36b), (59b), and (53) leads to the following: 𝜁^(𝑛+1)−𝑦^(𝑛+1)2_𝐺≤𝜁^(𝑛)−𝑦^(𝑛)2_𝐺+𝑝₄ℎ×𝑌^(𝑛)−𝑌^(𝑛)2+𝑌^{(𝑛−𝑚)}−𝑌(𝑡_𝑛−𝑚)²+𝑌^(𝑛)−𝑌^(𝑛)2. (66) Using Theorem 9, we have the following: 𝑌^(𝑛)−𝑌^(𝑛)≤𝑝₂𝐶₁₂⋅𝜁^(𝑛)−𝐻^ℎ𝑡_𝑛+ℎ𝑌^{(𝑛−𝑚)}−𝑌𝑡_𝑛−𝑚+ℎ𝑌^(𝑛)−𝑌^(𝑛), (67) which on substitution into (66) gives the following: 𝜁^(𝑛+1)−𝑦^(𝑛+1)_𝐺⁽²⁾≤1+3ℎ𝑝₄𝑝²₂𝐶₁₂²𝜆₂⋅𝜁^(𝑛)−𝐻^ℎ(𝑡_𝑛)²_𝐺+𝑝₄1+3ℎ²𝑝²₂ℎ𝑌^{(𝑛−𝑚)}−𝑌(𝑡_𝑛−𝑚)²+𝑝₄1+3ℎ²𝑝²₂ℎ2𝑠𝑑²ℎ^2𝜈, (68) where 𝜆₂ is the minimum characteristic value of 𝐺. On the other hand, 𝜁^(𝑛)−𝐻^ℎ(𝑡_𝑛)²_𝐺≤𝜁^(𝑛)−𝑦^(𝑛)2_𝐺+𝑦^(𝑛)−𝐻^ℎ(𝑡_𝑛)²_𝐺+2𝜁^(𝑛)−𝑦^(𝑛)_𝐺⋅𝑦^(𝑛)−𝐻^ℎ(𝑡_𝑛)_𝐺≤(1+ℎ)𝜁^(𝑛)−𝑦^(𝑛)2_𝐺+11+ℎ𝑦^(𝑛)−𝐻^ℎ𝑡_𝑛²_𝐺. (69) In view of (41a)–(41c) and (59a) and (59b), the application of Theorem 9 leads to the following: 𝑦^(𝑛)−𝐻^ℎ𝑡_𝑛≤𝑝₃Δ^ℎ𝑡_𝑛+𝛿^ℎ𝑡_𝑛,ℎ∈0,ℎ₀ (70) which gives 𝑦^(𝑛)−𝐻^ℎ(𝑡_𝑛)²_𝐺≤𝜆₁𝑝²₃Δ^ℎ𝑡_𝑛+𝛿^ℎ𝑡_𝑛²,ℎ∈0,ℎ₀, (71) where 𝜆₁ denotes the maximum eigenvalue of the matrix 𝐺. A combination of (40) and (66)–(71) leads to the following: 𝜁^(𝑛)−𝐻^ℎ(𝑡_𝑛)²_𝐺≤(1+ℎ)𝛾₁𝜁^(𝑛−1)−𝐻^ℎ𝑡_𝑛−1²_𝐺+𝛾₂𝑌^{(𝑛−𝑚−1)}−𝑌𝑡_{𝑛=𝑚−1}²+𝛾₃ℎ^2𝜈+11+41+ℎ𝜆₁𝑝²₃𝑝²₁ℎ^2𝑝+2,ℎ∈0,ℎ₃, (72) where ℎ₃=minℎ₀,ℎ₁,ℎ₂𝛾≤1,₁=1+3ℎ𝑝₄𝑝²₂𝐶₁₂²𝜆₂,𝛾₂=𝑝₄ℎ1+3𝑝²₂ℎ²,𝛾₃=2𝑠𝑑²𝑝₄1+3ℎ²𝑝²₂. (73) Considering Theorem 11, we further obtain the following: 𝜁^(𝑛)−𝐻^ℎ(𝑡_𝑛)²_𝐺≤(1+ℎ)^𝑛𝛾^𝑛₁𝜁⁽⁰⁾−𝐻^ℎ(𝑡₀)²_𝐺+(1+ℎ)𝛾₂×^𝑛_𝑖=1(1+ℎ)^𝑖−1𝛾₁^𝑖−1𝑌^{(𝑛−𝑚−𝑖)}−𝑌(𝑡_{𝑛=𝑚−𝑖})²_𝐺+(1+ℎ)𝛾₃ℎ^2𝜈+11+41+ℎ𝜆₁𝑝²₃𝑝²₁ℎ^2𝑝+2×^𝑛_𝑖=1(1+ℎ)^𝑖−1𝛾₁^𝑖−1. (74) Using discrete Bellman inequality, we have the following: 𝜁^(𝑛)−𝐻^ℎ(𝑡_𝑛)²_𝐺≤(1+ℎ)^𝑛𝛾^𝑛₁𝜁⁽⁰⁾−𝐻^ℎ(𝑡₀)²_𝐺+(1+ℎ)𝛾₂×^𝑛_𝑖=1(1+ℎ)^𝑖−1𝛾₁^0𝑖−1_{𝑗=1−𝑚}𝑌^(𝑗)−𝑌(𝑡_𝑗)²_𝐺+(1+ℎ)𝛾₃ℎ^2𝜈+11+41+ℎ𝜆₁𝑝²₃𝑝²₁ℎ^2𝑝+2×^𝑛_𝑖=1(1+ℎ)^𝑖−1𝛾₁^𝑖−1×exp(1+ℎ)𝛾₂^𝑛_𝑖=1(1+ℎ)^𝑖−1𝛾₁^𝑖−1, (75) where ℎ∈(0,ℎ₀], and ℎ₀=min{ℎ₀,ℎ₁,ℎ₂,ℎ₃}≤1.
Considering ‖𝐻(𝑡)−𝐻^ℎ(𝑡_𝑛)‖≤𝑝₀ℎ^𝑝, we obtain the following: 𝜁^(𝑛)−𝐻𝑡_𝑛≤𝜁^(𝑛)−𝐻^ℎ𝑡_𝑛+𝐻^ℎ𝑡_𝑛−𝐻𝑡_𝑛≤𝜁^(𝑛)−𝐻^ℎ𝑡_𝑛+𝑝₀ℎ^𝑝. (76) Considering (75) and (76), we can easily conclude that method (7a)–(7c) with quadrature formula (9) is 𝐷-convergent of order at least >min{𝑝,𝜈}. The proof is completed.

5. Numerical Experiments

Consider the following nonlinear Volterra delay integro-differential equations: 𝑦^′(𝑡)()𝑦(𝑡)𝜋=−6+sin𝑡+𝑦𝑡−4−^𝑡_{𝑡−𝜋/4}(𝑣)𝑦(𝑣)()𝑦(𝑡)()𝜋sin𝑑𝑣+5expcos𝑡,𝑡≥0,=expcos𝑡,−4≤𝑡≤0 (77) and its perturbed problem 𝑧^′(𝑡)()𝑧(𝑡)𝜋=−6+sin𝑡+𝑧𝑡−4−^𝑡_{𝑡−𝜋/4}(𝑣)𝑧(𝑣)()𝑧(𝑡)𝜋sin𝑑𝑣+5expcos𝑡,𝑡≥0,=2,−4≤𝑡≤0. (78) One may check that this system has the solution 𝑦(𝑡)=exp(cos𝑡) and 𝑎=−4, 𝛽=1/2, 𝜎=1/2, 𝛾=1. With Theorem 9, we conclude that system (77) satisfies stability and convergence properties.

Apply the two-step Runge-Kutta method induced by the GL method in [12] (79a)to (77) and (78), 𝑌_𝑗^(𝑛) is computed by the compound Simpson rule with an even integer 𝑚≥4 formula: 𝑌_𝑗^(𝑛)=ℎ3sin𝑡_𝑛+𝑐_𝑗ℎ𝑌_𝑗^(𝑛)+4^𝑚/2_𝑞=1sin𝑡_{𝑛−2𝑞+1}+𝑐_𝑗ℎ𝑌_𝑗^{(𝑛−2𝑞+1)}+2^𝑚−2/2_𝑞=1sin𝑡_𝑛−2𝑞+𝑐_𝑗ℎ𝑌_𝑗^{(𝑛−2𝑞)}+sin𝑡_𝑛−𝑚+𝑐_𝑗ℎ𝑌_𝑗^{(𝑛−𝑚)} (79b)Using the condition (28), we know this method is stable; the convergence order is shown in Table 1.

For Figures 1 and 2, it is obvious that the corresponding method for VDIDEs is stable and convergent, and the convergent order is min{4,4}.

Figure 1 

The numerical solution of (79a) and (79b) for (77) with ℎ=0.1.

Figure 2 

The numerical solution of (79a) and (79b) for (78) with ℎ=0.1.

6. Conclusions

This paper is devoted to the stability and convergence analysis of the two-step Runge-Kutta (TSRK) methods with compound quadrature formula for the numerical solution for a nonlinear Volterra delay integro-differential equations. Nonlinear stability and 𝐷-convergence are introduced and proved.

We believe that the results presented in this paper can be extended to other general DIDEs and NDIDEs. However, it is difficult to extend the results presented in this paper to more general delay differential equations; the discussion will be discussed later.

Acknowledgments

This work was supported by the Education Foundation of Heilongjiang Province of China (12523039) and the Doctor Foundation of Heilongjiang Institute of Technology (2012BJ27).