1 Introduction

It is well known that delay integro-differential equations (DIDEs) have a wide range of practical applications, such as control theory, epidemiology, and population dynamics. An extensive development in the theoretical analysis of numerical methods for DIDEs has been made in past several decades. What most attracts our attention is the theory of stability of discretization schemes. The stability of delay differential systems is called delay-dependent if the size of delay has an effect on stability, or it will be called delay-independent. Therefore, the stability of a numerical method is also divided into delay-dependent and delay-independent depending on which system the numerical method is applied to.

In recent years, a great deal of results about the stability of numerical methods for DIDEs have been obtained^[1-6], and linear multistep methods have been widely considered as a class of traditional methods. Baker and Ford^[7] first investigated the stability of linear multistep methods for scalar DIDEs. Luzyanina et al.^[8] further gave a series of stability results for linear multistep methods under a sufficient condition. Huang^[9] proved that every A-stable, strongly 0-stable linear multistep method of Pouzet type preserves the delay-independent stability of the linear multi-dimensional DIDEs. For other related topics, we refer the readers to Refs. [10]-^[15]. The delay-dependent stability of numerical methods for DIDEs has been discussed^[16-19], but few results of the delay-dependent stability have been presented for DIDEs of multi-dimensions.

The delay-dependent stability of numerical methods for the system

(1)

was proposed^{[1, 20]} and named D-stability. A numerical method for a delay differential system is D-stable on the condition that the resulting difference system from it is asymptotically stable for all asymptotically stable systems (1) and for all the step-sizes. It suggests that most of numerical methods cannot satisfy this restrictive condition.

In order to seek a weaker stability that can be used to a wider range of numerical methods for delay differential systems, Hu and Mitsui^[21] recently introduced weak delay-dependent stability. It only requires a difference system given by a numerical method to satisfy the conditions of asymptotically stable with a certain integer m, when the numerical method is applied to an asymptotically stable differential system and the integer m provides a step-size h=τ/m. In the present paper, we continue this study and consider the weak delay-dependent stability of linear multistep methods for the following differential system with distributed delays:

(2)

where the constant delay τ>0, , L and M are constant d×d matrices, K(s) is also a d×d matrix with piece-wise continuous elements on [-τ, 0], and (θ) is a continuous function on [-τ, 0].

Throughout the paper, I stands for the identity matrix, C is a complex plane, ||A|| means the F-norm of the matrix A, and the ith eigenvalue of A is denoted by λ_i(A). The real and imaginary parts of a complex number z are denoted by Re(z) and Im(z), respectively. The open left half-plane {s:Re(s) < 0} and the right half-plane {s:Re(s)≥0} are denoted by C^- and C⁺, respectively.

2 Preliminaries

In this paper, we discuss the stability of the system (2). Substituting the solution of exponential form x(t)=e^ztξ into the system (2), one gets the characteristic equation of system (2),

(3)

where z∈ C, and ξ∈^d is a nonzero vector. The roots of characteristic equation (3) are called characteristic roots. The asymptotic stability of the system (2) is completely determined by the distribution of the characteristic roots. That is to say, it is sufficient and necessary that all the characteristic roots lie in C^- for the asymptotic stability of system (2). As a special case of the conclusions in the literature^[22], we know that if there are unstable characteristic roots of (3), i.e., P(z)=0 for Re(z)≥0, they are bounded.

Lemma 1^[22]  Let z be an unstable characteristic root of (3) of system (2), i.e., Re(z)≥ 0. If there exists a nonnegative constant matrix with such entries that |K(s)|≤, i.e., k_ij(s) and , which are the entries of matrix K(s) and , respectively, and satisfy |k_ij(s)|≤ for all s∈[-τ, 0] and i, j=1, 2, …, d, then

(4)

The lemma means that there exists a bounded region in the right-half complex plane C⁺ which includes all the unstable characteristic roots of (3). According to Ref. [21], we obtain that the bounded region is a half-disk D_β={s: Re(s)≥ 0 and |s|≤ β}. The boundary of D_β is denoted by Γ_β.

A meromorphic function on an open subset D of the complex plane is a function that is analytic on D except for some poles. The argument principle will be stated in the following.

Lemma 2^[23] If P(s) is a meromorphic function in a region enclosed by a contour γ, and P(s) has no zeros and poles on γ, then

where Z and Y denote the numbers of zeros and poles of P(s) inside the contour γ, respectively, each zero and each pole are counted as many as their multiplicity, and Δ_γε_argP(s) represents the change of the argument of P(s) along γ.

Lemma 3^[22] Assume that the conditions of Lemma 1 hold, and the bound β is obtained. Then, the system (2) is asymptotically stable if and only if

(5)

(6)

hold, in which P(z) is defined by (3) and Δ_ΓβεargP(z) has the same meaning as in Lemma 2.

Now we give the definition of weak delay-dependent stability of numerical methods for (2).

Definition 1^[21] Suppose that the system (2) is asymptotically stable for given matrices L, M, K, and a fixed delay τ. A numerical method is called weakly delay-dependently stable for the system (2) if there exists a positive integer m such that the numerical solution x_n with h=τ/m satisfying x_n→0 as n→ ∞ for any initial function.

3 Delay-dependent stability of linear multistep methods

In this section, we present some weakly delay-dependently stable conditions of linear multistep methods for the system (2) based on the argument principle.

Considering the initial value problem of ordinary differential equations (ODEs),

(7)

a linear k-step method in the standard form is defined (e.g., Ref. [24]) by

(8)

where α_j and β_j are constants subject to the conditions α_k=1 and |α₀|+|β₀|≠0.

Investing the DIDEs with the initial condition,

(9a)

(9b)

we always assume that (9) has a (unique) solution y, where τ>0, y∈, and f and g are vector functions defined in the real domain.

The following schemes (e.g., Refs. [7] and [9]) are obtained by adopting the method (8) to the system (9):

(10)

where h=τ/m with m being a given positive integer, t_n=nh, n∈, and y_n+j and w_n+j are approximations to y(t_n+j) and , respectively.

As to the computation of the integral w(t_n+j), we use the compound quadrature formula,

(11)

with weights {v_r} independent of m and determined by those of a convergent quadrature rule

(12)

Here, h=τ/m, and ϕ is continuous on [0, τ]. In order to get a convergent formula, we make the same assumption as in Ref. [25] that there exists a finite constant v, independent of m, such that

(13)

Remark For the method (10)-(12), suppose y_l=ψ(t_l) for -m≤ l < 0 and the following assumptions are satisfied.

(ⅰ) The order of a zero-stable linear multistep method is p>0.

(ⅱ) The formula (12) satisfies .

(ⅲ) The additional starting values satisfy .

Then, according to Theorem 1.1 in Ref. [7], we can make clear that for any T < ∞,

where s=min{r, p, q}, and ξ_q is a constant.

Now applying the method (10)-(12) to the system (2), we get the following approximation scheme:

(14a)

(14b)

Lemma 4 The characteristic polynomial P_LMC(z) of the system (14) is given by

(15)

where I stands for the d-dimensional identity matrix.

Proof We derive the following difference equations by substituting (14b) into (14a):

(16)

Taking z-transform to (16) and introducing as z(x_n-m)=V(z), we obtain

Hence, the characteristic polynomial of the difference system (14) is given as desired.

For explicit linear multistep methods used to the system (2), we get the following result.

Theorem 1 Assume that the conditions of Lemma 3 hold. Then, an explicit linear multistep method for the system (2) is weakly delay-dependently stable if there exists a positive integer m such that the following conditions hold:

(ⅰ) the linear multistep method is k-step with the step-size h=τ/m};

(ⅱ) the characteristic polynomial P_LMC(z) never vanishes on the unit circle μ={z:|z|=1} and its change of argument satisfies

(17)

Proof The difference system (14) is asymptotically stable if and only if all the characteristic roots of P_LMC(z)=0 lie in the inner of the unit circle. Notice that the coefficient matrix of the term z^m+k in P_LMC(z)=0 is α_kI-hβ_kL-h²β_kv₀K(0). It is nonsingular because α_k=1 and β_k=0 hold in explicit methods. Therefore, the degree of P_LMC(z) is d(m+k), and it has d(m+k) roots in total by counting multiplicity. By the argument principle, the condition (ⅱ) implies that all the d(m+k) roots of P_LMC(z)=0 lie inside the unit circle μ={z:|z|=1}. Therefore, the difference system (14) is asymptotically stable. The proof completes.

When linear multistep methods are implicit, we derive the following result.

Theorem 2 Assume that the conditions of Lemma 3 hold. Then, an implicit linear multistep method for the system (2) is weakly delay-dependently stable if there exists a positive integer m such that the following conditions hold:

(ⅰ) the conditions in Theorem 1 hold;

(ⅱ) the eigenvalues λ_i of matrix hβ_kL+h²β_kv₀K(0) are never equal to unity for all i (1≤ i≤ d).

Proof The proof can be carried out similar to that of Theorem 1. Notice that the condition (ⅱ) ensures that the matrix I-hβ_kL-h²β_kv₀K(0), which is the coefficient matrix of z^m+k in P_LMC(z)=0, is nonsingular. Then, similar to the proof of Theorem 1, the implicit linear multistep methods for system (2) are weakly delay-dependently stable. The proof completes.

Now an algorithm is provided to check the conditions of Theorems 1-2.

Algorithm 1

Step 1 For a sufficiently large integer n∈, the nodes z₀, z₁, …, z_n-1 distributed upon the unit circle μ of the z-plane satisfy ε_argz_l=(2π)l/n. For each z_l, the characteristic polynomials of the numerical method are evaluated,

Also, we decompose P(z_l) into its real and imaginary parts.

Step 2 For each z_l, we check whether the determinant P(z_l) vanishes by examining its magnitude with the given small positive tolerance δ₁. If |P(z_l)|≤ δ₁ holds, z_l∈μ is a root of P(z). Therefore, the numerical method is not asymptotically stable and the algorithm stops. Otherwise, it proceeds to next step.

Step 3 Along the sequence {P(z_l)}, we examine whether holds with the preassigned tolerance δ₂. If holds, it shows that the change of the argument is d(m+k). Then, the numerical method for the system (2) is weakly delay-dependently stable. Otherwise, it is not weakly delay-dependently stable.

4 Numerical examples

In the following example, our underlying scheme of the method (10)-(12) is a 4-stage explicit linear multistep method with a compound trapezoid formula, and the coefficients of the linear multistep method are α₄=1, α₃=-1, α₂=α₁=α₀=0, β₄=0, β₃=55/24, β₂=-59/24, β₁= 37/24, β₀=-9/24. Moreover, the additional starting values are obtained by the classical fourth-order Runge-Kutta method.

We consider a three-dimensional delay differential system (2) with the parameter matrices,

where s ∈ [−τ, 0].

The initial vector function is defined as

From Lemma 2, we have

||L||=11.045 4, ||M||=8.660 3, and ||||=6.284 9. Now we examine the stability of the system by employing Algorithm 1.

Case 1 τ=0.5

The radius of unstable region β=22.848 2 is obtained by Lemma 1. Lemma 3 tells that the system with the known parameter matrices is asymptotically stable. When m=20, , the conditions of Theorem 1 are not satisfied, and the numerical solution shown in Fig. 1(a) is divergent. However, when m=352, , the assumptions of Theorem 1 are all satisfied, and the numerical solution shown in Fig. 1(b) is asymptotically stable.

Case 2 τ=2

The radius of unstable region β=32.275 5, and . Thus, we derive that the system is unstable according to Lemma 3. Theorem 1 does not work. In fact, its figures show a divergence for m=50, 350 given in Fig. 2. Several numerical experiments are taken for m>10, the results of which show a divergence of the numerical solutions.

5 Conclusions

The numerical example reveals that our main results work well in the actual computations. However, it is necessary to note that the criteria for weakly delay-dependent stability presented here are sufficient but not necessary. Not all linear multistep methods can use these criteria to check stability.