Stability of the Radau IA and Lobatto IIIC methods for neutral delay differential system

Lin QIU and Taketomo MITSUI^{∗} Graduate School of Human Informatics, Nagoya University Nagoya 464-8601, Japan

January 2000

Abstract

Numerical stability is considered for several Runge-Kutta methods to systems of neutral delay differential equations. The linear stability analysis is adopted to the system. Adapted with the equistage interpolation process as well as the continuous extension, the Runge-Kutta methods are shown to have the numerical stability similar to the analytical asymptotic stability with arbitrary stepsize, when certain assumptions hold for the logarithmic matrix norm on the coefficient matrices of the NDDE system.

1991 Mathematics Subject Classification: 65L05, 65L06, 65L20, 34K28 Keywords: Numerical stability; NGP-stability; Runge-Kutta with continuous extensions; Equistage interpolation process

1 Introduction

We are concerned with stability behaviours in numerical methods for the solution of initial value problems in neutral delay differential equations (NDDEs) of the form

u′(t) = f(t, u(t), u(t − τ), u′(t − τ)), t ≥ 0,  u(t) = g(t), −τ ≤ t ≤ 0,  (1.1)

where the delay τ is a given positive constant, f and g denote given vector-valued functions, and u(t) is the vector-valued unknown function to be solved for t ≥ 0.

In particular, stability of numerical methods will be studied based on the linear NDDEs of the form

u′(t) = Lu(t) + Mu(t − τ) + Nu′(t − τ), t ≥ 0,  (1.2)

where L, M and N are d × d complex matrices. Relevant works can be found in [1, 2, 9, 10, 11, 12, 13, 14, 16, 17, 23].

∗Supported by Grant-in-Aid for Scientific Research (A), Japan Society for the Promotion of Science (No.

11304004).

Bellen, Jackiewicz and Zennaro [1] investigated the stability properties of Runge-Kutta (RK) method for (1.2) with respect to the scalar complex parameters L, M and N, and have shown that every collocation RK method whose underlying RK is A-stable becomes NP-stable, i. e., the methods have a similar stability property to A-stability, when the delay is an integer-multiple of the stepsize.

For the system of NDDEs (1.2), the collocation RK method whose underlying RK is A-stable has been considered by Hu and Mitsui [9], Qiu, Yang and Kuang [18]. On the other hand, Kuang, Xiang and Tian [16] have shown that under the condition ∥ N ∥< 1, the system is asymptotically stable if and only if

(S1) every eigenvalue of the matrix L has negative real part,
(S2) for any pure imaginary λ, the spectral radius of (λI − L)^{−1}(M + λN) is

less than 1 in magnitude, and

(S3) the spectrum of L^{−1}M does not include − 1.

Note that the condition (S2) is readily verified to be equivalent to the condition that for the parametrized matrix Q(ξ) = (I − ξN)^{−1}(L + ξM)

(S2^{∗})  all the non-zero eigenvalues of the matrix Q(ξ) have negative real part
whenever |ξ| ≤ 1.

And we also have proved that a continuous collocation RK method applied to (1.2) is NGPstable if and only if the underlying method is A-stable (see [18]).

Moreover, for the case of underlying RK methods which are A-stable but its continuous extension is not natural, Koto [14] has showed that the Radau IA and Lobatto IIIC methods equipped with suitable continuous extensions are NP-stable for the scalar case of (1.2). This drives us to the question whether they also have a similar stability property for the system of NDDEs (1.2) with complex matrices L, M and N.

In this paper, we will show that the Radau IA and Lobatto IIIC methods have a similar stability property for the system (1.2), and that the property is also valid for arbitrary stepsize.

2 Continuous Runge-Kutta

For the initial-value problem (1.1), we consider an s-stage RK method. Let h denote the stepsize. To cope with arbitrary stepsize, we introduce the fractional part of the delay τ with respect to h. That is, for a certain natural number m, let δ satisfy (m − δ)h = τ and δ ∈ [0, 1).

The recurrence formula for RK at t_{n} = nh is expressed with

u_{n+1} = u_{n} + h ∑_{i=1}^{s} b_{i}U′_{n,i}  (2.1)

U′_{n,i} = f (t_{n} + c_{i}h, u_{n} + h ∑_{j=1}^{s} a_{ij}U′_{n,j},

u_{n−m+δ} + h ∑_{j=1}^{s} ω_{ij}U′_{n−m+δ,j}, ∑_{j=1}^{s} c_{ij}U′_{n−m+δ,j})  (2.2)

where c_{i} = ∑_{j=1}^{s}a_{ij} holds for 1 ≤ i ≤ s except for the cases of the 1-stage Radau IA and 2-stage

Lobatto IIIB methods. Here we introduce the vector b = (b_{1}, b_{2}, ⋯, b_{s})^{T} and the matrix A = (a_{ij})_{1≤i,j≤s} which define the underlying RK method for ordinary differential equations (ODEs).

The second argument in f of (2.2) can be interpreted as an approximation to u(t) at the intermediate point t_{n} + c_{i}h ([4]). Similarly, the third and the fourth arguments in f can stand for an approximation to u(t_{n} + c_{i}h − τ) and u′(t_{n} + c_{i}h − τ), respectively. However, we need a certain interpolation process to these values. Hereafter we adopt the so-called equistage interpolation procedure (in ’t Hout[7]). That is, for a pair of natural numbers (r, t), letting L_{p}(x) be the p-th Lagrange polynomial of degree (r + t) given by

L_{p}(x) = ∏_{k=−r,k≠p}^{t} \frac{x − k}{p − k} , (p = −r, −r + 1, …, t − 1, t)

we employ the approximations

u_{n−m+δ} = ∑_{p=−r}^{t} L_{p}(δ)u_{n−m+p},  (2.3)

U′_{n−m+δ,j} = ∑_{p=−r}^{t} L_{p}(δ)U′_{n−m+p,j}. (j = 1, …, s)  (2.4)

Note that we assume m ≥ t + 1.

Usually, the coefficients ω_{ij} and c_{ij} in (2.2) are obtained through the polynomial ω_{j}(θ), which defines the natural continuous extension (NCE, [21]) of the RK, by ω_{ij} = ω_{j}(c_{i}) and c_{ij} = ω′_{j}(c_{i}). In the two exceptional cases (1-stage Radau IA and 2-stage Lobatto IIB), we take the linear interpolation for their continuous extension. NCE will be discussed later.

By introducing the notations L̄ = hL, M̄ = hM, e = (1, 1, …, 1)^{T} ∈ IR^{s},

W = (ω_{j}(c_{i}))_{1≤i,j≤s},  Γ = (ω′_{j}(c_{i}))_{1≤i,j≤s}  (2.5)

and the Kronecker product symbol ⊗, the continuous RK method (2.1)-(2.2) applied to (1.2) can be expressed with the difference equation

U′_{n} − (Γ ⊗ N) ∑_{p=−r}^{t} L_{p}(δ)U′_{n−m+p} =

(e ⊗ L̄)u_{n} + (A ⊗ L̄)U′_{n} + (e ⊗ M̄) ∑_{p=−r}^{t} L_{p}(δ)u_{n−m+p}

+(W ⊗ M̄) ∑_{p=−r}^{t} L_{p}(δ)U′_{n−m+p},  (2.6)
u_{n+1} = u_{n} + (b^{T} ⊗ I_{d})U′_{n}.  (2.7)

Here U′_{l} stands for the (ds)-dimensional vector of the stage values at t_{l}. Eqs. (2.6) and (2.7) are also written in

( \begin{matrix}I − A ⊗ L̄ & 0 \\ −b^{T} ⊗ I_{d} & I_{d} \end{matrix} ) ( \begin{matrix}U′_{n} \\ u_{n+1} \end{matrix} ) =

( \begin{matrix}0 & e ⊗ L̄ \\ 0 & I_{d} \end{matrix} ) ( \begin{matrix}U′_{n−1} \\ u_{n} \end{matrix} ) + ( \begin{matrix}Γ ⊗ N + W ⊗ M̄ & 0 \\ 0 & 0 \end{matrix} ) ( \begin{matrix}∑_{p=−r}^{t} L_{p}(δ)U′_{n−m+p} \\ ∑_{p=−r}^{t} L_{p}(δ)u_{n−m+1+p} \end{matrix} )

+ ( \begin{matrix}0 & e ⊗ M̄ \\ 0 & 0 \end{matrix} ) ( \begin{matrix}∑_{p=−r}^{t} L_{p}(δ)U′_{n−m−1+p} \\ ∑_{p=−r}^{t} L_{p}(δ)u_{n−m+p} \end{matrix} ) .  (2.8)

Therefore the characteristic equation of the above difference equation is given by

det [z^{m}Q(z) + P (z)] = 0,

where the two (d(s + 1))-square matrices are

Q(z) = ( \begin{matrix}(I − A ⊗ L̄)z & −e ⊗ L̄ \\ −(b^{T} ⊗ I_{d})z & (z − 1)I_{d} \end{matrix} ) ,

P (z) = ( \begin{matrix}−(Γ ⊗ N + W ⊗ M̄) ∑_{p=−r}^{t} L_{p}(δ)z^{p+1} & −(e ⊗ M̄) ∑_{p=−r}^{t} L_{p}(δ)z^{p} \\ 0 & 0 \end{matrix} ) .

For notational convenience, we define the following matrix-to-matrix functions M_{X}(Y, Z) and R_{X}(Y, Z) associated with the continuous RK,

M_{X}(Y, Z) = I_{sd} − A ⊗ X − W ⊗ Y − Γ ⊗ Z,  (2.9)
R_{X}(Y, Z) = I_{sd} + (b^{T} ⊗ I_{d})M_{X}(Y, Z)^{−1}(e ⊗ I_{d})(X + Y).  (2.10)

Here X, Y and Z are d-square matrices, while M_{X} and R_{X} (sd)-square matrices.

In the sequel, we restrict ourselves in the case of r ≤ t ≤ r + 2 for the equistage interpolation procedure. As a matter of fact, the following lemma was shown in [19].

Lemma 2.1 Let γ(z, x) = ∑_{p=−r}^{t}L_{p}(x)z^{p+r} (z ∈ C I ). For any z and x satisfying |z| = 1 and

0 ≤ x < 1, the estimation |γ(z, x)| ≤ 1 holds if and only if the inequality r ≤ t ≤ r + 2 is valid.

Therefore, we can give a sufficient condition of the asymptotic stability of the continuous RK. Hereafter ρ[·] means the spectral radius of the referred square matrix.

Theorem 2.1 The equilibrium solution of the difference equation (2.8) is asymptotically stable if the condition

(C) det[ML̄(ζM̄, ζN)] ≠ 0 and ρ[RL̄(ζM̄, ζN)] < 1

is satisfied for any complex parameter ζ with |ζ| ≤ 1.

Proof. When ζ = 0, we have

ML̄(ζM̄, ζN) = I_{sd} − A ⊗ L̄, RL̄(ζM̄, ζN) = I_{sd} + (b^{T} ⊗ I_{d})ML̄(0, 0)^{−1}(e ⊗ I_{d})L̄.

Thus from the condition (C),

det[Q(z)] = z^{sd} det[I − A ⊗ L̄] det[zI − RL̄(0, 0)] ≠ 0, for any |z| ≥ 1.

Furthermore, we obtain

sup_{|z|=1} ρ[Q(z)^{−1}P (z)] = sup_{|z|=1}{|ζ̂| : det[ζ̂Q(z) − P (z)] = 0}

= ( inf_{|z|=1}{|ζ| : det[Q(z) + ζP (z)] = 0})^{−1} ,

and

det[Q(z) + ζP (z)]

= z^{sd} det[zI_{d} − RL̄(ζ ∑_{p=−r}^{t} L_{p}(δ)z^{p}M̄, ζ ∑_{p=−r}^{t} L_{p}(δ)z^{p}N)]

× det[ML̄(ζ ∑_{p=−r}^{t} L_{p}(δ)z^{p}M̄, ζ ∑_{p=−r}^{t} L_{p}(δ)z^{p}N)]

= z^{sd} det[zI_{d} − RL̄(ζ^{∗}M̄, ζ^{∗}N)] det[ML̄(ζ^{∗}M̄, ζ^{∗}N)],

where ζ^{∗} = ζ ∑_{p=−r}^{t} L_{p}(δ)z^{p}.

Lemma 2.1 yields |ζ^{∗}| ≤ 1 when |z| = 1 and 0 ≤ δ < 1. Then the condition (C) implies that sup_{|z|=1} ρ[Q(z)^{−1}P (z)] < 1. Thus the theorem is obtained by the theorem of in ’t Hout and Spijker [6].

Corollary 2.1 If det[ML̄(M̄, N)] ≠ 0 and ρ[RL̄(M̄, N)] < 1 hold for any L, M and N which satisfy the conditions (S1), (S2) and (S3), then RK method is asymptotically stable.

Proof. Since L, M and N satisfy the conditions (S1), (S2) and (S3) implies L, ζM and ζN satisfy (S1), (S2) and (S3) for any |ζ| ≤ 1.

Finally in this section, we refer to the proposition by [14].

Proposition 2.1 Let ω_{i}(θ), 1 ≤ i ≤ s, be polynomials which define an NCE of the RK method, then the following statements are true:

(i) The matrices W and Γ defined by (2.5) satisfy

W Γ = W, b^{T} Γ = b^{T} .  (2.11)

(ii) The inclusion σ[Γ] ⊂ {0, 1} holds. Here σ[Γ] denotes the spectrum of Γ.

3 Radau IA and Lobatto IIIC types

Let us discuss the Radau IA method applied to NDDEs. Zennaro (Theorems 5 and 7 in [21]) showed that the Radau IA method of order 2s − 1 has a unique NCE, which is called the collocation-based NCE, determined by

ω_{j}(θ) = ∫_{0}^{θ} l_{j}(ϑ)dϑ (1 ≤ j ≤ s),

where l_{j}(θ) are the Lagrange polynomials of degree (s − 1) associated with the RK abscissae {c_{i}}. That is to say,

l_{j}(θ) = ∏_{k=1,≠j}^{s} \frac{θ − c_{k}}{c_{j} − c_{k}} , (1 ≤ j ≤ s).

Due to the above expression the matrix W of the collocation-based NCE of Radau IA has its (i, j)-component given by

ω_{j}(c_{i}) = ∫_{0}^{c_{i}} l_{j}(ϑ)dϑ,  (3.1)

where the abscissa c_{i} is equal to that of Radau IA. We note that the matrix W coincides with the matrix A neither of Radau IA nor of Radau IIA underlying for ODEs. On the other hand, it is obvious

Γ = I.  (3.2)

For the case of the Lobatto IIIC method which is of order (2s − 2), it has at least one NCE called the projection-based NCE. As a matter of fact, the projection-based NCE of Lobatto IIIC employs the polynomial ω_{j}(θ) of degree (s − 1) defined by

ω_{j}(0) = 0, ω_{j}(1) = b_{j}, 1 ≤ j ≤ s,

∫_{0}^{1} θ^{q−1}ω_{j}(θ)dθ = \frac{1}{q}b_{j}(1 − c_{j}^{q}), 1 ≤ j ≤ s, 1 ≤ q ≤ s − 2.

Since the underlying s-stage Lobatto IIIC method satisfies the simplifying condition C(s−1) and ω_{j}(θ) are polynomials of degree (s − 1) fulfilling ω_{j}(0) = 0, we have the identity

AΓ = W.  (3.3)

As the abscissae of the Lobatto methods of all types (IIIA to IIIC) are same, we arrive at Proposition 3.1 The matrix W of the projection-based or the collocation-based NCE of Lobatto IIIC method is equal to the matrix A of the underlying Lobatto IIIA method.

The proof is in [14] (Proposition 2.5).

Next, we turn to the W-transformation. Define the matrix W by

W = (P_{j−1}(c_{i}))_{1≤i,j≤s},

where P_{k}(θ) denotes the normalized shifted Legendre polynomial of order k, i.e.,

P_{k}(θ) = \frac{\sqrt{2k + 1}}{k!} \frac{d^{k}}{dθ^{k}} (θ^{k}(θ − 1)^{k}).

Then by Eq. (3.1) and Proposition 3.1, we have

W^{−1}AW = ( \begin{matrix}\frac{1}{2} & −ξ_{1} \\ ξ_{1} & 0 & ... \\ & ... & ... & −ξ_{s−2} \\ & & ξ_{s−2} & 0 & −η_{s−1} \\ & & & ξ_{s−1} & k_{s} \end{matrix} ) ,  (3.4)

W^{−1}W W = ( \begin{matrix}\frac{1}{2} & −ξ_{1} \\ ξ_{1} & 0 & ... \\ & ... & ... & −ξ_{s−2} \\ & & ξ_{s−2} & 0 & −η̂_{s−1} \\ & & & ξ_{s−1} & k̂_{s} \end{matrix} ) ,  (3.5)

where the off-diagonal elements are ξ_{k} = 1/(2\sqrt{4k^{2} − 1}). However, in the case of the Lobatto IIIC method, we have further

η_{s−1} = \frac{2s − 1}{s − 1}ξ_{s−1} and η̂_{s−1} = 0,

and on the diagonal k_{s} = 1/(2s − 2) and k̂_{s} = 0, whereas, in the case of the Radau IA method

η_{s−1} = η̂_{s−1} = ξ_{s−1}, k_{s} = \frac{1}{4s − 2} and k̂_{s} = −\frac{1}{4s − 2}

(see [5]). Furthermore

b^{T} W = e_{1}^{T} , W^{−1}e = e_{1},  (3.6)

where e_{1} = (1, 0, ⋯, 0)^{T}. These properties of the W-transformation will be applied in our stability analysis.

4 Preliminary lemmata

A series of lemmata concerning with estimations for the matrices M_{X}(Y, Z) and R_{X}(Y, Z) will be shown in this section for our main theorems. Let μ_{2}[·] denote the logarithmic norm of a matrix induced by the Euclidean norm on C I ^{d}. For any d × d matrix X, it can be shown

μ_{2}[X] = max_{v∈C I ^{d},v≠0} \frac{ℜ(v^{∗}Xv)}{v^{∗}v},  (4.1)

where the symbol ^{∗} denotes the conjugate transpose of a vector (see [20]). Moreover, a positive definite matrix X is denoted by X > 0.

The following is Lemma 3.2 of [15].

Lemma 4.1 When H stands for a positive definite Hermitian matrix, the following three statements hold:

(i) If μ_{2}[HX] < 0, then σ[X] ⊂ C I ^{−}.
(ii) If the matrix X is invertible and μ_{2}[HX] ≤ 0 holds, then μ_{2}[HX^{−1}] ≤ 0.

(iii) If the matrix sum X + Y is invertible and X^{∗}H^{−1}X − Y^{∗}H^{−1}Y > 0 holds, then

μ_{2}[H(X + Y)^{−1}(Y − X)] < 0.

Lemma 4.2 Let A, b and W be those of the Radau IA with the collocation-type NCE. Assume that X, X + Y and I − Z are invertible and there is also a Hermitian matrix H > 0 satisfying the conditions μ_{2}[HX(I − Z)^{−1}] ≤ 0, μ_{2}[H(X + Y)(I − Z)^{−1}] ≤ 0 and

X^{∗}H^{−1}X − Y^{∗}H^{−1}Y > 0.  (4.2)

Then the matrix M_{X}(Y, Z) is non-singular, i.e.,

det[M_{X}(Y, Z)] ≠ 0,  (4.3)

the matrix R_{X}(Y, Z) can be expressed with

R_{X}(Y, Z) = I_{d} + 2 [−I_{d} + (I − Z)^{−1}Φ_{0}(X, Y, Z)(I − Z)]^{−1} ,  (4.4)

and its spectral radius is less than unity, i.e.,

ρ[R_{X}(Y, Z)] < 1.  (4.5)

Here the matrices {Φ_{k}(X, Y, Z)} is recursively defined from

Φ_{s}(X, Y, Z) = (I − Z)(Y − X)^{−1}(X + Y)(I − Z)^{−1}  (4.6)

through

Φ_{k}(X, Y, Z) = 2(2k + 1)(I − Z)(X + Y)^{−1} + Φ_{k+1}(X, Y, Z)^{−1}  (4.7)

for k = s − 1, s − 2, …, 0.

Proof. Put X̂ = X(I − Z)^{−1}, Ŷ = Y (I − Z)^{−1}. Due to (3.2), we have

M_{X}(Y, Z) = MX̂(Ŷ , 0)[I ⊗ (I − Z)],  (4.8)

and

R_{X}(Y, Z) = I_{sd} + b^{T} ⊗ I_{d}[I ⊗ (I − Z)]^{−1}[MX̂(Ŷ , 0)]^{−1}(e ⊗ I_{d})(X̂ + Ŷ)(I − Z)
= I_{sd} + b^{T} ⊗ (I − Z)^{−1}[MX̂(Ŷ , 0)]^{−1}(e ⊗ I_{d})(X̂ + Ŷ)(I − Z)
= (I − Z)^{−1}RX̂(Ŷ , 0)(I − Z).  (4.9)

Consider the linear equation

( \begin{matrix}I_{sd} − A ⊗ X − W ⊗ Y − I ⊗ Z & O \\ −b^{T} ⊗ I_{d} & I_{d} \end{matrix} ) ( \begin{matrix}V \\ v & 1 \end{matrix} ) = ( \begin{matrix}e ⊗ (X + Y)v_{0} \\ v & 0 \end{matrix} ) ,  (4.10)

where v_{1} ∈ C I ^{d} and V ∈ C I ^{sd} are unknown variables while v_{0} ∈ C I ^{d} a known variable arbitrarily given. Multiplying both sides of (4.10) by

( \begin{matrix}W^{−1} ⊗ I_{d} & O \\ O & I_{d} \end{matrix} )

and letting

[W^{−1} ⊗ (I − Z)]V = (V̂_{1}^{T} , V̂_{2}^{T} , ⋯, V̂_{s}^{T})^{T} = V̂ , v̂_{1} = (I − Z)v_{1}, v̂_{0} = (I − Z)v_{0} (4.11)

yield the following identities.

[I_{sd} − (W^{−1}AW) ⊗ X̂ − (W^{−1}W W) ⊗ Ŷ]V̂ = W^{−1}e ⊗ (X̂ + Ŷ)v̂_{0}, (4.12)
−(b^{T} ⊗ I_{d})(I ⊗ (I − Z)^{−1})(W ⊗ I_{d})V̂ + (I − Z)^{−1}v̂_{1} = (I − Z)^{−1}v̂_{0}.  (4.13)

Noting the identity (b^{T} ⊗ I_{d})(I ⊗ (I − Z)^{−1}) = (I − Z)^{−1}(b^{T} ⊗ I_{d}) and the equations (3.6), (3.4) and (3.5), we obtain the following decompositions.

[I_{d} − \frac{1}{2}(X̂ + Ŷ)]V̂_{1} + ξ_{1}(X̂ + Ŷ)V̂_{2} = (X̂ + Ŷ)v̂_{0},  (4.14)
−ξ_{k}(X̂ + Ŷ)V̂_{k} + V̂_{k+1} + ξ_{k+1}(X̂ + Ŷ)V̂_{k+2} = 0, (k = 1, …, s − 2), (4.15)
−ξ_{s−1}(X̂ + Ŷ)V̂_{s−1} + (I − \frac{1}{4s − 2}(X̂ − Ŷ)) V̂_{s} = 0,  (4.16)
−V̂_{1} + v̂_{1} = v̂_{0}.  (4.17)

Furthermore, substitution of ψ_{k} = V̂_{k}/\sqrt{2k − 1} (k = 1, …, s) transforms the equations (4.14) – (4.17) into

(−I_{d} + 2(X̂ + Ŷ)^{−1}) ψ_{1} + ψ_{2} = 2v̂_{0},  (4.18)
ψ_{k} = 2(2k + 1)(X̂ + Ŷ)^{−1}ψ_{k+1} + ψ_{k+2}, (k = 1, …, s − 2),  (4.19)
ψ_{s−1} = [2(2s − 1)(X̂ + Ŷ)^{−1} − (X̂ + Ŷ)^{−1}(X̂ − Ŷ)] ψ_{s},  (4.20)
v̂_{1} = v̂_{0} + ψ_{1}.  (4.21)

The second and third statements of Lemma 4.1 together with the equations (4.2) and (4.6) imply

μ_{2}[H(I − Z)^{−1}Φ_{s}(X, Y, Z)(I − Z)] = μ_{2}[H(Y − X)^{−1}(X + Y)] < 0,

which means σ[(I − Z)^{−1}Φ_{s}(X, Y, Z)(I − Z)] ⊂ C I ^{−} through the first statement of Lemma 4.1. By induction and Lemma 4.1, we can prove that σ[(I − Z)^{−1}Φ_{k}(X, Y, Z)(I − Z)] ⊂ C I ^{−} and Φ_{k}(X, Y, Z) is invertible for k = 0, 1, …, s.

Equation (4.20) implies ψ_{s−1} = Φ_{s−1}(X, Y, Z)ψ_{s}. Then by (4.19) we recursively obtain ψ_{k} = Φ_{k}(X, Y, Z)ψ_{k+1} (k = s − 2, …, 1). Henceforth Eqs. (4.18) and (4.21) yield

v̂_{1} = v̂_{0} + 2 (−I_{d} + 2(X̂ + Ŷ)^{−1} + Φ_{1}(X, Y, Z)^{−1})^{−1} v̂_{0},

which brings

v̂_{1} = v̂_{0} + 2 (−I_{d} + Φ_{0}(X, Y, Z))^{−1} v̂_{0},

then

v_{1} = v_{0} + 2 (−I_{d} + (I − Z)^{−1}Φ_{0}(X, Y, Z)(I − Z))^{−1} v_{0},  (4.22)

Since v_{0} = 0 implies that v_{1} as well as V vanishes, Eq. (4.3) holds. From (4.10), we can write v_{1} = R_{X}(Y, Z)v_{0}, which, together with (4.22), yields Eq. (4.4).

Furthermore, since the function 1 + 2/(−1 + z) = (z + 1)/(z − 1) maps the left-half plane ℜz < 0 onto the unit open disk, and the inclusion σ[(I − Z)^{−1}Φ_{0}(X, Y, Z)(I − Z)] ⊂ C I ^{−} holds, we can obtain (4.5) through the spectral mapping theorem.

As for the Lobatto IIIC method with the projection-based NCE, we have the following.

Lemma 4.3 Let A, b and W be those of the Lobatto IIIC method with the projection-based NCE. If X, X + Y and I − Z are invertible, and that there is a Hermitian matrix H > 0 which supports the inequalities μ_{2}[HX] ≤ 0 and μ_{2}[H(X + Y)(I − Z)^{−1}] ≤ 0, then we have

det[M_{X}(Y, Z)] ≠ 0,  (4.23)

and R_{X}(Y, Z) satisfies

R_{X}(Y, Z) = I_{d} + 2 (−I_{d} + (I − Z)^{−1}Φ_{0}(X, Y, Z)(I − Z))^{−1} ,  (4.24)

where the matrix Φ_{0}(X, Y, Z) is recursively defined as

Φ_{k}(X, Y, Z) = 2(2k + 1)(I − Z)(X + Y)^{−1} + Φ_{k+1}(X, Y, Z)^{−1},  (4.25)
for k = s − 2, s − 3, …, 0

with

Φ_{s−1}(X, Y, Z) = (2s − 2)(I − Z)(X + Y)^{−1}X^{−1}(X + Y)(I − Z)^{−1} − I_{d}.  (4.26)

and the inequality

ρ[R_{X}(Y, Z)] < 1  (4.27)

holds.

Proof. By Propositions 2.1 and (3.3), we have AΓ = W and W Γ = W. Thus, putting Ỹ = (XZ + Y)(I − Z)^{−1}, we obtain

M_{X}(Y, Z) = I − A ⊗ X − W ⊗ Y − Γ ⊗ Z
= (I − A ⊗ X − W ⊗ Ỹ)(I − Γ ⊗ Z)
= M_{X}(Ỹ , 0)(I − Γ ⊗ Z).

Since Proposition 2.1 gives b^{T} Γ = b^{T}, we can calculate as follows.

R_{X}(Y, Z) = I_{ds} + (b^{T} ⊗ I_{d})M_{X}(Y, Z)^{−1}(e ⊗ I_{d})(X + Y)
= I_{ds} + (b^{T} ⊗ I_{d})(I − Γ ⊗ Z)^{−1}M_{X}(Ỹ , 0)^{−1}(e ⊗ I_{d})(X + Ỹ)(I − Z)
= (I − Z)^{−1}[b^{T} ⊗ (I_{d} − Γ ⊗ Z)(I − Γ ⊗ Z)^{−1}M_{X}(Ỹ , 0)^{−1}(e ⊗ I_{d})(X + Ỹ)]
×(I − Z)
= (I − Z)^{−1}R_{X}(Ỹ , 0)(I − Z).

Consider the system of linear equations corresponding to (4.10) in this case. Note that the (1,1)-block of the coefficient matrix in the left-hand side is substituted with M_{X}(Y, Z) of this case. By putting

(W^{−1} ⊗ I_{d})(I − Γ ⊗ Z)V = (V̂_{1}^{T} , V̂_{2}^{T} , …, V̂_{s}^{T})^{T} , v̂_{1} = (I − Z)v_{1}, v̂_{0} = (I − Z)v_{0}.

the remaining part of the proof is analogous to that of the previous Lemma and we omit the repeating.

For the Lobatto IIIC method case with the collocation-based NCE we can have the following.

Lemma 4.4 Let A, b and W be those of the Lobatto IIIC method with the collocation-based NCE. If the three matrices X, X + Y and I − Z are invertible, and there is a Hermitian matrix H > 0 which supports the inequalities

μ_{2}[HX(I − Z)^{−1}] ≤ 0 and μ_{2}[H(X + Y)(I − Z)^{−1}] ≤ 0,

then the same results of Lemma 4.3 can be obtained for the matrices M_{X}(Y, Z) and R_{X}(X, Y), except that the expression for Φ_{s−1}(X, Y, Z) is replaced with

Φ_{s−1}(X, Y, Z) = (2s − 2)(I − Z)(X + Y)^{−1}(I − Z)X^{−1}(X + Y)(I − Z)^{−1} − I_{d}.

Taking into account that the matrix W coincides with the matrix A of the underlying Lobatto IIIA method, we can prove this lemma similarly to the above lemmata.

5 Main theorems and illustrative examples

Prior to the statement of our main theorems, we give the target concept of numerical stability precisely.

Definition 5.1 A numerical method is said to be NGP-stable if the numerical solution of the method when it is applied to any analytically stable system (1.2) diminishes to zero as n → ∞ with arbitrary stepsize.

Here the symbol NGP stands for the general P-stability for NDDEs.

Lemmata 4.2, 4.3 and 4.4 obviously imply the following main theorems, respectively.

Theorem 5.1 Assume that the conditions (S1), (S2) and (S3) hold and, moreover, that

(LCP) there is a d × d Hermitian matrix H > 0 which supports

μ_{2}[HL] ≤ 0 and μ_{2}[H(L + M)(I − N)^{−1})] ≤ 0.  (5.1)

Then the Lobatto IIIC method with the projection-based NCE is NGP-stable.

Theorem 5.2 Assume that the conditions (S1), (S2) and (S3) hold and, moreover, that

(LCI) there is a d × d Hermitian matrix H > 0 which supports

μ_{2}[HL(I − N)^{−1}] ≤ 0 and μ_{2}[H(L + M)(I − N)^{−1})] ≤ 0.  (5.2)

Then the Lobatto IIIC method with the collocation-based NCE is NGP-stable.

Theorem 5.3 Assume that the conditions (S1), (S2), (S3) and (LCI) are satisfied, and, furthermore, that the inequality

L^{∗}H^{−1}L − M^{∗}H^{−1}M > 0,

holds. Then the Radau IA method with the NCE is NGP-stable.

To illustrate our main theorems, we give a few simple examples of NDDE. Recall Eq. (1.2) given by

u′(t) = Lu(t) + Mu(t − τ) + Nu′(t − τ), t ≥ 0.

The first example is the case for

L = ( \begin{matrix}−10 & 0 \\ 0 & −10 \end{matrix} ) , M = ( \begin{matrix}−1 & 0 \\ −1 & 0 \end{matrix} ) , N = ( \begin{matrix}0 & 0 \\ \frac{1}{3} & 0 \end{matrix} ) .

It is easy to calculate that the spectrum of L is just −10, while that of L^{−1}M {0, 1/10}. Furthermore, we have

sup_{ℜz=0} ρ ((zI − L)^{−1}(M + zN)) = sup_{ℜz=0} |\frac{1}{10 + z} | < 1

Henceforth the conditions (S1), (S2) and (S3) are satisfied, and the system (1.2) is analytically asymptotically stable.

When we take the positive definite matrix H = ( \begin{matrix}3 & 0 \\ 0 & 1 \end{matrix} ), we can show

μ_{2}[HL] = −10,

μ_{2}[HL(I − N)^{−1}] = \frac{−60 + 5\sqrt{37}}{3} ≤ 0,

μ_{2}[H(L + M)(I − N)^{−1})] = \frac{−774 + 6\sqrt{4930}}{3}6 ≤ 0,

and L^{∗}H^{−1}L − M^{∗}H^{−1}M = ( \begin{matrix}32 & 0 \\ 0 & 100 \end{matrix} ) is positive definite.

Then by Theorems 5.1, 5.2 and 5.3, the numerical solutions of the Radau IA method, of the Lobatto IIIC method with the collocation-based NCE and of the Lobatto IIIC with the projection-based NCE are all NGP-stable.

The second example is the case of

L = ( \begin{matrix}−10 & 1 & 2 \\ 3 & −12 & 0 \\ 1 & 2 & −9 \end{matrix} ) , M = ( \begin{matrix}−1 & 0 & 3 \\ 1/2 & 1/2 & 1 \\ 1/2 & 3/2 & 0 \end{matrix} ) , N = \frac{1}{72} ( \begin{matrix}−1 & 5 & 2 \\ 4 & 0 & 3 \\ −2 & 4 & 1 \end{matrix} ) .

It has been shown to be analytically asymptotically stable in [8]. And we can obtain the estimation ∥ N ∥_{2}≈ 0.0982387 < 1, which implies (S1), (S2) and (S3).

When we take the identity matrix as H, we can obtain

μ_{2}[HL] ≈ −7.0761,
μ_{2}[HL(I − N)^{−1}] ≈ −7.46325,
μ_{2}[H(L + M)(I − N)^{−1})] ≈ −5.3965,

and L^{∗}H^{−1}L − M^{∗}H^{−1}M = ( \begin{matrix}217/2 & −45 & −53/2 \\ −45 & 293/2 & −33/2 \\ −53/2 & −33/2 & 75 \end{matrix} ) is positive definite. Hence the

numerical solutions by the Runge-Kutta methods considered in this paper are all NGPstable.

References

[1] A. Bellen, Z. Jackiewicz and M. Zennaro, Stability analysis of one-step methods for neutral delay-differential equations, Numer. Math. , 52(1988), 605-619.

[2] R. K. Brayton and R. A. Willoughby, On the numerical integration of a symmetric system of a difference-differential equations, J. Math. Anal. Appl. , 18(1967), 182-189.

[3] J. C. Butcher, The Numerical Analysis of Ordinary Differential Equations. Runge-Kutta and General Linear Methods, John Wiley & Sons, Chichester, 1987.

[4] K. Dekker and J. G. Verwer, Stability of Runge-Kutta Methods for Stiff Nonlinear Differential Equations. North-Holland, Amsterdam, 1984.

[5] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, Springer-Verlag, Berlin, Heidelberg, New York, 1991.

[6] K. J. in ’t Hout and M. N. Spijker, Stability analysis of numerical methods for delay differential equations, Numer. Math. , 59(1991), 807-814.

[7] K. J. in ’t Hout, A new interpolation procedure for adapting Runge-Kutta methods for delay differential equations, BIT, 32(1992), 634-649.

[8] G.-D. Hu, Stability and Numerical Analysis for Delay Differential Systems, Ph. D. thesis granted by Graduate School of Human Informatics, Nagoya University, Oct., 1996

[9] G. D. Hu and T. Mitsui, Stability of numerical methods for systems of neutral delay differential equations, BIT, 35(1995), 504-515.

[10] Z. Jackiewicz, One-step methods of any order for neutral functional -differential equations, SIAM J. Numer. Anal. , 21(1984), 486-511.

[11] Z. Jackiewicz, Quasilinear multistep methods and variable-step predictor-corrector methods for neutral functional differential equations. SIAM J. Numer. Anal. , 23(1986), 423-452.

[12] Z. Kamont and M. Kwapisz, On the Cauchy problem for differential delay equations in a Banach space, Math. Nachr. , 74(1976), 173-190.

[13] T. Koto, A stability property of A-stable collocation-based Runge-Kutta methods for neutral delay differential equations, BIT, 36(1996), 855-859.

[14] T. Koto, NP-stability of Runge-Kutta methods based on classical quadrature, BIT, 37(1997), 870-884.

[15] T. Koto, Stability of the Radau IA and Lobatto IIIC methods for delay differential equations, Numer. Math. , 79(1998), 569-580.

[16] J. X. Kuang, J. X. Xiang and H. J. Tian, The asymptotic stability of one-parameter methods for neutral differential equations, BIT, 34(1994), 400-408.

[17] Y. K. Liu, Stability analysis of θ-methods for neutral functional differential equations, Numer. Math. , 70(1995), 473-483.

[18] L. Qiu, B. Yang and J. X. Kuang, The NGP-stability of Runge-Kutta methods for systems of neutral delay differential equations, Numer. Math. , 81(1999), 451-459.

[19] G. Strang, Trigonometric polynomials and difference methods of maximum accuracy, J. Math. Phys., 41(1962), 147-154.

[20] T. Ström, On logarithmic norms, SIAM J. Numer. Anal. , 12(1975), 741-753.

[21] M. Zennaro, Natural continuous extensions of Runge-Kutta methods, Math. Comput. , 46, (1986), 119-133.

[22] M. Zennaro, Natural Runge-Kutta and projection methods, Numer. Math. , 53, (1988), 423-438.

[23] M. Zennaro, Delay differential equations: Theory and numerics, in M. Ainsworth et al. (ed.), Theory and Numerics of Ordinary and Partial Differential Equations, 291-333, Oxford University Press, Oxford, 1995.