# New Phase-Fitted and Amplification-Fitted Fourth-Order and Fifth-Order Runge-Kutta-Nyström Methods for Oscillatory Problems

Authors: K. W. Moo; N. Senu; F. Ismail; M. Suleiman
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

Two new Runge-Kutta-Nyström (RKN) methods are constructed for solving second-order differential equations with oscillatory solutions. These two new methods are constructed based on two existing RKN methods. Firstly, a three-stage fourth-order Garcia’s RKN method. Another method is Hairer’s RKN method of four-stage fifth-order. Both new derived methods have two variable coefficients with phase-lag of order infinity and zero amplification error (zero dissipative). Numerical tests are performed and the results show that the new methods are more accurate than the other methods in the literature.


---

## Body

## 1. Introduction

Throughout this paper, we are dealing with the initial value problems (IVP) related to second-order ODEs in the form:(1)y′′=f(t,y),y(t0)=y0,y′(t0)=y0′,
where f is independent of y′(t) explicitly and it is known that their solutions are periodic. This kind of problems often arise in many fields of applied sciences such as astronomy, quantum mechanics, physical chemistry, structural mechanics, and electronics.Over the last three decades, there are a number of numerical methods that have been derived by several authors based on different approaches such as minimal phase-lag, phase-fitted, and exponential-fitted for solving second-order oscillatory IVPs. See van der Houwen and Sommeijer [1], Senu et al. [2], van de Vyver [3], Vigo-Aguiar and Ramos [4], Ramos and Vigo-Aguiar [5], and Ismail [6]. The term phase-lag was first introduced by Brusa and Nigro [7] in 1980. In 1993, Simos derived a Runge-Kutta-Fehlberg method based on the idea of phase-lag of order infinity [8]. Recently, the idea of phase-lag of order infinity has been used to developed new numerical methods. Simos and Vigo-Aguiar presented in [9] a modified Runge-Kutta-Nyström method with phase-lag of order infinity for solving the Schrödinger equation and related problems. Papadopoulos et al. [10] presented a new phase-fitted RKN method for solving IVP with oscillating solutions.In this paper, we will combine the idea of phase-lag of order infinity and zero amplification error together. First, we will construct a three-stage phase-fitted and amplification-fitted RKN method which is based on Garcia's method of algebraic order four [11]. Then, we will construct a four-stage phase-fitted and amplification-fitted RKN method based on the coefficients of Hairer's method of algebraic order five [12].


## 2. Analysis Phase-Lag of the Methods

The generalm-stage RKN method for (1) is of the form
(2)yn=yn-1+hyn-1′+h2∑i=1m‍bif(tn-1+cih,Yi),yn′=yn-1′+h∑i=1mbi′‍f(tn-1+cih,Yi),
where
(3)Yi=yn-1+hciyn-1′+h2∑j=1i-1‍aijf(tn-1+cjh,Yj),i=1,…,m.In the investigation of phase-lag error of the method, we will use the test equation in the following:(4)y′′(t)=-v2y(t),v∈R.
By applying method (2) to test equation (4), it yields
(5)[ynhyn′]=Dn[y0hy0′],D=[A(z2)B(z2)A′(z2)B′(z2)],z=vh,
where A,B,A′, and B′ are polynomial in term of z2 and are totally determined by the parameters of method (2).It is given that the exact solution of (4) is
(6)y(tn)=σ1[exp(iz)]n+σ2[exp(-iz)]n,
where
(7)σ1,2=12[y0±(iy0′)v]orσ1,2=|σ|exp(±iχ).
By substituting (7) into (6), it yields
(8)y(tn)=2|σ|cos(χ+nz).
Then, we assume that the eigenvalues of D are ρ1, ρ2, and they will be called as the amplification factors of the RKN method. The consequent eigenvectors are [1,γ1]T, [1,γ2]T, where γi=A′/(ρi-B′), i=1,2. The numerical solution of (4) is
(9)yn=c1ρ1n+c2ρ2n,
where
(10)c1=-γ2y0-hy0′γ1-γ2,c2=γ1y0-hy0′γ1-γ2.
If ρ1 and ρ2 are complex conjugate, then c1,2=|c|exp(±iw) and ρ1,2=|ρ|exp(±ip¯). By substituting both into (9), we have
(11)yn=2|c||ρ|cos(w+np¯).
From the exact solution (8) and the numerical solution (11), we can forward to the following definition.Definition 1 (phase-lag, see [1]).
Apply RKN method (2) to test equation (4). Then we define the phase-lag Φ(z)=z-p¯. If Φ(z)=O(zq+1), then the RKN method is said to have phase-lag order q.Definition 2 (amplification error, see [1]).
Apply RKN method (2) to test equation (4). Then we define the quantity α(z)=1-|ρ| is the amplification error. If α(z)=O(zr+1), then the RKN method is said to have dissipation order r.Let us denote that(12)R(z2)=trace(D)=A(z2)+B′(z2),Q(z2)=det(D)=A(z2)B′(z2)-A′(z2)B(z2).
From Definition 1 it follows that
(13)Φ(z)=z-cos-1(R(z2)2Q(z2)),|ρ|=Q(z2).
From (13), we have the following definition for the case of phase-lag of order infinity.Definition 3 (phase-lag of order infinity, see [10]).
To obtain phase-lag of order infinity the relationΦ(z)=z-cos-1(R(z2)/2Q(z2))=0 must hold.When at a pointz, where α(z)=0, the method is said to have zero amplification error (zero-dissipative), and we have
(14)α(z)=1-Q(z2)=0⟶Q(z2)=1.
From the condition above, we have the following remark.Remark 4 (see [13]).
To achieve phase-lag of order infinity and zero amplification error, the relations below must hold(15)R(z2)-2cos(z)=0,Q(z2)-1=0.


## 3. Construction of the New Methods

In this section, we will present the construction of two new RKN methods. The first method is based on a three-stage RKN method with algebraic order four (see [11]). The second method is based on a four-stage RKN method with algebraic order five (see [12]).


### 3.1. New Fourth Order RKN Method

In this section, we want to derive an RKN method with phase-lag of order infinity and zero amplification error which is based on Garcia's RKN method of three-stage and algebraic order four as follows(16)01413245-412544125124411252641241633125264First of all, we have to compute the polynomialsA, B, A′, and B′ and substitute all of them into (12). Hence, we have R and Q in terms of parameters of RKN method. In order to construct the new method with phase-lag of order infinity and zero amplification error, we set a31 and a32 as free parameters but leave the rest of the coefficients. Then we substitute the value of the coefficients above into R and Q yielding
(17)R(z2)=-25a328448z6+(25a31264+75a32352+188)z4-z2+2,Q(z2)=(25a323168-51584+25a31792)z6+(-25a3296+788-25a3166)z4+1.
As mentioned in Remark 4, in order to achieve phase-lag of order infinity and zero amplification error, conditions (15) must be satisfied. Therefore, we substitute R and Q from (17) into (15) and solve the equations yielding
(18)a31=(-278784cos(z)+278784+14880z4-147840z2+8448z2cos(z)-534z6+5z8)×(50z4(600-76z2+z4))-1,a32=(176(z6-66z4-96z2cos(z)+672z2+1152cos(z)-1152)176(z6-66z4-96z2cos(z))×(25z4(600-76z2+z4))-1.For small value ofz, we use the Taylor series expansion as follows:
(19)a31=-4125+20975000z2-3196562500z4-2968923625000000z6-6207413543750000000z8-62556915132248125000000000z10+⋯,a32=44125-889375z2+4073281250z4-4511476562500z6+41287442968750000z8+82370983979492187500z10+⋯.
Hence, a new method is derived and we denote it by PFAFRKN4. This method has two variable coefficients a31 and a32 that depend on the product of the step-length h and the frequency of the method v.


### 3.2.  New Fifth Order RKN Method

In this section, we want to derive a new RKN method which is based on Hairer's RKN method of four-stage and algebraic order five as follows:(20)01515023-1277271310-235935143361003365433601433612533616233635336By using the same strategy, we can obtain a new fifth-order RKN method with phase-lag of order infinity and zero amplification error. In this case, we seta21 and a31 as free parameters and we found that R and Q are
(21)R(z2)=(-1720-a2124)z6+(112+9a3156+25a2184)z4-z2+2,Q(z2)=(-734560+3a21896-a31560)z8+(-19a212016+3a31112+1576)z6+(-196-9a3128-25a21336)z4+1.
Then we substitute (21) into (15) and solve the equations yielding
(22)a21=(495z8-2073600cos(z)+178920z4-11520z4cos(z)-13500z6+2073600-1209600z2-8z10+172800z2cos(z))×(50z4(-235z4+2460z2-9000+8z6))-1,a31=(-432000+1233z8-14460z6+37800z4+161280z2+432000cos(z)-19440z4cos(z)+54720z2cos(z)-38z10)×(27z4(-235z4+2460z2-9000+8z6))-1.
For small value of z, we use the Taylor series expansion as follows:
(23)a21=150-2625z2-14693937500z4-50910546875z6-3276761649687500000z8-616230539913302351562500000z10+⋯,a31=-127+1405z2+2891701000z4-13112150000z6-2648977841995000000z8-517048863711493231750000000z10+⋯.Thus, we have obtained a new method which has two variable coefficients,a21 and a31, that depend on the product of the step-size h and the frequency of the method v. We denote this method as PFAFRKN5. For linear problems, the frequency of the method is the same as the frequency of the exact solution. For nonlinear problems, however, the frequency of the method is usually different from the frequency of the exact solution but depends on the method itself (see [14]). For each specific product of vh, it helps to nullify the phase-lag error and the amplification error. For solving oscillatory problems, reducing the phase-lag error of the method is far more important than decreasing its algebraic error.


## 3.1. New Fourth Order RKN Method

In this section, we want to derive an RKN method with phase-lag of order infinity and zero amplification error which is based on Garcia's RKN method of three-stage and algebraic order four as follows(16)01413245-412544125124411252641241633125264First of all, we have to compute the polynomialsA, B, A′, and B′ and substitute all of them into (12). Hence, we have R and Q in terms of parameters of RKN method. In order to construct the new method with phase-lag of order infinity and zero amplification error, we set a31 and a32 as free parameters but leave the rest of the coefficients. Then we substitute the value of the coefficients above into R and Q yielding
(17)R(z2)=-25a328448z6+(25a31264+75a32352+188)z4-z2+2,Q(z2)=(25a323168-51584+25a31792)z6+(-25a3296+788-25a3166)z4+1.
As mentioned in Remark 4, in order to achieve phase-lag of order infinity and zero amplification error, conditions (15) must be satisfied. Therefore, we substitute R and Q from (17) into (15) and solve the equations yielding
(18)a31=(-278784cos(z)+278784+14880z4-147840z2+8448z2cos(z)-534z6+5z8)×(50z4(600-76z2+z4))-1,a32=(176(z6-66z4-96z2cos(z)+672z2+1152cos(z)-1152)176(z6-66z4-96z2cos(z))×(25z4(600-76z2+z4))-1.For small value ofz, we use the Taylor series expansion as follows:
(19)a31=-4125+20975000z2-3196562500z4-2968923625000000z6-6207413543750000000z8-62556915132248125000000000z10+⋯,a32=44125-889375z2+4073281250z4-4511476562500z6+41287442968750000z8+82370983979492187500z10+⋯.
Hence, a new method is derived and we denote it by PFAFRKN4. This method has two variable coefficients a31 and a32 that depend on the product of the step-length h and the frequency of the method v.


## 3.2.  New Fifth Order RKN Method

In this section, we want to derive a new RKN method which is based on Hairer's RKN method of four-stage and algebraic order five as follows:(20)01515023-1277271310-235935143361003365433601433612533616233635336By using the same strategy, we can obtain a new fifth-order RKN method with phase-lag of order infinity and zero amplification error. In this case, we seta21 and a31 as free parameters and we found that R and Q are
(21)R(z2)=(-1720-a2124)z6+(112+9a3156+25a2184)z4-z2+2,Q(z2)=(-734560+3a21896-a31560)z8+(-19a212016+3a31112+1576)z6+(-196-9a3128-25a21336)z4+1.
Then we substitute (21) into (15) and solve the equations yielding
(22)a21=(495z8-2073600cos(z)+178920z4-11520z4cos(z)-13500z6+2073600-1209600z2-8z10+172800z2cos(z))×(50z4(-235z4+2460z2-9000+8z6))-1,a31=(-432000+1233z8-14460z6+37800z4+161280z2+432000cos(z)-19440z4cos(z)+54720z2cos(z)-38z10)×(27z4(-235z4+2460z2-9000+8z6))-1.
For small value of z, we use the Taylor series expansion as follows:
(23)a21=150-2625z2-14693937500z4-50910546875z6-3276761649687500000z8-616230539913302351562500000z10+⋯,a31=-127+1405z2+2891701000z4-13112150000z6-2648977841995000000z8-517048863711493231750000000z10+⋯.Thus, we have obtained a new method which has two variable coefficients,a21 and a31, that depend on the product of the step-size h and the frequency of the method v. We denote this method as PFAFRKN5. For linear problems, the frequency of the method is the same as the frequency of the exact solution. For nonlinear problems, however, the frequency of the method is usually different from the frequency of the exact solution but depends on the method itself (see [14]). For each specific product of vh, it helps to nullify the phase-lag error and the amplification error. For solving oscillatory problems, reducing the phase-lag error of the method is far more important than decreasing its algebraic error.


## 4. Algebraic Order and Error Analysis

In this section, we will compute the local truncation error (LTE) of the new methods and verify the algebraic order of the methods.Firstly, we compute the Taylor expansions of the exact solutiony(xn+h), the first derivative of the exact solution y′(xn+h), the numerical solution yn+1, and the first derivative of the numerical solution yn+1′. Then, we compute the LTE of y and its first derivative y′ as follows:
(24)LTE=yn+1-y(xn+h),LTEder=yn+1′-y′(xn+h).The LTE andLTEder of method PFAFRKN4 are
(25)LTE=-h51440yn(5)+O(h6),LTEder=h55760(yn(6)+18v2yn(4))+O(h6).
From (25), we can verify that the algebraic order of PFAFRKN4 is 4 since all of the coefficients up to h4 vanished, whereas the LTE and LTEder of method PFAFRKN5 are
(26)LTE=h621600(yn(6)+12v2yn(4))+O(h7),LTEder=-h6108000(yn(7)+20v2yn(5))+O(h7).
From (26), method PFAFRKN5 is said to have algebraic order 5 since all of the coefficients up to h5 vanished. Table 1 shows a comparison of the properties of the methods derived.Table 1
Summary of the properties of the methods.
Method
q
r
∥
τ
(
p
+
1
)
∥
2
∥
τ
′
(
p
+
1
)
∥
2
DPC
DSC
S.I/P.I
SRKN3V
6
∞
2.750
(
-
3
)
9.007
(
-
4
)
2.126
(
-
4
)
—
(0, 7.57)
PFERKN4P
∞
5
—
—
—
—
—
OPTRKN4P
∞
∞
1
/
2160
—
—
—
—
RKN4G
4
5
1
/
1440
3.130
(
-
3
)
1
/
640
1
/
1440
(0, 8.77)
PFAFRKN4
∞
∞
1
/
1440
—
—
—
(0, 8.94)
OPTRKN5K1
∞
∞
—
1.598
(
-
4
)
—
—
—
OPTRKN5K2
∞
∞
—
1.598
(
-
4
)
—
—
—
RKN5H
6
5
1
/
21600
1.449
(
-
4
)
1
/
6720
1
/
3600
—
PFAFRKN5
∞
∞
—
—
—
—
(0, 8.01)
Note:2.750(-3)=2.750×10-3.
∥
τ
(
p
+
1
)
∥
2 is principal local truncation error constant for y.
∥
τ
′
(
p
+
1
)
∥
2 is principal local truncation error constant for y′.
DPC is dispersion constant.
DSC is dissipation constant.
P.I is periodicity interval.
S.I is stability interval.


## 5. Problems Tested and Numerical Results

In this section, we will apply the new methods to some second-order differential equation problems. The following explicit RKN methods are selected for the numerical comparisons.Fourth-order:(i)
PFAFRKN4: the new derived fourth-order RKN method;(ii)
RKN4G: the three-stage fourth-order RKN method derived by García et al. [11];(iii)
SRKN3V: the third-order symplectic RKN method with minimal phase-lag derived by van de Vyver [3];(iv)
PFERKN4P: the fourth-order RKN method with phase-lag of order infinity derived by Papadopoulos et al. [10];(v)
OPTRKN4P: the fourth-order optimized RKN method derived by Papadopoulos and Simos [13].Fifth-order:(i)
PFAFRKN5: the new derived fifth-order RKN method;(ii)
RKN5H: the four-stage fifth-order RKN method derived by Hairer et al. [12];(iii)
OPTRKN5K1: the fifth-order optimized RKN method with zero first derivative of phase-lag derived by Kosti et al. [15];(iv)
OPTRKN5K2: the fifth-order optimized RKN method with zero first derivative of amplification factor derived by Kosti et al. [16].Those methods are categorized into two categories according to algebraic order of each method for comparison purposes. The accuracy criteria taken is calculating thelog10 of the maximum absolute error,
(27)Accuracy=log10max∥y(tn)-yn∥,
where tn=t0+nh, n=1,2,…,(T-t0)/h. We test the problems in the following for the interval of [0,4000] with the step-size h=1/2N, N≥1.Problem 1 (non-linear, two body problem).
(28)
y
1
′
′
=
-
y
1
r
3
,
y
1
(
0
)
=
1
-
e
,
y
1
′
(
0
)
=
0
,
y
2
′
′
=
-
y
2
r
3
,
y
2
(
0
)
=
0
,
y
2
′
(
0
)
=
1
+
e
1
-
e
,
where r=y12+y22 and e is the eccentricity.
Exact solution:y1(t)=cos(u)-e, y2(t)=1-e2sin(u), where u can be found by solving the equation u-esin(u)-t=0.
Estimated frequency:v=r-3/2.
Source: Kosti et al. [15].Problem 2 (homogeneous).
(29)
d
2
y
(
t
)
d
t
2
=
-
64
y
(
t
)
,
y
(
0
)
=
1
,
y
′
(
0
)
=
-
2
.
Exact solution:y(t)=-(1/4)sin(8t)+cos(8t).
Estimated frequency:v=8.Problem 3 (inhomogeneous).
(30)
d
2
y
(
t
)
d
t
2
=
-
u
2
y
(
t
)
+
(
u
2
-
1
)
sin
(
t
)
,
y
(
0
)
=
1
,
y
′
(
0
)
=
u
+
1
.
This case is usingu=10.
Exact solution:y(t)=cos(ut)+sin(ut)+sin(t).
Estimated frequency:v=10.
Source: van der Houwen and Sommeijer [1].Problem 4 (inhomogeneous).
(31)
d
2
y
(
t
)
d
t
2
=
-
y
(
t
)
+
t
,
y
(
0
)
=
1
,
y
′
(
0
)
=
2
.Exact solution:y(t)=sin(t)+cos(t)+t.Estimated frequency:v=1.Source: Allen, Jr. and Wing [17].Problem 5 (inhomogeneous system).
(32)
d
2
y
1
(
t
)
d
t
2
=
-
u
2
y
1
(
t
)
+
u
2
f
(
t
)
+
f
′′
(
t
)
,
y
1
(
0
)
=
a
+
f
(
0
)
,
y
1
′
(
0
)
=
f
′
(
0
)
,
d
2
y
2
(
t
)
d
t
2
=
-
u
2
y
2
(
t
)
+
u
2
f
(
t
)
+
f
′′
(
t
)
,
y
2
(
0
)
=
f
(
0
)
,
y
2
′
(
0
)
=
u
a
+
f
′
(
0
)
.
Exact solution:y1(t)=acos(ut)+f(t), y2(t)=asin(ut)+f(t),f(t)=e-10t, and parameters u and a are 20 and 0.1, respectively.
Estimated frequency:v=20.
Source: Lambert and Watson [18].The numerical results are plotted in Figures1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, and 16. Figures 1–4 and Figures 9–12 display the efficiency curves where the common logarithm of the maximum global absolute error throughout the integration versus computational cost measured by number of function evaluations. Figures 5–8 and Figures 13–16 display the efficiency curves where 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. Figures 1–8 show the efficiency curves of the methods in the first category that consist of PFAFRKN4, RKN4G, SRKN3V, PFERKN4P, and OPTRKN4P, whereas Figures 9–16 show the efficiency curves of the fifth-order methods consist of PFAFRKN5, RKN5H, OPTRKN5K1, and OPTRKN5K2. These efficiency curves display a clear comparison among those methods.Figure 1
Efficiency curves for all methods for Problem1 with e=0 and h=1/2i, i=3,…,7.Figure 2
Efficiency curves for all methods for Problem1 with e=0.1 and h=1/2i, i=3,…,7.Figure 3
Efficiency curves for all methods for Problem1 with e=0.2 and h=1/2i, i=3,…,7.Figure 4
Efficiency curves for all methods for Problem1 with e=0.3 and h=1/2i, i=3,…,7.Figure 5
Efficiency curves for all methods for Problem2 with h=1/2i, i=4,…,8.Figure 6
Efficiency curves for all methods for Problem3 with h=1/2i, i=4,…,8.Figure 7
Efficiency curves for all methods for Problem4 with h=1/2i, i=1,…,5.Figure 8
Efficiency curves for all methods for Problem5 with tend=4000 and h=1/2i, i=5,…,9.Figure 9
Efficiency curves for all methods for Problem1 with e=0 and h=1/2i, i=3,…,7.Figure 10
Efficiency curves for all methods for Problem1 with e=0.1 and h=1/2i, i=3,…,7.Figure 11
Efficiency curves for all methods for Problem1 with e=0.2 and h=1/2i, i=3,…,7.Figure 12
Efficiency curves for all methods for Problem1 with e=0.3 and h=1/2i, i=3,…,7.Figure 13
Efficiency curves for all methods of order five for Problem2 with h=1/2i, i=4,…,8.Figure 14
Efficiency curves for all methods of order five for Problem3 with h=1/2i, i=4,…,8.Figure 15
Efficiency curves for all methods of order five for Problem4 with h=1/2i, i=1,…,5.Figure 16
Efficiency curves for all methods of order five for Problem5 with h=1/2i, i=5,…,9.From Figures1–8, we found that the new method PFAFRKN4 is the most accurate and efficient for solving Problems 1–5, followed by OPTRKN4P, SRKN3V, PFERKN4P, and the corresponding original method RKN4G, whereas for the fifth-order methods, from Figures 9–16, we observed that the new method PFAFRKN5 is the most efficient method, followed by optimized RKN methods OPTRKN5K1 and OPTRKN5K2 and the corresponding original method RKN5H. In addition, we can observe that the new methods are having the same computational cost with the corresponding original methods.


## 6. Conclusion

In this paper, we have derived two new phase-fitted and amplification-fitted RKN methods for solving second-order IVPs which are oscillatory in nature. First method is based on Garcia's fourth algebraic order RKN method and the second method is based on Hairer's four-stage fifth algebraic order RKN method. Numerical results show that both methods are more accurate and efficient for solving second-order differential equations with oscillating solutions.



---

*Source: 939367-2013-10-09.xml*