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

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 fourstage two-derivative Runge-Kutta (TDRK)method designed using the phase-fitted and amplification-fitted property.The stability of the newmethod is analyzed.Thenumerical experiments are carried out to show the efficiency of the derivedmethods in comparison with other existing Runge-Kutta (RK) methods.


Introduction
Consider the numerical solution of the IVPs for first-order Ordinary Differential Equations (ODEs) in the form of   =  (, ) ,  ( 0 ) =  0 , 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 phasefitted RK method with phase-lag of order infinity for the numerical solution of periodic IVPs based on the fifthalgebraic-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 phasefitting 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 phasefitted 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 amplificationfitted RK methods, respectively.
The explicit TDRK methods are presented with the coefficients in (3) using the Butcher tableau as follows: Explicit methods with minimal number of function evaluations can be developed by considering the methods in the form where  = 2, . . ., .
The above method is called special explicit TDRK methods.The unique part of this method is that it involves only one evaluation of  per step compared to many evaluations of  per step in traditional explicit RK methods.Its Butcher tableau is given as follows: The TDRK parameters â , b , and   are assumed to be real and  is the number of stages of the method.We introduce the -dimensional vectors b,  and  ×  matrix, Â where b = [ b1 , b2 , . . ., b ]  ,  = [ 1 ,  2 , . . .,   ]  , and Â = [â  ], respectively.
The order conditions for special explicit TDRK methods are given in Table 1.

Phase-Fitted and Amplification-Fitted Property
Consider the following linear scalar equation The exact solution of this equation with the initial value ( 0 ) =  0 satisfies where  0 () = exp (),  = V.The exact solution experiences a phase advance V = Δ and after a period of time Δ, the amplification remains constant.Applying the test equation (7) to the TDRK method yields where where  = [1, . . ., 1]  .() is the stability function of the TDRK method.Denote the real and imaginary part of () by (V) and (V), respectively.For small Δ, arg () = tan −1 ((V)/(V)) and |()| = √ 2 (V) +  2 (V).The following definition came from the analysis above.
Definition 1 (van der Houwen and Sommeijer [16]).The quantities are called the phase-lag (or dispersion) and the error of amplification factor (or dissipation) of the method, respectively.If then the method is called dispersive of order  and dissipative of order , respectively.If the method is called phase-fitted (or zero dispersive) and amplification-fitted (or zero dissipative), respectively.
where  (+1)  is the error coefficient of the method.The positive number is the error constant of the method.

Derivation of the New Phase-Fitted and Amplification-Fitted Method
A TDRK method is phase-fitted and amplification-fitted if and only if Theorem 2 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.Considering the stability function (10) for sixth-order four-stage TDRK method, by letting b1 and  4 as free parameters, we have Substitute matrices (18) into the stability function, (), which is (10) and separate the real part and the imaginary part of () and denote them as (V) and (V), respectively.For optimized value of maximum global error, the combination of b1 and  4 is chosen as free parameters.By using Theorem 2, solve (14) to find the coefficients of b1 and  4 and this leads to Solving (19) we will obtain As V → 0, the following Taylor expansions are obtained: The following expansions are obtained by direct calculation: 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, 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  (7)  1 = 17 720 ,  (7)  2 = 1 270 ,  (7)  3 = 1 480 ,  (7)  4 = 1 864 ,  (7)  5 = 4673 972000 . ( Therefore for TDRK4 (6), 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  (7)  1 = 13 720 − 1 10883911680 1 864 ,  (7)  5 = 4673 972000 . ( For PFAFTDRK4(6), we have (944784000000 ( 13 720 where 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).

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 where () is given as (10).
The stability polynomial of the PFAFTDRK4(6) method is given as follows: (30) Coefficient of 12  Coefficient of 14  Coefficient of 6  Original method The comparison of the stability region of the PFAFTDRK4(6) method up to V  , where  = 6, 12, 14, and its original method is plotted in Figure 1.
The stability interval of this method with the coefficients of V 6 , V 12 , and V 14 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 Δ the method can take so that it will always remain stable.It is known that V = Δ and the value of  comes from the test problems.Hence, dividing V with  will lead us to the value of Δ.The following stability test will demonstrate how the stability regions are used for practical purposes.We have where () is a smooth function.Take  = −1 and () = sin().The exact solution is () = ().
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 Δ, , and |()| by running the stability test.The method is stable when Δ = 4.15 where this is the biggest value Δ can take for the method to become stable in this stability test.The global errors are collected in Table 2 for a variety of Δ values.

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.
Exact solution is Total energy is given in [17] where Ψ depends on the initial conditions.
The following notations are used in Figures 2-14: (i) PFAFTDRK4 (6).New phase-fitted and amplificationfitted TDRK method of sixth order and four stages derived in this paper.
We represent the performance of these numerical results graphically in Figures 2-14.

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 amplificationfitted 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, Δ, 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.Figures 9-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 phasefitted and amplification-fitted properties.In Figure 9, as the value of Δ 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, Δ, and the frequency, .The derived method will converge to its original method as the value of Δ becomes smaller.

Table 1 :
Order conditions for special explicit TDRK methods.