# Phase-Fitted and Amplification-Fitted Higher Order Two-Derivative Runge-Kutta Method for the Numerical Solution of Orbital and Related Periodical IVPs

Authors: N. A. Ahmad; N. Senu; F. Ismail
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

A phase-fitted and amplification-fitted two-derivative Runge-Kutta (PFAFTDRK) method of high algebraic order for the numerical solution of first-order Initial Value Problems (IVPs) which possesses oscillatory solutions is derived. We present a sixth-order four-stage two-derivative Runge-Kutta (TDRK) method designed using the phase-fitted and amplification-fitted property. The stability of the new method is analyzed. The numerical experiments are carried out to show the efficiency of the derived methods in comparison with other existing Runge-Kutta (RK) methods.


---

## Body

## 1. Introduction

Consider the numerical solution of the IVPs for first-order Ordinary Differential Equations (ODEs) in the form of(1)q′=ft,q,qt0=q0,whose solutions show an observable oscillatory or periodical behavior. Such problems occur in several fields of applied sciences, for example, circuit simulation, molecular dynamics, orbital mechanics, mechanics, electronics, and astrophysics, which have attracted the concern of a number of researchers. In general, most problems with oscillatory or periodical behavior are second order or higher order. Hence, it is important to reduce the higher order problems to first-order problems in order to solve the ODEs (1).Several well-known authors in their papers have developed phase-fitted and amplification-fitted RK methods. Simos and Vigo-Aguiar [1] constructed a modified phase-fitted RK method with phase-lag of order infinity for the numerical solution of periodic IVPs based on the fifth-algebraic-order RK method of Dormand and Prince. Chen et al. [2] improved traditional RK methods by introducing frequency-depending weights in the update. With the phase-fitting and amplification-fitting conditions and algebraic order conditions, new practical RK integrators are obtained and two of the new methods have updates that are also phase-fitted and amplification-fitted.With the evolution of RK methods, Papadopoulos et al. [3] developed a new Runge-Kutta-Nyström (RKN) method for the numerical solution of the Schrödinger equation with phase-lag and amplification error of order infinity based on the fourth-order RKN method by Dormand, El-Mikkawy, and Prince. Meanwhile Moo et al. [4] derived two new RKN methods for solving second-order differential equations with oscillatory solutions based on two existing RKN methods, a fourth-order three-stage Garcias RKN method and fifth-order four-stage Hairers RKN method. The derived methods both have two variable coefficients with zero amplification error (zero dissipative) and phase-lag of order infinity.In the last few years, Senu et al. [5] constructed zero dissipative explicit RK method for solving second-order ODEs with periodical solutions which has algebraic order three with dissipation of order infinity. Recently, Fawzi et al. in their papers [6, 7] developed fourth-algebraic-order phase-fitted and amplification-fitted modified RK method and fourth-order seven-stage phase-fitted and amplification-fitted RK methods, respectively.Chan and Tsai [8] introduced special explicit TDRK methods by including the second derivative. It involves only one evaluation of f and a number of evaluations of g per step. They presented methods up to five stages and up to seventh order as well as some embedded pairs. Zhang et al. [9] proposed a new fifth-order trigonometrically fitted TDRK method for the numerical solution of the radial Schrödinger equation and oscillatory problems. Meanwhile, Fang et al. [10] and Chen et al. [11] derived two-fourth-order and three practical exponentially fitted TDRK methods, respectively. They compared the new methods with some well-known optimized codes and traditional exponentially fitted RK methods.Awoyemi and Idowu [12] proposed a class of hybrid collocation methods for the solution of general third-order ODEs. The efficiency of the methods is tested using two standard general problems of third-order ODEs. They found out that as the step length h decreases, the accuracy of the methods increases. A year later, Franco [13] presented a class of explicit two-step hybrid methods for the numerical solution of second-order IVPs with reduced number of stages per step. New hybrid methods which are up to orders five and six with optimized error constants are developed.A seventh-order three-step hybrid linear multistep method (HLMM) with three nonstep points is proposed by Jator [14] for the direct solution of the special second-order IVPs. They applied the method in block form as simultaneous numerical integrators over nonoverlapping intervals and then, expressing the method as a one-step method to analyze the stability property of the new method. Senu et al. [15] developed two-step optimized hybrid methods of fifth and sixth order for the integration of second-order oscillatory IVPs based on the existing nonzero dissipative hybrid methods with the requirement of phase-lag, dissipation, or amplification error and the differentiation of the phase-lag relations. They found out that the nonzero dissipative hybrid methods are more suitable to be optimized than phase-fitted methods.Up until now, there are no research findings related to phase-fitting in TDRK methods. Researchers have not yet explored the advantages or disadvantages of applying phase-fitted techniques to TDRK methods. Hence, in this paper, a new sixth-order four-stage phase-fitted and amplification-fitted TDRK methods is constructed. In Section2, an overview of TDRK method is given. In Section 3, phase-fitted and amplification-fitted conditions are considered. The new phase-fitted and amplification-fitted TDRK method is constructed in Section 4. Meanwhile in Section 5, the analysis of the stability property is discussed. The numerical results, discussion, and conclusion are dealt with in Sections 6, 7, and 8, respectively.


## 2. Two-Derivative Runge-Kutta Methods

Consider the scalar ODEs (1) with g:RN→RN. For this case, the second derivative is assumed to be known where(2)q′′=gq≔f′qfq,g:RN⟶RN.An explicit TDRK method for the numerical integration of IVPs (1) is given by(3)qn+1=qn+Δt∑i=1sbifqi+Δt2∑i=1sb^igQi,Qi=qn+Δt∑j=1saijfqj+Δt2∑j=1sa^ijgQj,where i=1,…,s.The explicit TDRK methods are presented with the coefficients in (3) using the Butcher tableau as follows:(4)
c
A
A
^
b
T
b
^
T
.Explicit methods with minimal number of function evaluations can be developed by considering the methods in the form(5)qn+1=qn+Δtftn,qn+Δt2∑i=1sb^igtn+Δtci,Qi,Qi=qn+Δtciftn,qn+Δt2∑j=1i-1a^ijgtn+Δtcj,Qj,where i=2,…,s.The above method is called special explicit TDRK methods. The unique part of this method is that it involves only one evaluation off per step compared to many evaluations of f per step in traditional explicit RK methods. Its Butcher tableau is given as follows: (6)
c
A
^
b
^
T
.The TDRK parametersa^ij, b^i, and ci are assumed to be real and s is the number of stages of the method. We introduce the s-dimensional vectors b^,c and s×s matrix, A^ where b^=b^1,b^2,…,b^sT,c=c1,c2,…,csT, and A^=a^ij, respectively.The order conditions for special explicit TDRK methods are given in Table1.Table 1
Order conditions for special explicit TDRK methods.
Order Conditions 1 b T e = 1 2 b ^ T e = 1 2 3 b ^ T c = 1 6 4 b ^ T c 2 = 1 12 5 b ^ T c 3 = 1 20 b ^ T A ^ c = 1 120 6 b ^ T c 4 = 1 30 b ^ T c A ^ c = 1 180 b ^ T A ^ c 2 = 1 360 7 b ^ T c 5 = 1 42 b ^ T c 2 A ^ c = 1 252 b ^ T c A ^ c 2 = 1 504 b ^ T A ^ c 3 = 1 840 b ^ T A ^ 2 c = 1 5040


## 3. Phase-Fitted and Amplification-Fitted Property

Consider the following linear scalar equation(7)q′=iλq.The exact solution of this equation with the initial value q(t0)=q0 satisfies(8)qt0+Δt=H0zq0,where H0(z)=exp⁡(z),z=iv. The exact solution experiences a phase advance v=λΔt and after a period of time Δt, the amplification remains constant.Applying the test equation (7) to the TDRK method yields(9)q1=Hzq0,where(10)Hz=1+v2b^I-v2A^-1e+iv+v3b^I-v2A^-1c,where e=1,…,1T.H ( z ) is the stability function of the TDRK method. Denote the real and imaginary part of H(z) by U(v) and V(v), respectively. For small Δt, arg⁡H(z)=tan-1⁡V(v)/U(v) and |H(z)|=U2(v)+V2(v). The following definition came from the analysis above.Definition 1 (van der Houwen and Sommeijer [16]).
The quantities(11)P~v=v-arg⁡Hz,D~v=1-Hzare called the phase-lag (or dispersion) and the error of amplification factor (or dissipation) of the method, respectively. If(12)P~v=cϕvα+1+Ovα+3,D~v=cdvβ+1+Ovβ+3,then the method is called dispersive of order α and dissipative of order β, respectively. If(13)P~v=0,D~v=0,the method is called phase-fitted (or zero dispersive) and amplification-fitted (or zero dissipative), respectively.Theorem 2.
The method is phase-fitted and amplification-fitted if and only if(14)Uv=cos⁡v,Vv=sin⁡v.When (5) are applied to the first-order ODEs (1), for any (ζ+1)th differentiable function g(q), the local truncation error LTE=q(t0+Δt)=O(Δtζ+1). Hence the method is said to have (algebraic) order ζ.Denote(15)ECζ+1v=∑i=1jτjζ+121/2,where τj(ζ+1) is the error coefficient of the method. The positive number(16)ECζ+1=limv→0⁡ECζ+1vis the error constant of the method.


## 4. Derivation of the New Phase-Fitted and Amplification-Fitted Method

A TDRK method is phase-fitted and amplification-fitted if and only if Theorem2 is satisfied. In deriving the new method, the phase-fitted and amplification-fitted property is combined to the existing TDRK method. Hence, the derivation of the new method will be discussed next.For this study, two-derivative sixth-algebraic-order method presented by Chan and Tsai [8] is considered. The coefficient of the method is given as follows.Butcher Tableau for Sixth-Order TDRK Method. (17)
​
0
1
4
1
32
2
3
-
2
81
20
81
1
5
4
-
6
5
9
20
3
40
64
225
27
200
1
180Considering the stability function (10) for sixth-order four-stage TDRK method, by letting b^1 and c4 as free parameters, we have(18)
I
=
1
0
0
0
0
1
0
0
0
0
1
0
0
0
0
1
,
e
=
1
1
1
1
,
b
^
=
b
^
1
64
225
27
200
1
180
,
c
=
0
1
4
2
3
c
4
,
A
^
=
0
1
32
0
-
2
81
20
81
0
5
4
-
6
5
9
20
0
.Substitute matrices (18) into the stability function, H(z), which is (10) and separate the real part and the imaginary part of H(z) and denote them as U(v) and V(v), respectively. For optimized value of maximum global error, the combination of b^1 and c4 is chosen as free parameters. By using Theorem 2, solve (14) to find the coefficients of b^1 and c4 and this leads to(19)1-1740v2-b^1v2+124v4-1720v6+151840v8=cos⁡v,v-1180c4v3-29180v3+1120v5-16480v7=sin⁡v.Solving (19) we will obtain(20)b^1=-151840v2-51840+51840cos⁡v+22032v2-2160v4+72v6-v8,c4=-136v36480sin⁡v-6480v+1044v3-54v5+v7.Asv→0, the following Taylor expansions are obtained:(21)b^1=340-1181440v6+13628800v8-1479001600v10+187178291200v12+⋯,c4=1+1126v4-12016v6+1221760v8-134594560v10+17264857600v12+⋯.The following expansions are obtained by direct calculation:(22)b^Te=1740-151840-v8+72v6-2160v4+22032v2-51840+51840cos⁡vv2=12+Ov6,b^Tc=29180-16480v7-54v5+1044v3-6480v+6480sin⁡vv3=16+Ov4,b^Tc2=790+1233280v7-54v5+1044v3-6480v+6480sin⁡v2v6=112+Ov4,b^Tc3=245-18398080v7-54v5+1044v3-6480v+6480sin⁡v3v9=120+Ov4,b^Tc4=136+1302330880v7-54v5+1044v3-6480v+6480sin⁡v4v12=130+Ov4,b^TA^c=1120,b^TcA^c=1180,b^TA^c2=1360.Hence, the coefficients given in (17) satisfy all the order conditions of order two to order six. But it failed to satisfy the order condition for order seven. For example,(23)b^Tc5=0≠142+Ov.Hence, it is a sixth-order method. The original method is obtained by Chan and Tsai [8] and it is denoted as TDRK4(6). The error coefficients of TDRK4(6) for order seven are given by(24)τ17=17720,τ27=1270,τ37=1480,τ47=1864,τ57=4673972000.Therefore for TDRK4(6),(25)EC7=1972000566865679.Since we have verified that this new method is order six, hence it is called PFAFTDRK4(6). The error coefficients of PFAFTDRK4(6) are given by(26)τ17=13720-110883911680v7-54v5+1044v3-6480v+6480sin⁡v5v15,τ27=1720-151840v7-54v5+1044v3-6480v+6480sin⁡vv3,τ37=1270,τ47=1864,τ57=4673972000.For PFAFTDRK4(6), we have(27)EC7v=197200094478400000013720-110883911680N5v152+36062554+9447840000001720-151840Nv321/2,where(28)N=v7-54v5+1044v3-6480v+6480sin⁡v.PFAFTDRK4(6) will reduce to its original method that is TDRK4(6) as v→0. Other than that, as v→0, PFAFTDRK4(6) will have the same error constant as TDRK4(6).


## 5. Stability of the New Method

In this section, the linear stability of the method developed is analyzed. Consider the test equation (7) where λ>0. Applying (7) to the special explicit TDRK method produces the difference equation(29)qn+1=Hzqn,z=iv,i2=-1,where H(z) is given as (10).Definition 3.
A TDRK method is said to be absolutely stable if|H(z)|<1 for all ∈(-v,0).The stability polynomial of the PFAFTDRK4(6) method is given as follows: (30)Hz=1+v+12v2+16v3+124v4+1120v5+1720v6+15040v7+172576v8-1362880v9+13628800v10+139916800v11-1479001600v12-16227020800v13+187178291200v14+⋯.The comparison of the stability region of the PFAFTDRK4(6) method up to vi, where i=6,12,14, and its original method is plotted in Figure 1.Figure 1
Stability region of PFAFTDRK4(6) method for different order.The stability interval of this method with the coefficients ofv6, v12, and v14 is (-4.270,0.000), (-4.157,0.000), and (-4.150,0.000), respectively. Observing from the stability regions plotted in Figure 1, when the order of the coefficients tend to infinity, the stability interval becomes closer to the stability interval of the original method. The stability interval of the original method is (-4.060,0.000).From the stability interval, we can actually find the biggest value ofΔt the method can take so that it will always remain stable. It is known that v=λΔt and the value of λ comes from the test problems. Hence, dividing v with λ will lead us to the value of Δt. The following stability test will demonstrate how the stability regions are used for practical purposes. We have(31)q′=λq-φ+φ′,q0=φ0,Reλ<0,t∈0,2000,where φ(t) is a smooth function. Take λ=-1 and φ(t)=sin(t). The exact solution is q(t)=φ(t).A method is stable if the maximum global error is small and converges to its exact solution. Otherwise, the method is unstable if it has a bigger maximum global error which means it is actually diverging from its exact solution. We will show the relationship betweenΔt,λ, and H(z) by running the stability test. The method is stable when Δt=4.15 where this is the biggest value Δt can take for the method to become stable in this stability test. The global errors are collected in Table 2 for a variety of Δt values.Table 2
Stability test for PFAFTDRK4(6) using coefficient of v12 with λ=-1 for variable Δt.
Δ t H ( z ) Global error 4.55 2.3772841250 1.627253 × 1 0 165 4.35 1.5452837990 5.363198 × 1 0 86 4.15 0.9838244901 1.032649 × 1 0 0 1.00 0.3678739259 8.849764 × 1 0 - 5 0.15 0.8607079765 1.334905 × 1 0 - 10 0.01 0.9900498337 6.514789 × 1 0 - 13


## 6. Problems Tested and Numerical Results

In this section, the performance of the proposed method PFAFTDRK4(6) is compared with existing RK methods by considering the following problems. All problems below are tested using C code for solving differential equations where the solutions are periodic.Problem 1 (harmonic oscillator).
(32) q 1 ′ t = q 2 t ,  q 1 0 = q 0 1 , q 2 ′ t = - ω 2 q 1 t ,  q 2 0 = q 0 2 , t ∈ 0 , t e n d .Exact solution is (33)q1t=c¯1sin⁡ωt+c¯2cos⁡ωt,q2t=c¯3ωcos⁡ωt-c¯4ωsin⁡ωt.Total energy is given in [17] (34)Eq1,q2=q122+q222=Ψ22,where Ψ depends on the initial conditions.Problem 2 (inhomogeneous problem [18]).
(35) q 1 ′ = q 2 ,  q 1 0 = 1 , q 2 ′ = - 100 q 1 + 99 sin ⁡ t ,  q 2 0 = 11 , t ∈ 0,1000 π .Exact solution is (36)q1t=cos⁡10t+sin⁡10t+sin⁡t,q2t=-10sin⁡10t+10cos⁡10t+cos⁡t.Problem 3 (an “almost” periodic orbit problem [19]).
(37) q 1 ′ = q 2 ,  q 1 0 = 1 , q 2 ′ = - q 1 + 0.001 cos ⁡ t ,  q 2 0 = 1 , q 3 ′ = q 4 ,  q 3 0 = 0 , q 4 ′ = - q 3 + 0.001 sin ⁡ t ,  q 4 0 = 0.995 , t ∈ 0,1000 π .Exact solution is (38)q1t=cos⁡t+0.0005tsin⁡t,q2t=-sin⁡t+0.0005tcos⁡t+0.0005tsin⁡t,q3t=sin⁡t-0.0005tcos⁡t,q4t=cos⁡t+0.0005tsin⁡t-0.0005cos⁡t.Problem 4 (two-body problem [20]).
(39) q 1 ′ = q 2 ,  q 1 0 = 1 , q 2 ′ = - q 1 q 1 2 + q 3 2 3 ,  q 2 0 = 0 , q 3 ′ = q 4 ,  q 3 0 = 0 , q 4 ′ = - q 3 q 1 2 + q 3 2 3 ,  q 4 0 = 1 , t ∈ 0,1000 π .Exact solution is (40)q1t=cos⁡t,q2t=-sin⁡t,q3t=sin⁡t,q4t=cos⁡t.Problem 5 (Duffing problem [21]).
(41) q 1 ′ = q 2 ,  q 1 0 = 0.200426728067 , q 2 ′ = - q 1 - q 1 3 + 0.002 cos ⁡ 1.01 t ,  q 2 0 = 0 , t ∈ 0,1000 π .Exact solution is (42)q1t=0.200179477536cos⁡1.01t+2.46946143×10-4cos⁡3.03t+3.04014×10-7cos⁡5.05t+3.74×10-10cos⁡7.07t,q2t=-0.2021812723sin⁡1.01t-7.482468133×10-4sin⁡3.03t-1.53527070×10-6sin⁡5.05t-2.64418×10-9sin⁡7.07t.Problem 6 (Prothero-Robinson problem [8]).
(43) q ′ = λ q - φ + φ ′ ,  q 0 = φ 0 ,   R e λ < 0 ,   t ∈ 0,1000 π ,where φ(t) is a smooth function. We take λ=-1 and φ(t)=sin⁡(t).
Exact solution isq(t)=φ(t).The following notations are used in Figures2–14:(i)
PFAFTDRK4(6). New phase-fitted and amplification-fitted TDRK method of sixth order and four stages derived in this paper.(ii)
TFTDRK3(5). Existing fifth-order three-stage trigonometrically fitted TDRK method developed by Zhang et al. [9].(iii)
TFRKS6(5). Existing fifth-order six-stage trigonometrically fitted RK method derived by Simos [22].(iv)
TFRKAS6(5). Existing fifth-order six-stage trigonometrically fitted RK method given in Anastassi and Simos [23].(v)
PFAFRKC7(5). Existing fifth-order seven-stage phase-fitted and amplification-fitted RK method developed by Chen et al. [2].(vi)
RKB7(6). Existing sixth-order seven-stage RK method developed by Butcher [24].(vii)
TDRK4(6). Existing sixth-order four-stage TDRK method given by Chan and Tsai [8].(viii)
ABMoulton. Existing fourth-order Adams-Bashforth-Moulton method (a predictor-corrector method) given in Lambert [25].Figure 2
Energy conservation. The logarithm error of energy (MAXERR) at each integration point when solving the harmonic oscillator (Problem1) for ω=8 with initial condition q01=1, q02=-2, and Δt=1/5 with tend=1000.Figure 3
The error at each integration point when solving the harmonic oscillator (Problem1) with ω=8 with initial condition q01=1, q02=-2, and Δt=1/20 with tend=1000.Figure 4
The global error at each integration point when solving the inhomogeneous problem (Problem2) with Δt=1/10.Figure 5
The global error at each integration point when solving the “almost” periodic problem (Problem3) with Δt=1/2.Figure 6
The global error at each integration point when solving the two-body problem (Problem4) with Δt=1/2.Figure 7
The global error at each integration point when solving the Duffing problem (Problem5) with Δt=1/10.Figure 8
The global error at each integration point when solving the Prothero-Robinson problem (Problem6) with Δt=1/20.Figure 9
The efficiency curve for the harmonic oscillator (Problem1) with λ=8 and time step Δt=1/10-1/50(i), i=0,…,4 with tend=1000.Figure 10
The efficiency curve for the inhomogeneous problem (Problem2) with time step Δt=1/2i, i=4,…,8.Figure 11
The efficiency curve for the “almost” periodic problem (Problem3) with time step Δt=1/2i, i=1,…,5.Figure 12
The efficiency curve for the two-body problem (Problem4) with time step Δt=1/2i, i=4,…,8.Figure 13
The efficiency curve for the Duffing problem (Problem5) with time step Δt=1/2i, i=1,…,5.Figure 14
The efficiency curve for the Prothero-Robinson problem (Problem6) with time step Δt=1/2i, i=1,…,5.We represent the performance of these numerical results graphically in Figures2–14.


## 7. Discussion

The results show the typical properties of the new phase-fitted and amplification-fitted TDRK method, PFAFTDRK4(6), which have been derived earlier. The derived method is compared with some well-known existing RK methods. Figure 2 shows the error of the energy at each integration point. From the figures, the phase-fitted and amplification-fitted TDRK method conserved the energy by having smaller magnitude of energy error compared to TDRK4(6), RKB7(6), and ABMoulton. In Figures 3–8, the logarithm number of global error versus the time of integration for different time step, Δt, is plotted for various physical problems. From Figures 3–5, it is observed that global error produced by the PFAFTDRK4(6) method is smaller compared to TDRK4(6), RKB7(6), and ABMoulton. Meanwhile in Figures 6–8, the global errors between PFAFTDRK4(6) and TDRK4(6) are quite close to each other but still the derived method has the smallest global error compared to the others.Next, the global error and the efficiency of the method are plotted over a long period of integration. Figures9–14 represent the efficiency and accuracy of the method developed by plotting the graph of the logarithm of the maximum global error against the logarithm number of function evaluations for longer periods of computations. From the graphs plotted, it can be seen that the PFAFTDRK4(6) method has the smallest maximum global error compared to other existing RK methods which have trigonometrically fitted and phase-fitted and amplification-fitted properties. In Figure 9, as the value of Δt becomes smaller, the maximum global error of the PFAFTDRK4(6) method seems to flatten at the end of the curve. The accuracy of the method depends on the step-size, Δt, and the frequency, λ. The derived method will converge to its original method as the value of Δt becomes smaller.


## 8. Conclusion

In this research, a new phase-fitted and amplification-fitted higher order TDRK method is developed. Based on the numerical results obtained, it can be concluded that the new PFAFTDRK4(6) method is more promising compared to other well-known existing explicit RK methods in terms of accuracy and the number of function evaluations per step.



---

*Source: 1871278-2017-02-14.xml*