# Efficient Two-Derivative Runge-Kutta-Nyström Methods for Solving General Second-Order Ordinary Differential Equationsy

Authors: T. S. Mohamed; N. Senu; Z. B. Ibrahim; N. M. A. Nik Long
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

This paper proposes and investigates a special class of explicit Runge-Kutta-Nyström (RKN) methods for problems in the formy
′
′
(
x
)
=
f
(
x
,
y
,
y
′
) including third derivatives and denoted as STDRKN. The methods involve one evaluation of second derivative and many evaluations of third derivative per step. In this study, methods with two and three stages of orders four and five, respectively, are presented. The stability property of the methods is discussed. Numerical experiments have clearly shown the accuracy and the efficiency of the new methods.


---

## Body

## 1. Introduction

In this article, we are interested in initial value problems (IVPs) of second-order ordinary differential equations (ODEs):(1)
y
′
′
x
=
f
x
,
y
x
,
y
′
x
,
y
x
0
=
α
,
y
′
x
0
=
β
,
x
∈
x
0
,
x
e
n
d
,where y
∈
R
N
,
f
:
R
×
R
N
×
R
N
→
R
N are continuous vector valued functions. This type of problems arises naturally in many applied science fields such as the Kepler problems in celestial mechanics, quantum physics, and Newton’s second law in classical mechanics (see Dormand [1], Hairer et al. [2], and Kristensson [3]).Problems (1) in which the first derivative does not appear explicitly are an important subclass of second order (ODEs). Thus, several numerical methods for directly solving this subclass have been presented (see Dormand [1], Hairer et al. [2], Butcher [4], Lambert [5], and Senu [6]). In the case of direct solutions for the general second order (IVPs), some numerical methods have been proposed (see Chen et al. [7], Franco [8], Jator [9], Awoyemi [10], Wu et al. [11], Wu and Wang [12], and Chawla and Sharma [13]). The objective of this paper is to design STDRKN methods with a minimal number of function evaluation. This paper is organized as follows: In Section 3, we construct STDRKN methods; the stability analysis of STDRKN methods is discussed in Section 4; and numerical results are given in Section 5.


## 2. The Formulation of STDRKN Methods

In many problems in applications the third derivative(2)
y
′
′
′
x
=
g
x
,
y
,
y
′
=
f
x
x
,
y
,
y
′
+
f
y
x
,
y
,
y
′
y
′
+
f
y
′
x
,
y
,
y
′
f
x
,
y
,
y
′is available and easy to obtain. This derivative can be computed but the evaluation of g
(
x
,
y
,
y
′
) requires the evaluation of f
(
x
,
y
,
y
′
), f
x
(
x
,
y
,
y
′
), f
y
′
(
x
,
y
,
y
′
)
f
y
(
x
,
y
,
y
′
)
y
′. Therefore, in the scalar case (differential systems of dimension one), an evaluation of the third derivative g
(
x
,
y
,
y
′
) can be as expensive as four evaluations of the second derivative f
(
x
,
y
,
y
′
) and at least as two f-evaluations. An s-stage two-derivative Runge-Kutta-Nyström (TDRKN) method for (1) is defined by the formula (see Chen et al. [7])(3)
y
n
+
1
=
y
n
+
h
y
n
′
+
h
2
∑
i
=
1
s
b
i
f
x
n
+
c
i
h
,
Y
i
,
Y
i
′
+
h
3
∑
i
=
1
s
d
i
g
x
n
+
c
i
h
,
Y
i
,
Y
i
′
,
y
n
+
1
′
=
y
n
′
+
h
∑
i
=
1
s
b
i
′
f
x
n
+
c
i
h
,
Y
i
,
Y
i
′
+
h
2
∑
i
=
1
s
d
i
′
g
x
n
+
c
i
h
,
Y
i
,
Y
i
′
,where(4)
Y
i
=
y
n
+
c
i
h
y
n
′
+
h
2
∑
j
=
1
s
a
i
,
j
f
x
n
+
c
i
h
,
Y
j
,
Y
j
′
+
h
3
∑
j
=
1
s
r
i
,
j
g
x
n
+
c
i
h
,
Y
j
,
Y
j
′
,
Y
i
′
=
y
n
′
+
h
∑
j
=
1
s
a
i
,
j
′
f
x
n
+
c
i
h
,
Y
j
,
Y
j
′
+
h
2
∑
j
=
1
s
r
i
,
j
′
g
x
n
+
c
i
h
,
Y
j
,
Y
j
′
,where c
i
,
b
i
,
d
i
,
b
i
′
,
d
i
′
,
a
i
,
j
,
r
i
,
j
,
a
i
,
j
′
,
r
i
,
j
′
,
i
,
j
=
1
,
…
,
s, are real numbers. This method can also be written in Butcher’s tableau of coefficients as given in Table 1.Table 1
Butcher tableau for TDRKN methods.
C
A
R
A
′
R
′
b
T
d
T
b
′
T
d
′
TIn this paper, a special part of TDRKN method is studied that has the form(5)
y
n
+
1
=
y
n
+
h
y
n
′
+
h
2
2
f
x
n
,
y
n
,
y
n
′
+
h
3
∑
i
=
1
s
b
i
g
x
n
+
c
i
h
,
Y
i
,
Y
i
′
,
y
n
+
1
′
=
y
n
′
+
h
f
x
n
,
y
n
,
y
n
′
+
h
2
∑
i
=
1
s
d
i
g
x
n
+
c
i
h
,
Y
i
,
Y
i
′
,where(6)
Y
i
=
y
n
+
c
i
h
y
n
′
+
1
2
c
i
2
h
2
f
x
n
,
y
n
,
y
n
′
+
h
3
∑
j
=
1
s
a
i
,
j
g
x
n
+
c
i
h
,
Y
j
,
Y
j
′
,
Y
i
′
=
y
n
′
+
c
i
h
f
x
n
,
y
n
,
y
n
′
+
h
2
∑
j
=
1
s
r
i
,
j
g
x
n
+
c
i
h
,
Y
j
,
Y
j
′
. An alternative expression of formula (5) is given as follows:(7)
y
n
+
1
=
y
n
+
h
y
n
′
+
h
2
2
f
x
n
,
y
n
,
y
n
′
+
h
3
∑
i
=
1
s
b
i
k
i
,
y
n
+
1
′
=
y
n
′
+
h
f
x
n
,
y
n
,
y
n
′
+
h
2
∑
i
=
1
s
d
i
k
i
,where(8)
k
i
=
g
x
n
+
c
i
,
y
n
+
c
i
h
y
n
′
+
1
2
c
i
2
h
2
f
x
n
,
y
n
,
y
n
′
+
h
3
∑
j
=
1
s
a
i
,
j
k
j
,
y
n
′
+
c
i
h
f
x
n
,
y
n
,
y
n
′
+
h
2
∑
j
=
1
s
r
i
,
j
k
j
.This STDRKN method can be written in Butcher’s tableau as shown in Table2. The STDRKN methods are explicit methods if a
i
,
j
=
0
,
r
i
,
j
=
0 for i
≤
j and are implicit method if a
i
,
j
≠
0
,
r
i
,
j
≠
0 for i
≤
j
. STDRKN methods involve only one evaluation of f and many evaluations of g per step.Table 2
Butcher tableau for STDRKN methods.
C
A
R
b
T
d
T


## 3. Construction of STDRKN Methods

In this section, our effort is to determine the coefficients of the STDRKN methods as given in (7). Hence, using the Taylor series expansion in (3) with the Taylor series expansion of y
n and y
n
′ and by comparing the coefficients of the power of h, we obtained the order conditions of STDRKN methods for y and y
′ as in (10)–(18), while the rooted trees for STDRKN methods up to order five based on [7] are given in Table 3.Table 3
Root trees for STDRKN methods up to order five.
Order
ρ
(
t
)
Tree
t
α
(
t
)
Density
ϒ
(
t
)
Elementary weight
Φ
(
t
)
Elementary differential
F
(
t
)
(
y
,
y
′
)
0
ϕ
1
1
y
1
1
1
y
′
2
1
2
e
f
3
1
6
c
f
y
′
y
′
3
1
6
c
f
y
′
′
f
4
1
12
c
2
f
y
y
′
′
(
y
′
,
y
′
)
4
2
12
c
2
f
y
y
′
′
′
(
y
′
,
f
)
4
1
12
c
2
/
2
f
y
′
y
′
′
′
(
f
,
f
)
4
1
24
c
f
y
′
′
f
y
′
y
′
4
1
24
1
/
2
c
2
f
y
′
y
′
4
1
24
c
f
y
′
′
f
y
′
′
f
5
1
20
c
3
f
y
y
y
3
(
y
′
,
y
′
,
y
′
)
5
3
20
c
3
f
y
y
y
′
3
(
y
′
,
y
′
,
f
)
5
3
20
c
3
f
y
y
′
y
′
3
(
y
′
,
f
,
f
)
5
1
20
1
/
2
c
2
f
y
′
y
′
y
′
3
(
f
,
f
,
f
)
5
3
40
R
c
f
y
y
′
′
(
y
′
,
f
)
5
1
40
1
/
2
c
3
f
y
y
′
′
′
(
y
′
,
f
y
′
y
′
)
5
3
40
R
c
f
y
y
′
′
′
(
y
′
,
f
y
′
′
f
)
5
1
40
R
c
f
y
′
y
′
′
(
f
,
f
)
5
1
40
R
c
f
y
′
y
′
′
′
(
f
,
f
y
′
y
′
)
5
1
60
c
2
f
y
′
′
f
y
y
′
′
(
y
′
,
y
′
)
5
1
60
c
2
f
y
′
′
f
y
y
′
′
′
(
y
′
,
f
)
5
1
60
c
2
f
y
′
′
f
y
′
y
′
′
′
(
f
,
f
)
5
2
120
c
f
y
′
f
y
′
y
′
5
1
120
c
f
y
′
f
y
′
′
f
5
1
120
c
f
y
′
′
f
y
′
f
5
2
120
c
f
y
′
′
f
y
′
′
f
y
′
y
′
5
2
120
c
f
y
′
′
f
y
′
′
f
y
′
′
y
′The following simplifying assumption is suggested in practice:(9)
∑
j
=
1
i
-
1
r
i
,
j
=
c
i
2
2
,
i
=
2
,
…
,
s
.The following are the order conditions for explicit STDRKN.The order conditions fory:Third Order
(10)
∑
i
=
1
s
b
i
=
1
6
.
Fourth Order
(11)
∑
i
=
1
s
b
i
c
i
=
1
24
.
Fifth Order
(12)
1
2
∑
i
=
1
s
b
i
c
i
2
=
1
120
.
Sixth Order
(13)
∑
i
=
1
s
b
i
c
i
3
=
1
120
,
∑
i
=
1
s
∑
j
=
1
i
-
1
b
i
a
i
,
j
=
1
720
,
∑
i
=
1
s
∑
j
=
2
i
-
1
b
i
r
i
,
j
c
j
=
1
720
.The order conditions for y
′:Second Order
(14)
∑
i
=
1
s
d
i
=
1
2
.
Third Order
(15)
∑
i
=
1
s
d
i
c
i
=
1
6
.
Fourth Order
(16)
∑
i
=
1
s
d
i
c
i
2
=
1
12
.
Fifth Order
(17)
∑
i
=
1
s
d
i
c
i
3
=
1
20
,
∑
i
=
1
s
∑
j
=
1
i
-
1
d
i
a
i
,
j
=
1
120
,
∑
i
=
1
s
∑
j
=
1
i
-
1
d
i
r
i
,
j
c
j
=
1
120
.
Sixth Order
(18)
∑
i
=
1
s
d
i
c
i
4
=
1
30
,
∑
i
=
1
n
∑
j
=
1
i
-
1
d
i
c
i
a
i
,
j
=
1
180
,
∑
i
=
1
s
∑
j
=
1
i
-
1
d
i
a
i
,
j
c
j
=
1
720
,
∑
i
=
1
s
∑
j
=
1
i
-
1
d
i
r
i
,
j
c
j
2
=
1
360
,
∑
i
=
1
s
∑
j
=
1
i
-
1
d
i
c
i
r
i
,
j
c
j
=
1
180
.


### 3.1. Two-Stage STDRKN Method of Order Four

To derive the fourth-order STDRKN method, we use the algebraic order conditions up to order four in the equationsy and y
′, that is, (10)-(11) and (14)–(16), respectively. However, we get a system of equations which consists of 5 nonlinear equations with 6 unknown variables. Solving these equations with the simplifying assumption r
2,1
=
c
2
2
/
2, the resulting system has one free parameter a
2,1; this free parameter can be selected by minimizing error equations. The truncation error norms and global error of the fifth-order condition are defined as follows:(19)
τ
5
2
=
∑
i
=
1
n
1
τ
i
5
2
=
1
480
,
τ
′
5
2
=
∑
i
=
1
n
1
′
τ
i
′
5
2
=
1
120
2
+
1600
a
2,1
2
-
80
a
2,1
,
τ
g
5
2
=
∑
i
=
1
n
1
τ
i
5
2
+
∑
i
=
1
n
1
′
τ
i
′
5
2
=
1
480
33
+
25600
a
2,1
2
-
1280
a
2,1
.The error equations τ
g
5
2 have a minimum value at a
2,1
=
1
/
40 which produces(20)
τ
5
2
=
2.08
×
1
0
-
3
,
τ
′
5
2
=
8.3
×
1
0
-
3
,
τ
g
5
2
=
8.5
×
1
0
-
3
,where τ
(
5
) and τ
′
(
5
) are error terms of the fifth-order conditions for y
n and y
n
′, respectively. The coefficients of the new method can be given in Butcher tableau and denoted by STDRKN4(
2
) as seen in Table 4.Table 4
The STDRKN4(
2
) method.
0
0
0
1/2
1
/
40
0
1
/
8
0
1
/
12
1
/
12
1
/
6
1
/
3


### 3.2. Three-Stage STDRKN Method of Order Five

According to the order conditions in (10)–(12) and (14)–(17), we obtain system of equations that consists of 9 nonlinear equations with 14 unknowns variables that need to be solved. Solving the system simultaneously with simplifying assumption conditions r
2,1
=
c
2
2
/
2
,
r
3,1
=
c
3
2
/
2
-
r
3,2 yields a solution with three free parameters a
2,1
,
a
3,2, and c
3 as follows: (21)
a
3,1
=
-
1
10
-
3
+
5
c
3
1500
a
2,1
c
3
2
-
3000
a
2,1
c
3
3
+
2000
a
2,1
c
3
4
-
250
a
2,1
c
3
-
30
a
3,2
+
50
a
3,2
c
3
+
9
c
3
-
45
c
3
2
+
80
c
3
3
-
50
c
3
4
,
b
1
=
10
c
3
2
-
8
c
3
+
1
24
-
3
+
5
c
3
c
3
,
b
2
=
5
2
-
9
c
3
+
10
c
3
2
2
c
3
-
1
24
3
-
10
c
3
+
10
c
3
2
-
3
+
5
c
3
,
b
3
=
-
c
3
-
1
24
c
3
3
-
10
c
3
+
10
c
3
2
,
c
2
=
-
3
+
5
c
3
5
2
c
3
-
1
,
d
1
=
10
c
3
2
-
8
c
3
+
1
12
-
3
+
5
c
3
c
3
,
d
2
=
25
1
-
4
c
3
+
4
c
3
2
2
c
3
-
1
12
3
-
10
c
3
+
10
c
3
2
-
3
+
5
c
3
,
d
3
=
1
12
c
3
3
-
10
c
3
+
10
c
3
2
,
r
3,2
=
3
-
10
c
3
+
10
c
3
2
c
3
2
c
3
-
1
2
-
3
+
5
c
3
.Substituting the above solution into τ
(
6
)
2, τ
′
(
6
)
2, and τ
g
(
6
)
2 gives us (22)
τ
6
2
=
1
3600
-
3
+
5
c
3
2
c
3
-
1
2
-
900000
a
2,1
c
3
6
+
18000000
a
2,1
2
c
3
6
+
33750
c
3
6
-
115500
c
3
5
-
54000000
a
2,1
2
c
3
5
+
2940000
a
2,1
c
3
5
+
67500000
a
2,1
2
c
3
4
+
163525
c
3
4
-
3990000
a
2,1
c
3
4
-
11460
c
3
+
2880000
a
2,1
c
3
3
-
122650
c
3
3
+
51445
c
3
2
+
16875000
a
2,1
2
c
3
2
-
1166250
a
2,1
c
3
2
+
251250
a
2,1
c
3
-
3375000
a
2,1
2
c
3
+
1062
+
281250
a
2,1
2
-
22500
a
2,1
-
45000000
a
2,1
2
c
3
3
1
/
2
.
τ
′
6
2
=
1
3600
-
3
+
5
c
3
c
3
3
-
10
c
3
+
10
c
3
2
2
c
3
-
1
291600
a
3,2
2
+
1124760000
a
2,1
c
3
5
+
20250000
a
2,1
2
c
3
2
-
16110000000
a
2,1
2
c
3
5
+
502200
a
3,2
c
3
2
-
378000000
a
2,1
2
c
3
3
+
3195000000
a
2,1
2
c
3
4
+
53685000000
a
2,1
2
c
3
6
-
3582480000
a
2,1
c
3
6
+
27000000
c
3
12
-
123588000000
a
2,1
2
c
3
7
+
7880280000
a
2,1
c
3
7
-
51676200
c
3
5
+
199116000000
a
2,1
2
c
3
8
-
12124200000
a
2,1
c
3
8
-
48600
a
3,2
c
3
+
12883200000
a
2,1
c
3
9
-
221760000000
a
2,1
2
c
3
9
+
83673
c
3
2
+
163440000000
a
2,1
2
c
3
10
-
9048000000
a
2,1
c
3
10
-
1428030
c
3
3
+
359320000
c
3
10
+
14400000000
a
2,1
2
c
3
12
-
1500000
c
3
7
a
3,2
-
72000000000
a
2,1
2
c
3
11
+
3792000000
a
2,1
c
3
11
-
146400000
c
3
11
+
11099385
c
3
4
+
159450100
c
3
6
-
340446000
c
3
7
+
509060500
c
3
8
-
526050000
c
3
9
-
720000000
a
2,1
c
3
12
+
5139000
c
3
4
a
3,2
-
6870000
c
3
5
a
3,2
-
1944000
a
3,2
2
c
3
+
4860000
a
3,2
2
c
3
2
-
2187000
c
3
3
a
3,2
+
4950000
a
3,2
c
3
6
-
5400000
a
3,2
2
c
3
3
+
2250000
a
3,2
2
c
3
4
-
1620000
a
2,1
c
3
2
-
233370000
a
2,1
c
3
4
+
28890000
a
2,1
c
3
3
1
/
2
,
τ
g
6
2
=
∑
i
=
1
n
1
τ
i
6
2
+
∑
i
=
1
n
1
′
τ
i
′
6
2
,Since the expression error norm (22) is so complicated and not suitable to write and minimizing (22) with respect to the free parameters a
2,1
,
a
3,2
,
c
3, and r
3,2, we get a
2,1
=
3
/
1000
,
a
3,2
=
1
/
25, c
3
=
18
/
25, and r
3,2
=
4059
/
15625. These values give us(23)
τ
6
2
=
2.2
×
1
0
-
4
,
τ
′
6
2
=
6.8
×
1
0
-
4
,
τ
g
6
2
=
7.2
×
1
0
-
4
.The coefficients of the new method can be given in Butcher tableau and denoted by STDRKN5(
3
) as seen in Table 5.Table 5
The STDRKN5(
3
) method.
0
0
0
3/11
3
/
1000
0
9
/
242
0
18/25
36221
/
1562500
1
/
25
0
-
9
/
15625
4059
/
15625
0
53
/
1296
121
/
1107
875
/
53136
53
/
648
1331
/
4428
3125
/
26568


## 3.1. Two-Stage STDRKN Method of Order Four

To derive the fourth-order STDRKN method, we use the algebraic order conditions up to order four in the equationsy and y
′, that is, (10)-(11) and (14)–(16), respectively. However, we get a system of equations which consists of 5 nonlinear equations with 6 unknown variables. Solving these equations with the simplifying assumption r
2,1
=
c
2
2
/
2, the resulting system has one free parameter a
2,1; this free parameter can be selected by minimizing error equations. The truncation error norms and global error of the fifth-order condition are defined as follows:(19)
τ
5
2
=
∑
i
=
1
n
1
τ
i
5
2
=
1
480
,
τ
′
5
2
=
∑
i
=
1
n
1
′
τ
i
′
5
2
=
1
120
2
+
1600
a
2,1
2
-
80
a
2,1
,
τ
g
5
2
=
∑
i
=
1
n
1
τ
i
5
2
+
∑
i
=
1
n
1
′
τ
i
′
5
2
=
1
480
33
+
25600
a
2,1
2
-
1280
a
2,1
.The error equations τ
g
5
2 have a minimum value at a
2,1
=
1
/
40 which produces(20)
τ
5
2
=
2.08
×
1
0
-
3
,
τ
′
5
2
=
8.3
×
1
0
-
3
,
τ
g
5
2
=
8.5
×
1
0
-
3
,where τ
(
5
) and τ
′
(
5
) are error terms of the fifth-order conditions for y
n and y
n
′, respectively. The coefficients of the new method can be given in Butcher tableau and denoted by STDRKN4(
2
) as seen in Table 4.Table 4
The STDRKN4(
2
) method.
0
0
0
1/2
1
/
40
0
1
/
8
0
1
/
12
1
/
12
1
/
6
1
/
3


## 3.2. Three-Stage STDRKN Method of Order Five

According to the order conditions in (10)–(12) and (14)–(17), we obtain system of equations that consists of 9 nonlinear equations with 14 unknowns variables that need to be solved. Solving the system simultaneously with simplifying assumption conditions r
2,1
=
c
2
2
/
2
,
r
3,1
=
c
3
2
/
2
-
r
3,2 yields a solution with three free parameters a
2,1
,
a
3,2, and c
3 as follows: (21)
a
3,1
=
-
1
10
-
3
+
5
c
3
1500
a
2,1
c
3
2
-
3000
a
2,1
c
3
3
+
2000
a
2,1
c
3
4
-
250
a
2,1
c
3
-
30
a
3,2
+
50
a
3,2
c
3
+
9
c
3
-
45
c
3
2
+
80
c
3
3
-
50
c
3
4
,
b
1
=
10
c
3
2
-
8
c
3
+
1
24
-
3
+
5
c
3
c
3
,
b
2
=
5
2
-
9
c
3
+
10
c
3
2
2
c
3
-
1
24
3
-
10
c
3
+
10
c
3
2
-
3
+
5
c
3
,
b
3
=
-
c
3
-
1
24
c
3
3
-
10
c
3
+
10
c
3
2
,
c
2
=
-
3
+
5
c
3
5
2
c
3
-
1
,
d
1
=
10
c
3
2
-
8
c
3
+
1
12
-
3
+
5
c
3
c
3
,
d
2
=
25
1
-
4
c
3
+
4
c
3
2
2
c
3
-
1
12
3
-
10
c
3
+
10
c
3
2
-
3
+
5
c
3
,
d
3
=
1
12
c
3
3
-
10
c
3
+
10
c
3
2
,
r
3,2
=
3
-
10
c
3
+
10
c
3
2
c
3
2
c
3
-
1
2
-
3
+
5
c
3
.Substituting the above solution into τ
(
6
)
2, τ
′
(
6
)
2, and τ
g
(
6
)
2 gives us (22)
τ
6
2
=
1
3600
-
3
+
5
c
3
2
c
3
-
1
2
-
900000
a
2,1
c
3
6
+
18000000
a
2,1
2
c
3
6
+
33750
c
3
6
-
115500
c
3
5
-
54000000
a
2,1
2
c
3
5
+
2940000
a
2,1
c
3
5
+
67500000
a
2,1
2
c
3
4
+
163525
c
3
4
-
3990000
a
2,1
c
3
4
-
11460
c
3
+
2880000
a
2,1
c
3
3
-
122650
c
3
3
+
51445
c
3
2
+
16875000
a
2,1
2
c
3
2
-
1166250
a
2,1
c
3
2
+
251250
a
2,1
c
3
-
3375000
a
2,1
2
c
3
+
1062
+
281250
a
2,1
2
-
22500
a
2,1
-
45000000
a
2,1
2
c
3
3
1
/
2
.
τ
′
6
2
=
1
3600
-
3
+
5
c
3
c
3
3
-
10
c
3
+
10
c
3
2
2
c
3
-
1
291600
a
3,2
2
+
1124760000
a
2,1
c
3
5
+
20250000
a
2,1
2
c
3
2
-
16110000000
a
2,1
2
c
3
5
+
502200
a
3,2
c
3
2
-
378000000
a
2,1
2
c
3
3
+
3195000000
a
2,1
2
c
3
4
+
53685000000
a
2,1
2
c
3
6
-
3582480000
a
2,1
c
3
6
+
27000000
c
3
12
-
123588000000
a
2,1
2
c
3
7
+
7880280000
a
2,1
c
3
7
-
51676200
c
3
5
+
199116000000
a
2,1
2
c
3
8
-
12124200000
a
2,1
c
3
8
-
48600
a
3,2
c
3
+
12883200000
a
2,1
c
3
9
-
221760000000
a
2,1
2
c
3
9
+
83673
c
3
2
+
163440000000
a
2,1
2
c
3
10
-
9048000000
a
2,1
c
3
10
-
1428030
c
3
3
+
359320000
c
3
10
+
14400000000
a
2,1
2
c
3
12
-
1500000
c
3
7
a
3,2
-
72000000000
a
2,1
2
c
3
11
+
3792000000
a
2,1
c
3
11
-
146400000
c
3
11
+
11099385
c
3
4
+
159450100
c
3
6
-
340446000
c
3
7
+
509060500
c
3
8
-
526050000
c
3
9
-
720000000
a
2,1
c
3
12
+
5139000
c
3
4
a
3,2
-
6870000
c
3
5
a
3,2
-
1944000
a
3,2
2
c
3
+
4860000
a
3,2
2
c
3
2
-
2187000
c
3
3
a
3,2
+
4950000
a
3,2
c
3
6
-
5400000
a
3,2
2
c
3
3
+
2250000
a
3,2
2
c
3
4
-
1620000
a
2,1
c
3
2
-
233370000
a
2,1
c
3
4
+
28890000
a
2,1
c
3
3
1
/
2
,
τ
g
6
2
=
∑
i
=
1
n
1
τ
i
6
2
+
∑
i
=
1
n
1
′
τ
i
′
6
2
,Since the expression error norm (22) is so complicated and not suitable to write and minimizing (22) with respect to the free parameters a
2,1
,
a
3,2
,
c
3, and r
3,2, we get a
2,1
=
3
/
1000
,
a
3,2
=
1
/
25, c
3
=
18
/
25, and r
3,2
=
4059
/
15625. These values give us(23)
τ
6
2
=
2.2
×
1
0
-
4
,
τ
′
6
2
=
6.8
×
1
0
-
4
,
τ
g
6
2
=
7.2
×
1
0
-
4
.The coefficients of the new method can be given in Butcher tableau and denoted by STDRKN5(
3
) as seen in Table 5.Table 5
The STDRKN5(
3
) method.
0
0
0
3/11
3
/
1000
0
9
/
242
0
18/25
36221
/
1562500
1
/
25
0
-
9
/
15625
4059
/
15625
0
53
/
1296
121
/
1107
875
/
53136
53
/
648
1331
/
4428
3125
/
26568


## 4. Stability of the STDRKN Methods

In this part, we study the linear stability of the STDRKN methods. We use the test problem (see [7, 14])(24)
y
′
′
=
-
ω
2
y
+
ε
y
′
.Applying STDRKN method (5) to test problem (24), we obtain (25)
y
n
+
1
h
y
n
+
1
′
=
M
v
,
z
y
n
h
y
n
′
,where (26)
M
v
,
z
=
m
1,1
m
1,2
m
2,1
m
2,2with(27)
m
1,1
=
1
-
v
2
2
+
z
2
-
v
2
b
T
N
1
-
1
e
-
v
2
z
z
2
-
v
2
N
1
-
1
R
N
2
-
1
A
N
1
-
1
e
N
3
-
v
2
z
b
T
z
2
-
v
2
N
1
-
1
R
N
2
-
1
-
N
2
-
1
N
4
,
m
1,2
=
1
+
z
2
+
z
2
-
v
2
b
T
N
1
-
1
e
-
v
2
z
z
2
-
v
2
N
1
-
1
R
N
2
-
1
A
N
1
-
1
e
N
3
-
v
2
z
b
T
z
2
-
v
2
N
1
-
1
R
N
2
-
1
-
N
2
-
1
N
4
,
m
2,1
=
-
v
2
+
z
2
-
v
2
d
T
N
1
-
1
e
-
v
2
z
z
2
-
v
2
N
1
-
1
R
N
2
-
1
A
N
1
-
1
e
N
3
-
v
2
z
d
T
z
2
-
v
2
N
1
-
1
R
N
2
-
1
-
N
2
-
1
N
4
,
m
2,2
=
1
+
z
+
z
2
-
v
2
d
T
N
1
-
1
e
-
v
2
z
z
2
-
v
2
N
1
-
1
R
N
2
-
1
A
N
1
-
1
e
N
3
-
v
2
z
d
T
z
2
-
v
2
N
1
-
1
R
N
2
-
1
-
N
2
-
1
N
4
,where N
1
=
I
-
(
z
2
-
v
2
)
R
(28)
N
2
=
I
+
v
2
z
A
+
v
2
z
z
2
-
v
2
A
N
1
-
1
R
,
N
3
=
-
c
1
v
2
1
+
c
1
z
⋯
⋯
⋯
⋯
-
c
s
v
2
1
+
c
s
z
,
N
4
=
1
-
1
2
c
1
2
z
c
1
+
1
2
c
1
2
z
⋯
⋯
⋯
⋯
1
-
1
2
c
s
2
z
c
s
+
1
2
c
s
2
z
,where v
=
ω
h and z
=
ε
h. The matrix M
(
v
,
z
) is called stability matrix, while the stability region of STDRKN method is defined by (29)
S
R
=
v
,
z
:
λ
i
M
<
1
,
i
=
1,2
.
λ
i are eigenvalues of M
(
v
,
z
)
. The stability regions for STDRKN4(
2
) and STDRKN5(
3
) are shown in Figures 1 and 2, respectively.Figure 1
The stability region for STDRKN4(
2
) method.Figure 2
The stability region for STDRKN5(
3
) method.


## 5. Numerical Experiments

In this section, we test the effectiveness of the new methods of orders four and five on the same problems for comparison. The numerical methods used for comparison are given as follows:(i)
STDRKN5(
3): new special explicit two-derivative RKN method of fifth order derived in this paper
(ii)
STDRKN4(
2): new special explicit two-derivative RKN method of fourth order derived in this paper
(iii)
TDRKN5(
3): three-stage fifth-order two-derivative RKN method derived in [7]
(iv)
TDRKN4(
2): two-stage fourth-order two-derivative RKN method derived in [7]
(v)
RKNG5(
6): the classical six-stage fifth-order RKN method which is the limit method of ARKNGV5 as the frequency matrix M
→
0 derived in [8]
(vi)
RKNG4: the classical four-stage fourth-order RKN method derived in [2].Problem 1.
Consider the linear homogeneous problem(30)
y
′
′
=
x
y
-
y
′
,
y
0
=
1
,
y
′
0
=
-
1
,
x
e
n
d
=
12 whose analytic solution is y
x
=
e
-
x
-
1
/
2
x
2.Problem 2.
Consider the linear nonhomogeneous problem(31)
y
′
′
=
-
y
′
+
cos
⁡
x
,
y
0
=
-
1
2
,
y
′
0
=
1
2
,
x
e
n
d
=
10 whose analytic solution is y
x
=
1
/
2
sin
⁡
x
-
cos
⁡
x
.Problem 3.
Consider the nonlinear problem(32)
y
′
′
=
1
40
10
-
y
y
′
,
y
0
=
1
,
y
′
0
=
19
80
,
x
e
n
d
=
10 whose analytic solution is y
(
x
)
=
20
/
1
+
19
e
-
x
/
4
.Problem 4.
Consider the famous Van der Pol equation (see [15]) (33)
y
′
′
=
-
y
+
δ
1
-
y
2
y
′ with the initial values(34)
y
0
=
2
+
1
96
δ
2
+
1033
552960
δ
4
+
1019689
55738368000
δ
6
,
y
′
0
=
0
.
This is a nonlinear equation. Here we takeδ
=
0.8
×
1
0
-
4. The problem is integrated on the interval [
0,100
] with stepsizes h
=
1
/
2
i
,
i
=
1,2
,
3,4. Since the exact solution of the problem is not available, when estimating the error of each method, we use RKNG4 method given in [2] as a reference numerical solution with a very small stepsize.Problem 5.
Consider the linear homogeneous system (see [7]) (35)
y
′
′
=
-
M
y
+
K
y
′
, where (36)
M
=
1
-
9
40
27
40
0
9
2
-
3
2
0
3
8
63
8
,
K
=
0
21
40
3
10
0
1
2
1
0
-
7
8
-
1
2 with the initial values (37)
y
0
=
-
1
2
,
1
6
,
-
1
6
T
,
y
′
0
=
1,2
,
1
T
.
The problem is integrated on the interval[
0,12
] with stepsizes h
=
0.2
/
2
i
,
i
=
0,1
,
2,3.
The exact solution of the problem is given by(38)
y
x
=
sin
⁡
x
-
1
2
cos
⁡
2
x
sin
⁡
2
x
+
1
2
cos
⁡
2
x
-
1
3
cos
⁡
3
x
1
3
sin
⁡
3
x
-
1
2
cos
⁡
2
x
+
1
3
cos
⁡
3
x
.Problem 6.
Consider the linear nonhomogeneous system (see [8]) (39)
y
′
′
=
-
M
y
+
12
ϵ
5
K
y
′
+
ϵ
2
L
x
,
M
=
13
-
12
-
12
13
,
K
=
3
2
-
2
-
3
,
L
=
36
5
sin
⁡
x
+
24
sin
⁡
5
x
-
24
5
sin
⁡
x
-
36
sin
⁡
5
xwith the initial values(40)
y
0
=
ϵ
,
ϵ
T
,
y
′
0
=
-
4,6
T
,
ϵ
=
1
0
-
3
.
The problem is integrated on the interval[
0,5
] with stepsizes h
=
0.1
/
2
i
,
i
=
1,2
,
3,4.
The exact solution is(41)
y
x
=
sin
⁡
x
-
sin
⁡
5
x
+
ϵ
cos
⁡
x
sin
⁡
x
+
sin
⁡
5
x
+
ϵ
cos
⁡
5
x
.Problem 7.
Consider the damped wave equation with periodic boundary conditions (see [15]) (42)
∂
u
2
∂
t
2
+
δ
∂
u
∂
t
=
∂
u
2
∂
x
2
+
f
u
,
-
1
<
x
<
1
,
t
>
0
,
u
-
1
,
t
=
u
1
,
t
.
A semidiscretization in the spatial variable by second-order symmetric differences leads to the following system of second-order ODEs in time:(43)
U
¨
+
K
U
=
F
U
,
U
˙
,
0
<
t
<
t
e
n
d
, where U
(
t
)
=
(
u
1
(
t
)
,
…
,
u
N
(
t
)
)
T with u
i
(
t
)
≈
u
(
x
i
,
t
)
,
i
=
1,2
,
…
,
N, and (44)
K
=
1
Δ
x
2
2
-
1
-
1
-
1
2
-
1
·
·
·
·
·
·
·
·
·
-
1
2
-
1
-
1
-
1
2 with Δ
x
=
2
/
N and x
i
=
-
1
+
i
Δ
x and F
(
U
,
U
˙
)
=
f
(
u
1
)
-
δ
u
˙
1
,
…
,
f
(
u
N
)
-
δ
u
˙
N. In this experiment, we take f
(
u
)
=
sin
⁡
u and δ
=
0.08 and the initial conditions as (45)
U
0
=
π
i
N
,
U
t
0
=
N
0.01
+
sin
⁡
2
π
i
N
i
=
1
N
.
We chooseN
=
40 and integrate the problem in the interval [
0,100
] with the stepsizes h
=
1
/
2
i
,
i
=
5,6
,
7,8. The reference numerical solution is obtained by method RKNG4 with a very small stepsize.


## 6. Conclusion

In this study, the special class of explicit two-derivative Runge-Kutta-Nyström methods of order up to five that involve onef-evaluation and minimal number of g-evaluations was derived. Figures 3–9 display the efficiency curves showing the common logarithm of the maximum global absolute error throughout the integration versus computational cost measured by time used by each method in the same computation machine. An advantage of the STDRKN methods over the general classical Runge-Kutta-Nyström methods and TDRKN methods is that they can reach higher order with fewer functions evaluations per step and also give us higher stage order than RKN. Some tested problems were performed. From Figures 3, 4, 5, 6, 7, 8, and 9, the numerical results showed that the new methods agreed very well with the existing methods in the literature and required less time compared to the existing methods.Figure 3
The time curves for all methods for Problem1.Figure 4
The time curves for all methods for Problem2.Figure 5
The time curves for all methods for Problem3.Figure 6
The time curves for all methods for Problem4.Figure 7
The time curves for all methods for Problem5.Figure 8
The time curves for all methods for Problem6.Figure 9
The time curves for all methods for Problem7.



---

*Source: 2393015-2018-03-20.xml*