# Stability of Runge-Kutta Methods for Neutral Delay Differential Equations

Authors: Liping Wen; Xiong Liu; Yuexin Yu
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

This paper is concerned with the numerical stability of a class of nonlinear neutral delay differential equations. The numerical stability results are obtained for(
k
,
l
)-algebraically stable Runge-Kutta methods when they are applied to this type of problem. Numerical examples are given to confirm our theoretical results.


---

## Body

## 1. Introduction

Let〈
·
,
·
〉 be an inner product on C
N and let · be the corresponding norm. Consider the initial value problems (IVPs) of nonlinear neutral delay differential equations (NDDEs) as follows (cf. [1]): (1)
y
′
t
=
f
t
,
y
t
,
y
t
-
τ
1
,
y
′
t
-
τ
2
,
t
≥
t
0
,
y
t
=
φ
t
,
t
0
-
τ
≤
t
≤
t
0
,where τ
1, τ
2 are real positive constants, τ
=
max
⁡
{
τ
1
,
τ
2
}, φ
:
[
t
0
-
τ
,
t
0
]
→
C
N is a continuously differentiable function, and f
:
[
t
0
,
+
∞
)
×
C
N
×
C
N
×
C
N
→
C
N is a continuous function satisfying(2)
R
e
u
1
-
u
2
,
f
t
,
u
1
,
v
,
w
-
f
t
,
u
2
,
v
,
w
≤
α
u
1
-
u
2
2
,
∀
t
≥
t
0
,
u
1
,
u
2
,
v
,
w
∈
C
N
,
(3)
f
t
,
u
,
v
1
,
w
-
f
t
,
u
,
v
2
,
w
≤
β
1
v
1
-
v
2
,
∀
t
≥
t
0
,
u
,
v
1
,
v
2
,
w
∈
C
N
,
(4)
f
t
,
u
,
v
,
w
1
-
f
t
,
u
,
v
,
w
2
≤
β
2
w
1
-
w
2
,
∀
t
≥
t
0
,
u
,
v
,
w
1
,
w
2
∈
C
N
,
(5)
F
t
,
u
,
v
1
,
x
1
,
y
,
z
-
F
t
,
u
,
v
2
,
x
2
,
y
,
z
≤
γ
1
v
1
-
v
2
+
γ
2
x
1
-
x
2
,
∀
t
≥
t
0
,
u
,
v
1
,
v
2
,
x
1
,
x
2
,
y
,
z
∈
C
N
,where (6)
F
t
,
u
,
v
,
x
,
y
,
z
≔
f
t
,
u
,
v
,
f
t
-
τ
2
,
x
,
y
,
zand α, β
1, β
2, γ
1, and γ
2 are real constants and β
1
≥
0, 0
≤
β
2
<
1, γ
1
>
0, and γ
2
>
0.The existence and uniqueness of solution to (1) are discussed in [2, 3]. Throughout this paper we assume that problem (1) has a unique true solution y
(
t
), and we will still use the symbol D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
) presented in literature [1] to denote the problem class consisting of all problem (1) satisfying conditions (2)–(5).In order to investigate the stability of NDDEs, we have to introduce the perturbed problem of (1), which is defined by the same function f
(
t
,
u
,
v
,
w
), but with another initial condition:(7)
z
′
t
=
f
t
,
z
t
,
z
t
-
τ
1
,
z
′
t
-
τ
2
,
t
≥
t
0
,
z
t
=
ψ
t
,
t
0
-
τ
≤
t
≤
t
0
,where ψ
:
[
t
0
-
τ
,
t
0
]
→
C
N is a continuous and differentiable function. The unique exact solution of problem (7) is denoted by z
(
t
).Remark 1.
In 2004, Wang and Li [4] considered nonlinear NDDEs of form as(8)
y
′
t
=
f
t
,
y
t
,
y
t
-
τ
,
y
′
t
-
τ
,
t
≥
0
,
y
t
=
φ
t
,
-
τ
≤
t
≤
0
,and they obtained a series of stability results of theoretical solution and numerical solution which was given by backward Euler methods. Recently, they further investigated the stability of one-leg methods [5], Runge-Kutta methods [6], and continuous Runge-Kutta-type methods [7] for the solution to problem (8), respectively. More research in this field can be found in literatures [2, 3, 8–15]. The research of numerical methods for more extensive problem class can be seen in [16–19].The difference between (1) and (8) is that the quantities τ
1, τ
2 in (1) can be unequal. However, many real-world phenomena can be described by (1) but not with (8). For the application of this type in the real world, one is referred to [1], because, in these problems, the change of the state variable depends on the state of the past some time and, in addition, a state change of another some time (cf. [20, 21]). It is easy to see that the existent analytic and numerical stability results for (8) in the abovementioned literatures can not be applied to problem (1).Remark 2.
In 2012, the authors of the present paper considered problem (1), and a sufficient condition for the stability of the problem itself is given [1]. This result is described below.Theorem 3 (cf. [1]).
Assume that (1) belongs to D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
) with α
<
0 and (9)
γ
1
+
γ
2
-
β
2
α
-
α
≤
1
.Then one has (10)
y
t
-
z
t
≤
ϕ
1
+
β
1
ϕ
1
+
β
2
ϕ
2
-
α
∑
i
=
0
τ
2
/
τ
1
1
+
β
1
-
α
i
,where (11)
ϕ
1
=
s
u
p
t
0
-
τ
≤
s
≤
t
0
φ
s
-
ψ
s
,
ϕ
2
=
s
u
p
t
0
-
τ
≤
s
≤
t
0
φ
′
s
-
ψ
′
s
,and the notation x denotes the largest integer smaller than or equal to x.Remark 4.
In [1], the numerical stability results are obtained for A-stable one-leg methods when they are applied to problem (1). This paper pursues this and further investigates the stability of Runge-Kutta methods for problem (1). It is well known that the Runge-Kutta-type methods are a class of important and common numerical methods for solving differential equations. Therefore, it is important to analyse whether or not Runge-Kutta methods inherit the stability of the underlying problem when they are applied to (1). This also is the motivation of this paper.The organization of this paper is as follows. In Section2, the numerical stability results of (
k
,
l
)-algebraically stable Runge-Kutta methods for (1) are given. Finally, in Section 3, some numerical experiments are given which confirm the theoretical results obtained in this paper.


## 2. The Stability of Runge-Kutta Methods

The adaptation of thes stage Runge-Kutta method(12)
A
,
b
,
c
=
c
1
a
11
a
12
⋯
a
1
s
c
2
a
21
a
22
⋯
a
2
s
⋮
⋯
⋯
⋯
⋯
c
s
a
s
1
a
s
2
⋯
a
s
s
b
1
b
2
⋯
b
sfor ODEs can generally lead to an s stage Runge-Kutta method(13)
Y
i
n
=
y
n
+
h
∑
j
=
1
s
a
i
j
f
t
n
+
c
j
h
,
Y
j
n
,
Y
¯
j
n
,
Y
~
j
n
,
i
=
1,2
,
…
,
s
,
y
n
+
1
=
y
n
+
h
∑
j
=
1
s
b
j
f
t
n
+
c
j
h
,
Y
j
n
,
Y
¯
j
n
,
Y
~
j
n
,
n
=
0,1
,
…for solving problem (1) in NDDEs; here, h
>
0 is the step size and t
n
=
t
0
+
n
h. The values y
n, Y
i
(
n
), Y
¯
i
(
n
), and Y
~
i
(
n
) are approximations to y
(
t
n
), y
(
t
n
+
c
i
h
), y
(
t
n
+
c
i
h
-
τ
1
), and y
′
(
t
n
+
c
i
h
-
τ
2
), respectively. y
0
=
φ
(
t
0
). The symbols A and b are linear mappings corresponding to the matrices A
=
[
a
i
j
]
∈
R
s
×
s and b
T
=
[
b
1
,
b
2
,
…
,
b
s
]
∈
R
s, respectively. We always assume that ∑
i
=
1
s
b
i
=
1; each c
i
=
∑
j
=
1
s
a
i
j
∈
[
0,1
].In this paper, we always leth
=
τ
2
/
m
2 for some positive integer m
2 and τ
1
=
(
m
1
-
δ
)
h with integer m
1
≥
1 and δ
∈
[
0,1
). Let Y
j
(
i
)
=
φ
(
t
i
+
c
j
h
) for i
<
0, 1
≤
j
≤
s. Define (cf. [22]) (14)
Y
¯
j
n
=
∑
i
=
-
μ
ν
L
i
δ
Y
j
n
-
m
1
+
i
,
t
n
+
c
j
h
-
τ
1
>
t
0
,
φ
t
n
+
c
j
h
-
τ
1
,
t
n
+
c
j
h
-
τ
1
≤
t
0
,
j
=
1,2
,
…
,
s
,where (15)
L
i
x
=
∏
k
=
-
μ
,
k
≠
i
ν
x
-
k
i
-
k
,
x
∈
0,1
,and we assume integers μ
,
ν
≥
0 and m
1
≥
ν
+
1 so as to guarantee that, in the interpolation procedure for Y
¯
j
(
n
), no unknown values Y
j
(
k
) with k
≥
n are used.The argumentsY
~
i
(
n
) are given by(16)
Y
~
i
n
=
f
t
n
+
c
i
h
-
τ
2
,
Y
i
n
-
m
2
,
Y
¯
i
n
-
m
2
,
Y
~
i
n
-
m
2
,
t
n
+
c
i
h
-
τ
2
>
t
0
,
φ
′
t
n
+
c
i
h
-
τ
2
,
t
n
+
c
i
h
-
τ
2
≤
t
0
,
i
=
1,2
,
…
,
s
.Similarly, the adaptation of the Runge-Kutta method(
A
,
b
,
c
) with the same interpolation procedures to problem (7) leads to the following process: (17)
Z
i
n
=
z
n
+
h
∑
j
=
1
s
a
i
j
f
t
n
+
c
j
h
,
Z
j
n
,
Z
¯
j
n
,
Z
~
j
n
,
i
=
1,2
,
…
,
s
,
z
n
+
1
=
z
n
+
h
∑
j
=
1
s
b
j
f
t
n
+
c
j
h
,
Z
j
n
,
Z
¯
j
n
,
Z
~
j
n
,
n
=
0,1
,
…
,where z
0
=
ψ
(
t
0
).Definition 5 (cf. [23]).
Letk, l be real constants. A method (
A
,
b
,
c
) is said to be (
k
,
l
)-algebraically stable if there exists a diagonal nonnegative matrix D
=
diag
⁡
(
d
1
,
d
2
,
…
,
d
s
) such that M
=
[
M
i
j
] is nonnegative definite, where (18)
M
=
k
-
1
-
2
l
e
T
D
e
e
T
D
-
b
T
-
2
l
e
T
D
A
D
e
-
b
-
2
l
A
T
D
e
D
A
+
A
T
D
-
b
b
T
-
2
l
A
T
D
A
.As an important special case, a (
1,0
)-algebraically stable method is called algebraically stable for short.In the following, we give the main result of this paper and its proof.Theorem 6.
Assume that method (12) is (
k
,
l
)-algebraically stable with 0
<
k
≤
1 and that problems (1) and (7) belong to D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
) with α
<
0, β
1
≤
γ
1. Then when (19)
2
α
+
γ
1
+
L
0
2
μ
+
ν
+
1
2
γ
1
+
2
γ
2
1
-
β
2
h
≤
2
l
,method (13) with (14) and (16) is stable; that is, (20)
y
n
-
z
n
≤
C
max
⁡
m
a
x
t
0
-
τ
≤
t
≤
t
0
φ
t
-
ψ
t
,
m
a
x
t
0
-
τ
≤
t
≤
t
0
φ
′
t
-
ψ
′
t
,where the constant C
>
0 depends only on α, β
1, β
2, γ
1, γ
2, τ
1, and τ
2 and the sequences {
y
n
} and {
z
n
} are two approximation solutions produced by (13) and (17), respectively.Proof.
For anyn
>
0, let q
i
(
n
)
=
(
n
+
c
i
)
/
m
2, i
=
1,2
,
…
,
s; then we have (21)
t
0
+
q
i
n
τ
2
≤
t
n
+
c
i
h
≤
t
0
+
q
i
n
+
1
τ
2
.
Let(22)
w
n
=
y
n
-
z
n
,
W
i
n
=
Y
i
n
-
Z
i
n
,
W
¯
i
n
=
Y
¯
i
n
-
Z
¯
i
n
,
Q
i
n
=
f
t
n
+
c
i
h
,
Y
i
n
,
Y
¯
i
n
,
Y
~
i
n
-
f
t
n
+
c
i
h
,
Z
i
n
,
Z
¯
i
n
,
Z
~
i
n
.Then from (13) and (17) we obtain that (23)
W
i
n
=
w
n
+
h
∑
j
=
1
s
a
i
j
Q
j
n
,
i
=
1,2
,
…
,
s
,
w
n
+
1
=
w
n
+
h
∑
j
=
1
s
b
j
Q
j
n
.As in [23], we can easily obtain that(24)
w
n
+
1
2
-
k
w
n
2
-
2
∑
j
=
1
s
d
j
R
e
W
j
n
,
h
Q
j
n
-
l
W
j
n
=
-
∑
i
=
1
s
+
1
∑
j
=
1
s
+
1
M
i
j
u
i
,
u
j
,where (25)
u
1
=
w
n
,
u
i
=
h
Q
i
-
1
n
,
i
=
2,3
,
…
,
s
+
1
.Noting that the method is (
k
,
l
)-algebraically stable; then (24) reads (26)
w
n
+
1
2
≤
k
w
n
2
+
2
∑
j
=
1
s
d
j
R
e
W
j
n
,
h
Q
j
n
-
l
W
j
n
.It follows from (2) that(27)
2
R
e
W
i
n
,
h
Q
i
n
=
2
h
R
e
W
i
n
,
f
t
n
+
c
i
h
,
Y
i
n
,
Y
¯
i
n
,
Y
~
i
n
-
f
t
n
+
c
i
h
,
Z
i
n
,
Z
¯
i
n
,
Z
~
i
n
≤
2
h
α
W
i
n
2
+
2
h
W
i
n
Γ
i
n
,where (28)
Γ
i
n
=
f
t
n
+
c
i
h
,
Z
i
n
,
Y
¯
i
n
,
Y
~
i
n
-
f
t
n
+
c
i
h
,
Z
i
n
,
Z
¯
i
n
,
Z
~
i
n
.When t
n
+
c
i
h
≤
t
0
+
τ
2, from (3), (4), and (16), we have (29)
Γ
i
n
≤
β
1
W
¯
i
n
+
β
2
ϕ
2
.Otherwise, when t
n
+
c
i
h
>
t
0
+
τ
2, from (4), (5), and (16), we have (30)
Γ
i
n
=
f
t
n
+
c
i
h
,
Z
i
n
,
Y
¯
i
n
,
f
t
n
+
c
i
h
-
τ
2
,
Y
i
n
-
m
2
,
Y
¯
i
n
-
m
2
,
Y
~
i
n
-
m
2
-
f
t
n
+
c
i
h
,
Z
i
n
,
Z
¯
i
n
,
f
t
n
+
c
i
h
-
τ
2
,
Z
i
n
-
m
2
,
Y
¯
i
n
-
m
2
,
Y
~
i
n
-
m
2
+
f
t
n
+
c
i
h
,
Z
i
n
,
Z
¯
i
n
,
f
t
n
+
c
i
h
-
τ
2
,
Z
i
n
-
m
2
,
Y
¯
i
n
-
m
2
,
Y
~
i
n
-
m
2
-
f
t
n
+
c
i
h
,
Z
i
n
,
Z
¯
i
n
,
f
t
n
+
c
i
h
-
τ
2
,
Z
i
n
-
m
2
,
Z
¯
i
n
-
m
2
,
Z
~
i
n
-
m
2
≤
γ
1
W
¯
i
n
+
γ
2
W
i
n
-
m
2
+
β
2
Γ
i
n
-
m
2
.From this and by induction we have (31)
Γ
i
n
≤
γ
1
∑
j
=
1
q
i
n
β
2
j
-
1
W
¯
i
n
-
j
-
1
m
2
+
γ
2
∑
j
=
1
q
i
n
β
2
j
-
1
W
i
n
-
j
m
2
+
β
2
q
i
n
Γ
i
n
-
q
i
n
m
2
≤
γ
1
∑
j
=
1
q
i
n
+
1
β
2
j
-
1
W
¯
i
n
-
j
-
1
m
2
+
γ
2
∑
j
=
1
q
i
n
β
2
j
-
1
W
i
n
-
j
m
2
+
β
2
q
i
n
β
2
ϕ
2
,where we have used (29) and β
1
≤
γ
1. Apparently, inequality (29) can be integrated to (31) by q
i
(
n
)
=
0. Note that, here and later, when q
<
p, the value of sum ∑
p
q
=
0.
Inequality (27) together with (31) yields (32)
2
R
e
W
i
n
,
h
Q
i
n
≤
h
2
α
+
γ
1
+
γ
2
1
-
β
2
W
i
n
2
+
h
γ
1
∑
j
=
1
q
i
n
+
1
β
2
j
-
1
W
¯
i
n
-
j
-
1
m
2
2
+
h
γ
2
∑
j
=
1
q
i
n
β
2
j
-
1
W
i
n
-
j
m
2
2
+
h
β
2
q
i
n
β
2
2
γ
2
ϕ
2
2
.Denote (33)
μ
1
=
h
2
α
+
γ
1
+
γ
2
1
-
β
2
-
2
l
.Then, from (26) and noting that 0
<
k
≤
1, we can easily obtain (34)
w
n
+
1
2
≤
k
w
n
2
+
μ
1
∑
i
=
1
s
d
i
W
i
n
2
+
h
∑
i
=
1
s
d
i
γ
1
∑
j
=
1
q
i
n
+
1
β
2
j
-
1
W
¯
i
n
-
j
-
1
m
2
2
+
γ
2
∑
j
=
1
q
i
n
β
2
j
-
1
W
i
n
-
j
m
2
2
+
β
2
q
i
n
β
2
2
γ
2
ϕ
2
2
≤
w
0
2
+
μ
1
∑
j
=
0
n
∑
i
=
1
s
d
i
W
i
j
2
+
h
γ
1
∑
μ
=
0
n
∑
i
=
1
s
d
i
∑
j
=
1
q
i
μ
+
1
β
2
j
-
1
W
¯
i
μ
-
j
-
1
m
2
2
+
h
γ
2
∑
μ
=
0
n
∑
i
=
1
s
d
i
∑
j
=
1
q
i
μ
β
2
j
-
1
W
i
μ
-
j
m
2
2
+
h
∑
μ
=
0
n
∑
i
=
1
s
d
i
β
2
q
i
μ
β
2
2
γ
2
ϕ
2
2
≤
w
0
2
+
μ
1
∑
i
=
1
s
d
i
∑
j
=
0
n
W
i
j
2
+
h
γ
2
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
n
W
i
j
2
+
h
γ
1
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
n
W
¯
i
j
2
+
h
d
m
2
1
-
β
2
β
2
2
γ
2
ϕ
2
2
≤
w
0
2
+
μ
1
+
h
γ
2
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
n
W
i
j
2
+
h
γ
1
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
n
W
¯
i
j
2
+
d
τ
2
1
-
β
2
β
2
2
γ
2
ϕ
2
2
,where d
=
∑
i
=
1
s
d
i.
In order to estimate the termW
¯
i
(
j
) in (34), we denote J
0
=
{
1,2
,
…
,
s
}, J
1
=
{
j
∈
J
0
∣
c
j
+
δ
≤
1
}. Then when n
≥
m
1 or n
=
m
1
-
1 but j
∈
J
0
-
J
1, from (14), we have (35)
W
¯
j
n
=
Y
¯
j
n
-
Z
¯
j
n
=
∑
i
=
-
μ
ν
L
i
δ
Y
j
n
-
m
1
+
i
-
Z
j
n
-
m
1
+
i
≤
L
0
∑
i
=
-
μ
ν
W
j
n
-
m
1
+
i
,where (36)
L
0
=
m
a
x
-
μ
≤
i
≤
ν
s
u
p
x
∈
0,1
L
i
x
.Therefore, we obtain that (37)
W
¯
j
n
2
≤
L
0
2
μ
+
ν
+
1
∑
i
=
-
μ
ν
W
j
n
-
m
1
+
i
2
.
Otherwise, whenn
≤
m
1
-
2 or n
=
m
1
-
1 but j
∈
J
1, we have (38)
W
¯
j
n
≤
ϕ
1
.
Therefore, whenn
≤
m
1
-
2, (34) with (38) gives (39)
w
n
+
1
2
≤
w
0
2
+
μ
1
+
h
γ
2
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
n
W
i
j
2
+
h
γ
1
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
n
ϕ
1
2
+
d
τ
2
1
-
β
2
β
2
2
γ
2
ϕ
2
2
≤
w
0
2
+
μ
1
+
h
γ
2
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
n
W
i
j
2
+
d
τ
1
γ
1
1
-
β
2
ϕ
1
2
+
d
τ
2
1
-
β
2
β
2
2
γ
2
ϕ
2
2
.
Whenn
=
m
1
-
1, (34) gives (40)
w
m
1
2
≤
w
0
2
+
μ
1
+
h
γ
2
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
m
1
-
1
W
i
j
2
+
h
γ
1
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
m
1
-
1
W
¯
i
j
2
+
d
τ
2
1
-
β
2
β
2
2
γ
2
ϕ
2
2
≤
w
0
2
+
μ
1
+
h
γ
2
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
m
1
-
1
W
i
j
2
+
2
d
τ
1
γ
1
1
-
β
2
ϕ
1
2
+
d
τ
2
1
-
β
2
β
2
2
γ
2
ϕ
2
2
+
h
γ
1
1
-
β
2
∑
i
∈
J
0
-
J
1
d
i
W
¯
i
m
1
-
1
2
≤
w
0
2
+
μ
1
+
h
γ
2
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
m
1
-
1
W
i
j
2
+
2
d
τ
1
γ
1
1
-
β
2
ϕ
1
2
+
d
τ
2
1
-
β
2
β
2
2
γ
2
ϕ
2
2
+
L
0
2
μ
+
ν
+
1
h
γ
1
1
-
β
2
∑
i
∈
J
0
-
J
1
d
i
∑
j
=
-
μ
ν
W
i
-
1
+
j
2
≤
w
0
2
+
μ
1
+
h
γ
2
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
m
1
-
1
W
i
j
2
+
2
d
τ
1
γ
1
1
-
β
2
ϕ
1
2
+
d
τ
2
1
-
β
2
β
2
2
γ
2
ϕ
2
2
+
L
0
2
μ
+
ν
+
1
h
γ
1
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
-
μ
ν
W
i
-
1
+
j
2
.
Whenn
≥
m
1, from (34) and (37) we have (41)
w
n
+
1
2
≤
w
0
2
+
μ
1
+
h
γ
2
1
-
β
2
∑
i
=
1
s
d
i
∑
j
=
0
n
W
i
j
2
+
2
d
τ
1
γ
1
1
-
β
2
ϕ
1
2
+
d
τ
2
1
-
β
2
β
2
2
γ
2
ϕ
2
2
+
L
0
2
μ
+
ν
+
1
h
γ
1
1
-
β
2
∑
k
=
m
1
-
1
n
∑
i
=
1
s
d
i
∑
j
=
-
μ
ν
W
i
k
-
m
1
+
j
2
≤
w
0
2
+
h
2
α
+
γ
1
+
L
0
2
μ
+
ν
+
1
2
γ
1
+
2
γ
2
1
-
β
2
-
2
l
∑
i
=
1
s
d
i
∑
j
=
0
n
W
i
j
2
+
2
+
L
0
2
μ
+
ν
+
1
μ
+
1
μ
+
2
2
d
τ
1
γ
1
1
-
β
2
ϕ
1
2
+
d
τ
2
1
-
β
2
β
2
2
γ
2
ϕ
2
2
.By condition (19), a combination of (39), (40), and (41) shows the method is stable, which completes the proof of Theorem 6.For the case where the Runge-Kutta method is algebraically stable, from Theorem6 we can obtain following result.Corollary 7.
Assume that method (12) is algebraically stable and that problems (1) and (7) belong to D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
) with α
<
0, β
1
≤
γ
1. Then when (42)
2
α
+
γ
1
1
+
L
0
2
μ
+
ν
+
1
2
+
2
γ
2
1
-
β
2
≤
0
,method (13) with (14) and (16) is unconditionally stable for step size h.Remark 8.
It is well known that the Gauss, Radau IA, Radau IIA, and Lobatto IIIC Runge-Kutta methods are all algebraically stable; therefore, from Corollary7, these methods with (14) and (16) are stable when applied to (1) ∈
D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
) with α
<
0, β
1
≤
γ
1 provided (42) satisfied. Paper [24] also shows that the Gauss and Radau methods have some good stability properties for delay equations. This is in harmony with the results in this paper.Remark 9.
Whenτ
1
=
τ
2 in right-hand side function of problem (1), from Theorem 6 we can obtain the corresponding stability result which is similar to that obtained by Wang in [6], but the result of Theorem 6 is more extensive than that in [6], because, in [6], for processing the delay items, only piecewise linear interpolation is considered, which will lead to higher order method appear order reduction phenomenon for the problem which possesses sufficiently smooth solution.Remark 10.
For the method which has order at most 2, we can chooseν
=
1, μ
=
0 in (15) to process Y
¯
j
(
n
). Here (14) can be reduced to (43)
Y
¯
j
n
=
1
-
δ
Y
j
n
-
m
1
+
δ
Y
j
n
-
m
1
+
1
,
t
n
+
c
j
h
-
τ
1
>
t
0
,
φ
t
n
+
c
j
h
-
τ
1
,
t
n
+
c
j
h
-
τ
1
≤
t
0
,
j
=
1,2
,
…
,
s
,which is widely used by many authors for numerically solving differential equations with delay term, such as [6, 8, 25]. Theorem 6 and Corollary 7 are suitable for this situation certainly. Notable is, in this case, the conclusion of the theorem can be more perfect, if we deduce it directly based on (43).


## 3. Numerical Experiment

We consider the following nonlinear NDDE:(44)
y
1
′
t
=
-
6
y
1
t
+
0.1
sin
⁡
y
2
t
-
1
sin
⁡
y
1
′
t
-
0.33
+
ϕ
1
t
,
t
≥
0
,
y
2
′
t
=
-
4
y
2
t
-
0.2
cos
⁡
y
1
t
-
1
cos
⁡
y
2
′
t
-
0.33
+
ϕ
2
t
,
t
≥
0
,
y
1
t
=
sin
⁡
5
t
,
y
2
t
=
cos
⁡
5
t
,
-
1
≤
t
≤
0
,where (45)
ϕ
1
t
=
5
cos
⁡
5
t
+
6
sin
⁡
5
t
-
0.1
sin
⁡
cos
⁡
5
t
-
1
sin
⁡
5
cos
⁡
5
t
-
0.33
,
ϕ
2
t
=
-
5
sin
⁡
5
t
+
4
cos
⁡
5
t
+
0.2
cos
⁡
sin
⁡
5
t
-
1
cos
⁡
5
sin
⁡
5
t
-
0.33
.Letα
=
-
4, β
1
=
0.2, β
2
=
0.2, γ
1
=
0.2, and γ
2
=
1.2; then problem (44) belongs to D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
) and condition (9) is satisfied, where the inner product is standard inner product. Theorem 3 implies that system (44) is stable.In order to solve problem (44), we consider the 2-stage Radau IIA Runge-Kutta method: (46)
Y
1
n
=
y
n
+
5
12
h
f
t
n
+
1
3
h
,
Y
1
n
-
1
12
h
f
t
n
+
h
,
Y
2
n
,
Y
2
n
=
y
n
+
3
4
h
f
t
n
+
1
3
h
,
Y
1
n
+
1
4
h
f
t
n
+
h
,
Y
2
n
,
y
n
+
1
=
Y
2
n
.Because method (46) is A-stable, then from Theorem 6 the numerical solutions will preserve the stability of problem (44) itself. Let h
=
0.033 and h
=
0.0033; the numerical solutions y
1
n
,
y
2
n and z
1
n
,
z
2
n are computed, and the errors y
n
-
z
n are shown in Figures 1 and 2, respectively, where y
n
=
[
y
1
n
,
y
2
n
]
T and z
n
=
[
z
1
n
,
z
2
n
]
T are given by applying (46) to (44) and its perturbed problem, which takes initial functions(47)
y
1
t
=
sin
⁡
6
t
+
0.1
,
y
2
t
=
cos
⁡
4
t
-
0.2
,
-
1
≤
t
≤
0
.Figure 1
The norm of differencey
n
-
z
n with step size h
=
0.033.Figure 2
The norm of differencey
n
-
z
n with step size h
=
0.0033.These results further confirm our conclusions.



---

*Source: 127027-2015-11-19.xml*