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.

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:

An -stage Runge-Kutta-Nyström (RKN) method for the numerical integration of the IVP is given by where

The RKN parameters , and are assumed to be real and is the number of stages of the method. Introduce the -dimensional vectors , and and matrix , where , and , respectively. The latter contains the class of Diagonally Implicit RKN (DIRKN) methods for which all the entries in the diagonal of are equal. An embedded pair of DIRKN methods is based on the method of order and the other DIRKN method of order and can be expressed in Butcher notation by the table of coefficients (see Table 1).

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., [1221]).

Most of the numerical methods developed for solving (1.1) are in constant stepsize (see [13, 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

By applying the general method (1.2)-(1.3) to the test equation (2.1) yields where . Here is the stability matrix of the RKN method and its characteristic polynomial is the stability polynomial of the RKN method. Solving the difference system (2.2), the computed solution is given by The exact solution of (2.1) is given by 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 . If , then the RKN method is said to have phase-lag order . Additionally, the quantity is called amplification error. If , then the RKN method is said to have dissipation order .

Let us denote From Definition 2.1, it follows that Let us denote and in the following form: where is diagonal element in the Butcher tableau.

Based on the functions and defined as (2.8), a few properties of the functions and 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 functions and are consistent, dispersive, and dissipative of orders , and , respectively,
(2) An RKN method of algebraic order , dispersion of order , and dissipation order possess functions R and S that are consistent, dispersive, and dissipative of orders , and .
(3) If , then the order of consistency and dispersion of and is equal.

Proof (see Van der Houwen and Sommeijer [2]). Based on the above theorem, the dispersion relations are developed. For the dispersion relation of order six () in terms of and is and for the dispersion relations up to order eight () for are given by
The following quantity is used to determine the dissipation constant of the formula.(i) for (ii) for

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 and in (2.2), we obtain a difference equation of the form 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, . 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 tends to infinity. For convenience, we note the following definition adopted for method (2.2).

Definition 2.3. An interval is called the interval of absolute stability of the method (2.2) if for all , .

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 Section 2. The RKN parameters must satisfy the following algebraic conditions to find fourth-order accuracy as given in Hairer and Wanner [27]: For most methods, the need to satisfy

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), (b), where are called error coefficients for and , 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 is determined by the expressions where and are solutions using the third-order formula. These local error estimations can be used to control the step size by the standard formula [28, 29] where 0.9 is a safety factor, 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 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 ( with dispersive order six ( and dissipative order five ( 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 and be a free parameter. Therefore, the following solution of one-parameter family is obtained:

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: and solving for yields 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 it will give the method with PLTE for and , respectively. The dissipation constant and the stability interval are and , respectively.

3.1.2. The Third-Order Formula

Based on the values of and 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 and in terms of , and is the same as the fourth-order formula in Section 3.1.1. With this solution, and are functions in terms of . Next, we plot the graph for and against . We only consider with and lying between , and , respectively. From numerical experiments, the optimal pair and giving and , respectively, and . We denote this pair as DIRKN4(3)6New method (see Table 2).

3.2. The Four-Stage DIRKN Formula
3.2.1. The Fourth-Order Formula

To derive four-stage fourth-order () with dispersive order eight (), (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:

By substituting the above solution to the dispersion relation of order eight (), (2.12) gives us expression in terms of and solving for will give us the values , , , , , , and . Numerical results show that choosing will give us smallest dissipation constant hence more accurate scheme compared to the others. The dissipation constant and the stability interval are and , respectively.

3.2.2. The Third-Order Formula

Based on the values of and in Section 3.2.1, here we solve third-order embedded formula for the values of and obtaining where and are free parameters. The and are functions in terms of and while and are in terms of . By setting , we have and in terms of .

Similarly, we plot the graph for and against . Consider with and lying between and , respectively, while with and lying between and , respectively. From numerical experiments, the optimal pair chosen are and . With these values will give and . The pair we denote by DIRKN4(3)8New method (see Table 3).

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].

Problem 1 (homogenous). One has Exact solution

Problem 2 (inhomogeneous problem studied by Allen and Wing [31]). One has Exact solution

Problem 3 (problem studied by van der Houwen and Sommeijer [1]). One has where .
Exact solution is . Numerical result is for the case .

Problem 4 (inhomogeneous system studied by Franco [8]). One has Exact solution

Problem 5 (The Duffing's equations as given in [32]). One has The exact solution computed by the Galerkin method and given by with and .

Problem 6 (Inhomogeneous problem studied by Lambert and Watson [33]). One has Exact solution is is chosen to be and parameters and are 20 and 0.1, respectively.

Problem 7 (an almost periodic orbit problem given in stiefel and bettis [34]). One has Exact solution , .

Problem 8 (linear Strehmel-Weiner problem studied in Cong [35]). One has Exact solution

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, Table 4 shows the absolute maximum error for fourth-order DIRKN4(3)6New, DIRKN4(3)8New, PFRKN, and RKND methods when solving Problems 18 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, Figures 18 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 18, 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 Table 4 and Figures 18, 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.

Acknowledgments

This work is partially supported by IPTA Fundamental Research Grant, Universiti Putra Malaysia (Project no. 05-10-07-385FR) and UPM Research University Grant Scheme (RUGS) (Project no. 05-01-10-0900RU).