A FOURTH ORDER IMPLICIT SYMMETRIC AND SYMPLECTIC EXPONENTIALLY FITTED RUNGE-KUTTA-NYSTR¨OM METHOD FOR SOLVING OSCILLATORY PROBLEMS

. In this paper, we derive an implicit symmetric, symplectic and exponentially ﬁtted Runge-Kutta-Nystr¨om (ISSEFRKN) method. The new integrator ISSEFRKN2 is of fourth order and integrates exactly diﬀerential systems whose solutions can be expressed as linear combinations of functions from the set { exp( λt ) , exp( − λt ) | λ ∈ C } , or equivalently { sin( ωt ) , cos( ωt ) | λ = iω, ω ∈ R } . We analysis the periodicity stability of the derived method ISSEFRKN2. Some the existing implicit RKN methods in the literature are used to compare with ISSEFRKN2 for several oscillatory problems. Numerical results show that the method ISSEFRKN2 possess a more accuracy among them.

1. Introduction. In this paper we focus on the initial value problems (IVP) related to a system of second-order ODEs of the form y = f (x, y), y(x 0 ) = y 0 , y (x 0 ) = y 0 , x ∈ [0, x end ], (1) whose solutions exhibit an oscillatory character. Problems of this type are of great interest in applied sciences such as molecular dynamics, orbital mechanics, and electronics. High accuracy of integration is often required in these areas. Until now there are broadly two categories of approaches to numerical integration of the IVP (1): indirect and direct. On one hand, if a new variable u is introduced to represent the first derivative y , then IVP (1) is turned into the partitioned system of first order equations y = u, u = f (x, y), y(x 0 ) = y 0 , u(x 0 ) = y 0 , and the problem can be solved by the general Runge-Kutta (RK) methods or partitioned Runge-Kutta (PRK) methods (see Refs. [16,4,17,5,22,23,6]). On the a ij f (t 0 + c j h, Y j ), i = 1, · · · , s, which can be expressed in the Butcher tableau as c s 1 1 a s1 · · · a ss 1 1b 1 · · ·b s 1 b 1 · · · b s The objective of this section is to specify when the RKN method (3) is symmetric, symplectic and exponentially fitted. This is the cornerstone of our paper. In the following subsections, we will put forward these three important properties step by step.
2.1. Symmetry conditions. The key to understanding symmetry is the concept of the adjoint method. We denote a one-step method for second-order ODEs (1) as Φ h : (y 0 , y 0 ) T → (y 1 , y 1 ) T . Here, from y 0 to y 1 , the variable goes forward with a step h. Then the symmetry of Φ h is defined as follows.
Definition 2.1. The adjoint method Φ * h of a one-step method Φ h is the inverse map of the original method with reversed time step −h, i.e., Φ * h := Φ −1 −h . In other words, In the case of the s-stage RKN method (3), a set of sufficient conditions for the symmetric method are given by In this paper we consider the method (3) whose coefficients are z-dependent, as we do for exponentially fitted type methods (see Refs. [22,23,25]). Then the method (3) is symmetric if its coefficients satisfy the following conditions where z = iωh, ω is the principal frequency of the problem. We assume that the coefficients of the method (3) are even functions of h, as we frequently encounter in the case of EFRKN methods, so that these conditions reduce to the conditions (4).

2.2.
Symplectic conditions. Now, we turn to the symplectic conditions for the scheme (3). Symplecticity is defined for a Hamiltonian system. On many occasions, the problem under consideration takes the form of a Hamiltonian systeṁ where S is a symmetric positive definite constant matrix. This system is equivalent to the second-order equation (1) with f (x, q) = −S −1 ∂ ∂q U (x, q). The following definition can be found in [8]. Here we only list it without explanation.
Definition 2.2. A one-step method is symplectic if for every smooth Hamiltonian function H and for every step size h, the corresponding flow preserves the differential 2-form Accordingly, the scheme (3) for the problem (1) is symplectic if and only if For the left side of this equation, we have Eliminating dy 0 in the second term of this equation by inserting Eq.(3) we obtain Therefore, Eq. (5) holds if the following conditions are satisfied 2.3. Exponential fitting conditions. Following Albrecht's approach (see Refs. [1,2]), each stage of the scheme (3) can be viewed as a linear multistep method on a non-equidistant grid. With each stage one can associate a linear functions as follows: • for the internal stages, a ij y (x + c j h), i = 1, 2, · · · , s; • for the final stages, By requiring the internal and final stages vanish for the functions from the set {exp(±iωx)} leads to the following equations Note that cosh(z) = (e z + e −z )/2 and sinh(z) = (e z − e −z )/2, then the equations (7) imply that In this paper, we call the method (3) satisfied the exponentially fitted (EF) conditions (8) and (9) as exponentially fitted RKN (EFRKN) method.
3. Algebraic order conditions. In this section, we will present algebraic order conditions for exponentially fitted Runge-Kutta-Nystöm (EFRKN) methods. For an EFRKN method, the local truncation errors in the approximations of its solution and its derivative can be expressed as where F (j) (y 0 ) denotes an elementary differential and the terms d or equivalently, Since our RKN method is a particular case of the RKN method considered in [7], following the approach in [7] we consider the following assumptions The order conditions up to fifth order for the RKN method (3) are the following ones: Order 1 requires: Order 2 requires in addition: Order 3 requires in addition: Order 4 requires in addition: Order 5 requires in addition: From Theorem 2.1 in [7], we know that the EFRKN method (3) has algebraic order at least 2.
4. Construction of implicit symmetric symplectic EFRKN method. In this section we construct implicit EFRKN method under the symmetry, symplecticity and exponential fitting conditions obtained in the previous section.

WENJUAN ZHAI AND BINGZHEN CHEN
Until now, we obtain an implicit symmetric and symplectic exponentially fitted Runge-Kutta-Nyström method which coefficients are given by We denote this method as ISSEFRKN2. In order to specify the algebra order of ISSEFRKN2, we give the Taylor expansions of the coefficients.
From the taylor expansions, we can verify that our method ISSEFRKN2 satisfies algebraic conditions up to fourth order, but doesn't satisfy the fifth order condition d i c 4 i − 1 5 = 0. So, the method ISSEFRKN2 is of order 4. The method ISSEFRKN2 is exponentially fitted. So when the solution of (1) can be expressed as linear combinations of functions from the set {exp(±iωx)}, ISSEFRKN2 has higher efficiency and competence than other integrators which are not exponentially fitted. This will be shown in the numerical studies. 5. Periodicity region of the new method. Now we start to analyze the stability property of our new method. Stability means that the numerical solutions remain bounded as we move further away from the starting point. For classical RKN methods, the stability properties are checked using the second order linear test model Recall that the new symmetric and symplectic exponentially fitted implicit RKN method derived in the previous section is dependent on the complex number λ = iω, where ω > 0 is an estimate of the dominant frequency. Applying an s-stage ISSEFRKN method (3) to the test model (19) yields where The stability behavior of the numerical solution depends on the eigenvalues or the spectrum of the stability matrix M = M (H 2 , ν 2 ). Eliminating y 0 and y 1 from (20) and the equation that is obtained from (20) by replacing the subscript 0 by 1 gives the difference equation Accordingly, the characteristic equation is given by (iii) If R s = (0, ∞) × (0, ∞) except possibly for a discrete set of curves, the method is A-stable; (iv) If R p = (0, ∞) × (0, ∞) except possibly for a discrete set of curves, the method is P-stable.
The periodicity region of the method ISSEFRKN2 is depicted in Figure 1.  . This method is not exponentially fitted. • ISSEFRKN2: The symmetric and symplectic exponentially fitted two-stage fourth-order RKN method (18) proposed in this paper.
Compared with our method ISSEFRKN2, the method DIRKNRaed or DIRKNNora is neither symmetric, symplectic nor EF, ISSRKN2 is not EF. In our numerical experiments we have solved the non-linear equations with the Newton iteration method and taking initial values Y . The iteration is carried out until the difference between the Euclidean norm of two successive iterations attains 10 −8 . The maximum number of iterations is 1000.
The criterion used in the numerical comparisons is the usual test based on computing the maximum global error in the solution over the whole integration interval. In Figures 2-6 we show the decimal logarithm of the maximum global error (log10(err)) versus the number of steps required by each code on a logarithmic scale (log10(nsteps)). All computations are carried out in double precision arithmetic (16 significant digits of accuracy). Problem 1. We consider the linear problem with variable coefficients y + 4x 2 y = (4x 2 − ω 2 ) sin(ωx) − 2 sin(x 2 ), x ∈ [0, x end ] y(0) = 1, y (0) = ω, whose analytic solution is given by y(x) = sin(ωx) + cos(x 2 ). This solution represents a periodic motion that involves a constant frequency and a variable frequency. In our test we choose the parameter values ω = 10, λ = 10i, x end = 10, h = 1/2 m , m = 4, 5, 6, 7, and the numerical results are stated in Figure 2.
As we can see, its solution is independent on µ.
In this problem, the parameters are chosen as µ = 0.25, λ = √ µi = 0.5i, x end = 10, and the numerical results presented in Figure 6 have been computed with the integration steps h = 1/2 m , m = 1, 2, 3, 4. From Figures 2-6, we can find that symmetric and symplectic method ISSRKN2 is more efficient than the nonsymmetric or nonsymplectic methods. But, all of them are not better than the exponentially fitted RKN method ISSEFRKN2 when the exact solutions can be expressed in term of triangular functions.  7. Conclusions. In this paper a two-stage IEFRKN integrator which is symmetric and symplectic have been derived. Like the existing EFRKN integrators (see [25] for example), the coefficients of the new method depend on the product of the dominant frequency ω and the step size h. When the parameter z(= ωh) approaches to zero, the ISSEFRKN method reduces to the classical RKN method. The numerical experiments carried out show that the new method is more efficient than the twostage classical symmetric and symplectic RKN integrator and other RKN methods used in the numerical studies.