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.

1. Introduction

Consider the numerical solution of the IVPs for first-order Ordinary Differential Equations (ODEs) in the form ofwhose 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 and a number of evaluations of 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 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 Section 2, 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 . For this case, the second derivative is assumed to be known where

An explicit TDRK method for the numerical integration of IVPs (1) is given bywhere

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 formwhere

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 , , and are assumed to be real and is the number of stages of the method. We introduce the -dimensional vectors and matrix, where , and , respectively.

The order conditions for special explicit TDRK methods are given in Table 1.

3. Phase-Fitted and Amplification-Fitted Property

Consider the following linear scalar equationThe exact solution of this equation with the initial value satisfieswhere . The exact solution experiences a phase advance and after a period of time , the amplification remains constant.

Applying the test equation (7) to the TDRK method yieldswherewhere .

is the stability function of the TDRK method. Denote the real and imaginary part of by and , respectively. For small , and . The following definition came from the analysis above.

Definition 1 (van der Houwen and Sommeijer [16]). The quantitiesare called the phase-lag (or dispersion) and the error of amplification factor (or dissipation) of the method, respectively. Ifthen the method is called dispersive of order and dissipative of order , respectively. Ifthe 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

When (5) are applied to the first-order ODEs (1), for any th differentiable function , the local truncation error . Hence the method is said to have (algebraic) order .

Denotewhere is the error coefficient of the method. The positive numberis 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 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 and 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 and , respectively. For optimized value of maximum global error, the combination of and is chosen as free parameters. By using Theorem 2, solve (14) to find the coefficients of and and this leads to

Solving (19) we will obtain

As , 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(). The error coefficients of TDRK4() for order seven are given byTherefore for TDRK4(),

Since we have verified that this new method is order six, hence it is called PFAFTDRK4(). The error coefficients of PFAFTDRK4() are given byFor PFAFTDRK4(), we havewhere

PFAFTDRK4() will reduce to its original method that is TDRK4() as . Other than that, as , PFAFTDRK4() will have the same error constant as TDRK4().

5. Stability of the New Method

In this section, the linear stability of the method developed is analyzed. Consider the test equation (7) where . Applying (7) to the special explicit TDRK method produces the difference equationwhere is given as (10).

Definition 3. A TDRK method is said to be absolutely stable if for all .

The stability polynomial of the PFAFTDRK4() method is given as follows:

The comparison of the stability region of the PFAFTDRK4() method up to , where , and its original method is plotted in Figure 1.

The stability interval of this method with the coefficients of , , and is , , and , 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 .

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 and the value of comes from the test problems. Hence, dividing with will lead us to the value of . The following stability test will demonstrate how the stability regions are used for practical purposes. We havewhere is a smooth function. Take and . 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 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.

6. Problems Tested and Numerical Results

In this section, the performance of the proposed method PFAFTDRK4() 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). Exact solution is Total energy is given in [17] where depends on the initial conditions.

Problem 2 (inhomogeneous problem [18]). Exact solution is

Problem 3 (an “almost” periodic orbit problem [19]). Exact solution is

Problem 4 (two-body problem [20]). Exact solution is

Problem 5 (Duffing problem [21]). Exact solution is

Problem 6 (Prothero-Robinson problem [8]). where is a smooth function. We take and .
Exact solution is .

The following notations are used in Figures 214:(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].

We represent the performance of these numerical results graphically in Figures 214.

7. Discussion

The results show the typical properties of the new phase-fitted and amplification-fitted TDRK method, PFAFTDRK4(), 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(), RKB7(), and ABMoulton. In Figures 38, the logarithm number of global error versus the time of integration for different time step, , is plotted for various physical problems. From Figures 35, it is observed that global error produced by the PFAFTDRK4() method is smaller compared to TDRK4(), RKB7(), and ABMoulton. Meanwhile in Figures 68, the global errors between PFAFTDRK4() and TDRK4() 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 914 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() 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 becomes smaller, the maximum global error of the PFAFTDRK4() 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.

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() 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.

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.