Delay-dependent stability of numerical methods for delay differential systems of neutral type

We are concerned with stability of numerical methods for delay differential systems of neutral type. In particular, delay-dependent stability of numerical methods is investigated. By means of the H-matrix norm, a necessary and sufficient condition for the asymptotic stability of analytic solution of linear neutral differential systems is derived. Then, based on the argument principle, sufficient conditions for delay-dependent stability of Runge–Kutta and linear multi-step methods are presented, respectively. Furthermore, two algorithms are provided for checking delay-dependent stability of analytical and numerical solutions, respectively. Numerical examples are given to illustrate the main results.


Introduction
We are concerned with linear delay differential system of neutral type with a single delay described bẏ with the condition ρ(N ) < 1, (1.2) where u(t) ∈ R d , parameter matrices L , M, N ∈ R d×d and delay τ > 0. Throughout the present paper, the NDDE system refers to Eq. (1.1) with (1.2). Stability of delay and neutral systems can be divided into two categories according to its dependence upon the size of delays. The stability which does not depend on delays is called delay-independent, otherwise it is referred to as delay-dependent (they are discussed in [1,2,4,[9][10][11]19]). Each case is extended to the neutral type and the reader can refer to [2,10] for the delay-independent and to [2,9,11,19] for the delaydependent case. Then, the stability of numerical methods is also divided into delayindependent and delay-dependent according to which system the method is applied to. The present paper is devoted to study the delay-dependent stability of numerical methods for the systems of neutral type, for we believe it gives more precise insight of the methods.
Although numerical delay-independent stability has been discussed in [1,2,10], only a few works have been reported for the delay-dependent case [2,7,13,18]. The literature [2,7,18] proposed the delay-dependent stability of numerical methods for the systemu (t) = Lu(t) + Mu(t − τ ), (1.3) and called it as D-stability. However, as pointed out in [2], delay-dependent stability of numerical methods for the system is a real challenge. In fact, in [18] it is proved that for any A-stable natural Runge-Kutta method for system (1.3) and for any step-size h = τ/m (m is a positive integer), there exists an asymptotically stable linear system of dimension d = 4 for which u n → 0 as n → ∞ does not hold. Later, the same negative result for d = 2 is obtained in [7]. In short, no A-stable natural Runge-Kutta method for system (

1.3) is D-stable.
This suggests that the definition of D-stability is too restricted. Its definition requires the resulting difference system from a Runge-Kutta method is asymptotically stable for all asymptotically stable system (1.3) and for all the step-sizes. Hence almost all the standard Runge-Kutta methods are excluded in the D-stability sense. However, practical computations shows it is not the case. Our experiences are as follows. When we numerically solve an asymptotically stable NDDE system, it is possible that by carefully choosing a small step-size h = τ/m numerical solutions by a standard Runge-Kutta or linear multi-step method exhibit an asymptotically stable behaviour. Thus we are motivated to present a relaxed definition for delay-dependent stability of numerical methods. We will give a new definition of delay-dependent stability, which only requires the conditions of asymptotically stable difference system derived from a numerical method for an asymptotically stable NDDE system with a certain integer giving the step-size h = τ/m. This paper is organized as follows. A necessary and sufficient condition for the asymptotic stability of analytic solution of linear neutral differential systems is derived in Sect. 2. In Sect. 3, the new delay-dependent stability of Runge-Kutta and linear multi-step methods are discussed, respectively. Numerical examples in Sect. 4 are provided to illustrate the main results. Conclusions are given in Sect. 5.
Throughout the paper, for a complex matrix F, F means any matrix norm. The Hermitian conjugate of a complex vector or matrix is denoted by x * or F * . The j th eigenvalue of F is denoted by λ j (F). The symbol ρ(F) represents its spectral radius. Re z and Im z stand for the real and the imaginary parts of a complex number z, respectively. The open left half-plane {s : Re s < 0} is denoted by C − and the right half-plane {s : Re s ≥ 0} by C + .

Preliminaries
In this section, several lemmas and definitions are reviewed which will be used in Sects. 3

and 4.
A function P(s) is said to be meromorphic in a domain D if it is analytic throughout D except poles. The following argument principle is well-known. Lemma 2.1 (e.g. [5]) Suppose that (i) a function P(s) is meromorphic in the domain interior to a positively oriented simple closed contour γ ; (ii) P(s) is analytic and nonzero on γ ; (iii) counting multiplicities, Z is the number of zeros and Y is the number of poles of P(s) inside γ.
where γ arg P(s) represents the change of the argument of P(s) along γ .
The following result is useful to check the existence of inverse matrix. Then, the Lemma guarantees the unique existence of H in our case of Lemma 2.3, for we assume a stable matrix F with respect to the unit circle. Furthermore, we note that the matrix H is computable in our case and computing algorithms have been developed (e.g. [6]). In fact, the software Matlab has the functions slstst (for stable Stein equations) and slgest (for generalized Stein equations) ( [3]). Thus an introduction of H -norm by the following definition is meaningful in conjunction with the stable matrix with respect to the unit circle.
Definition 2.2 (e.g. [12]) For any vector x, any matrix F and any positive definite matrix H , the H -norm of x and F are defined, respectively, by By the above definition, we have the following formula for estimating the H -norm of F.
where λ max (W ) stands for the maximal eigenvalue of a symmetric matrix W.
Proof By Definition 2.2, we can calculate which is the desired result.
For a stable matrix, the following Lemma asserts an evaluation of its H -norm by ρ(H ).
which implies the assertion of (2.5).

Stability criterion
In this section, by means of the H -matrix norm, a necessary and sufficient condition for asymptotic stability of the NDDE system is presented to extend the results in [11]. We will consider the stability of the system (1.1). Its characteristic equation reads whose root is called a characteristic root. The asymptotic stability of the NDDE system is determined by the location of the characteristic roots. That is, the system is asymptotically stable if and only if all the characteristic roots lie in the open left complex half-plane C − [9]. If there exists unstable characteristic roots of (3.1), i.e., P(s) = 0 for Re s ≥ 0, they are bounded due to the following result. .
rewrites P(s) as Since ρ(N exp(−τ s)) < 1, the equation must hold. This implies the s is an eigenvalue of the matrix W (s) and there exists an Since Re s ≥ 0, we have According to (3.5) and (3.6), we can estimate Thus the proof is completed.
The theorem means that there is a bounded region in the right-half complex plane C + which includes all the unstable characteristic roots of (3.1).

Remark 3.2
Since ρ(N ) ≤ N holds generally, the condition ρ(N ) < 1 in Theorem 3.1 is weaker than the condition N < 1 in [11]. Hence the theorem is an extension of the results in [11]. Moreover, the boundary of D β is denoted by Γ β . See Fig.1 By means of the argument principle, a computable criterion for system (1.1) can be derived. Now we present the main result in this section. The following theorem will exclude all the unstable characteristic root of Eq. (3.1) from the set D β .

Algorithm 1
Step 0 Solve Stein equation Then as the boundary of D β we have the curve Γ β , which consists of two parts, i.e., the segment {s = it; −β ≤ t ≤ β} and the half-circle {s; |s| = β and − π/2 ≤ arg s ≤ π/2}. Step 1 Take a sufficiently large integer n ∈ N and distribute n node points {s j } ( j = 1, 2, . . . , n) on Γ β as uniformly as possible. For each s j , we evaluate P(s j ) by computing the determinant as Also we decompose P(s j ) into its real and imaginary parts.
Step 2 We examine whether P(s j ) = 0 holds for each s j ( j = 1, . . . , n) by checking its magnitude satisfies |P(s j )| ≤ δ 1 with the preassigned tolerance δ 1 . If it holds, i.e., s j ∈ Γ β is a root of P(s), then the NDDE system is not asymptotically stable and stop the algorithm. Otherwise, to go to the next step.
Step 3 We examine whether Γ β arg P(s) = 0 holds along the sequence {P(s j )} by checking | Γ β arg P(s)| ≤ δ 2 with the preassigned tolerance δ 2 . If it holds, it means that the change of the argument is 0 along the sequence {P(s j )}, then the system is asymptotically stable, otherwise not stable.

Delay-dependent stability of numerical methods
In this section, a delay-dependent stability of both Runge-Kutta (RK) and linear multistep (LM) methods for linear delay differential systems of neutral type is discussed. Based on the argument principle, sufficient conditions for delay-dependent stability of the numerical methods are presented. First, we assume the numerical solution we are now discussing gives a sequence of approximate values {u 0 , u 1 , . . . , u n , . . .} of {u(t 0 ), u(t 1 ), . . . , u(t n ), . . .} of (1.1) on certain equidistant step-values {t n (= nh)} with the step-size h. We introduce the following new definition of delay-dependent stability of numerical methods for system (1.1).

Definition 4.1
Assume that the NDDE system is asymptotically stable for given matrices L , M, N and a delay τ . A numerical method is called weakly delay-dependently stable for system (1.1) if there exists a positive integer m such that the step-size h = τ/m and the numerical solution u n with h satisfies u n → 0 as n → ∞ for any initial function.

Remark 4.1
In [2,7,18], the definition of delay-dependent stability, which is called D-stability there, is too restricted since it requires for all asymptotically stable delay differential system (1.3) and for all the natural numbers m the resulting numerical solution is asymptotically stable. Therefore, almost all the standard RK methods are excluded in the D-stability sense. It is obvious that Definition 4.1 is weaker than D-stability in the literature.

Runge-Kutta case
For the initial value problem of ordinary differential equations (ODEs) f (t, y(t)), for t ≥ 0 (4.1) an s-stage (implicit) Runge-Kutta method for ODEs (4.1) is defined (e.g., [16]) by We associate the s-square matrix A and the s-dimensional vector b as We analyse RK method which is extended to apply to the NDDE system. Taking advantage of the constant step-size h = τ/m, we can employ the so-called natural Runge-Kutta scheme [1,2] for (1.1) given by and where the symbol K ,i means the i-th stage value of the RK at the -th step-point.

Lemma 4.1
The characteristic polynomial P RK (z) of the linear difference system (4.4) and (4.5) is given by where the s-dimensional vector e is defined by e = (1, 1, . . . , 1) T . Proof The difference system (4.4) and (4.5) can be expressed with the Kronecker product as follows: where the (ds)-dimensional vector K n means Hence the dimension of the vector K n u n+1 becomes (s + 1)d.
Taking z-transform to (4.7) and introducing as Z Hence, the characteristic polynomial of difference system is given as desired.
For an explicit RK method, that is, the one with a i j = 0 for i ≤ j, we have the following result.

Theorem 4.1 For an explicit RK method, assume that
(i) the NDDE system is asymptotically stable for given matrices L , M, N and delay τ (therefore, Theorem 3.2 holds); (ii) the RK method is of s-stage and natural with the step-size h = τ/m; (iii) The characteristic polynomial P RK (z) never vanishes on the unit circle μ = {z : |z| = 1} and its change of argument satisfies Then the RK method for the NDDE system is weakly delay-dependently stable. Proof The difference system (4.7) is asymptotically stable if and only if all the characteristic root of P RK (z) = 0 lie in the inner of the unit circle. Notice that the coefficient matrix of the term z m+1 in P RK (z) = 0 is Since is also nonsingular. Since the degree of P RK (z) is d(s+1)(m+1), it has d(s+1)(m+1) roots in total by counting their multiplicity. By the argument principle, the condition (iii) implies that the condition |z| < 1 holds for all the d(s + 1)(m + 1) roots of P RK (z) = 0. The proof completes.
For an implicit RK method applied to the NDDE system, we derive the following result.

Theorem 4.2 For an implicit RK method, assume that (i) the conditions (i),(ii) and (iii) in Theorem 4.1 hold; (ii) the product hλ i (A)λ j (L) never equals to unity for all i (1 ≤ i ≤ s) and j (1 ≤ j ≤ d).
Then the RK method for the NDDE system is weakly delay-dependently stable.
Proof The proof can be carried out similarly to that of Theorem 4.1. Notice that the condition (ii) ensures the matrix is nonsingular since the matrix I sd − h(A L) is nonsingular. Thus the degree of the polynomial P RK (z) becomes d(s + 1)(m + 1).

Linear multi-step case
When applied to ODEs (4.1), a linear k-step method is given as follows (e.g., [16]): k j=0 α j y n+ j = h k j=0 β j f n+ j , (4.9) where h stands for the step-size, α j , β j are the formula parameters satisfying α k = 1 and |α 0 | + |β 0 | = 0. and f j = f (t j , y j ). As in the case of RK method, in order to apply the LM method to the delay differential systems of neutral type (1.1), we employ the step-size h of an integral fraction of τ , i.e., h = τ/m for a certain m ∈ N.
For the application to the linear system (1.1), we introduce the same reformulation given by [10]. That is, with an auxiliary variable v(t) the system (1.1) is rewritten as the simultaneous system:u (4.10) and Then, an application of the method (4.9) to (4.10) and (4.11) yields and v n+ j = Lu n+ j + Mu n−m+ j + N v n−m+ j . (4.13) This reformulation can ease the derivation of its characteristic polynomial.

Lemma 4.2
The characteristic polynomial of (4.12) and (4.13) is given by The above is the difference equation to be solved after eliminating v's. Taking its Z -transform and putting Z {u n−m } = U (z), we obtain ⎡ Hence, the characteristic polynomial is given as desired.
For delay-dependent stability of an LM method, we have the following result.

Algorithm 2
Step 1 Taking a sufficiently big n ∈ N, we compute n nodes {z 0 , z 1 , . . . , z n−1 } upon the unit circle μ of z-plane so as arg z = (2π) /n. For each z ( = 0, 1 . . . n − 1), we evaluate the characteristic polynomials of the numerical scheme. That is, in RK case, we evaluate it by computing the determinant 0 , while in LM case by computing the determinant Also we decompose P(z ) into its real and imaginary parts.
Step 2 We examine whether P(z ) = 0 holds for each z ( = 0, 1, . . . , n − 1) by checking its magnitude satisfies |P(z )| ≤ δ 1 with the preassigned small positive tolerance δ 1 . If it holds, i.e., z ∈ μ is a root of P(z), then the numerical scheme for the NDDE system is not asymptotically stable and stop the algorithm. Otherwise, to go to the next step.
Step 3 We examine whether 1 2π μ arg P(z ) = d(s + 1)(m + 1) in RK case (or d(k + m) in LM case ) holds along the sequence {P(z )} by checking If it holds, it means that the change of the argument along the sequence {P(z )} is d(s + 1)(m + 1) in RK case (or d(k + m) in LM case ), then the numerical scheme for the NDDE system is asymptotically stable, otherwise not stable.

Remark 4.3
From the above three theorems, in order to solve numerically an asymptotically stable delay differential system of neutral type by an RK or LM method, it is enough for us to choose a positive integer m such that the resulting difference system is asymptotically stable. When the positive integer m is determined, the step-size h = τ/m is obtained such that the RK or LM method can work. Theorems 4.1, 4.2 and 4.3 do not contradict the results reported in [7,18] which assert for any A-stable natural RK method for system (1.3) and for any step-size h = τ/m, (where m is a positive integer) there exists an asymptotically stable linear system for which the resulting system is not asymptotically stable, i.e., u n → 0 as n → ∞ does not hold. [15] need information of all the coefficients of the characteristic polynomial P RK (z) (or P L M (z)). It is a well-known ill-posed problem to compute all the coefficients of the characteristic polynomial for a high dimensional matrix [8]. Although Schur-Cohn and Jury stability criteria can be applied to the resulting difference systems from RK and LM methods in theoretical sense, they can not work well in practice when m or d are big. Algorithm 2 avoids computation of the coefficients of the characteristic polynomial P RK (z) (or P L M (z)). It evaluates the determinant P RK (z ) or P L M (z ) through the elementary row (or column) operations which are relatively efficient ways [14].

Numerical examples
In this section, two numerical examples are given to demonstrate the main results in Sects. 2 and 3. The matrix norm F 2 is defined by which is equal to the special case of the H -norm for H = I , i.e., F I = F 2 holds. The number of significant digits is 4, and all the equality signs hold in this sense. The classical fourth-order RK formulas for ODEs (4.1), which is given as is our underlying scheme of the natural RK for (1.1).

Example 1
The two-dimensional linear neutral system (1.1) with the parameter matrices We have that ρ(N ) = 0. Let the initial vector function be The case of τ = 1.1. Algorithm 1 is employed to check stability of the system. As we have Γ β arg P(s) = 0 along the curve Γ β , Theorem 3.2 tells that the system with the given parameter matrices is asymptotically stable. Now using Algorithm 2 with n = 3.2 × 10 5 to check delay-dependent stability of the RK method for the system. When m = 100, we obtain that μ arg P RK (z) = d(s + 1)(m + 1) = 2(4 + 1)(100 + 1) = 1010. Theorem 4.1 asserts the RK method is weakly delay-dependently stable. The numerical solution, which is converging to 0, is depicted in Fig. 2. Conversely, when m = 10, we obtain μ arg P RK (z) = 97 = d(s + 1)(m + 1) = 2(4 + 1)(10 + 1) = 110 and the theorem does not hold. The numerical solution is divergent and its behaviour is shown in Fig. 3. The case of τ = 3.0. Again we employ Algorithm 1 to check stability of the system. As we have Γ β arg P(s) = 2 along the curve Γ β , Theorem 3.2 tells that the system with the given parameter matrices is not asymptotically stable. Then the assumptions of Theorem 4.1 do not hold and the numerical solution is not guaranteed to be asymptotically stable. In fact, its figure given in Fig. 4 shows a divergence for m = 100. We also carried out several numerical experiments for m > 100, whose numerical solutions are still divergent.

Example 2
We take the following four-dimensional linear neutral system with the parameter matrices given by Since N 2 < 1, we can take H = I and we have Let the initial vector function be The case of τ = 0.1. The procedure goes similarly to Example 1. Since our computation gives Γ β arg P(s) = 0 along the curve Γ β , according to Theorem 3.2, we can confirm that the system with the given parameter matrices is asymptotically stable. Now using Algorithm 2 with n = 3.2 × 10 5 to check delay-dependent stability of the RK method for the system. In the case of m = 100, our computation gives that μ arg P RK (z) = d(s + 1)(m + 1) = 4(4 + 1)(100 + 1) = 2020. Thus, by Theorem 4.1, the RK method is weakly delay-dependently stable. In fact, the numerical solution is convergent as shown in Fig. 5. Even when m = 1, we can compute μ arg P RK (z) = d(s + 1)(m + 1) = 4(4 + 1)(1 + 1) = 40. Hence again the RK method is weakly delay-dependently stable and the numerical solution is convergent as shown in Fig. 6. The case of τ = 0.3. Since the computation in Algorithm 1 gives Γ β arg P(s) = 2 along the curve Γ β , according to Theorem 3.2, we can confirm that the system with the given parameter matrices is not asymptotically stable. Then, the assumptions of Theorem 4.1 do not hold and the numerical solution is divergent for m = 100, whose computed results are shown in Fig. 7. We also carried out several other numerical experiments for m > 100, but the numerical solutions are still divergent.
The two numerical examples show that our main results work well in the actual computations. Therefore, we can employ them to seek a step-size h so as to give numerically stable solution by Runge-Kutta or linear multi-step method. That is, assume that we are given an asymptotically stable system (1.1), then we can select an appropriate natural number m for the step-size h = τ/m which yields an asymptotically stable numerical solution by utilizing our criterion.

Conclusions
For the asymptotic stability of analytical solution of the NDDE system (1.1) with (1.2), Theorem 3.2 provides a novel necessary and sufficient condition, which extends the results in the literature. Then, we present a relaxed definition, Definition 4.1, for delaydependent stability of numerical methods. By applying this, Theorems 4.1, 4.2 and 4.3 provide sufficient conditions for delay-dependent stability of RK and LM methods, respectively. They show that it is possible to seek a step-size h = τ/m such that the resulting difference system by a standard RK or a LM method is asymptotically stable for an asymptotically stable NDDE system. Theorems can not guarantee an existence of such a step-size h = τ/m for all asymptotically stable NDDE systems. However, for many cases of the system the theorems can hold. Furthermore, two algorithms are provided for checking delay-dependent stability of analytical and numerical solutions, respectively. and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.