A 5 ( 4 ) Embedded Pair of Explicit Trigonometrically-Fitted Runge – Kutta – Nyström Methods for the Numerical Solution of Oscillatory Initial Value Problems

Musa A. Demba 1,2,†, Norazak Senu 1,† and Fudziah Ismail 2,* 1 Department of Mathematics, Universiti Putra Malaysia, Serdang 43400, Malaysia; musdem2004@gmail.com (M.A.D.); Norazak@upm.edu.my (N.S.) 2 Department of Mathematics and Institute for Mathematical Research, Universiti Putra Malaysia, Serdang 43400, Malaysia * Correspondence: fudziah@upm.edu.my; Tel.: +601-6971-5654 or +234-803-551-8557 † These authors contributed equally to this work.


Introduction
In this study, our focus is on the numerical solution of the initial value problems (IVPs) of second-order differential equations; whose first derivative does not appear explicitly of the form : y = f (x, y), x ∈ [x 0 , X], y(x 0 ) = y 0 , y (x 0 ) = y 0 , whose solutions have a noticeable oscillatory character, where y ∈ d and f : d → d are sufficiently differentiable.Such problems frequently occur in several areas of applied sciences such as: theoretical physics, celestial mechanics, nuclear chemistry, nuclear physics, electronics, molecular dynamics and elsewhere.Many numerical methods have been developed for the numerical integration of Equation ( 1), among them are Runge-Kutta-Nyström (RKN) methods.Exponentially/trigonometrically-fitted RKN methods have been studied by Simos in [1],Van de Vyver in [2], Kalogiratou and Simos in [3] and Shiwei Liu et al. in [4].Franco in [5] proposed a 5(3) pair of explicit adapted Runge-Kutta-Nyström methods for the numerical integration of perturbed oscillators, Van de Vyver in [6,7] proposed a Runge-Kutta-Nyström pair for the numerical integration of perturbed oscillators, a 5(3) pair of explicit RKN methods for oscillatory problems.Senu et al. in [8] proposed an embedded explicit RKN method for solving oscillatory problems.Recently, Tsitouras in [9] proposed fitted modifications of RKN pairs, Franco et al. in [10] proposed two new embedded pair of explicit Runge-Kutta methods adapted to the numerical solution of oscillatory problems, and Anastassi and Kosti in [11] proposed a 6(4) optimized embedded Runge-Kutta-Nyström pair for the numerical solution of periodic problems.
In this paper, we develop a new efficient embedded explicit trigonometrically-fitted RKN method based on the technique proposed by Simos in [12] for Runge-Kutta (RK) methods.The constructed method can exactly integrate the test equation y = −w 2 y and the numerical results show the efficacy of the new method.The remaining part of this paper is arranged as follows: in Section 2, we give the basic theory of an explicit Runge-Kutta-Nyström pair, the basic definition of trigonometrically-fitted RKN method and the derivation of an explicit trigonometrically-fitted RKN method.Section 3 deals with the derivation of the proposed method.In Section 4, we analyze the algebraic order of the new pair from their local truncation errors and give the interval of absolute stability of the new pair.In Section 5, we present the numerical results.In Section 6, we give a brief discussion about the graphs and the last section deals with the conclusions.

Basic Theory
The general form of an explicit k-stage RKN method is given by: where y n+1 and y n+1 represent the approximation of y(x n+1 ) and y (x n+1 ), respectively, where x n+1 = x n + h, n = 0, 1, ..., or in Butcher Tableau as : An embedded p(q) pair of RKN methods is based on the method (c, A, b, d) of order p and the other RKN method (c, A, b, d) of order q(q < p).The higher order produces the solution (y n+1 , y n+1 ), while the lower order method produces the solution ŷn+1 and ŷ n+1 , which is used for the estimation of local truncation error only.An embedded pair is characterized by the Butcher tableau below: In this work, a variable step size algorithm using the embedded technique is used because it provides cheap local error estimation.Local error estimation at the point x n+1 = x n + h is determined by δ n+1 = ŷn+1 − y n+1 and δ n+1 = ŷ n+1 − y n+1 .
Let Est n+1 = max( δ n+1 ∞ , δ n+1 ∞ ) represent local error estimation to control the step size h.For the numerical integration of the oscillatory problems, we use the step-size control procedure given by Shiwei in [4]: where Tol is the tolerance (requested local error).It should be noted that the N th order approximation y n is used as the initial value for the (n+1)th step, that is to say, the embedded pair is applied in local extrapolation mode or higher order mode.Definition 1. Runge-Kutta-Nyström method (2)-( 4) is said to be trigonometrically-fitted if it integrates exactly the function e iwx and e −iwx or equivalently sin(wx) and cos(wx) with w > 0 the principal frequency of the problem when applied to the test equation y = −w 2 y; Leading to a system of equations as derived below: When an explicit Runge-Kutta-Nyström method (2)-( 4) is applied to the test equation y = −w 2 y, the method becomes: where Let y n = e Iwx .By computing the value of y n+1 , y n and y n+1 and substituting in the Equations ( 5)- (8) and by using e Iv = cos(v) + I sin(v) and comparing the real and imaginary part, we obtain the following system of equations: where v = wh.

Derivation of the Proposed Method
In this section, we will construct a new efficient embedded explicit trigonometrically-fitted RKN method.
In this study, Embedded Runge-Kutta-Nyström 5(4)M pair (ERKN5(4)M) is used as given by Senu in [13].The coefficients of the method are given in Table 1.Solving the above system of Equations ( 9)-( 12) using the coefficients of the lower order method (order 4) for b 2 , b 3 , d 2 , d 3 , we obtain the solution as given in Equation ( 13).

Algebraic Order and Error Analysis
In this section, we will find the local truncation error of the new methods and verify their algebraic order.We first find the Taylor series expansion of the actual solution y(x n + h), the first derivative of the actual solution y (x n + h), the approximate solution y n+1 , and the first derivative of the approximate solution y n+1 .The local truncation error (LTE) of y and its first derivative y is given as: The LTE and LTE der of the lower order method (order 4) are: ( From Equation (18), we can see that the algebraic order of the lower order method is 4 because all of the coefficients up to h 4 vanished.Similarly, the LTE and LTE der of the higher order method (order 5) are : From Equation (19), the higher order method has order 5 because all of the coeffients up to h 5 vanished.

Absolute Stability Analysis of the New Embedded Pair
The linear stability of the RKN method (2)-( 4) is derived by applying it to the test equation y = −w 2 y.Setting H = −(wh) 2 , the numerical solution satisfies the following: It is assumed that P(H) has complex conjugate eigenvalues for sufficiently small values of v [14].With this assumption, a periodic numerical solution is obtained whose characteristic depends on the eigenvalues of P(H), which is called the stability matrix and its characteristic equation can be written as λ 2 − tr(P(H))λ + det(P(H)) = 0.
Hence, we obtain the approximate interval of absolute stability of the higher order method (order 5) of our new embedded pair as (−36.99,0) and the lower order method (order 4) has no interval of absolute stability.
The exact solution is The exact solution is The exact solution is where = 0.001, Ψ = 0.1 and x end = 100.
The numerical results are shown in Tables 2-5.

Discussion
For all of the problems tested and from Figures 1-4, we can deduce that our new method has a lower number of function evaluations per step, which signifies that our new method has less computational costs than the other existing methods, and is therefore more suitable for solving second order ordinary differential equations with oscillatory solutions than the other existing methods in the literature.

Conclusions
In this study, we have derived a new efficient 5(4) embedded pair of explicit trigonometrically-fitted Runge-Kutta-Nyström methods for the solution of oscillatory initial value problems.The numerical results obtained indicate that the function evaluations per step of the new method are less when compared with the other existing embedded pairs.Hence, the new method has less computational costs than the other existing methods, and, therefore, the efficiency of the new method is higher than the other existing methods.

( 14 )
As v → 0, the coefficients b 2 , b 3 , d 2 and d 3 of the lower order method reduces to the coefficients of the original method (lower order ).That is to say, b 2 (0), b 3 (0), d 2 (0) and d 3 (0) are identical to b 2 , b 3 , d 2 and d 3 of the lower order method in ERKN5(4)M.
a 11 a 12 a 13 a 14 a 21 a 22 a 23 a 24 a 31 a 32 a 33 a 34 a 41 a 42 a 43 a 44

Table 2 .
Numerical results for Problem 1.

Table 3 .
Numerical results for Problem 2.

Table 4 .
Numerical results for Problem 3.

Table 5 .
Numerical results for Problem 4.