# Stability Analysis of Runge-Kutta Methods for Nonlinear Functional Differential and Functional Equations

Authors: Yuexin Yu; Zhongyan Liu; Liping Wen
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

This paper is concerned with the numerical stability of Runge-Kutta methods for a class of nonlinear functional differential and functional equations. The sufficient conditions for the stability and asymptotic stability of(
k
,
l
)-algebraically stable Runge-Kutta methods are derived. A numerical test is given to confirm the theoretical results.


---

## Body

## 1. Introduction

This paper is concerned with the numerical solution of the following nonlinear functional differential and functional equations (FDFEs):(1)
y
′
(
t
)
=
f
(
t
,
y
(
t
)
,
y
(
t
-
τ
)
,
z
(
t
-
τ
)
)
,
z
(
t
)
=
g
(
t
,
y
(
t
)
,
y
(
t
-
τ
)
,
z
(
t
-
τ
)
)
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
t
≥
0
,
with the initial conditions
(2)
y
(
t
)
=
φ
(
t
)
,
z
(
t
)
=
ψ
(
t
)
,
t
≤
0
,
where τ
>
0 is a real constant, y and z are unknown vectors of complex functions, f and g are given vectors of complex functions with appropriate domains of definition, and φ and ψ are given vectors of complex functions which satisfy the consistency relation
(3)
ψ
(
0
)
=
g
(
0
,
φ
(
0
)
,
φ
(
-
τ
)
,
ψ
(
-
τ
)
)
.Systems of the form (1) are sometimes called hybrid systems [1] or coupled delay differential and difference equations [2, 3]. They arise widely in the fields of science and technology, such as control systems, physics, and biology (see [1–6] and the references therein). In particular, they include neutral delay differential equations (NDDEs) as special cases. In fact, the explicit NDDEs
(4)
y
′
(
t
)
=
f
(
t
,
y
(
t
)
,
y
(
t
-
τ
)
,
y
′
(
t
-
τ
)
)
are equivalent to
(5)
y
′
(
t
)
=
f
(
t
,
y
(
t
)
,
y
(
t
-
τ
)
,
z
(
t
-
τ
)
)
,
z
(
t
)
=
f
(
t
,
y
(
t
)
,
y
(
t
-
τ
)
,
z
(
t
-
τ
)
)
,
while the implicit NDDEs
(6)
[
z
(
t
)
-
G
(
t
,
z
(
t
-
τ
)
)
]
′
=
f
(
t
,
z
(
t
)
,
z
(
t
-
τ
)
)
are equivalent to
(7)
y
′
(
t
)
=
f
(
t
,
y
(
t
)
+
G
(
t
,
z
(
t
-
τ
)
)
,
z
(
t
-
τ
)
)
,
z
(
t
)
=
y
(
t
)
+
G
(
t
,
z
(
t
-
τ
)
)
.In recent years, numerical methods for explicit NDDEs and implicit NDDEs have been studied extensively and a significant number of numerical stability results have been found (see [7–20]). However, the above results of numerical stability cannot be applied to the more general problem (1). In 1999, Liu [6] discussed the numerical stability of Runge-Kutta collocation methods with a constrained grid and linear θ-methods with a uniformed grid for linear systems of FDFEs:
(8)
y
′
(
t
)
+
A
1
y
(
t
)
+
A
2
y
(
t
-
τ
)
+
B
1
z
(
t
-
τ
)
=
0
,
z
(
t
)
+
A
3
y
(
t
)
+
A
4
y
(
t
-
τ
)
+
B
2
z
(
t
-
τ
)
=
0
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
t
≥
0
,
where τ is a positive constant and A
1, A
2, A
3, A
4, B
1, and B
2 are the coefficient matrices. Then, the asymptotic stability of linear multistep methods, one-leg methods, Runge-Kutta methods, multistep Runge-Kutta methods, and Rosenbrock methods for linear systems of FDFEs (8) was investigated in papers [21–23], respectively. Recently, Yu and Li [24] and Yu and Wen [25] dealt with the stability and asymptotic stability of the analytical and numerical solutions (obtained by one-leg methods) of nonlinear FDFEs (1), respectively. In the present paper, we further discuss the numerical stability of Runge-Kutta methods for the nonlinear FDFEs. The sufficient conditions for the stability and asymptotic stability of (
k
,
l
)-algebraically stable Runge-Kutta methods are derived.


## 2. Stability of the Problem ClassD
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
)

Let〈
·
,
·
〉 be an inner product and ∥
·
∥ the corresponding norm in complex N-dimensional space ℂ
N; assume that the mappings f and g in (1) satisfy the following conditions:
(9)
Re
〈
u
1
-
u
2
,
f
(
t
,
u
1
,
v
,
w
)
-
f
(
t
,
u
2
,
v
,
w
)
〉
≤
α
∥
u
1
-
u
2
∥
2
,
∀
t
≥
0
,
u
1
,
u
2
,
v
,
w
∈
ℂ
N
,
(10)
∥
f
(
t
,
u
,
v
1
,
w
1
)
-
f
(
t
,
u
,
v
2
,
w
2
)
∥
≤
β
1
∥
v
1
-
v
2
∥
+
β
2
∥
w
1
-
w
2
∥
,
h
h
∀
t
≥
0
,
u
,
v
1
,
v
2
,
w
1
,
w
2
∈
ℂ
N
,
(11)
∥
g
(
t
,
u
1
,
v
1
,
w
1
)
-
g
(
t
,
u
2
,
v
2
,
w
2
)
∥
≤
γ
1
∥
u
1
-
u
2
∥
+
γ
2
∥
v
1
-
v
2
∥
+
δ
∥
w
1
-
w
2
∥
,
h
h
h
h
h
h
h
h
h
∀
t
≥
0
,
u
1
,
u
2
,
v
1
,
v
2
,
w
1
,
w
2
∈
ℂ
N
,
where α, β
1, β
2, γ
1, γ
2, and δ are real constants and δ
<
1.Throughout this paper, we assume that the problem (1) has unique exact solution y
(
t
), z
(
t
) and denote the problem class consisting of all problems (1) with (9)–(11) by class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
).Remark 1.
Inequality 2.1 means that we admit of stiffness of the problem, that is, admitting large value for the classical Lipschitz constant off
(
t
,
u
,
v
,
w
) with respect to the second argument u (for the concept of stiffness we refer to [26, 27]).Remark 2.
Linear systems of FDFEs (8) belong to the class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
), where α
=
μ
(
-
A
1
), μ
(
·
) is the logarithmic matrix norm corresponding to the inner product norm in ℂ
N, and β
1
=
∥
A
2
∥, β
2
=
∥
B
1
∥, γ
1
=
∥
A
3
∥, γ
2
=
∥
A
4
∥, δ
=
∥
B
2
∥.For problems of the classD
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
), we have the following stability results (see [24]).Theorem 3.
Suppose the problem (1) belongs to the class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
) and α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
≤
0. Then one has the following inequalities:
(12)
∥
y
(
t
)
-
y
~
(
t
)
∥
≤
max
⁡
-
τ
≤
t
≤
0
{
∥
φ
(
t
)
-
φ
~
(
t
)
∥
,
∥
ψ
(
t
)
-
ψ
~
(
t
)
∥
}
,
t
>
0
,
∥
z
(
t
)
-
z
~
(
t
)
∥
≤
(
γ
1
+
γ
2
1
-
δ
+
δ
)
max
⁡
-
τ
≤
t
≤
0
{
∥
φ
(
t
)
-
φ
~
(
t
)
∥
,
∥
ψ
(
t
)
-
ψ
~
(
t
)
∥
}
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
t
>
0
.Here and later,y
~
(
t
), z
~
(
t
) denote the solution of any given perturbed problem of (1):
(13)
y
~
′
(
t
)
=
f
(
t
,
y
~
(
t
)
,
y
~
(
t
-
τ
)
,
z
~
(
t
-
τ
)
)
,
z
~
(
t
)
=
g
(
t
,
y
~
(
t
)
,
y
~
(
t
-
τ
)
,
z
~
(
t
-
τ
)
)
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
t
≥
0
,
with the initial conditions
(14)
y
~
(
t
)
=
φ
~
(
t
)
,
z
~
(
t
)
=
ψ
~
(
t
)
,
t
≤
0
,
which satisfy the consistency relation
(15)
ψ
~
(
0
)
=
g
(
0
,
φ
~
(
0
)
,
φ
~
(
-
τ
)
,
ψ
~
(
-
τ
)
)
.Theorem 4.
Suppose the problem (1) belongs to the class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
) and α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
<
0. Then one has
(16)
lim
⁡
t
→
+
∞
∥
y
(
t
)
-
y
~
(
t
)
∥
=
0
,
lim
⁡
t
→
+
∞
∥
z
(
t
)
-
z
~
(
t
)
∥
=
0
,
which characterizes the asymptotic stability property of the problem (1).


## 3. Stability Analysis of Runge-Kutta Methods for FDFEs

Ans-stage Runge-Kutta method for ordinary differential equations (ODEs) can be expressed as
(17)
c
A
b
T
=
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
s
,
where A
=
(
a
i
j
)
∈
ℝ
s
×
s, b
=
(
b
1
,
b
2
,
…
,
b
s
)
T
∈
ℝ
s, and c
=
(
c
1
,
c
2
,
…
,
c
s
)
T
∈
ℝ
s. In this paper we always assume that 0
≤
c
i
≤
1 (i
=
1,2
,
…
,
s) and ∑
j
=
1
s
‍
b
j
=
1.The adaptation of the Runge-Kutta method (17) for solving the problem (1) leads to
(18)
Y
i
(
n
)
=
y
n
+
h
∑
j
=
1
s
‍
a
i
j
f
(
t
n
+
c
j
h
,
Y
j
(
n
)
,
Y
j
(
n
-
m
)
,
Z
j
(
n
-
m
)
)
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
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
-
m
)
,
Z
j
(
n
-
m
)
)
,
h
z
n
=
g
(
t
n
,
y
n
,
y
n
-
m
,
z
n
-
m
)
,
where the integration step size h
=
τ
/
m, m is an arbitrarily given positive integer, t
n
=
n
h, Y
i
(
n
), y
n, and z
n denote approximations to y
(
t
n
+
c
i
h
), y
(
t
n
), and z
(
t
n
), respectively, y
n
=
φ
(
t
n
) and z
n
=
ψ
(
t
n
) for n
≤
0, and Z
i
(
n
) is an approximation to z
(
t
n
+
c
i
h
) which is obtained by using the following formula:
(19)
Z
i
(
n
)
=
g
(
t
n
+
c
i
h
,
Y
i
(
n
)
,
Y
i
(
n
-
m
)
,
Z
i
(
n
-
m
)
)
,
i
=
1,2
,
…
,
s
,
where Y
i
(
n
)
=
φ
(
t
n
+
c
i
h
) and Z
i
(
n
)
=
ψ
(
t
n
+
c
i
h
) for t
n
+
c
i
h
≤
0.Similarly, applying the same method to the perturbed problem (13), we have
(20)
Y
~
i
(
n
)
=
y
~
n
+
h
∑
j
=
1
s
‍
a
i
j
f
(
t
n
+
c
j
h
,
Y
~
j
(
n
)
,
Y
~
j
(
n
-
m
)
,
Z
~
j
(
n
-
m
)
)
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
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
-
m
)
,
Z
~
j
(
n
-
m
)
)
,
h
z
~
n
=
g
(
t
n
,
y
~
n
,
y
~
n
-
m
,
z
~
n
-
m
)
,
where Y
~
i
(
n
), y
~
n, and z
~
n denote approximations to y
~
(
t
n
+
c
i
h
), y
~
(
t
n
), and z
~
(
t
n
), respectively, y
~
n
=
φ
~
(
t
n
) and z
~
n
=
ψ
~
(
t
n
) for n
≤
0, and Z
~
i
(
n
) is an approximation to z
~
(
t
n
+
c
i
h
) which is obtained by using the following formula:
(21)
Z
~
i
(
n
)
=
g
(
t
n
+
c
i
h
,
Y
~
i
(
n
)
,
Y
~
i
(
n
-
m
)
,
Z
~
i
(
n
-
m
)
)
,
i
=
1,2
,
…
,
s
,
where Y
~
i
(
n
)
=
φ
~
(
t
n
+
c
i
h
) and Z
~
i
(
n
)
=
ψ
~
(
t
n
+
c
i
h
) for t
n
+
c
i
h
≤
0.Definition 5 (see [28]).
Letk
,
l be real constants with k
>
0. A Runge-Kutta method (17) 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
]
∈
ℝ
(
s
+
1
)
×
(
s
+
1
) is nonnegative definite, where
(22)
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
)
,
and e
=
(
1,1
,
…
,
1
)
T
∈
ℝ
s. Particularly, the (
1,0
)-algebraically stable method is called algebraically stable for short.Theorem 6.
Assume that the Runge-Kutta method (17) is (
k
,
l
)-algebraically stable with 0
<
k
≤
1. Then the numerical solutions y
n, z
n and y
~
n, z
~
n, obtained by applying the corresponding method (18) to the problems (1) and (13) which belong to the class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
) with (
α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
)
h
≤
l, respectively, satisfy the global stability inequalities
(23)
∥
y
n
-
y
~
n
∥
≤
C
max
⁡
-
τ
≤
t
≤
0
{
∥
φ
(
t
)
-
φ
~
(
t
)
∥
,
∥
ψ
(
t
)
-
ψ
~
(
t
)
∥
}
,
(24)
∥
z
n
-
z
~
n
∥
≤
(
C
(
γ
1
+
γ
2
)
1
-
δ
+
δ
)
max
⁡
-
τ
≤
t
≤
0
{
∥
φ
(
t
)
-
φ
~
(
t
)
∥
,
∥
ψ
~
(
t
)
∥
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
∥
ψ
(
t
)
-
ψ
~
(
t
)
∥
}
,
where C depends only on the method, β
1, β
2, γ
1, γ
2, δ, and τ.Proof.
Let(25)
w
n
=
y
n
-
y
~
n
,
W
i
(
n
)
=
Y
i
(
n
)
-
Y
~
i
(
n
)
,
i
=
1,2
,
…
,
s
,
Q
i
(
n
)
=
f
(
t
n
+
c
i
h
,
Y
i
(
n
)
,
Y
i
(
n
-
m
)
,
Z
i
(
n
-
m
)
)
-
f
(
t
n
+
c
i
h
,
Y
~
i
(
n
)
,
Y
~
i
(
n
-
m
)
,
Z
~
i
(
n
-
m
)
)
,
S
=
max
⁡
-
τ
≤
t
≤
0
{
∥
φ
(
t
)
-
φ
~
(
t
)
∥
,
∥
ψ
(
t
)
-
ψ
~
(
t
)
∥
}
,
and q
j
=
⌊
(
n
+
c
j
)
/
m
⌋ (j
=
1,2
,
…
,
s
,
n
>
0), where ⌊
·
⌋ denotes the integer part; then q
j
τ
≤
t
n
+
c
j
h
<
(
q
j
+
1
)
τ.
It follows from (18) and (20) that
(26)
W
i
(
n
)
=
w
n
+
h
∑
j
=
1
s
‍
a
i
j
Q
j
(
n
)
,
i
=
1,2
,
…
,
s
,
(27)
w
n
+
1
=
w
n
+
h
∑
j
=
1
s
‍
b
j
Q
j
(
n
)
,
n
=
0,1
,
…
.
Thus, it is easily obtained that (see [28])
(28)
∥
w
n
+
1
∥
2
-
k
∥
w
n
∥
2
-
2
∑
j
=
1
s
‍
d
j
Re
〈
W
j
(
n
)
,
h
Q
j
(
n
)
-
l
W
j
(
n
)
〉
=
-
∑
i
=
1
s
+
1
‍
∑
j
=
1
s
+
1
‍
M
i
j
〈
ξ
i
,
ξ
j
〉
,
where ξ
1
=
w
n, ξ
i
=
h
Q
i
-
1
(
n
), i
=
2,3
,
…
,
s
+
1. In view of (
k
,
l
)-algebraic stability of the method and 0
<
k
≤
1, we get
(29)
∥
w
n
+
1
∥
2
≤
∥
w
n
∥
2
+
2
∑
j
=
1
s
‍
d
j
Re
〈
W
j
(
n
)
,
h
Q
j
(
n
)
-
l
W
j
(
n
)
〉
.
By using conditions (9)–(11), we have
(30)
2
Re
〈
W
j
(
n
)
,
h
Q
j
(
n
)
〉
=
2
h
Re
〈
W
j
(
n
)
,
f
(
t
n
+
c
j
h
,
Y
j
(
n
)
,
Y
j
(
n
-
m
)
,
Z
j
(
n
-
m
)
)
h
h
h
h
h
h
h
h
h
h
-
f
(
t
n
+
c
j
h
,
Y
~
j
(
n
)
,
Y
~
j
(
n
-
m
)
,
Z
~
j
(
n
-
m
)
)
〉
≤
2
h
α
∥
W
j
(
n
)
∥
2
+
2
h
∥
W
j
(
n
)
∥
×
∥
f
(
t
n
+
c
j
h
,
Y
~
j
(
n
)
,
Y
j
(
n
-
m
)
,
Z
j
(
n
-
m
)
)
h
h
h
h
h
h
-
f
(
t
n
+
c
j
h
,
Y
~
j
(
n
)
,
Y
~
j
(
n
-
m
)
,
Z
~
j
(
n
-
m
)
)
∥
.
When t
n
+
c
j
h
-
τ
<
0, that is, q
j
=
0, (30) leads to
(31)
2
Re
〈
W
j
(
n
)
,
h
Q
j
(
n
)
〉
≤
2
h
α
∥
W
j
(
n
)
∥
2
+
2
h
∥
W
j
(
n
)
∥
(
β
1
+
β
2
)
S
≤
(
2
α
+
β
1
+
β
2
)
h
∥
W
j
(
n
)
∥
2
+
h
(
β
1
+
β
2
)
S
2
.
On the other hand, when t
n
+
c
j
h
-
τ
≥
0, that is, q
j
≥
1, using conditions (9)–(11) and δ
<
1, (30) leads to
(32)
2
Re
〈
W
j
(
n
)
,
h
Q
j
(
n
)
〉
≤
2
h
α
∥
W
j
(
n
)
∥
2
+
2
h
∥
W
j
(
n
)
∥
×
(
β
1
∥
W
j
(
n
-
m
)
∥
+
β
2
∥
Z
j
(
n
-
m
)
-
Z
~
j
(
n
-
m
)
∥
)
≤
2
h
α
∥
W
j
(
n
)
∥
2
+
2
h
∥
W
j
(
n
)
∥
×
[
β
1
∥
W
j
(
n
-
m
)
∥
+
β
2
(
γ
1
∥
W
j
(
n
-
m
)
∥
+
γ
2
∥
W
j
(
n
-
2
m
)
∥
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
δ
∥
Z
j
(
n
-
2
m
)
-
Z
~
j
(
n
-
2
m
)
∥
)
]
≤
2
h
α
∥
W
j
(
n
)
∥
2
+
2
h
∥
W
j
(
n
)
∥
×
(
β
1
∥
W
j
(
n
-
m
)
∥
+
β
2
γ
1
∑
i
=
1
q
j
‍
δ
i
-
1
∥
W
j
(
n
-
i
m
)
∥
h
h
h
h
h
h
h
+
β
2
γ
2
∑
i
=
2
q
j
‍
δ
i
-
2
∥
W
j
(
n
-
i
m
)
∥
+
β
2
γ
2
δ
q
j
-
1
S
+
β
2
δ
q
j
S
)
≤
2
h
α
∥
W
j
(
n
)
∥
2
+
h
β
1
(
∥
W
j
(
n
)
∥
2
+
∥
W
j
(
n
-
m
)
∥
2
)
+
h
β
2
γ
1
∑
i
=
1
q
j
‍
δ
i
-
1
(
∥
W
j
(
n
)
∥
2
+
∥
W
j
(
n
-
i
m
)
∥
2
)
+
h
β
2
γ
2
∑
i
=
2
q
j
‍
δ
i
-
2
(
∥
W
j
(
n
)
∥
2
+
∥
W
j
(
n
-
i
m
)
∥
2
)
+
h
β
2
δ
q
j
-
1
(
γ
2
+
δ
)
(
∥
W
j
(
n
)
∥
2
+
S
2
)
≤
(
2
α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
h
∥
W
j
(
n
)
∥
2
+
h
β
1
∥
W
j
(
n
-
m
)
∥
2
+
h
β
2
γ
1
∑
i
=
1
q
j
‍
δ
i
-
1
∥
W
j
(
n
-
i
m
)
∥
2
+
h
β
2
γ
2
∑
i
=
2
q
j
‍
δ
i
-
2
∥
W
j
(
n
-
i
m
)
∥
2
+
h
β
2
δ
q
j
-
1
(
γ
2
+
δ
)
S
2
.
Here and below, we define ∑
k
=
s
t
‍ equal to 0 for t
<
s. Combining (31) and (32) yields
(33)
2
Re
〈
W
j
(
n
)
,
h
Q
j
(
n
)
〉
≤
(
2
α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
h
∥
W
j
(
n
)
∥
2
+
h
β
1
∥
W
j
(
n
-
m
)
∥
2
+
h
β
2
γ
1
∑
i
=
1
q
j
‍
δ
i
-
1
∥
W
j
(
n
-
i
m
)
∥
2
+
h
β
2
γ
2
∑
i
=
2
q
j
‍
δ
i
-
2
∥
W
j
(
n
-
i
m
)
∥
2
+
h
β
2
δ
q
j
-
1
(
γ
2
+
δ
)
S
2
.
Substituting (33) into (29) and using condition (
α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
)
h
≤
l, we obtain
(34)
∥
w
n
+
1
∥
2
≤
∥
w
n
∥
2
+
∑
j
=
1
s
‍
d
j
[
(
(
β
2
(
γ
1
+
γ
2
)
1
-
δ
2
α
+
β
1
+
β
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
h
-
2
l
)
∥
W
j
(
n
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
1
∥
W
j
(
n
-
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
γ
1
∑
i
=
1
q
j
‍
δ
i
-
1
∥
W
j
(
n
-
i
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
γ
2
∑
i
=
2
q
j
‍
δ
i
-
2
∥
W
j
(
n
-
i
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
δ
q
j
-
1
(
γ
2
+
δ
)
S
2
(
(
2
α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
h
-
2
l
)
]
≤
∥
w
n
∥
2
+
∑
j
=
1
s
‍
d
j
[
-
(
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
h
∥
W
j
(
n
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
1
∥
W
j
(
n
-
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
γ
1
∑
i
=
1
q
j
δ
i
-
1
∥
W
j
(
n
-
i
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
γ
2
∑
i
=
2
q
j
‍
δ
i
-
2
∥
W
j
(
n
-
i
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
δ
q
j
-
1
(
γ
2
+
δ
)
S
2
(
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
]
.
By induction, (34) gives
(35)
∥
w
n
+
1
∥
2
≤
∥
w
0
∥
2
+
∑
k
=
0
n
‍
∑
j
=
1
s
‍
d
j
[
-
(
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
×
h
∥
W
j
(
k
)
∥
2
+
h
β
1
∥
W
j
(
k
-
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
γ
1
∑
i
=
1
⌊
(
k
+
c
j
)
/
m
⌋
‍
δ
i
-
1
∥
W
j
(
k
-
i
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
γ
2
∑
i
=
2
⌊
(
k
+
c
j
)
/
m
⌋
‍
δ
i
-
2
∥
W
j
(
k
-
i
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
(
γ
2
+
δ
)
S
2
δ
δ
⌊
(
k
+
c
j
)
/
m
⌋
]
≤
∥
w
0
∥
2
+
∑
j
=
1
s
‍
d
j
[
-
(
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
×
h
∑
k
=
0
n
‍
∥
W
j
(
k
)
∥
2
+
h
β
1
∑
k
=
-
m
n
‍
∥
W
j
(
k
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
γ
1
1
-
δ
∑
k
=
-
m
n
‍
∥
W
j
(
k
-
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
γ
2
1
-
δ
∑
k
=
-
m
n
‍
∥
W
j
(
k
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
β
2
(
γ
2
+
δ
)
τ
S
2
δ
(
1
-
δ
)
(
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
h
∑
k
=
0
n
‍
]
≤
∥
w
0
∥
2
+
∑
j
=
1
s
‍
d
j
[
(
-
(
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
1
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
h
)
∑
k
=
0
n
‍
∥
W
j
(
k
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
τ
(
β
1
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
max
⁡
-
m
≤
k
≤
-
1
∥
W
j
(
k
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
β
2
(
γ
2
+
δ
)
τ
S
2
δ
(
1
-
δ
)
(
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
]
≤
∥
w
0
∥
2
+
∑
j
=
1
s
‍
d
j
(
β
1
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
+
β
2
(
γ
2
+
δ
)
δ
(
1
-
δ
)
)
τ
S
2
.
Therefore, there is a real constant C depending only on the method, β
1, β
2, γ
1, γ
2, δ, and τ such that the inequality (23) holds. On the other hand, using condition (11) and δ
<
1, we have
(36)
∥
z
n
-
z
~
n
∥
=
∥
g
(
t
n
,
y
n
,
y
n
-
m
,
z
n
-
m
)
-
g
(
t
n
,
y
~
n
,
y
~
n
-
m
,
z
~
n
-
m
)
∥
≤
γ
1
∥
w
n
∥
+
γ
2
∥
w
n
-
m
∥
+
δ
∥
z
n
-
m
-
z
~
n
-
m
∥
≤
(
C
(
γ
1
+
γ
2
)
1
-
δ
+
δ
)
×
max
⁡
-
τ
≤
t
≤
0
{
∥
φ
(
t
)
-
φ
~
(
t
)
∥
,
∥
ψ
(
t
)
-
ψ
~
(
t
)
∥
}
,
and this completes the proof of Theorem 6.Particularly, for the algebraically stable Runge-Kutta method, we have the following.Corollary 7.
Assume that the Runge-Kutta method (17) is algebraically stable. Then the numerical solutions y
n, z
n and y
~
n, z
~
n, obtained by applying the corresponding method (18) to the problems (1) and (13) which belong to the class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
) with α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
≤
0, respectively, satisfy the global stability inequalities
(37)
∥
y
n
-
y
~
n
∥
≤
C
max
⁡
-
τ
≤
t
≤
0
{
∥
φ
(
t
)
-
φ
~
(
t
)
∥
,
∥
ψ
(
t
)
-
ψ
~
(
t
)
∥
}
,
∥
z
n
-
z
~
n
∥
≤
(
C
(
γ
1
+
γ
2
)
1
-
δ
+
δ
)
max
⁡
-
τ
≤
t
≤
0
{
∥
φ
(
t
)
-
φ
~
(
t
)
∥
,
∥
ψ
~
(
t
)
∥
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
∥
ψ
(
t
)
-
ψ
~
(
t
)
∥
}
,
where C depends only on the method, β
1, β
2, γ
1, γ
2, δ, and τ.Remark 8.
It is well known that the formulae Gauss, Radau IA, Radau IIA, and Lobatto IIIC (for ODEs) are all algebraically stable. Therefore, in terms of Corollary7, the corresponding methods are globally stable for solving the nonlinear FDFEs of the class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
) which satisfy the condition α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
≤
0.In the following, we further discuss the asymptotic stability of the Runge-Kutta method.Theorem 9.
Assume that the Runge-Kutta method (17) is (
k
,
l
)-algebraically stable with 0
<
k
≤
1, D
>
0 (i.e., the matrix D is positive definite), det
⁡
A
≠
0, and |
1
-
b
T
A
-
1
e
|
<
1. Then the numerical solutions y
n, z
n and y
~
n, z
~
n, obtained by applying the corresponding method (18) to the problems (1) and (13) which belong to the class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
) with (
α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
)
h
<
l, respectively, satisfy
(38)
lim
⁡
n
→
+
∞
∥
y
n
-
y
~
n
∥
=
0
,
lim
⁡
n
→
+
∞
∥
z
n
-
z
~
n
∥
=
0
.
The relations (38) characterize the asymptotic stability property of the method.Proof.
Letσ
=
(
α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
)
h
-
l; thus σ
<
0. In terms of the proof of Theorem 6, we have
(39)
∥
w
n
+
1
∥
2
≤
∥
w
n
∥
2
+
∑
j
=
1
d
‍
d
j
[
(
2
σ
-
(
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
h
)
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
×
∥
W
j
(
n
)
∥
2
+
h
β
1
∥
W
j
(
n
-
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
γ
1
∑
i
=
1
q
j
‍
δ
i
-
1
∥
W
j
(
n
-
i
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
γ
2
∑
i
=
2
q
j
‍
δ
i
-
2
∥
W
j
(
n
-
i
m
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
h
β
2
δ
q
j
-
1
(
γ
2
+
δ
)
S
2
(
(
β
2
(
γ
1
+
γ
2
)
1
-
δ
)
)
]
.
By induction, (39) gives
(40)
∥
w
n
+
1
∥
2
≤
∥
w
0
∥
2
+
∑
j
=
1
s
‍
d
j
[
(
β
2
(
γ
2
+
δ
)
δ
(
1
-
δ
)
)
(
2
σ
-
β
2
h
)
∑
k
=
0
n
‍
∥
W
j
(
k
)
∥
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
(
β
1
+
β
2
(
γ
1
+
γ
2
)
1
-
δ
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
β
2
(
γ
2
+
δ
)
δ
(
1
-
δ
)
)
τ
S
2
∑
k
=
0
n
‍
]
.
Since (
2
σ
-
β
2
h
)
<
0 and D
>
0, we easily obtain that
(41)
lim
⁡
n
→
+
∞
∥
W
j
(
n
)
∥
=
0
,
j
=
1,2
,
…
,
s
.
On the other hand, while det
⁡
A
≠
0, denote G
=
[
g
i
j
]
=
A
-
1; thus (26) yields
(42)
Q
i
(
n
)
=
1
h
∑
j
=
1
s
‍
g
i
j
(
W
j
(
n
)
-
w
n
)
,
i
=
1,2
,
…
,
s
.
Substituting (42) into (27) leads to
(43)
w
n
+
1
=
(
1
-
b
T
A
-
1
e
)
w
n
+
∑
i
=
1
s
‍
∑
j
=
1
s
‍
b
i
g
i
j
W
j
(
n
)
.
Noting that |
1
-
b
T
A
-
1
e
|
<
1 and (41), we easily obtain
(44)
lim
⁡
n
→
+
∞
∥
w
n
∥
=
lim
⁡
n
→
+
∞
∥
y
n
-
y
~
n
∥
=
0
.
Furthermore, using condition (11), we have
(45)
∥
z
n
-
z
~
n
∥
=
∥
g
(
t
n
,
y
n
,
y
n
-
m
,
z
n
-
m
)
-
g
(
t
n
,
y
~
n
,
y
~
n
-
m
,
z
~
n
-
m
)
∥
≤
γ
1
∥
w
n
∥
+
γ
2
∥
w
n
-
m
∥
+
δ
∥
z
n
-
m
-
z
~
n
-
m
∥
.
Considering that δ
<
1 and (44), (45) leads to
(46)
lim
⁡
n
→
+
∞
∥
z
n
-
z
~
n
∥
=
0
,
and this completes the proof of Theorem 9.Particularly, for the algebraically stable Runge-Kutta method, we have the following.Corollary 10.
Assume that the Runge-Kutta method (17) is algebraically stable with D
>
0, det
⁡
A
≠
0, and |
1
-
b
T
A
-
1
e
|
<
1. Then the numerical solutions y
n, z
n and y
~
n, z
~
n, obtained by applying the corresponding method (18) to the problems (1) and (13) which belong to the class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
) with α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
<
0, respectively, satisfy
(47)
lim
⁡
n
→
+
∞
∥
y
n
-
y
~
n
∥
=
0
,
lim
⁡
n
→
+
∞
∥
z
n
-
z
~
n
∥
=
0
.Remark 11.
It is well known that the formulae Radau IA, Radau IIA, and Lobatto IIIC (for ODEs) are algebraically stable withD
>
0, det
⁡
A
≠
0 and |
1
-
b
T
A
-
1
e
|
<
1. Therefore, in terms of Corollary 10, the corresponding methods are asymptotically stable for solving the nonlinear FDFEs of the class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
) which satisfy the condition α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
<
0.Remark 12.
In the paper [25], it is proved that an A-stable one-leg method is globally stable and a strongly A-stable one-leg method is asymptotically stable for FDFEs. However, any A-stable one-leg method has order at most two. In the present paper, the stability results are based on (
k
,
l
)-algebraic stability of Runge-Kutta methods, which, in general, can be of high order.


## 4. Numerical Experiments

Consider the following initial value problem:(48)
∂
u
(
x
,
t
)
∂
t
=
∂
2
u
(
x
,
t
)
∂
x
2
+
sin
u
(
x
,
t
-
π
)
cos
⁡
v
(
x
,
t
-
π
)
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
0
<
x
<
1
,
t
≥
0
,
v
(
x
,
t
)
=
sin
u
(
x
,
t
)
sin
u
(
x
,
t
-
π
)
+
0.2
v
(
x
,
t
-
π
)
,
u
(
x
,
t
)
=
(
x
-
x
2
)
sin
t
,
v
(
x
,
t
)
=
(
x
2
-
x
)
sin
t
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
0
<
x
<
1
,
t
≤
0
,
u
(
0
,
t
)
=
u
(
1
,
t
)
=
0
,
v
(
0
,
t
)
=
v
(
1
,
t
)
=
0
.
After application of the numerical method of lines, we obtain the following FDFEs:
(49)
u
i
′
(
t
)
=
1
Δ
x
2
(
u
i
-
1
(
t
)
-
2
u
i
(
t
)
+
u
i
+
1
(
t
)
)
+
sin
u
i
(
t
-
π
)
cos
⁡
v
i
(
t
-
π
)
,
t
≥
0
,
v
i
(
t
)
=
sin
u
i
(
t
)
sin
u
i
(
t
-
π
)
+
0.2
v
i
(
t
-
π
)
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
=
1,2
,
…
,
N
-
1
,
u
i
(
t
)
=
i
Δ
x
(
1
-
i
Δ
x
)
sin
t
,
v
i
(
t
)
=
i
Δ
x
(
i
Δ
x
-
1
)
sin
t
,
t
≤
0
,
u
0
(
t
)
=
u
N
(
t
)
≡
0
,
v
0
(
t
)
=
v
N
(
t
)
≡
0
,
where Δ
x is the spatial step, N is a natural number such that N
Δ
x
=
1, x
i
=
i
Δ
x, i
=
1,2
,
…
,
N
-
1, and u
i
(
t
)
=
u
(
x
i
,
t
), v
i
(
t
)
=
v
(
x
i
,
t
). Then, the problem (49) belongs to the class D
(
α
,
β
1
,
β
2
,
γ
1
,
γ
2
,
δ
) with
(50)
α
=
-
4
N
2
sin
2
π
2
N
,
β
1
=
β
2
=
γ
1
=
γ
2
=
1
,
δ
=
0.2
,
where the inner product is standard inner product. We take Δ
x
=
0.01 (i.e., N
=
100) for the numerical method of lines; thus the condition α
+
β
1
+
β
2
+
β
2
(
γ
1
+
γ
2
)
/
(
1
-
δ
)
≤
0 (<
0) is satisfied, which means the analytical solution of the problem (49) is stable and asymptotically stable.As an example, we consider the 2-stage Radau IIA method:(51)
c
A
b
T
=
1
/
3
5
/
12
-
1
/
12
1
3
/
4
1
/
4
3
/
4
1
/
4
,
for solving the problem (49) and its perturbed problem, where the initial conditions of the perturbed problem are
(52)
u
~
i
(
t
)
=
5
i
Δ
x
(
1
-
i
Δ
x
)
sin
t
,
v
~
i
(
t
)
=
5
i
Δ
x
(
i
Δ
x
-
1
)
sin
t
,
2222222222222222
t
≤
0
.According to the results of Corollaries7 and 10, the corresponding method (for FDFEs) will be stable and asymptotically stable. We denote the numerical solutions of problem (49) and its perturbed problem u
n, v
n and u
~
n, v
~
n, respectively, where u
n and v
n are approximations to [
u
1
(
t
n
)
,
u
2
(
t
n
)
,
…
,
u
N
-
1
(
t
n
)
]
T and [
v
1
(
t
n
)
,
v
2
(
t
n
)
,
…
,
v
N
-
1
(
t
n
)
]
T, respectively. The values ∥
u
n
-
u
~
n
∥ and ∥
v
n
-
v
~
n
∥ are listed in Figure 1 (where the abscissa denotes variable t).Values∥
u
n
-
u
~
n
∥ (a) and values ∥
v
n
-
v
~
n
∥ (b) with h
=
π
/
100.
(a)
(b)As a comparison, we consider the explicit 3-stage Runge-Kutta method(53)
c
A
b
T
=
0
0
0
0
1
/
2
1
/
2
0
0
1
-
1
2
0
1
/
6
4
/
6
1
/
6
,
for solving the problem (49) and its perturbed problem, and list the values ∥
u
n
-
u
~
n
∥ and ∥
v
n
-
v
~
n
∥ in Figure 2.Values∥
u
n
-
u
~
n
∥ (a) and values ∥
v
n
-
v
~
n
∥ (b) with h
=
π
/
100.
(a)
(b)From Figure1, one can see that the values ∥
u
n
-
u
~
n
∥ and ∥
v
n
-
v
~
n
∥ are bounded and tend to zero. This coincides with the results of Corollaries 7 and 10. However, for the explicit 3-stage Runge-Kutta method (53), which is not algebraically stable, the situation is inverse as one can see that the values ∥
u
n
-
u
~
n
∥ are divergent as n
→
+
∞.



---

*Source: 607827-2014-05-06.xml*