# A Zero-Dissipative Phase-Fitted Fourth Order Diagonally Implicit Runge-Kutta-Nyström Method for Solving Oscillatory Problems

Authors: K. W. Moo; N. Senu; F. Ismail; M. Suleiman
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A new diagonally implicit Runge-Kutta-Nyström (DIRKN) method is constructed for solving second order differential equations with oscillatory solutions. The method is originally based on existing DIRKN method derived by Senu et al. which is three-stage and fourth algebraic order. The new derived method has a variable coefficient with phase-lag of order infinity. The numerical experiments are carried out and the results show the efficiency and accuracy of the new method in comparison with the other DIRKN methods in the literature.


---

## Body

## 1. Introduction

In this paper, we are dealing with the initial value problems (IVPs) related to second order ordinary differential equations (ODEs) in the form:(1)y′′=f(t,y),y(t0)=y0,y′(t0)=y0′,
where y′(t) does not appear explicitly and know that their solutions are periodic. These problems arise in many fields of applied sciences such as astronomy, quantum mechanics, physical chemistry, structural mechanics, and electronics.For solving this kind of problems, one of the numerical methods is Runge-Kutta-Nyström (RKN) method. RKN method is frequently used due to its computational advantage (see [1]). RKN method consists of two main classes; they are explicit method and implicit method. Generally, it is quite difficult to handle oscillatory problems of the form (1) with classical RKN methods. The term phase-lag was first introduced by Brusa and Nigro [2]. Phase-lag is the angle between the analytical solution and the numerical solution. For solving oscillatory problems, phase-lag is an important property. In [3], van der Houwen and Sommeijer implemented the phase-lag theory to Runge-Kutta (RK) and RKN methods. They presented a few explicit RK and RKN methods with reduced phase errors. Based on the minimal phase-lag theory, Senu et al. constructed a zero-dissipative RKN method in [4]. In [5], Simos derived a Runge-Kutta-Fehlberg method with phase-lag of order infinity. In [6], Papadopoulos et al. extended the idea of phase-lag of order infinity to RKN method and presented a phase-fitted RKN method. Based on idea of phase-lag of order infinity, Papadopoulos and Simos introduced a new methodology to construct optimized RKN methods in [7]. Since then, Kosti et al. constructed two optimized RKN methods with fifth algebraic order in [8, 9]. Notice that all methods mentioned above are explicit methods.In [10], Sommeijer presented two DIRKN methods with nonempty interval of periodicity. van der Houwen and Sommeijer extended the phase-lag theory to DIRKN methods and presented a few DIRKN methods with high phase-lag order for solving oscillatory problems in [11]. Sharp et al. also presented a few two-stage and three-stage DIRKN methods with high phase-lag order in [12]. Senu et al. extended the DIRKN methods to phase-lag order up to order eight and higher dissipative order in [13, 14]. However, there is no DIRKN method with phase-lag of order infinity being done yet. This motivates us to develop the DIRKN method with phase-lag of order infinity. In this paper, we will construct a three-stage phase-fitted DIRKN method which is based on a three-stage method of algebraic order four derived by Senu et al. [15].


## 2. Phase Properties of RKN Method

The general form of am-stage implicit RKN method for (1) is given by
(2)yn=yn-1+hyn-1′+h2∑i=1m‍bif(tn-1+cih,Yi),yn′=yn-1′+h∑i=1m‍bi′f(tn-1+cih,Yi),
where
(3)Yi=yn-1+hciyn-1′+h2∑j=1maijf(tn-1+cjh,Yj),WWWWWWWWWWWWi=1,…,m.The DIRKN method above can be expressed nicely in a Butcher table as shown as follows:(4)c1κc2a21κ⋮⋮⋮⋱cmam,1am,1⋯κb1b2⋯bmb1′b2′⋯bm′.For diagonally implicit RKN methods, the diagonal elements are equal. We denote the diagonal elements asκ so that a11=a22=⋯=amm=κ.The phase-lag error of method (2) is investigated by using the homogeneous test equation in the following:
(5)y′′(t)=-v2y(t),v∈R.Applying method (2) to the test equation (5) yields
(6)[ynhyn′]=Dn[y0hy0′],D=[A(z2)B(z2)A′(z2)B′(z2)],z=vh,
where D is the stability matrix of the RKN method and A,  B,  A′, and B′ are polynomials in terms of z2 and totally determined by the parameters of method (2). The characteristic equation of D can be written as
(7)ξ2-trace(D)ξ+det⁡(D)=0
which is the stability polynomial of the RKN method.It is given that the exact solution of (5) is
(8)y(tn)=σ1[exp⁡⁡(iz)]n+σ2[exp⁡⁡(-iz)]n,
where
(9)σ1,2=12[y0±(iy0′)v],orσ1,2=|σ|exp⁡(±iχ).
Substituting (9) into (8) yields
(10)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 (5) is
(11)yn=c1ρ1n+c2ρ2n,
where
(12)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 (10), we have
(13)yn=2|c||ρ|cos⁡(w+np¯).Hence, we have the exact solution (10) and the numerical solution (13) of (5) in the similar form. From (10) and (13) we have the following definition.Definition 1 (phase-lag and amplification error; see [3]).
Apply the RKN method (2) to the test equation (5). Then we define the phase-lag as Φ(z)=z-p¯. If Φ(z)=O(zq+1), then the RKN method is said to have phase-lag order q. In addition, the quantity α(z)=1-|ρ| is called amplification error. If α(z)=O(zr+1), then the RKN method is said to have dissipation order r.Then, we denote that(14)R(z2)=trace(D),Q(z2)=det⁡(D).
For DIRKN method, let us denote R and Q in the following form:
(15)R(z2)=2+α1z2+⋯+αmz2m(1+κz2)m,Q(z2)=1+β1z2+⋯+βmz2m(1+κz2)m.
From Definition 1 it follows that
(16)Φ(z)=z-cos⁡-1(R(z2)2Q(z2)),|ρ|=Q(z2).
From (12), it leads us to the definition of phase-lag of order infinity.Definition 2 (phase-lag of order infinity; see [6]).
To obtain phase-lag of order infinity the following relation must hold(17)Φ(z)=z-cos⁡-1(R(z2)2Q(z2))=0.In addition, whenα(z)=0 at a point z, the method is said to have zero amplification error (zero-dissipative). Hence, we have
(18)α(z)=1-Q(z2)=0.Let us denoteξ1,2 as the roots of (7); then we have the following definitions.Definition 3 (absolute stability interval; see [16]).
An interval(0,Ha) is called the interval of absolute stability of the method (2) if, for all H∈(0,Ha),  ξ1,2<1.Definition 4 (periodicity interval; see [10]).
An interval(0,Hp) is called the interval of periodicity of method (2) if, for all H∈(0,Hp), ξ1,2 are complex conjugate and of modulus one.


## 3. Construction of the New Method

In this section, we will present the construction of a zero-dissipative DIRKN method with phase-lag of order infinity. The method is originally based on a three-stage DIRKN method with algebraic order four (see Senu et al. [15]) as shown as follows:(19)12-3616-31212-36016-31212-3603616-312014+31214-31201212.Since we want to derive a new method with phase-lag order of order infinity, there must be a variable coefficient. Therefore, we seta31 as free parameter first but let the rest of the coefficients remain the same. Then, we have to compute the stability matrix D which is determined by the coefficients in the Butcher table. From matrix D, we can compute R and Q in form of (11):
(20)R(z2)=(-288+2(11-36a31+12a313-83)z4WWWiW+2(24+243)z2)WWwW×(-144+(43-7)z4+(243-48)z2)-1,Q(z2)=-144+(-7+24a313+43)z4+(243-48)z2-144+(43-7)z4+(243-48)z2.Now, we haveR and Q in terms of z2 and a31. From Definition 2, condition (13) must be satisfied in order to have phase-lag of order infinity. Therefore, we substitute R and Q from (15) into (13) and solve the equation for a31 which yields
(21)a31=[1368576+7879683-3248643z2w+624z63+720z6-559872z2-18z8w+1866240(cos⁡(z))2-6(cos⁡(z))2z8w+5184(cos⁡(z))2z4+165888(cos⁡(z))2z2w+1078272(cos⁡(z))23+26z83w+3456(cos⁡(z))2z43+96(cos⁡(z))2z63w+4(cos⁡(z))2z83+96768(cos⁡(z))2z23w+(-9cos⁡(z)z4-72cos⁡(z)z2-144cos⁡(z)3www-432cos⁡(z)+5cos⁡(z)z43wiw+24cos⁡(z)3z2)T]U,
where
(22)T=30025728+173352963-47969283z2WW+8544z63+13824z43+14880z6WW+24192z4-8308224z2+224z8WW+15012864(cos⁡(z))2+4(cos⁡(z))2z8WW+672(cos⁡(z))2z6+44928(cos⁡(z))2z4WW+124z83+1340928(cos⁡(z))2z2WW+8667648(cos⁡(z))23+25920(cos⁡(z))2z43WW+384(cos⁡(z))2z63+2(cos⁡(z))2z83WW+774144(cos⁡(z))2z23,U=-1144z4(243z2+5763+z4+48z2+1008).For small value ofz, we use the Taylor series expansion of a31 as follows:
(23)a31=V1z2+V2z4+V3z6+V4z8+V5z10+⋯,
where
(24)V1=847066076+4890538273360(2093+362)3(43+7),V2=77050517871906+44485137234545315120(2093+362)4(43+7)2,V3=232556319057343003+134266453409505479325200(2093+362)5(43+7)3,V4=(4869865288111386704234053WiW+2811618035008323543967243)WW×(29937600(2093+362)6(43+7)4)-1,V5=(3129408889209269105724330330613WiW+1806765064589379268561365285913)WW×(10897286400(2093+362)7(43+7)5)-1.Hence, a new method is derived and we denote it by PFDIRKN4. This method has one variable coefficienta31 that depends on the product of the step-length h and the frequency v. For each specific product of vh, it helps to nullify the phase-lag error. Hence, the accuracy of the method always depends on the nature of the problem and the properties of the method. For solving oscillatory problems, reducing the phase-lag error of the method is far more important than decreasing its algebraic error. The behavior of a31 value for different value of z is illustrated in Figure 1. For small value of z, α(z) tends to zero; therefore the new method is zero-dissipative too. The behavior of ∥τ(p+1)∥2 and ∥τ′(p+1)∥2 values of PFDIRKN4 for different value of z is illustrated in Figure 2. Table 1 shows a comparison of the properties of the method 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
PFDIRKN4
∞
∞
—
—
—
—
(0, 4.02)
DIRKN4Senu
4
∞
2.05
(
-
3
)
1.48
(
-
3
)
2.23
(
-
4
)
—
(0, 8.20)
DIRKN4Senu2
6
5
1.88
(
-
3
)
1.70
(
-
3
)
4.71
(
-
5
)
1.19
(
-
4
)
(0, 8.10)
DIRKN4Franco
4
5
1.04
(
-
3
)
2.57
(
-
3
)
1.74
(
-
3
)
7.76
(
-
3
)
(0,∞)
DIRKN4Raed
4
5
3.13
(
-
2
)
1.71
(
-
2
)
2.33
(
-
2
)
1.80
(
-
2
)
—
DIRKN4Sharp
6
5
1.85
(
-
3
)
6.26
(
-
4
)
2.82
(
-
2
)
1.03
(
-
2
)
(0, 12.00)
DIRKN4Som
4
5
6.35
(
-
4
)
1.59
(
-
4
)
3.18
(
-
4
)
1.43
(
-
4
)
(0, 9.51)
Note:∥τ(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.
1.23
(
-
3
)
=
1.23
×
1
0
-
3.Figure 1
Graph fora31 against z value.Figure 2
Graph for-log⁡10∥τ(p+1)∥2 and -log⁡10∥τ′(p+1)∥2 against z value.


## 4. Problems Tested and Numerical Results

In this section, we will apply the new method to some second order differential equation problems. The following DIRKN methods are selected for the numerical comparisons.(i)
PFDIRKN4: the new derived fourth order zero-dissipative phase-fitted DIRKN method.(ii)
DIRKN4Senu: the three-stage fourth order DIRKN method derived by Senu et al. [15].(iii)
DIRKN4Senu2: the three-stage fourth order DIRKN method with phase-lag order six derived by Senu et al. [15].(iv)
DIRKN4Franco: the A-stable fourth order DIRKN method derived by Franco and Gómez [16].(v)
DIRKN4Raed: the fourth order DIRKN method derived by Al-Khasawneh et al. [17].(vi)
DIRKN4Sharp: the three-stage fourth order with phase-lag of order six DIRKN method by Sharp et al. [12].(vii)
DIRKN4Som: the fourth order with phase-lag of order four by Sommeijer [10].The accuracy criteria taken are calculating thelog⁡10 of the maximum global absolute error:
(25)Accuracy=log⁡10max⁡∥y(tn)-yn∥,
where tn=t0+nh,n=1,2,…,(T-t0)/h. We test the problems in the following for each method with the step-size h=1/2N,  N≥1.Problem 1 (homogeneous).
Consider(26)d2y(t)dt2=-64y(t),y(0)=1,y′(0)=-2.Exact solution:y(t)=-(1/4)sin(8t)+cos⁡(8t).Estimated frequency:v=8.0.Problem 2 (inhomogeneous).
Consider(27)d2y(t)dt2=-u2y(t)+(u2-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.0.Source: van der Houwen and Sommeijer [3].Problem 3 (an almost periodic orbital problem).
Consider(28)d2y1(t)dt2+y1(t)=0.001cos⁡(t),y1(0)=1,y1′(0)=0,d2y2(t)dt2+y2(t)=0.001sin(t),y2(0)=0,y2′(0)=0.9995.Exact solution:y1(t)=cos⁡(t)+0.0005tsin(t),  y2(t)=sin(t)-0.0005tcos⁡(t).Estimated frequency:v=1.0.Source: Stiefel and Bettis [18].Problem 4 (two-body problem).
Consider(29)y1′′=-y1(y12+y22)3,y1(0)=1,y1′(0)=0,y2′′=-y2(y12+y22)3,y2(0)=0,y2′(0)=1.Exact solution:y1(t)=cos⁡(t),  y2(t)=sin(t).Estimated frequency:v=1.0.Source: Dormand et al. [19].Problem 5 (inhomogeneous system).
Consider(30)d2y1(t)dt2=-u2y1(t)+u2f(t)+f′′(t),y1(0)=a+f(0),y1′(0)=f′(0),d2y2(t)dt2=-u2y2(t)+u2f(t)+f′′(t),y2(0)=f(0),y2′(0)=ua+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.0.Source: Lambert and Watson [20].Problem 6 (Duffing problem).
Consider(31)y′′+y+y3=Fcos⁡(Ωt),y(0)=0.200426728067,y′(0)=0,
where F=0.002 and Ω=1.01.Exact solution:y(t)=∑i=04θ2i+1cos⁡[(2i+1)Ωt], where θ1=0.200179477536, θ3=0.246946143×10-3, θ5=0.304014×10-6, θ7=0.374×10-9, and θ9<10-12.Estimated frequency:v=0.9.Source: van de Vyver [21].Problem 7 (linear Strehmel-Weiner problem).
Consider(32)y′′(t)=(-20.20-9.67989.6-10000-6004.2-9.60-5.8)y(t)WWW+(150cos⁡(10t)75cos⁡(10t)75cos⁡(10t)),y(0)=(12-2),y′(0)=(000).Exact solution:(33)y(t)=(cos⁡(t)+2cos⁡(5t)-2cos⁡(10t)2cos⁡(t)+cos⁡(5t)-cos⁡(10t)-2cos⁡(t)+cos⁡(5t)-cos⁡(10t)).Estimated frequency:v=5.0.Source: Cong [22].The numerical results are plotted in Figures3, 4, 5, 6, 7, 8, and 9. Figures 3–9 display the efficiency curves where the accuracy versus the computational cost measured by the time used by each method in a same computation machine. Figures 3–9 show the efficiency curves of the methods that consist of PFDIRKN4, DIRKN4Senu, DIRKN4Senu2, DIRKN4Franco, DIRKN4Raed, DIRKN4Sharp, and DIRKN4Som. These efficiency curves display a clear comparison among those methods.Figure 3
The efficiency curves for all methods for Problem1 with tend=4000 and h=1/2i,  i=5,…,9.Figure 4
The efficiency curves for all methods for Problem2 with tend=4000 and h=1/2i,  i=5,…,9.Figure 5
The efficiency curves for all methods for Problem3 with tend=4000 and h=1/2i,  i=2,…,6.Figure 6
The efficiency curves for all methods for Problem4 with tend=4000 and h=1/2i,  i=3,…,7.Figure 7
The efficiency curves for all methods for Problem5 with tend=4000 and h=1/2i,  i=6,…,10.Figure 8
The efficiency curves for all methods for Problem6 with tend=4000 and h=1/2i,  i=2,…,6.Figure 9
The efficiency curves for all methods for Problem7 with tend=1000 and h=1/2i,  i=8,…,12.From Figures3 to 9, we found that the new method PFDIRKN4 is the most accurate and efficient for solving Problems 1–7, followed by DIRKN4Senu2, DIRKN4Senu, DIRKN4Som, DIRKN4Sharp, DIRKN4Franco, and DIRKN4Raed. From Figures 3 to 9, we can conclude that high phase-lag order of the method is very important for solving oscillating problems. Since both methods PFDIRKN4 and DIRKN4Senu2 are far more accurate than other methods which are having higher phase-lag order. In an addition, we can observe that the new method is having the same computational cost with the corresponding original method but higher accuracy.


## 5. Conclusion

In this paper, we have derived a new DIRKN method for solving second order IVPs which are oscillatory in nature. The new method is derived based on a fourth algebraic order DIRKN method [15] with zero dissipation. The new method has a variable coefficient that will help to nullify the phase-lag error. Numerical results show that the new method is more accurate and efficient for solving second order differential equations with oscillating solutions in nature.



---

*Source: 985120-2014-05-25.xml*