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

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. TwoDIRKNpairs which are threeand fourstage 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.


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: y f t, y , y t 0 y 0 , y t 0 y 0 .

1.1
An m-stage Runge-Kutta-Nystr öm RKN method for the numerical integration of the IVP is given by 3 where The RKN parameters a ij , b j , b j , and c j are assumed to be real and m is the number of stages of the method.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.

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 y iν 2 y t .

2.1
By applying the general method 1. where Here D H 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 y t n 2|σ| cos χ nz .

2.5
Equations 2.4 and 2.5 led to the following definition.Let us denote R z 2 and S z 2 in the following form: where λ 2λ 2 is diagonal element in the Butcher tableau.
Based on the functions R z 2 and S z 2 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 functions R z 2 and S z 2 are consistent, dispersive, and dissipative of orders r, q, and v, respectively, 2.9 (2) An RKN method of algebraic order r, 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) If S z 2 ≡ 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.For m 3, r 4 the dispersion relation of order six q 6 in terms of α i and β i is and for the dispersion relations up to order eight q 8 for m 4, r 4 are given by The following quantity is used to determine the dissipation constant of the formula.

2.13
ii for m 4

2.14
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 y n and y n 1 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, 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 .

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 c i 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: 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 t n 1 is determined by the expressions δ n 1 y n 1 − y n 1 , δ n 1 y n 1 − y n 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 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 .

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 b 1 0 and λ be a free parameter.Therefore, the following solution of one-parameter family is obtained:

3.9
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.11
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 fourthorder diagonally implicit RKN methods with dispersive order six.For λ −0.1015757589 it will give the method with PLTE

The Third-Order Formula
Based on the values of A and c in Section 3.1.1,we now derive the three-stage third-order embedded formula.Solving equations 3.1 -3.0.6367, respectively, and τ 4 1.624880 × 10 −3 .We denote this pair as DIRKN4 3 6New method see Table 2 .

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:  By substituting the above solution to the dispersion relation of order eight q 8 , 2.12 gives us expression in terms of λ

The Third-Order Formula
Based on the values of A and c in Section 3.     and τ 4 1.294485 × 10 −3 and τ 4  1.645005 × 10 −3 .The pair we denote by DIRKN4 3 8New method see Table 3 .

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

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 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, Figures 1-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.

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 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.
Introduce the m-dimensional vectors c, b, and b and m × m matrix A, where c c 1 , c 2 , . . ., c m T , b b 1 , b 2 , . . ., b m T , b b 1 , b 2 , . . ., b m T , and A a ij , 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
y n and y n , respectively.The dissipation constant and the stability interval are 1.19×10 −4 z 6 O z 8 and −8.10, 0 , respectively.

3 simultaneously yields a solution for b 1 and b 2 in terms of b 3 ,b 1 −0.147183860011593 1 .39296300725792 b 3 b 2 0.647183860011593 − 2 .39296300725792 b 3 3. 13 and
b i is the same as the fourth-order formula in Section 3.1.1.With this solution, B5 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 C5 λ 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 −5 z 6 O z 8 and −8.188, 0 , respectively.

2 . 1 , 4 b 1 where b 3 , b 4 ,
here we solve third-order embedded formula for the values of b i and b i obtaining b 1 −0.159774247344685 1.51211971235225 b 3 b 2 0.659774247344687 − 2.51211971235225 b 3 − b 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 for B 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 B5 1.369, C 5 0.460, B5 1.277, C 5 0.819 3.17
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., 12-21 .

Table 4 :
Numerical results for Problem 1 to Problem 8 with t end 10 4 .