The Stability of Two-Step Runge-Kutta Methods for Neutral Delay Integro Differential-Algebraic Equations with Many Delays

Abstract

This paper studies the asymptotic stability of the two-step Runge-Kutta methods for neutral delay integro differential-algebraic equations with many delays. It proves that A-stable two-step Runge-Kutta methods are asymptotically stable for neutral delay integro differential-algebraic equations with many delays.

1. Introduction

The stability of numerical methods for delay differential equations has been intensively studied in [1–3] for many years. These equations appeared in a wide variety of scientific and engineering fields, such as circuit analysis, computer-aided design power systems, and optimal control. The structure for these, the order of convergence, and the asymptotic stability of numerical methods have been studied in [4–6]. Zhu and Petzold investigated the asymptotic stability of neutral delay differential equations with 𝜃-methods, Runge-Kutta methods, BDF methods, and linear multistep methods [7]. Zhao et al. studied the stability of neutral delay differential equations with Rosenbrock methods [8]. Yu et al. studied the general neutral delay differential equations with multistep methods [9]. More recently, there is a growing interest in the analysis of delay integro differential equations. Baker and Ford [10] studied the asymptotic stability of a class of linear multistep (LM) methods for scalar linear delay integro differential equations; Koto [11] dealt with the linear stability of Runge-Kutta (R-K) methods for systems of delay integro differential equations; Huang and Vandewalle [12] gave sufficient and necessary stability conditions for exact and discrete solutions of linear scalar delay integro differential equations, and Luzyanina et al. [13] developed computational procedures for determining the stability of delay integro differential equations. Zhang and Vandewalle [14] gave the stability criteria for exact and discrete solution of neutral multidelay integro differential equations. Although the stability of numerical methods for delay integro differential equations has been very intensively studied, the stability of delay integro differential equations with many delays has not been studied so far.

In this paper, we focus on the asymptotic stability of numerical methods for neutral delay integro differential-algebraic equations with many delays. This paper is structured as follows. In Section 2 we give asymptotic stability of the analytical solution and introduce two-step Runge-Kutta methods and the stability region. In Section 3, we deal with the asymptotic stability of two-step Runge-Kutta method for neutral delay integro differential-algebraic equations with many delays; the theoretical results are proved. In Section 4, an example is given to illustrate the theoretical results.

2. Asymptotic Stability of the Analytical Solution

2.1. Asymptotic Stability of the Analytical Solution of Neutral Delay Integro Differential-Algebraic Equation with Many Delays

In this section, we consider the following linear system: 𝐴𝑢^′(𝑡)+𝐵𝑢(𝑡)+^𝑚_𝑞=1𝐶_𝑞𝑢^′𝑡−𝜏_𝑞+^𝑚_𝑞=1𝐷_𝑞𝑢𝑡−𝜏_𝑞+^𝑚_𝑞=1𝐺_𝑞^𝑡_{𝑡−𝜏𝑞}𝑢(𝛿)𝑢(𝑡)(𝑡)[),𝑑𝛿=0,𝑡≥0,=𝜑,𝑡∈−𝜏,0 (1) where 𝐴, 𝐵, 𝐶_𝑞, 𝐷_𝑞, 𝐺_𝑞∈𝑅^𝑑×𝑑, 𝐴 is a singular matrix, 𝜏_𝑞 is a given positive delay constant (𝑞=1,2,…,𝑚), and 0<𝜏₁≤𝜏₂≤⋯≤𝜏_𝑚=𝜏. 𝜑(𝑡) denotes a given vector-valued function and 𝑢(𝑡) is a vector-valued unknown function to be solved for 𝑡≥0.

In order to obtain the characteristic equation of system (1), we focus on the exponential solutions 𝑢(𝑡)=𝑒^𝑠𝑡𝑥 of (1); here 𝑥=(𝑥₁,𝑥₂,…,𝑥_𝑑)^Τ∈𝐶^𝑑 denotes the unknown vector. Then we have 𝑢(𝑡)=𝑒^𝑠𝑡𝑥₁,𝑒^𝑠𝑡𝑥₂,…,𝑒^𝑠𝑡𝑥_𝑑^Τ,𝑢^′(𝑡)=𝑠𝑒^𝑠𝑡𝑥,𝑢𝑡−𝜏_𝑞=𝑒^{𝑠(𝑡−𝜏𝑞)}𝑥₁,𝑒^{𝑠(𝑡−𝜏𝑞)}𝑥₂,…,𝑒^{𝑠(𝑡−𝜏𝑞)}𝑥_𝑑^Τ,𝑞=1,2,…,𝑚,𝑢^′𝑡−𝜏_𝑞=𝑠𝑒^{𝑠(𝑡−𝜏𝑞)}𝑥,𝑞=1,2,…,𝑚. (2a)(2b)(2c)(2d)

Substituting the above results into (1), we have the following equation: 𝑠𝐴+𝐵+𝑠^𝑚_𝑞=1𝐶_𝑞𝑒^{−𝑠𝜏𝑞}+^𝑚_𝑞=1𝐷_𝑞𝑒^{−𝑠𝜏𝑞}+𝑠^𝑚−1_𝑞=1𝐺_𝑞1−𝑒^{−𝑠𝜏𝑞}𝑥=0. (2e)

The existence of a nonzero 𝑥 in (2e) implies the characteristic equation of system (1) holds; that is, the following equation holds: det𝑠𝐴+𝐵+𝑠^𝑚_𝑞=1𝐶_𝑞𝑒^{−𝑠𝜏𝑞}+^𝑚_𝑞=1𝐷_𝑞𝑒^{−𝑠𝜏𝑞}+𝑠^𝑚−1_𝑞=1𝐺_𝑞1−𝑒^{−𝑠𝜏𝑞}=0. (3)

Definition 1 (see [13]). Equation (1) is said to be asymptotically stable, if for any continuous differential initial function and for any delay 𝜏_𝑞>0, 𝑞=1,2,…,𝑚 the analytical solution to (1) satisfies lim_𝑡→∞𝑢(𝑡)=0.

We know that the stability of analytical solution can be studied via the characteristic equation, so we give a criterion for the asymptotic stability of (1), which is based on the following lemmas.

Lemma 2 (see [14]). Assume SupRe(𝜆)∶𝑝(𝜆)=0<0, (4) where 𝑝(𝜆)=det{𝜆𝐴+𝐵+𝜆∑^𝑚_𝑞=1𝐶_𝑞𝑒^{−𝜆𝜏𝑞}+∑^𝑚_𝑞=1𝐷_𝑞𝑒^{−𝜆𝜏𝑞}+∑^𝑚_𝑞=1𝐺_𝑞𝜏_𝑞𝜂(𝑒^{−𝜆𝜏𝑞})} is the characteristic polynomial of (1). Then, system (1) is asymptotically stable.

Where 𝜂(𝑧) is a complex function defined by 𝜂(𝑧)=1−𝑧{}ln𝑧,𝑧∈𝐶⧵0,10,𝑧=0−1,𝑧=1. (5) And ln𝑧=ln|𝑧|+𝑖arg𝑧 (𝑧=0,1;−𝜋<arg𝑧≤𝜋) is the principal branch of the multivalued complex natural logarithm.

Lemma 3 (see [14]). Function 𝜂(𝑧) is analytic in 𝐶⧵𝑅⁻₀ and satisfies |𝜂(𝑧)|≤1 for |𝑧|≤1, where 𝑅⁻₀={𝑥∈𝑅∶𝑥≤0}.

Lemma 4. If the matrix (𝐴+∑^𝑚_𝑞=1𝐶_𝑞𝑒^{−𝜆𝜏𝑞}) is invertible for Re(𝜆)≥𝑟, where 𝑟∈𝑅, then the function 𝑝(𝜆)=det𝜆²𝐼_𝑑+𝐴+^𝑚_𝑞=1𝐶_𝑞𝑒^{−𝜆𝜏𝑞−1}×𝜆𝐵+𝜆^𝑚_𝑞=1𝐷_𝑞𝑒^{−𝜆𝜏𝑞}+^𝑚_𝑞=1𝐺_𝑞1−𝑒^{−𝜆𝜏𝑞} (6) has at most a finite number of zeros for Re(𝜆)≥𝑟.

Proof. When Re(𝜆)≥𝑟, the function 𝑝(𝜆) can be expanded into the following form: 𝑝(𝜆)=𝜆^2𝑑+𝜓_2𝑑−1𝑒^−𝜆𝜏1,𝑒^−𝜆𝜏2,…,𝑒^{−𝜆𝜏𝑚1}𝜆^2𝑑−1+⋯+𝜓₀𝑒^−𝜆𝜏1,𝑒^−𝜆𝜏2,…,𝑒^{−𝜆𝜏𝑚1}, (7) where 𝜓_𝑖(𝑒^−𝜆𝜏1,𝑒^−𝜆𝜏2,…,𝑒^{−𝜆𝜏𝑚1}), 𝑖=0,1,…,2𝑑−1, are rational functions for the expressions 𝑒^−𝜆𝜏1,𝑒^−𝜆𝜏2,…,𝑒^{−𝜆𝜏𝑚1}, and they have no poles for Re(𝜆)≥𝑟.
Since 𝜏_𝑖>0, we have that 𝑒^{−𝜆𝜏𝑖}=𝑒^{−𝜏𝑖Re(𝜆)}≤𝑒^{−𝜏𝑖𝑟},forRe(𝜆)≥𝑟. (8)
Hence, there exist constants 𝐾_𝑖>0 such that 𝜓_𝑖𝑒^−𝜆𝜏1,𝑒^−𝜆𝜏2,…,𝑒^{−𝜆𝜏𝑚1}≤𝐾_𝑖,𝑖=0,1,…,2𝑑−1. (9) Let 𝑀 be a positive number large enough such that 𝐾_2𝑑−1𝑀+𝐾_2𝑑−2𝑀²𝐾+⋯+⁰𝑀^2𝑑<1, (10)
which implies that, for Re(𝜆)≥𝑟 and |𝜆|≥𝑀, 𝑝(𝜆)≥|𝜆|^2𝑑𝐾1−^2𝑑−1𝑀−𝐾_2𝑑−2𝑀²𝐾−⋯−⁰𝑀^2𝑑>0. (11) That is, 𝑝(𝜆)≠0 in the set {𝜆∶Re𝜆≥𝑟,|𝜆|≥𝑀}.
By the isolation property of the zeros for analytic functions, 𝑝(𝜆) has at most a finite number of zeros in the set {𝜆∶Re𝜆≥𝑟,|𝜆|<𝑀}; this proves the lemma.

In the following, we denote the spectrum of a square matrix 𝐴 by 𝜎(𝐴) and introduce the set 𝐶⁻={𝑧∈𝐶∶Re(𝑧)<0}. (12)

Theorem 5. System (1) is asymptotically stable if the following conditions are satisfied:(a)det(𝐴+∑^𝑚_𝑞=1𝜉_𝑞𝐶_𝑞)≠0 for |𝜉_𝑞|≤1,(b)𝜎(𝐺(𝜉))⊆𝐶⁻ for 𝜉=(𝜉₁,𝜉₂,…,𝜉_𝑚)^Τwith |𝜉_𝑞|≤1, where 𝐺(𝜉)=𝐴+^𝑚_𝑞=1𝜉_𝑞𝐶_𝑞⁻¹×−𝐵−^𝑚_𝑞=1𝜉_𝑞𝐷_𝑞−^𝑚_𝑞=1𝜂𝜉_𝑞𝐺_𝑞𝜏_𝑞. (13)

Proof. When |𝜉_𝑞|≤1, 𝑞=1,2,…,𝑚, condition (a) leads to 𝑃𝜆,𝜉₁,𝜉₂,…,𝜉_𝑚=det𝜆𝐴+𝜆^𝑚_𝑞=1𝜉_𝑞𝐶_𝑞++𝐵_𝑚𝑞=1𝜉_𝑞𝐷_𝑞+^𝑚_𝑞=1𝜂𝜉_𝑞𝐺_𝑞𝜏_𝑞=det𝐴+^𝑚_𝑞=1𝜉_𝑞𝐶_𝑞det𝜆𝐼_𝑑−𝐺(𝜉). (14)
Condition (b) leads to 𝑃(𝜆)=𝑃𝜆,𝑒^−𝜆𝜏1,𝑒^−𝜆𝜏2,…,𝑒^{−𝜆𝜏𝑚}(𝜆)≠0forRe≥0. (15)
Hence Sup{Re(𝜆)∶𝑃(𝜆)=0}≤0. (16)
Now we will show that the strict inequality in (16) holds.
Define 𝐹𝜉₁,𝜉₂,…,𝜉_𝑚=det𝐴+^𝑚_𝑞=1𝜉_𝑞𝐶_𝑞, (17) and then 𝐹(𝜉₁,𝜉₂,…,𝜉_𝑚) is a multivariate polynomial and is nonzero on the compact domain defined by |𝜉_𝑞|≤1, 𝑞=1,2,…,𝑚, and equal to 1 at the origin. Hence, its modulus is bounded; that is, 𝐹𝜉₁,𝜉₂,…,𝜉_𝑚𝜉≥𝜀>0,when_𝑞≤1,for𝑞=1,2,…,𝑚. (18)
By the continuity of 𝐹, there exists a 𝛿>0 such that 𝐹𝜉₁,𝜉₂,…,𝜉_𝑚𝜉>0,when_𝑞≤1+𝛿,for𝑞=1,2,…,𝑚. (19)
It follows from this that det𝐴+^𝑚_𝑞=1𝜉_𝑞𝐶_𝑞𝑒≠0,when^{−𝜆𝜏𝑞}≤1+𝛿,for𝑞=1,2,…,𝑚. (20)
Let 𝑟 be the strictly positive number 𝑟=ln(1+𝛿)/𝜏; then det𝐴+^𝑚_𝑞=1𝑒^{−𝜆𝜏𝑞}𝐶_𝑞≠0forRe(𝜆)≥−𝑟. (21)
Thus, the equation 𝑃(𝜆)=0 has only a finite number of roots when Re(𝜆)≥−𝑟, and it holds true for the equation 𝑝(𝜆)=0 by condition (a). Combined with (16) we get that the characteristic equation has at most a finite number of roots in the region {𝜆∶−𝑟≤Re(𝜆)<0}.
Let −𝛾=max_{(𝜆)−𝑟≤Re<0}{Re(𝜆)}; (22) then 𝛾>0.
When Re(𝜆)>−𝛾, the characteristic equation 𝑝(𝜆)=0 has no root. Hence, a strict inequality holds in (16). By Lemma 2, the proof is completed.

2.2. The Two-Step Runge-Kutta Methods and the Stability Region

Consider the two-step Runge-Kutta method: 𝑌^(𝑛)=ℎ𝐶₁₁𝐹𝑡_𝑛,𝑌^(𝑛)+𝐶₁₂𝑦^(𝑛−1)𝑦^(𝑛)=ℎ𝐶₂₁𝐹𝑡_𝑛,𝑌^(𝑛)+𝐶₂₂𝑦^(𝑛−1) (23) for solving the initial value problem (1).

In order to simplify the analysis, we consider two-step Runge-Kutta method (TSRK) of the form𝑢_𝑛+1=(1−𝜃)𝑢_𝑛+𝜃𝑢_𝑛−1+ℎ^𝑠_𝑗=1𝑏̂_𝑗𝑓𝑡_𝑗,𝑈^𝑗_𝑛+ℎ^𝑠_𝑗=1𝑏̃_𝑗𝑓𝑡_𝑗̃,𝑈^𝑗_𝑛−1,𝑈^𝑖_𝑛=𝑢_𝑛+ℎ^𝑠_𝑗=1𝑎_𝑖,𝑗𝑓𝑡_𝑗,𝑈^𝑗_𝑛𝑖=1,2,…,𝑠, (24a)(24b)where 𝑡_𝑗=𝑡_𝑛+𝑐_𝑗ℎ, 𝑡_𝑗̃=𝑡_𝑛−1+𝑐_𝑗ℎ, 𝑢_𝑖 is an approximation to 𝑢(𝑡_𝑖), ℎ is a fixed step-size, 𝜃, 𝑏̂_𝑗, 𝑏̃_𝑗, 𝑎_𝑖,𝑗, and 𝑐_𝑗 are coefficients of the method, 0≤𝜃≤1.

These methods are a subclass of general linear methods introduced by Butcher [15] and could be possibly also referred to as two-step hybrid methods. They generalize 𝑘-step collocation methods (with 𝑘=2) for ordinary differential equations (ODEs) studied by Lie and Nørsett [16] and Lie [17] and two-step Runge-Kutta methods for ODEs investigated by Byrne and Lambert [18]. The variable stepsize continuous two-step Runge-Kutta methods for ODEs were investigated by Jackiewicz and Tracogna [19]. Here we will represent (24a) and (24b) by the following table of the coefficients: 𝐶𝐴𝑏̂^T𝜃𝑏̃^T=𝑐₁𝑎₁₁𝑎₁₂⋯𝑎_1𝑠𝑐₂𝑎₂₁𝑎₂₂⋯𝑎_2𝑠⋮𝑐⋮⋮⋱⋮_𝑠𝑎_𝑠1𝑎_𝑠2⋯𝑎_𝑠𝑠𝑏̂₁𝑏̂₂⋯𝑏̂_𝑠𝜃𝑏̃₁𝑏̃₂⋯𝑏̃_𝑠, (25) where 𝑐_𝑖=∑^𝑠_𝑗=1𝑎_𝑖𝑗 and ∑^𝑠_𝑗=1(𝑏̂_𝑗+𝑏̃_𝑗)=1+𝜃.

Apply (24a) and (24b) to the basic test equation 𝑢^′(𝑡)=𝑎𝑢(𝑡)𝑡≥0,Re𝑎<0, (26) which gives the following equations: 𝑢_𝑛+1=(1−𝜃)𝑢_𝑛+𝜃𝑢_𝑛−1+𝛼^𝑠_𝑗=1𝑏̂_𝑗𝑈^𝑗_𝑛+𝑏̃_𝑗𝑈^𝑗_𝑛−1𝑈^𝑖_𝑛=𝑢_𝑛+𝛼^𝑠_𝑗=1𝑎_𝑖,𝑗𝑈^𝑗_𝑛. (27)

Rewriting (27) we obtain 𝑢_𝑖+1=𝑅(𝛼,𝜃)𝑢_𝑖+𝑆(𝛼,𝜃)𝑢_𝑖−1, (28) where 𝑅(𝛼,𝜃)=1−𝜃+𝛼𝑏̂^Τ𝐼−𝛼𝐴⁻¹𝑆()𝑏̃𝑒,𝛼,𝜃=𝜃+𝛼^Τ𝐼−𝛼𝐴⁻¹[]𝑒,𝛼=𝑎ℎ,𝑒=1,1,…,1_𝑠^Τ. (29) To investigate the stability properties of (24a) and (24b) with (26), we must investigate the asymptotic behaviors of the solution to (28). This is determined by the location of roots of the characteristic polynomial 𝜑(𝑧)=𝑧²−𝑅(𝑎,𝜃)𝑧−𝑆(𝑎,𝜃). (30)

The stability region of the two-step Runge-Kutta methods (24a) and (24b) is the set of all points 𝛼 for which the roots of 𝜑(𝑧) are inside or on the unit circle with those on the unit circle being simple. If 𝜑(𝑧) is a Schur polynomial for any 𝛼 with Re𝛼<0, the stability of the two-step Runge-Kutta method contains the negative half plane; the method is said to be A-stable for ODEs.

3. Asymptotic Stability of TSRK Methods for Neutral Delay Integro Differential-Algebraic Equation with Many Delays

In this section, we will confine our discussion to neutral delay integro differential-algebraic equation with commensurate delays, that is, systems of the form (1) with 𝜏_𝑞=𝑞𝜏, 𝜏=𝑀ℎ, 𝑀 is a positive integer, 𝑞=1,2,…,𝑚.

Definition 6 (see [20]). A numerical method for asymptotically stable system (1) is called asymptotically stable if the numerical solution satisfies lim_𝑛→∞𝑢_𝑛=0. (31)

Applying the two-step method (24a) and (24b) to (1), we have 𝑢_𝑛+1=(1−𝜃)𝑢_𝑛+𝜃𝑢_𝑛−1+^𝑠_𝑗=1𝑏̂_𝑗𝐾_𝑛,𝑗+^𝑠_𝑗=1𝑏̃_𝑗𝐾_{𝑛−1,𝑗},𝐴𝐾_𝑛,𝑖+ℎ𝐵𝑢_𝑛+^𝑠_𝑗=1𝑎_𝑖𝑗𝑘_𝑛,𝑗+^𝑚_𝑞=1𝐶_𝑞𝐾_{𝑛−𝑞𝑀,𝑖}+ℎ^𝑚_𝑞=1𝐷_𝑞𝑢_{𝑛−𝑞𝑀}+^𝑠_𝑗=1𝑎_𝑖𝑗𝐾_{𝑛−𝑞𝑀,𝑗}+ℎ^𝑚_𝑞=1^𝑞𝑀_𝛾=0𝛾_𝑟𝐺_𝑞⋅𝑢_𝑛−𝑟+^𝑠_𝑗=1𝑎_𝑖𝑗𝑘_{𝑛−𝑟,𝑗}=0for𝑖=1,2,…,𝑠, (32)(33) where 𝐾_𝑛,𝑖=[𝐾¹_𝑛,𝑖,𝐾²_𝑛,𝑖,…,𝐾^𝑑_𝑛,𝑖]^Τ, 𝑖=1,2,…,𝑠, are stage derivatives multiplied by ℎ.

Let 𝑏̂^T=𝑏̂₁,𝑏̂₂,…,𝑏̂_𝑠,𝑏̃^T=𝑏̃₁,𝑏̃₂,…,𝑏̃_𝑠,𝐴=𝑎_𝑖𝑗. (34)

We assume that all the eigenvalues of 𝐴 have positive real part. Rearrange the variables of the stage derivatives as 𝐾_𝑛=𝐾¹_𝑛,1,𝐾¹_𝑛,2,…,𝐾¹_𝑛,𝑠,𝐾²_𝑛,1,𝐾²_𝑛,2,…,𝐾²_𝑛,𝑠,…,𝐾^𝑑_𝑛,1,𝐾^𝑑_𝑛,2,…,𝐾^𝑑_𝑛,𝑠^Τ. (35) Define 𝑌_𝑛=𝐾^Τ_𝑛,𝑢^Τ_𝑛+1^Τ,𝐵=ℎ𝐵,𝐷_𝑞=ℎ𝐷_𝑞,𝐺_𝑞=ℎ²𝐺_𝑞. (36) Rewrite (32) and (33) in the form 𝐴⊗𝐼_𝑠+𝐵⊗𝐴0−𝐼_𝑑⊗𝑏̂^Τ𝐼_𝑑𝑌_𝑛+0𝐵⊗𝑒−𝐼_𝑑⊗𝑏̃^Τ−(1−𝜃)𝐼_𝑑𝑌_𝑛−1+000−𝜃𝐼_𝑑𝑌_𝑛−2+^𝑚_𝑞=1𝐶_𝑞⊗𝐼_𝑠+𝐷_𝑞⊗𝑌𝐴000^{𝑛−𝑞𝑀}+^𝑚_𝑞=10𝐷_𝑞𝑌⊗𝑒00^{𝑛−𝑞𝑀−1}+^𝑚_𝑞=1^𝑀𝑞_𝑟=0𝐶_𝑞⊗𝐼_𝑠+𝐷_𝑞⊗𝑌𝐴000^𝑛−𝑟+^𝑚_𝑞=1^𝑀𝑞_𝑟=00𝛾_𝑟𝐺_𝑞𝑌⊗𝑒00^{𝑛−𝑟−1}=0. (37)

The characteristic polynomial of (37) is given by 𝑝(𝑧)𝑇=det¹(𝑧)𝑇₂(𝑧)𝑇₃(𝑧)𝑇₄(𝑧),𝑧∈𝐶, (38) where 𝑇₁(𝑧)=𝑧²𝐴⊗𝐼_𝑠+𝐴+𝐵⊗_𝑚𝑞=1𝐶_𝑞⊗𝐼_𝑠+𝐷_𝑞⊗𝐴𝑧^−𝑞𝑀+^𝑚_𝑞=1^𝑀𝑞_𝑟=0𝛾_𝑟𝐺_𝑞⊗𝐴𝑧^−𝑟,𝑇₂(𝑧)=𝑧𝐵⊗𝑒+^𝑚_𝑞=1𝐷_𝑞⊗𝑒𝑧^−𝑞𝑀+^𝑚_𝑞=1^𝑀𝑞_𝑟=0𝛾_𝑟𝐺_𝑞⊗𝑒𝑟^−𝑟,𝑇₃(𝑧)=−𝑧²𝐼_𝑑⊗𝑏̂^Τ−𝑧𝐼_𝑑⊗𝑏̃^Τ,𝑇₄(𝑧)=𝑧²𝐼_𝑑−𝑧(1−𝜃)𝐼_𝑑−𝜃𝐼_𝑑. (39) Following from the theorem on difference equations, we get that if all the zeros 𝑧 of (38) satisfy |𝑧|<1, then lim_𝑛→∞𝑌_𝑛=0. (40) Hence, we formulate the following lemmas.

Lemma 7 (see [21]). If all the zeros 𝑧 of (38) satisfy |𝑧|<1, then numerical method ((32) and (33)) satisfies lim_𝑛→∞𝑢_𝑛=0. (41)

Lemma 8. Assume that condition (a) of Theorem 5 holds and assume that [𝐼_𝑠−𝜆_l(𝑟(𝑧))𝐴] are invertible for |𝑧|≥1, where 𝑟(𝑧)=𝐴+^𝑚_𝑞=1𝑧^−𝑞𝑀𝐶_𝑞⁻¹×−𝐵−^𝑚_𝑞=1𝑧^−𝑞𝑀𝐷_𝑞−^𝑚_𝑞=1^𝑞𝑀_𝑟=0𝛾_𝑟𝑧^−𝑟𝐺_𝑞; (42) then, det[𝑇₁(𝑧)]≠0, for |𝑧|≥1.

Proof. Condition (a) in Theorem 5 implies that the matrix (𝐴+∑^𝑚_𝑞=1𝑧^−𝑞𝑀𝐶_𝑞) is invertible for |𝑧|≥1; then 𝑇₁(𝑧)=𝑧²[(𝐴+∑^𝑚_𝑞=1𝑧^−𝑞𝑀𝐶_𝑞)⊗𝐼_𝑠][𝐼_𝑑⊗𝐼_𝑠−𝑟(𝑧)⊗𝐴].
We have that det𝑇₁(𝑧)=𝑧^2𝑠𝑑det𝐴+^𝑚_𝑞=1𝑧^−𝑞𝑀𝐶_𝑞^𝑠×det𝐼_𝑑⊗𝐼_𝑠−𝑟(𝑧)⊗𝐴=𝑧^2𝑠𝑑det𝐴+^𝑚_𝑞=1𝑧^−𝑞𝑀𝐶_𝑞^𝑠×^𝑑𝑠_𝑙=1𝑗=11−𝜆_𝑙(𝑟(𝑧))𝜆_𝑗𝐴. (43) The matrix 𝐼_𝑠−𝜆_𝑙(𝑟(𝑧))𝐴 is invertible meaning that 𝜆_𝑙(𝑟(𝑧))𝜆_𝑗(𝐴)≠1 for all 𝑙,𝑗.
Hence, det[𝑇₁(𝑧)]≠0, for |𝑧|≥1.

Theorem 9. If the system ((32) and (33)) satisfies Lemma 8 and the following conditions,(ã)det(𝐴+∑^𝑚_𝑞=1𝜉_𝑞𝐶_𝑞)≠0, for |𝜉_𝑞|≤1;(b)|1−𝜃−𝑧⁻¹𝜃−(𝑏̂^Τ+𝑧⁻¹𝑏̃^Τ)[𝐼_𝑑⊗𝐼_𝑠−𝑟(𝑧)⊗𝐴]⁻¹[𝑒_𝑠⊗𝑟(𝑧)]|<1, for |𝑧|≥1,then the solution of the TSRK methods for (1) is asymptotically stable.

Proof. By Lemma 7, we need to prove that all the zeros of (38) satisfy |𝑧|<1.
If these were not true, there would exist a 𝑧₀∈𝐶 with |𝑧₀|≥1, such that 𝑇det¹𝑧₀𝑇₂𝑧₀𝑇₃𝑧₀𝑇₄𝑧₀=0. (44) By Lemma 8, we have that det[𝑇₁(𝑧₀)]≠0.
Hence, (44) is equivalent to det𝑇₄𝑧₀−𝑇₃𝑧₀𝑇₁⁻¹𝑧₀𝑇₂𝑧₀=0. (45)
Using the Kronecker product [5, chapter 4], we have that det𝑇₄𝑧₀−𝑇₃𝑧₀𝑇₁⁻¹𝑧₀𝑇₂𝑧₀=det𝑧²₀𝐼_𝑑−𝑧₀(1−𝜃)𝐼_𝑑−𝜃𝐼_𝑑−𝑧²₀𝐼_𝑑⊗𝑏̂^Τ+𝑧₀𝐼_𝑑⊗𝑏̃^Τ×𝑧₀⁻²𝐴+^𝑚_𝑞=1𝑧₀^−𝑞𝑀𝐶_𝑞⊗𝐼_𝑠⁻¹⋅𝐼_𝑑⊗𝐼_𝑠−𝑟𝑧₀⊗𝐴⁻¹⋅𝑧₀𝐵⊗𝑒+^𝑚_𝑞=1𝐷_𝑞⊗𝑒𝑧₀^−𝑞𝑀+^𝑚_𝑞=1^𝑀𝑞_𝑟=0𝛾_𝑟𝐺_𝑞⊗𝑒𝑧₀^−𝑟=𝑧^𝑑₀det𝑧₀𝐼_𝑑−(1−𝜃)𝐼_𝑑−𝑧₀⁻¹𝜃𝐼_𝑑−𝐼_𝑑⊗𝑏̂^Τ+𝑧₀⁻¹𝐼_𝑑⊗𝑏̃^Τ×𝐼_𝑑⊗𝐼_𝑠−𝑟𝑧₀⊗𝐴⁻¹×𝑒⊗𝑟𝑧₀. (46) Combining (45) and (46) gives that det𝑧₀𝐼_𝑑−(1−𝜃)𝐼_𝑑−𝑧₀⁻¹𝜃𝐼_𝑑−𝐼_𝑑⊗𝑏̂^Τ+𝑧₀⁻¹𝐼_𝑑⊗𝑏̃^Τ⋅𝐼_𝑑⊗𝐼_𝑠−𝑟𝑧₀⊗𝐴⁻¹𝑒⊗𝑟𝑧₀=0, (47) which implies |1−𝜃−𝑧₀⁻¹𝜃−(𝑏̂^Τ+𝑧₀⁻¹𝑏̃^Τ)[𝐼_𝑑⊗𝐼_𝑠−𝑟(𝑧₀)⊗𝐴]⁻¹[𝑒⊗𝑟(𝑧₀)]|≥1.
This contradicts the assumption that (𝑏̃) for |𝑧₀|≥1. Hence, the theorem is proved.

4. Numerical Experiments

Example 1. Consider the following linear system: 𝐴𝑢^′(𝑡)+𝐵𝑢(𝑡)+^𝑚_𝑞=1𝐶_𝑞𝑢^′𝑡−𝜏_𝑞+^𝑚_𝑞=1𝐷_𝑞𝑢𝑡−𝜏_𝑞+^𝑚_𝑞=1𝐺_𝑞^𝑡_{𝑡−𝜏𝑞}𝑢(𝛿)𝑢(𝑡)(𝑡)[),𝑑𝛿=0,𝑡≥0,=𝜑,𝑡∈−𝜏,0 (48) with 𝑚=2, 𝜏₁=1, 𝜏₂=2, and 𝜑(𝑡)=(cos(𝑡),sin(𝑡),cos(𝑡))^Τ for 𝑡∈[−2,0],
where ,𝐷𝐴=100000010,𝐵=1.50000−0.400.80.5₁=0.20.40000.100.2−0.2,𝐷₂=,𝐺0.100000.100.150.8₁=0.20000−0.200.10.25,𝐺₂=,0.10000−0.200.050.25 (49) and 𝐶₁=−0.3𝐴 and 𝐶₂=−0.5𝐴.
Here the matrix coefficients satisfy Theorem 9. Hence, the system is asymptotically stable.
We choose the A-stable TSRK methods as follows [22]: 𝐶𝐴𝑏̂^T𝜃𝑏̃^T=51325151976026012440103256−1060915616073439156160636886846889107451673728061448158637134314592160169774493669797621872982199107451673728052658918227,134314592160 (50) where 𝑐_𝑖=^𝑠_𝑗=1𝑎_𝑖𝑗,^𝑠_𝑗=1𝑏̂_𝑗+𝑏̃_𝑗=1+𝜃. (51)
It can be easily seen that the A-stable TSRK method is asymptotically stable, which illuminates the conclusion of Theorem 9.

5. Conclusions

This paper develops the asymptotic stability of the two-step Runge-Kutta methods for neutral delay integro differential-algebraic equations with many delays. It studies the asymptotic stability of the analytical solution and introduces two step Runge-Kutta methods and the stability region. It also deals with the asymptotic stability of two-step Runge-Kutta method for neutral delay integro differential-algebraic equations with many delays and proves that the A-stable two-step Runge-Kutta methods are asymptotically stable for neutral delay integro differential-algebraic equations with many delays.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work was supported by the Doctor Foundation of Heilongjiang Institute of Technology (2012BJ27).