# New 4(3) Pairs Diagonally Implicit Runge-Kutta-Nyström Method for Periodic IVPs

Authors: Norazak Senu; Mohamed Suleiman; Fudziah Ismail; Norihan Md Arifin
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

New 4(3) pairs Diagonally Implicit Runge-Kutta-Nyström (DIRKN) methods with reduced phase-lag are developed for the integration of initial value problems for second-order ordinary differential equations possessing oscillating solutions. Two DIRKN pairs which are three- and four-stage with high order of dispersion embedded with the third-order formula for the estimation of the local truncation error. These new methods are more efficient when compared with current methods of similar type and with the L-stable Runge-Kutta pair derived by Butcher and Chen (2000) for the numerical integration of second-order differential equations with periodic solutions.


---

## Body

## 1. Introduction

In many scientific areas of engineering and applied sciences such as celestial mechanics, quantum mechanics, elastodynamics, theoretical physics and chemistry, and electronics, oscillatory problems of second-order ordinary differential equations (ODEs) can be found. An oscillatory problems of second-order ODEs have the following form:(1.1)y′′=f(t,y),y(t0)=y0,y′(t0)=y0′.Anm-stage Runge-Kutta-Nyström (RKN) method for the numerical integration of the IVP is given by
(1.2)yn+1=yn+hyn′+h2∑i=1m‍bif(tn+cih,Yi)(1.3)yn+1′=yn′+h∑i=1m‍bi′f(tn+cih,Yi),
where
(1.4)Yi=yn+cihyn′+h2∑j=1m‍aijf(tn+cih,Yj).The RKN parametersaij,bj,bj′, and cj are assumed to be real and m is the number of stages of the method. Introduce the m-dimensional vectors c,b, and b′ and m×m matrix A, where c=[c1,c2,…,cm]T,b=[b1,b2,…,bm]T,b′=[b1′,b2′,…,bm′]T, and A=[aij], respectively. The latter contains the class of Diagonally Implicit RKN (DIRKN) methods for which all the entries in the diagonal of A are equal. An embedded r(s) pair of DIRKN methods is based on the method (c,A,b,b′) of order r and the other DIRKN method (c,A,b^,b^′) of order s(s<r) and can be expressed in Butcher notation by the table of coefficients (see Table 1).Table 1
m-stage DIRKN pair.
c
1
λ
^
c
2
a
21
λ
^
c
3
a
31
a
32
λ
^
⋮
⋮
⋮
⋮
λ
^
c
m
a
m
,
1
a
m
,
2
…
a
m
,
m
-
1
λ
^
b
1
b
2
…
b
m
-
1
b
m
b
1
'
b
2
'
…
b
m
-
1
'
b
m
'
b
^
1
b
^
2
…
b
^
m
-
1
b
^
m
b
^
1
'
b
^
2
'
…
b
^
m
-
1
'
b
^
m
'Several authors in their papers have developed numerical methods for this class of problems, for example, van der Houwen and Sommeijer [1, 2], Sideridis and Simos [3], García et al. [4], and Senu et al. [5]. Next, Van de Vyver [6] and Senu et al. [7] obtained explicit RKN method with minimal phase lag. Franco [8] have developed explicit hybrid method for periodic IVPs. For implicit RKN methods, see for example, van der Houwen and Sommeijer [2], Sharp et al. [9]. Imoni et al. [10] and Al-Khasawneh et al. [11] have developed general purpose of DIRKN methods with variable stepsize which is not related to dispersion property. Another classes of numerical methods for solving (1.1) are exponentially fitted or phase fitted in which the period or frequency is known in advance (see e.g., [12–21]).Most of the numerical methods developed for solving (1.1) are in constant stepsize (see [1–3, 6, 22, 23]). In this paper, the development of efficient DIRKN methods with reduced phase-lag in variable stepsize is studied. The strategies introduced in Dormand et al. [24] and Simos [18] for finding the optimized pair is used and also new implementation code is discussed in this paper.In this paper, dispersion relations are developed and used together with algebraic conditions to be solved explicitly. In the following section, the construction of the new 4(3) pairs of DIRKN method is described. Its coefficients are given using the Butcher tableau notation. Finally, numerical tests on second-order differential equation problems possessing an oscillatory solutions are performed.


## 2. Analysis of Phase-Lag and Stability

In this section, we will discuss the analysis of phase-lag for RKN method. The first analysis of phase-lag was carried out by Bursa and Nigro [25], then followed by Gladwell and Thomas [26] for the linear multistep method, and for explicit and implicit Runge-Kutta(-Nyström) methods by van der Houwen and Sommeijer [1, 2].The phase-lag analysis of the method (1.2)-(1.3) is investigated using the homogeneous test equation
(2.1)y′′=(iν)2y(t).By applying the general method (1.2)-(1.3) to the test equation (2.1) yields
(2.2)[yn+1hyn+1′]=D[ynhyn′],z=νh,D(H)=[1-HbT(I+HA)-1e1-HbT(I+HA)-1c-Hb'T(I+HA)-1e1-Hb'T(I+HA)-1c],
where H=z2,e=[1⋯1]T,c=[c1⋯cm]T. Here D(H) is the stability matrix of the RKN method and its characteristic polynomial
(2.3)ξ2-tr(D(z2))ξ+det(D(z2))=0
is the stability polynomial of the RKN method. Solving the difference system (2.2), the computed solution is given by
(2.4)yn=2|c||ρ|ncos(ω+nϕ).
The exact solution of (2.1) is given by
(2.5)y(tn)=2|σ|cos(χ+nz).
Equations (2.4) and (2.5) led to the following definition.Definition 2.1 (phase-lag).
Apply the RKN method (1.2)-(1.3) to (1.1). Next, we define the phase-lag φ(z)=z-ϕ. If φ(z)=O(zq+1), then the RKN method is said to have phase-lag order q. Additionally, the quantity α(z)=1-|ρ| is called amplification error. If α(z)=O(zv+1), then the RKN method is said to have dissipation order v.Let us denote(2.6)R(z2)=trace(D),S(z2)=det(D).
From Definition 2.1, it follows that
(2.7)φ(z)=z-cos-1(R(z2)2S(z2)),|ρ|=S(z2).
Let us denote R(z2) and S(z2) in the following form:
(2.8)R(z2)=2+α1z2+⋯+αmz2m(1+λ^z2)m,S(z2)=1+β1z2+⋯+βmz2m(1+λ^z2)m,
where λ^=2λ2 is diagonal element in the Butcher tableau.Based on the functionsR(z2) and S(z2) defined as (2.8), a few properties of the functions R and S are summarized in the following theorem which is introduced by Van der Houwen and Sommeijer [2]. The development of dispersion relations is according to the following theorem.Theorem 2.2.
(1) The functionsR(z2) and S(z2) are consistent, dispersive, and dissipative of orders r,q, and v, respectively,
(2.9)eiz[2cos(z)-R(z2)]+S(z2)-1=O(zr+2)R(z2)-2S(z2)cos(z)=O(zq+2)S(z2)-1=O(zv+1).
(2) An RKN method of algebraic orderr, dispersion of order q, and dissipation order v possess functions R and S that are consistent, dispersive, and dissipative of orders r,q, and v.
(3) IfS(z2)≡1, then the order of consistency and dispersion of R and S is equal.Proof (see Van der Houwen and Sommeijer [2]).
Based on the above theorem, the dispersion relations are developed. Form=3,r=4 the dispersion relation of order six (q=6) in terms of αi and βi is
(2.10)order6β3-α3=-8λ6+12λ4+1360-λ22
and for the dispersion relations up to order eight (q=8) for m=4,r=4 are given by
(2.11)order6α3-β3=32λ6-24λ4+2λ23-1360,(2.12)order812α3-β4+α4=16λ8-10λ4+14λ245-32240.
The following quantity is used to determine the dissipation constant of the formula.(i)
form=3(2.13)1-|ρ|=(3λ2-12β1)z2-(152λ4+12β2-32β1λ2-18β12)z4-(-352λ6-32β2λ2+154β1λ4-14β1β2+38β12λ2+12β3+116β13)z6+O(z8).(ii)
form=4(2.14)1-|ρ|=(4λ2-12β1)z2+(-12λ4-12β2+18β12+2β1λ2)z4+(14β1β2-12β12λ2+2β2λ2-12β3-116β13+32λ6-6β1λ4)z6+O(z8).Notice that the fourth-order method is already dispersive order four and dissipative order five due to consistency of the method. Furthermore, dispersive order is even and dissipative order is odd.We next discuss the stability properties of method for solving (1.1) by considering the scalar test problem (2.1) applied to the method (1.2)-(1.3), from which the expression given in (2.2) is obtained. Eliminating yn′ and yn+1′ in (2.2), we obtain a difference equation of the form
(2.15)yn+2-R(H)yn+1+S(H)yn=0.
The characteristic equation associated with (2.15) is given as in (2.3). Since our concerned here is with the analysis of high-order dispersive RKN method, we therefore drop the necessary condition of periodicity interval that is, S(H)≡1. Hence, the method derived will be with empty interval of periodicity. We now consider the interval of absolute stability of RKN method. We therefore need the characteristic equation (2.3) to have roots with modulus less than one so that approximate solution will converge to zero as n tends to infinity. For convenience, we note the following definition adopted for method (2.2).Definition 2.3.
An interval(-Ha,0) is called the interval of absolute stability of the method (2.2) if for all H∈(-Ha,0), |ξ1,2|<1.


## 3. Construction of the Method

In the following, we will derive a three-stage fourth-order and a four-stage fourth-order accurate DIRKN method with dispersive order six and eight, respectively, by taking into account the dispersion relation in Section2. The RKN parameters must satisfy the following algebraic conditions to find fourth-order accuracy as given in Hairer and Wanner [27]:
(3.1)order1∑‍bi′=1(3.2)order2∑‍bi=12,∑‍bi′ci=12(3.3)order3∑‍bici=16,12∑‍bi′ci2=16(3.4)order412∑‍bici2=124,16∑‍bi′ci3=124,∑‍bi′aijcj=124(3.5)order516∑‍bici3=1120,∑‍biaijcj=1120,124∑‍bi′ci4=1120,(3.6)14∑‍bi′ciaijcj=1120,12∑‍bi′aijcj2=1120.
For most methods, the ci need to satisfy
(3.7)12ci2=∑j=1m‍aij(i=1,…,m).The following strategies are used for developing our new efficient pairs.(1)
The high-order DIRKN formula with high order of dispersion. Our aim is to find the ratioκ (phase-lag order/algebraic order) as high as possible and the dissipation constant is “small”.(2)
The following quantities as in [24] should be as small as possible:
(a)
C(s+2)=∥τ^(s+2)-τ(s+2)∥2/∥τ^(s+1)∥2andC'(s+2)=∥τ^'(s+2)-τ'(s+2)∥2/∥τ^'(s+1)∥2,(b)
B(s+2)=∥τ^(s+2)∥2/∥τ^(s+1)∥2andB'(s+2)=∥τ^'(s+2)∥2/∥τ^'(s+1)∥2,whereτ^(s+1)andτ^'(s+1) are called error coefficients for y^n+1 and y^n+1′, respectively.(3)
The strategy to control the error is based on the phase-lag procedure first introduced by Simos [18] and also see [19, 20]. A local error estimation at the point tn+1 is determined by the expressions δn+1=y^n+1-yn+1,δn+1′=y^'n+1-yn+1′ where y^n+1 and y^'n+1 are solutions using the third-order formula. These local error estimations can be used to control the step size h by the standard formula [28, 29]
(3.8)hn+1=0.9hn(TolEst)1/(s+1),where 0.9 is a safety factor,Est=max{∥δn+1∥∞,∥δn+1′∥∞} represents the local error estimation at each step and Tol is the accuracy required which is the maximum allowable local error. If Est < Tol, then the step is accepted, and we applied the accepted procedure of performing local extrapolation (or higher-order mode) meaning that the more accurate approximation will be used to advance the integration. If Est ≥ Tol, then the step is rejected and the h will be updated using formula (3.8).


### 3.1. The Three-Stage DIRKN Formula

#### 3.1.1. The Fourth-Order Formula

In this section, we derive the fourth-order formula with dispersive order six and dissipative order five. Sharp et al. [9] stated that fourth-order method with dispersive order eight does not exist. Therefore, the method of algebraic order four (r=4) with dispersive order six (q=6) and dissipative order five (v=5) is now considered. From algebraic conditions (3.1)–(3.4) and (3.7), it formed eleven equations with thirteen unknowns to be solved. We let b1=0 and λ be a free parameter. Therefore, the following solution of one-parameter family is obtained:
(3.9)a31=[288λ3-24λ-72λ2-24λ23+3-3+12λ3][12(12λ-3+3)],a21=-2λ2+16-312,a32=-1+96λ3-8λ-24λ22(12λ-3+3),b1=0,b2=14+312,b3=14-312,b1′=0,b2′=12,b3′=12,c1=2λ,c2=12-36,c3=12+36.From the above solution, we are going to derive a method with dispersion of order-six. The six order dispersion relation (2.10) needs to be satisfied and this can be written in terms of RKN parameters which correspond to the above family of solution yields the following equation:
(3.10)(28803λ4+(960-14403)λ3+(120-403)λ2+(1203-192)λ-113+18)240(12λ-3+3)=0,
and solving for λ yields
(3.11)-0.1015757589,0.09374433416,0.2097189023,and0.1056624327.
The first two values will give us a nonempty stability interval while the others will produce the methods with empty stability interval. Taking the first two values of λ gives us two fourth-order diagonally implicit RKN methods with dispersive order six. For λ=-0.1015757589 it will give the method with PLTE
(3.12)∥τ(5)∥=1.875825×10-3,∥τ′(5)∥2=1.697439×10-3
for yn and yn′, respectively. The dissipation constant and the stability interval are 1.19×10-4z6+O(z8) and (-8.10,0), respectively.


#### 3.1.2. The Third-Order Formula

Based on the values ofA and c in Section 3.1.1, we now derive the three-stage third-order embedded formula. Solving equations (3.1)–(3.3) simultaneously yields a solution for b^1 and b^2 in terms of b^3,
(3.13)b^1=-0.147183860011593+1.39296300725792b^3b^2=0.647183860011593-2.39296300725792b^3
and b^'i is the same as the fourth-order formula in Section 3.1.1. With this solution, B(5) and C(5) are functions in terms of b^3. Next, we plot the graph for B(5) and C(5) against b^3. We only consider b^3=[0.107,0.3] with B(5) and C(5) lying between [2.612,0.189], and [0.909,0.189], respectively. From numerical experiments, the optimal pair b^3=0.1085 and giving B(5)=1.3184 and C(5)=0.6367, respectively, and ∥τ^(4)∥=1.624880×10-3. We denote this pair as DIRKN4(3)6New method (see Table 2).Table 2
The DIRKN4(3)6New method.
-
0.2031515178
0.02063526960
1
/
2
-
√
3
/
6
0.001693829777
0.02063526960
1
/
2
+
√
3
/
6
−0.0040532720
0.2944222365
0.02063526960
0
1
/
4
+
3
/
12
1
/
4
-
3
/
12
0
1
/
2
1
/
2
0.0039526263
0.3875473737
0.1085
0
1
/
2
1
/
2


### 3.2. The Four-Stage DIRKN Formula

#### 3.2.1. The Fourth-Order Formula

To derive four-stage fourth-order (r=4) with dispersive order eight (q=8), (3.1)–(3.4) and (3.7) together with the dispersion relation of order six equation (2.11) are solved simultaneously and will yield the following solution:
(3.14)c1=2λ,c2=12-36,c3=12+36,c4=12-36,a21=16-312-2λ2,a32=16+312-2λ2,a43=16-312-2λ2,a31=a41=a42=0,a11=a22=a33=a44=2λ2,b1=b1′=b2′=0,b3′=14-312,b2=3(80λ2-1)10(3-3+24λ23+24λ-12λ3-288λ3+72λ2),b4=-1-60λ23-15λ+5λ3+360λ3+1203λ35(3-3+24λ23+24λ-12λ3-288λ3+72λ2).By substituting the above solution to the dispersion relation of order eight (q=8), (2.12) gives us expression in terms of λ(3.15)(5806080λ7-1451520λ6-1451520λ63+241920λ53-967680λ5+60480λ4+181440λ43-806403λ3+147168λ3+44856λ2-29736λ23+924λ3-1752λ-585+3493)/[120960(12λ-3+3)]=0.0
and solving for λ will give us the values -0.2752157925, -0.08524516029, 0.04719733276, 0.1682412065, 0.2490198846, 0.6846776634, and -0.1056624327. Numerical results show that choosing λ=-0.08524516029 will give us smallest dissipation constant hence more accurate scheme compared to the others. The dissipation constant and the stability interval are 4.84×10-5z6+O(z8) and (-8.188,0), respectively.


#### 3.2.2. The Third-Order Formula

Based on the values ofA and c in Section 3.2.1, here we solve third-order embedded formula for the values of b^i and b^i′ obtaining
(3.16)b^1=-0.159774247344685+1.51211971235225b^3b^2=0.659774247344687-2.51211971235225b^3-b^4b^1′=0,b^2′=12-b^4′,b^3′=12,
where b^3,b^4, and b^4′ are free parameters. The B(5) and C(5) are functions in terms of b^3 and b^4 while B'(5) and C'(5) are in terms of b^4′. By setting b^3=0.108, we have B(5) and C(5) in terms of b^4.Similarly, we plot the graph forB(5) and C(5) against b^4. Consider b^4∈[0.13,0.3] with B(5) and C(5) lying between [1.324,2.170] and [0.302,1.745], respectively, while b^4′∈[0.0,0.4] with B'(5) and C'(5) lying between [0.732,2.45] and [0.592,1.165], respectively. From numerical experiments, the optimal pair chosen are b^4=0.14 and b^4′=0.28. With these values will give
(3.17)B(5)=1.369,C(5)=0.460,B'(5)=1.277,C'(5)=0.819
and ∥τ^(4)∥=1.294485×10-3and∥τ^'(4)∥=1.645005×10-3. The pair we denote by DIRKN4(3)8New method (see Table 3).Table 3
The DIRKN4(3)8New method.
c
1
λ
^
1
/
2
-
3
/
6
(
1
/
6
-
3
/
12
-
λ
^
)
λ
^
(
1
/
2
+
3
/
6
)
0
(
1
/
6
+
3
/
12
-
λ
^
)
λ
^
1
/
2
-
3
/
6
0
0
(
1
/
6
-
3
/
12
-
λ
^
)
λ
^
0
b
2
1
/
4
-
3
/
12
b
4
0
0
1
/
2
1
/
2
b
^
1
'
b
^
2
'
0.108
0.14
0
0.22
0.5
0.28
Wherec1=-0.1704903206, b2=0.2332957499, b4=0.1610418175, b^1'=0.00353468159, b^2'=0.24846531841 and λ^=2λ2=0.01453347471.


## 3.1. The Three-Stage DIRKN Formula

### 3.1.1. The Fourth-Order Formula

In this section, we derive the fourth-order formula with dispersive order six and dissipative order five. Sharp et al. [9] stated that fourth-order method with dispersive order eight does not exist. Therefore, the method of algebraic order four (r=4) with dispersive order six (q=6) and dissipative order five (v=5) is now considered. From algebraic conditions (3.1)–(3.4) and (3.7), it formed eleven equations with thirteen unknowns to be solved. We let b1=0 and λ be a free parameter. Therefore, the following solution of one-parameter family is obtained:
(3.9)a31=[288λ3-24λ-72λ2-24λ23+3-3+12λ3][12(12λ-3+3)],a21=-2λ2+16-312,a32=-1+96λ3-8λ-24λ22(12λ-3+3),b1=0,b2=14+312,b3=14-312,b1′=0,b2′=12,b3′=12,c1=2λ,c2=12-36,c3=12+36.From the above solution, we are going to derive a method with dispersion of order-six. The six order dispersion relation (2.10) needs to be satisfied and this can be written in terms of RKN parameters which correspond to the above family of solution yields the following equation:
(3.10)(28803λ4+(960-14403)λ3+(120-403)λ2+(1203-192)λ-113+18)240(12λ-3+3)=0,
and solving for λ yields
(3.11)-0.1015757589,0.09374433416,0.2097189023,and0.1056624327.
The first two values will give us a nonempty stability interval while the others will produce the methods with empty stability interval. Taking the first two values of λ gives us two fourth-order diagonally implicit RKN methods with dispersive order six. For λ=-0.1015757589 it will give the method with PLTE
(3.12)∥τ(5)∥=1.875825×10-3,∥τ′(5)∥2=1.697439×10-3
for yn and yn′, respectively. The dissipation constant and the stability interval are 1.19×10-4z6+O(z8) and (-8.10,0), respectively.


### 3.1.2. The Third-Order Formula

Based on the values ofA and c in Section 3.1.1, we now derive the three-stage third-order embedded formula. Solving equations (3.1)–(3.3) simultaneously yields a solution for b^1 and b^2 in terms of b^3,
(3.13)b^1=-0.147183860011593+1.39296300725792b^3b^2=0.647183860011593-2.39296300725792b^3
and b^'i is the same as the fourth-order formula in Section 3.1.1. With this solution, B(5) and C(5) are functions in terms of b^3. Next, we plot the graph for B(5) and C(5) against b^3. We only consider b^3=[0.107,0.3] with B(5) and C(5) lying between [2.612,0.189], and [0.909,0.189], respectively. From numerical experiments, the optimal pair b^3=0.1085 and giving B(5)=1.3184 and C(5)=0.6367, respectively, and ∥τ^(4)∥=1.624880×10-3. We denote this pair as DIRKN4(3)6New method (see Table 2).Table 2
The DIRKN4(3)6New method.
-
0.2031515178
0.02063526960
1
/
2
-
√
3
/
6
0.001693829777
0.02063526960
1
/
2
+
√
3
/
6
−0.0040532720
0.2944222365
0.02063526960
0
1
/
4
+
3
/
12
1
/
4
-
3
/
12
0
1
/
2
1
/
2
0.0039526263
0.3875473737
0.1085
0
1
/
2
1
/
2


## 3.1.1. The Fourth-Order Formula

In this section, we derive the fourth-order formula with dispersive order six and dissipative order five. Sharp et al. [9] stated that fourth-order method with dispersive order eight does not exist. Therefore, the method of algebraic order four (r=4) with dispersive order six (q=6) and dissipative order five (v=5) is now considered. From algebraic conditions (3.1)–(3.4) and (3.7), it formed eleven equations with thirteen unknowns to be solved. We let b1=0 and λ be a free parameter. Therefore, the following solution of one-parameter family is obtained:
(3.9)a31=[288λ3-24λ-72λ2-24λ23+3-3+12λ3][12(12λ-3+3)],a21=-2λ2+16-312,a32=-1+96λ3-8λ-24λ22(12λ-3+3),b1=0,b2=14+312,b3=14-312,b1′=0,b2′=12,b3′=12,c1=2λ,c2=12-36,c3=12+36.From the above solution, we are going to derive a method with dispersion of order-six. The six order dispersion relation (2.10) needs to be satisfied and this can be written in terms of RKN parameters which correspond to the above family of solution yields the following equation:
(3.10)(28803λ4+(960-14403)λ3+(120-403)λ2+(1203-192)λ-113+18)240(12λ-3+3)=0,
and solving for λ yields
(3.11)-0.1015757589,0.09374433416,0.2097189023,and0.1056624327.
The first two values will give us a nonempty stability interval while the others will produce the methods with empty stability interval. Taking the first two values of λ gives us two fourth-order diagonally implicit RKN methods with dispersive order six. For λ=-0.1015757589 it will give the method with PLTE
(3.12)∥τ(5)∥=1.875825×10-3,∥τ′(5)∥2=1.697439×10-3
for yn and yn′, respectively. The dissipation constant and the stability interval are 1.19×10-4z6+O(z8) and (-8.10,0), respectively.


## 3.1.2. The Third-Order Formula

Based on the values ofA and c in Section 3.1.1, we now derive the three-stage third-order embedded formula. Solving equations (3.1)–(3.3) simultaneously yields a solution for b^1 and b^2 in terms of b^3,
(3.13)b^1=-0.147183860011593+1.39296300725792b^3b^2=0.647183860011593-2.39296300725792b^3
and b^'i is the same as the fourth-order formula in Section 3.1.1. With this solution, B(5) and C(5) are functions in terms of b^3. Next, we plot the graph for B(5) and C(5) against b^3. We only consider b^3=[0.107,0.3] with B(5) and C(5) lying between [2.612,0.189], and [0.909,0.189], respectively. From numerical experiments, the optimal pair b^3=0.1085 and giving B(5)=1.3184 and C(5)=0.6367, respectively, and ∥τ^(4)∥=1.624880×10-3. We denote this pair as DIRKN4(3)6New method (see Table 2).Table 2
The DIRKN4(3)6New method.
-
0.2031515178
0.02063526960
1
/
2
-
√
3
/
6
0.001693829777
0.02063526960
1
/
2
+
√
3
/
6
−0.0040532720
0.2944222365
0.02063526960
0
1
/
4
+
3
/
12
1
/
4
-
3
/
12
0
1
/
2
1
/
2
0.0039526263
0.3875473737
0.1085
0
1
/
2
1
/
2


## 3.2. The Four-Stage DIRKN Formula

### 3.2.1. The Fourth-Order Formula

To derive four-stage fourth-order (r=4) with dispersive order eight (q=8), (3.1)–(3.4) and (3.7) together with the dispersion relation of order six equation (2.11) are solved simultaneously and will yield the following solution:
(3.14)c1=2λ,c2=12-36,c3=12+36,c4=12-36,a21=16-312-2λ2,a32=16+312-2λ2,a43=16-312-2λ2,a31=a41=a42=0,a11=a22=a33=a44=2λ2,b1=b1′=b2′=0,b3′=14-312,b2=3(80λ2-1)10(3-3+24λ23+24λ-12λ3-288λ3+72λ2),b4=-1-60λ23-15λ+5λ3+360λ3+1203λ35(3-3+24λ23+24λ-12λ3-288λ3+72λ2).By substituting the above solution to the dispersion relation of order eight (q=8), (2.12) gives us expression in terms of λ(3.15)(5806080λ7-1451520λ6-1451520λ63+241920λ53-967680λ5+60480λ4+181440λ43-806403λ3+147168λ3+44856λ2-29736λ23+924λ3-1752λ-585+3493)/[120960(12λ-3+3)]=0.0
and solving for λ will give us the values -0.2752157925, -0.08524516029, 0.04719733276, 0.1682412065, 0.2490198846, 0.6846776634, and -0.1056624327. Numerical results show that choosing λ=-0.08524516029 will give us smallest dissipation constant hence more accurate scheme compared to the others. The dissipation constant and the stability interval are 4.84×10-5z6+O(z8) and (-8.188,0), respectively.


### 3.2.2. The Third-Order Formula

Based on the values ofA and c in Section 3.2.1, here we solve third-order embedded formula for the values of b^i and b^i′ obtaining
(3.16)b^1=-0.159774247344685+1.51211971235225b^3b^2=0.659774247344687-2.51211971235225b^3-b^4b^1′=0,b^2′=12-b^4′,b^3′=12,
where b^3,b^4, and b^4′ are free parameters. The B(5) and C(5) are functions in terms of b^3 and b^4 while B'(5) and C'(5) are in terms of b^4′. By setting b^3=0.108, we have B(5) and C(5) in terms of b^4.Similarly, we plot the graph forB(5) and C(5) against b^4. Consider b^4∈[0.13,0.3] with B(5) and C(5) lying between [1.324,2.170] and [0.302,1.745], respectively, while b^4′∈[0.0,0.4] with B'(5) and C'(5) lying between [0.732,2.45] and [0.592,1.165], respectively. From numerical experiments, the optimal pair chosen are b^4=0.14 and b^4′=0.28. With these values will give
(3.17)B(5)=1.369,C(5)=0.460,B'(5)=1.277,C'(5)=0.819
and ∥τ^(4)∥=1.294485×10-3and∥τ^'(4)∥=1.645005×10-3. The pair we denote by DIRKN4(3)8New method (see Table 3).Table 3
The DIRKN4(3)8New method.
c
1
λ
^
1
/
2
-
3
/
6
(
1
/
6
-
3
/
12
-
λ
^
)
λ
^
(
1
/
2
+
3
/
6
)
0
(
1
/
6
+
3
/
12
-
λ
^
)
λ
^
1
/
2
-
3
/
6
0
0
(
1
/
6
-
3
/
12
-
λ
^
)
λ
^
0
b
2
1
/
4
-
3
/
12
b
4
0
0
1
/
2
1
/
2
b
^
1
'
b
^
2
'
0.108
0.14
0
0.22
0.5
0.28
Wherec1=-0.1704903206, b2=0.2332957499, b4=0.1610418175, b^1'=0.00353468159, b^2'=0.24846531841 and λ^=2λ2=0.01453347471.


## 3.2.1. The Fourth-Order Formula

To derive four-stage fourth-order (r=4) with dispersive order eight (q=8), (3.1)–(3.4) and (3.7) together with the dispersion relation of order six equation (2.11) are solved simultaneously and will yield the following solution:
(3.14)c1=2λ,c2=12-36,c3=12+36,c4=12-36,a21=16-312-2λ2,a32=16+312-2λ2,a43=16-312-2λ2,a31=a41=a42=0,a11=a22=a33=a44=2λ2,b1=b1′=b2′=0,b3′=14-312,b2=3(80λ2-1)10(3-3+24λ23+24λ-12λ3-288λ3+72λ2),b4=-1-60λ23-15λ+5λ3+360λ3+1203λ35(3-3+24λ23+24λ-12λ3-288λ3+72λ2).By substituting the above solution to the dispersion relation of order eight (q=8), (2.12) gives us expression in terms of λ(3.15)(5806080λ7-1451520λ6-1451520λ63+241920λ53-967680λ5+60480λ4+181440λ43-806403λ3+147168λ3+44856λ2-29736λ23+924λ3-1752λ-585+3493)/[120960(12λ-3+3)]=0.0
and solving for λ will give us the values -0.2752157925, -0.08524516029, 0.04719733276, 0.1682412065, 0.2490198846, 0.6846776634, and -0.1056624327. Numerical results show that choosing λ=-0.08524516029 will give us smallest dissipation constant hence more accurate scheme compared to the others. The dissipation constant and the stability interval are 4.84×10-5z6+O(z8) and (-8.188,0), respectively.


## 3.2.2. The Third-Order Formula

Based on the values ofA and c in Section 3.2.1, here we solve third-order embedded formula for the values of b^i and b^i′ obtaining
(3.16)b^1=-0.159774247344685+1.51211971235225b^3b^2=0.659774247344687-2.51211971235225b^3-b^4b^1′=0,b^2′=12-b^4′,b^3′=12,
where b^3,b^4, and b^4′ are free parameters. The B(5) and C(5) are functions in terms of b^3 and b^4 while B'(5) and C'(5) are in terms of b^4′. By setting b^3=0.108, we have B(5) and C(5) in terms of b^4.Similarly, we plot the graph forB(5) and C(5) against b^4. Consider b^4∈[0.13,0.3] with B(5) and C(5) lying between [1.324,2.170] and [0.302,1.745], respectively, while b^4′∈[0.0,0.4] with B'(5) and C'(5) lying between [0.732,2.45] and [0.592,1.165], respectively. From numerical experiments, the optimal pair chosen are b^4=0.14 and b^4′=0.28. With these values will give
(3.17)B(5)=1.369,C(5)=0.460,B'(5)=1.277,C'(5)=0.819
and ∥τ^(4)∥=1.294485×10-3and∥τ^'(4)∥=1.645005×10-3. The pair we denote by DIRKN4(3)8New method (see Table 3).Table 3
The DIRKN4(3)8New method.
c
1
λ
^
1
/
2
-
3
/
6
(
1
/
6
-
3
/
12
-
λ
^
)
λ
^
(
1
/
2
+
3
/
6
)
0
(
1
/
6
+
3
/
12
-
λ
^
)
λ
^
1
/
2
-
3
/
6
0
0
(
1
/
6
-
3
/
12
-
λ
^
)
λ
^
0
b
2
1
/
4
-
3
/
12
b
4
0
0
1
/
2
1
/
2
b
^
1
'
b
^
2
'
0.108
0.14
0
0.22
0.5
0.28
Wherec1=-0.1704903206, b2=0.2332957499, b4=0.1610418175, b^1'=0.00353468159, b^2'=0.24846531841 and λ^=2λ2=0.01453347471.


## 4. Problems Tested

In order to evaluate the effectiveness of the new embedded methods, we solved several problems which have oscillatory solutions. The code developed uses constant and variable step size mode and the results obtained are compared with the methods proposed in [10, 11, 21, 29, 30]. Table 4 and Figures 1, 2, 3, 4, 5, 6, 7, and 8 show the numerical results for all methods used. These codes have been denoted by the following:(i)
DIRKN4(3)8New: a new 4(3) pair with phase-lag order 8 derived in this paper.(ii)
DIRKN4(3)6New: a new 4(3) pair with phase-lag order 6 derived in this paper.(iii)
RKND: a fourth-order explicit RKN derived by Dormand [29].(iv)
FPRKN: a phase-fitted fourth-order explicit RKN derived by Papadopoulos et al. [21].(v)
DIRKN4(3)Imoni: a 4(3) RKN pair of three-stage derived by Imoni et al. [10].(vi)
DIRKN4(3)Raed: a 4(3) RKN pair derived by Al-Khasawneh et al. [11].(vii)
DIRK4(3)Butcher: a 4(3) Runge-Kutta pair with six-stage, L-stable, and FSAL property derived by Butcher and Chen [30].Table 4
Numerical results for Problem1 to Problem 8 with tend=104.
Problem
h
DIRKN4(3)6
DIRKN4(3)8
PFRKN
RKND
1
0.025
3.641739(−2)
5.515151(−4)
7.685603(−3)
3.246585(−1)
0.0125
1.121169(−3)
1.412513(−5)
2.381676(−4)
2.064614(−2)
0.00625
3.522474(−5)
3.162919(−6)
8.044719(−6)
1.293771(−3)
2
0.25
4.968941(−3)
8.176026(−5)
1.059320(−3)
4.581100(−2)
0.125
1.553957(−4)
2.266041(−6)
3.290638(−5)
2.860912(−3)
0.0625
4.858102(−6)
3.894456(−8)
1.025317(−6)
1.789827(−4)
3
0.025
5.050886(−2)
7.713254(−4)
1.065804(−2)
4.502978(−1)
0.0125
1.554479(−3)
1.993991(−5)
3.302186(−4)
2.863241(−2)
0.00625
4.884725(−5)
4.407458(−6)
1.114967(−5)
1.794225(−3)
4
0.025
5.050886(−2)
7.713235(−4)
1.065799(−2)
4.502978(−1)
0.0125
1.554479(−3)
1.993979(−5)
3.302153(−4)
2.863241(−2)
0.00625
4.884726(−5)
4.407458(−6)
1.114963(−5)
1.794225(−3)
5
0.25
3.142511(−5)
2.225689(−6)
5.906399(−6)
7.329390(−5)
0.125
9.872922(−7)
1.311623(−7)
2.060768(−7)
4.155967(−6)
0.0625
3.122323(−8)
7.739299(−9)
8.477345(−9)
4.155967(−6)
6
0.025
8.130019(−1)
5.143138(−3)
6.323161(−2)
1.581409(−1)
0.0125
7.270152(−3)
1.088848(−4)
1.513299(−3)
6.194539(−2)
0.00625
2.200129(−4)
3.162925(−6)
4.672316(−5)
4.048100(−3)
7
0.2
3.103741(−3)
4.927970(−5)
6.483674(−4)
1.326941(−2)
0.1
9.720586(−5)
1.472241(−6)
1.980997(−5)
8.293575(−4)
0.05
3.068984(−6)
7.146919(−8)
6.187652(−7)
5.183626(−5)
8
0.004
2.227290(−6)
2.224132(−6)
1.495525(−5)
1.548512(−5)
0.002
1.231234(−5)
1.231204(−5)
1.310806(−5)
1.314136(−5)
0.001
1.735729(−5)
1.735723(−5)
1.740702(−5)
1.740910(−5)Figure 1
Efficiency curves for Problem1.Figure 2
Efficiency curves for Problem2.Figure 3
Efficiency curves for Problem3.Figure 4
Efficiency curves for Problem4.Figure 5
Efficiency curves for Problem5.Figure 6
Efficiency curves for Problem6.Figure 7
Efficiency curves for Problem7.Figure 8
Efficiency curves for Problem8.Problem 1 (homogenous).
One has(4.1)d2y(t)dt2=-100y(t),y(0)=1,y′(0)=-2,0≤t≤10.
Exact solution y(t)=-(1/5)sin(10t)+cos(10t)Problem 2 (inhomogeneous problem studied by Allen and Wing [31]).
One has(4.2)d2y(t)dt2=-y(t)+t,y(0)=1,y′(0)=2,0≤t≤15π.
Exact solution y(t)=sin(t)+cos(t)+tProblem 3 (problem studied by van der Houwen and Sommeijer [1]).
One has(4.3)y′′(t)=-v2y(t)+(v2-1)sin(t),y(0)=1,y′(0)=v+1,
where v≫1,t∈[0,50].
Exact solution isy(t)=cos(vt)+sin(vt)+sin(t). Numerical result is for the case v=10.Problem 4 (inhomogeneous system studied by Franco [8]).
One has(4.4)y′′(t)+(1012-992-9921012)y(t)=ε(932cos(2t)-992sin(2t)932sin(2t)-992cos(2t)),y(0)=(-1+ε1),y′(0)=(-1010+2ε),t∈[0,10].
Exact solution
(4.5)y(t)=(-cos(10t)-sin(10t)+εcos(2t)cos(10t)+sin(10t)+εsin(2t))Problem 5 (The Duffing's equations as given in [32]).
One has(4.6)y′′(t)+y(t)+(y(t))3=0.002cos(1.01t),y(0)=0.200426728067,y′(0)=0,t∈[0,10].
The exact solution computed by the Galerkin method and given by
(4.7)y(t)=∑i=04‍a2i+1cos[1.01(2i+1)t],
with a1=0.200179477536,a3=0.246946143×10-3,a5=0.304014×10-6,a7=0.374×10-9, and a9<10-12.Problem 6 (Inhomogeneous problem studied by Lambert and Watson [33]).
One has(4.8)d2y1(t)dt2=-v2y1(t)+v2f(t)+f′′(t),y1(0)=a+f(0),y1′(0)=f′(0),d2y2(t)dt2=-v2y2(t)+v2f(t)+f′′(t),y2(0)=f(0),y2′(0)=va+f′(0),0≤t≤20
Exact solution is y1(t)=acos(vt)+f(t),y2(t)=asin(vt)+f(t),f(t) is chosen to be e-0.05t and parameters v and a are 20 and 0.1, respectively.Problem 7 (an almost periodic orbit problem given in stiefel and bettis [34]).
One has(4.9)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,0≤t≤1000.
Exact solution y1(t)=cos(t)+0.0005tsin(t), y2(t)=sin(t)-0.0005tcos(t).Problem 8 (linear Strehmel-Weiner problem studied in Cong [35]).
One has(4.10)y′′(t)=(-20.20-9.67989.6-10000-6004.2-9.60-5.8)y(t)+(150cos(10t)75cos(10t)75cos(10t)),(4.11)y(0)=(12-2),y′(0)=(000),t=[0,10]
Exact solution (4.12)y(t)=(cos(t)+2cos(5t)-2cos(10t)2cos(t)+cos(5t)-cos(10t)-2cos(t)+cos(5t)-cos(10t)).


## 5. Numerical Results

In this section, we evaluate the effectiveness of the new DIRKN pairs derived in the previous section when they are applied to the numerical solution of eight problems which are model and nonmodel problems for constant and variable step size.For constant step size, Table4 shows the absolute maximum error for fourth-order DIRKN4(3)6New, DIRKN4(3)8New, PFRKN, and RKND methods when solving Problems 1–8 with three different step values. From numerical results in Table 4, we observed that the new DIRKN4(3)8New method is more accurate compared with PFRKN and RKND method which is not related to the dispersion of the method. Also the new DIRKN4(3)6New method is more accurate compared with RKND method while comparable with PFRKN method for certain problem. Notice that all the methods are of the same algebraic order.For variable step size, Figures1–8 show the decimal logarithm of the maximum global error for the solution (MAXE) versus the function evaluations and also the decimal logarithm of the maximum global error for the solution (MAXE) versus the time taken. From Figures 1–8, we observed that DIRKN4(3)6New and DIRKN4(3)8New performed better compared to DIRKN4(3)Raed, DIRKN4(3)Imoni, and DIRK4(3)Butcher for integrating second-order differential equations possessing an oscillatory solution in terms of function evaluations. This is due to the fact that when using DIRK4(3)Butcher, the second-order system of ODEs needs to be transformed to a first-order system and hence the dimension is doubled. Furthermore, the DIRK4(3)Butcher method has five effective stages per step. This means that the function evaluations per step for DIRK4(3)Butcher are more than the function evaluations per step used in the DIRKN4(3)6New and DIRKN4(3)8New methods which have three and four function evaluations per step, respectively. Meanwhile the DIRKN4(3)Raed, DIRKN4(3)Imoni methods are less accurate compared with the DIRKN4(3)6New and DIRKN4(3)8New when phase lag and dissipation is considered. Therefore the new methods converges faster and consequently less steps are needed for a specified value of tolerance even though DIRKN4(3)Raed and DIRKN4(3)Imoni have the same number of stages per step with DIRKN4(3)8New and DIRKN4(3)6New method respectively.


## 6. Conclusion

In this paper, we have derived two 4(3) pairs, namely, DIRKN4(3)6New and DIRKN4(3)8New which have dispersive order six and eight, respectively, with “small” dissipation constant which is suitable for problems with oscillating solutions and moderate frequency. From the results shown in Table4 and Figures 1–8, we conclude that the new methods are more efficient for integrating second-order equations possessing an oscillatory solution when compared to the RKND derived in [29], PFRKN derived in [21], DIRK4(3)Butcher as in [30] and also with the others the same type of pairs for example, DIRKN4(3)Imoni pair derived in [10] and DIRKN4(3)Raed derived in [11]. The DIRKN4(3)6New method is the most efficient method since it has three function evaluations per step.



---

*Source: 324989-2012-09-20.xml*