Implicit symmetric and symplectic exponentially fitted modified Runge–Kutta–Nyström methods for solving oscillatory problems

Symplectic exponentially fitted RK and RKN methods have been extensively studied by many researchers. Due to their good property, they have been applied to many problems such as pendulum, Morse oscillator, harmonic oscillator, Lennard–Jones oscillator, Kepler’s orbit problem, and so on. In this paper, we construct an implicit symmetric and symplectic exponentially fitted modified Runge–Kutta–Nyström (ISSEFMRKN) method. The new integrator integrates exactly differential systems whose solutions can be expressed as linear combinations of functions from the set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\{\exp(\lambda t),\exp(-\lambda t)\}$\end{document}{exp(λt),exp(−λt)}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda\in\mathbb{C}$\end{document}λ∈C, or equivalently \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\{\sin(\omega t),\cos(\omega t)\}$\end{document}{sin(ωt),cos(ωt)} when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda=i\omega$\end{document}λ=iω, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\omega \in\mathbb{R}$\end{document}ω∈R. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$z=\lambda h$\end{document}z=λh approaches zero, the ISSEFMRKN method reduces to the classical symplectic, symmetric RKN integrator. Numerical experiments are accompanied to show the efficiency and competence of the new method compared with some efficient codes in the literature.


Introduction
In plenty of applied sciences such as celestial mechanics, astrophysics, chemistry, electronics, molecular dynamics, and so forth, the following second-order ODEs initial value problems (IVP) often arise: whose solutions exhibit an oscillatory character. Such problems are of great interest and have been studied extensively. Roughly speaking, there are two categories of approaches to numerical integration of the IVP (1): indirect and direct. In the first place, if a new variable v is introduced to represent the first derivative y , then (1) can be transformed into a partitioned system of first-order equations y = v, v = f (t, y), y(x 0 ) = y 0 , v(t 0 ) = y 0 .
This new reformulated problem can be solved by the general Runge-Kutta (RK) methods, partitioned Runge-Kutta (PRK) methods, or two-step methods (see Refs. [5,6,10,19,20,25,26,29,32,33]). In the second place, the Runge-Kutta-Nyström (RKN) method, which was introduced by Nyström in 1925, is designed to handle the second-order problem (1) directly. The detailed information can be seen in [4]. From then on, there have been many researchers focused on the RKN method. Subsequently, a lot of variations of RKN methods were given, such as [27,28,31,34] and others. The research on RK, RKN is tremendous, but it mainly focuses on explicit methods due to easier coding and less time consuming in comparison to implicit methods. The implicit methods are more suitable for solving stiff ODEs than the explicit methods. There are some researchers who work on the implicit RKN methods, such as [16-18, 22, 24]. If the solutions of ODEs satisfy a conservation law, such as dynamical systems for which total energy is conserved, the symplectic methods [8,9,30] should be considered. The term symplectic essentially means area preserving in a phase space. Approximate solutions generated by symplectic methods are conservative even at finite resolution, in contrast with numerical methods that generate approximate solutions that are conservative only in the limit as the time step size approaches zero. Symplectic methods have been applied to many problems such as pendulum, Morse oscillator, harmonic oscillator, Lennard-Jones oscillator, Kepler's orbit problem, and so on. As pointed out in Chap. V and Chap. XI of [13], symmetric methods show a better long time behavior than non-symmetric ones when solving the reversible differential system. So, some symmetric and symplectic RKN methods are proposed such as [24,34].
During the last thirty years, many researchers have been working on exponentially fitted RK or RKN methods. This technique was first analyzed in theory by Gautschi [12] and Lyche [21]. Exponentially fitted methods which intend to integrate exactly differential systems whose solutions can be expressed as linear combinations of functions from {exp(λt), exp(-λt), λ ∈ C}, or equivalently, from {sin(ωt), cos(ωt), ω ∈ R} with λ = iω, i 2 = -1, share better behaviors when applied to oscillatory problems than non-symplectic methods. Therefore, it has become an indispensable tool for solving oscillatory problems. The construction of exponentially fitted RK(N) methods is originally due to Paternoster [23], and a detailed exposition of exponentially fitted methods with an extensive bibliography on this subject can be found in Ixaru and Vanden Berghe [15].
Recently, the authors [35] have given a two-stage implicit symmetric and symplectic exponentially fitted Runge-Kutta-Nyström method (ISSEFRKN). It shows a good behavior compared with some existing methods. Exactly, this method is not a complete exponential fitting method. It can be seen from the process of derivation that there are two different expressions of b 1 . So the authors make them as close as possible by choosing a parameter θ = ± √ 3 6 . To avoid this happening, we investigate the construction of two-stage implicit symmetric and symplectic exponentially fitted modified Runge-Kutta-Nyström (ISSEFMRKN). Compared with the ISSEFRKN method, we add modified parameters for the term h in the internal stages. Consequently, we can obtain unique expression of every coefficient which is not true for ISSEFRKN. The new method ISSEFMRKN also reduces to the classical symplectic, symmetric RKN integrator when the parameter z approaches zero.
This paper is organized as follows: In Sect. 2 we present the notations and definitions to be used in the rest of the paper as well as some previous results on symmetric and symplectic methods. In Sect. 3 we make a study of the local truncation error, obtaining the order conditions (up to fifth order) for this class of methods. In Sect. 4, we derive the new two-stage implicit symplectic and symmetric EFMRKN integrator based on the former conditions. In Sect. 5, we devote to some numerical experiments. The numerical results show that the new method is more accurate and efficient compared with some other implicit RKN integrators. Finally, Sect. 6 involves in some conclusions.

Conditions for symmetry, symplecticity, exponential fitting of modified RKN methods
In this paper, we deduce a class of exponentially fitted RKN methods which integrate exactly second-order differential systems whose solutions can be expressed as linear combinations of the set of functions {exp(λt), exp(-λt), λ ∈ C}, or equivalently, from {sin(ωt), cos(ωt), ω ∈ R} with λ = iω, i 2 = -1. This means that the internal stage and the final stage have to integrate exactly these sets of functions. In order to do so, we must introduce some modifications to the ordinary RKN scheme. Here we consider the following s-stage modified implicit RKN method for the second-order ODEs (1): which can be expressed in the Butcher tableau as follows: It should be noted that scheme (2) coincides with the classical s-stage RKN formula when the coefficients γ i = 1, i = 1, . . . , s, and the remaining coefficients are constant. Now, we set out to derive the three corner stones to construct our method. In the following subsections, 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.

Symmetry conditions
The key to understanding symmetry is the concept of the adjoint method.
In this paper we consider scheme (2) whose coefficients are functions of z, as we do for exponentially fitted type methods [32,33]. Then the conditions for methods (2) to be symmetric are given by the following lemma. (2) is symmetric if its coefficients satisfy the following conditions:

Lemma 1 The modified RKN method
where z = iωh, ω is the principal frequency of the problem.
Proof Following the idea of Hairer et al. [14], (2), then we obtain the following equations: Note that, we do not change z to -z, because we will derive exponentially fitted methods for which its coefficients are even functions of z. From equation (4), one can easily obtain Inserting equations (5) and (6) into Y s+1-i , we have Comparing equations (7), (6), (5) with the counterpart in (2) respectively, we can obtain the symmetric conditions (3).
Omitting the variable z, i.e., all the coefficients are constants, they reduce to symmetric conditions for the traditional RKN methods.

Symplecticity conditions
Now, we turn to the symplecticity conditions for scheme (2). Symplecticity is defined for a Hamiltonian system. On many occasions, the problem under consideration takes the form of a Hamiltonian systeṁ where M is a symmetric positive definite constant matrix. This system is equivalent to the second-order equation (1) To the end of obtaining the symplectic conditions for (2), the following definition, which can be found in [14], is essential. 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 where the one-forms dp i , respectively dq i , map a tangent vector ξ to its ith, respectively (n + i)th, component. Here, we assume p and q all have n components. Furthermore, dp i ∧ dq i is a bilinear map acting on a pair of vectors and satisfies Grassmann's rules for exterior multiplication dp i ∧ dp j = -dp j ∧ dp i , dp i ∧ dp i = 0.
Hamiltonian systems have been seen to possess two remarkable properties: (a) the solutions preserve the Hamiltonian H(p, q); (b) the corresponding flow is symplectic, i.e., preserves the differential 2-form dp ∧ dq.
As is pointed out by Feng [7], "It is natural to look forward to those discrete systems which preserve as much as possible the intrinsic properties of the continuous system. " Based on this definition, we can easily obtain the symplectic conditions for RKN formula (2). (2) is symplectic if the following conditions are satisfied:
For the left-hand side of this equation, we have Inserting y 0 in (2) to the second term of this equation, we obtain Therefore, considering property (8), we can obtain that (10) holds if conditions (9) are satisfied.

Exponential fitting conditions
Following Albrecht's approach [2,3], each stage of scheme (2) can be viewed as a linear multi-step method on a non-equidistant grid. With each stage one can associate a linear function as follows: • for the internal stages, a ij y (t + c j h), i = 1, 2, . . . , s; • for the final stages, The exponentially fitted conditions can be obtained by requiring that these functions vanish for the functions from the set {exp(±λt)|λ ∈ R or λ ∈ iR}. The exponentially fitted conditions are stated in the following lemma. (2) is exponentially fitted if the following conditions are satisfied:

Lemma 3 The modified RKN method
Proof Requiring the internal stages vanishing when y(t) taken as e λt , e -λt respectively, we have Denoting z = λh, (12) leads to ⎧ ⎨ ⎩ e c i z -1c i γ i zz 2 s j=1 a ij e c j z = 0, e -c i z -1 + c i γ i zz 2 s j=1 a ij e -c j z = 0.
Similarly, for the final stages, we have It follows that Together with (14) and (15), we obtain the desired results.

Algebraic order conditions
As we all know, a numerical method having higher algebraic order will have higher accuracy. To the end of specifying the order of the new method, we will present algebraic order conditions for exponentially fitted modified Runge-Kutta-Nystöm (EFMRKN) methods. As it occurs in the case of a classical RKN method, for an EFMRKN method, the local truncation errors in the approximations of the solution and its derivative can be expressed as where F (j) (y 0 ) denotes an elementary differential and the terms d  For the EFRKN method, the coefficients are dependent on z. Specifically, they are even functions of z. Thus, we use the following Taylor expansions: Using these assumptions and following the way given in [14] (pp. 143-148), we obtain the terms of the local truncation error. Roughly speaking, comparing the Taylor expansions of y 1 , y 1 in (2) with and the Taylor expansions of y(t 0 + h), y (t 0 + h) at t 0 , the order k condition can be obtained by making the coefficients of h i , i = 1, . . . , k, equal. Following this procedure, the order conditions (up to the fifth order) for the EFMRKN methods considered in this paper can be obtained. Order 1 requires: Order 2 requires in addition: Order 3 requires in addition: Order 4 requires in addition: Order 5 requires in addition:

Numerical experiments
To test the numerical performance of the method ISSEFMRKN2, we carry out experiments on four problems to illustrate the effectiveness and efficiency. The codes used for comparison are • DIRKNRaed: The embedded diagonally implicit RKN 4(3) pair method proposed by Al-Khasawneh et al. in [1]. • DIRKNNora: The three-stage fourth-order diagonally implicit RKN method proposed by Senu et al. in [27]. • ISSRKN2: The symmetric and symplectic two-stage fourth-order implicit RKN method proposed by MENG-ZHAO QIN et al. in [24] with a 11 = -13+160θ 2 +720θ 4 2880θ 2 and θ = ± √ 3 6 , i.e., a 11 = 1 45 . • ISSEFRKN2: The symmetric and symplectic exponentially fitted two-stage RKN method proposed in [35] which is of order 4. • ISSEFMRKN2: The symmetric and symplectic exponentially fitted two-stage modified RKN method (25) proposed in this paper which is of order 4. As we can see, the proposed method ISSEFMRKN2 is implicit. The main computation is in computing the non-linear equations In this paper, we use the Newton iteration method with initial values Y (0) 1 = Y (0) 2 = y(0). Here we use two stopping criteria for iteration. First, iterations are carried out until the difference between the Euclidean norm of two successive iterations attains 10 -8 . Second, if the first criterion fails, then the iteration will be forced to terminate after running 1000 times.
This solution represents a periodic motion with different frequencies. In our test we choose the parameter values λ = 10i, t end = 10, ε = 10 -3 , and the numerical results stated in Fig. 5 have been computed with steps h = 1/2 m , m = 3, 4, 5, 6.
From Figs. 1-5, we can find that the implicit modified EFRKN method ISSEFMRKN2 is more efficient than ISSEFRKN2 and the symmetric and symplectic method ISSRK2. ISSEFRKN2 does not possess higher accuracy than ISSRKN2 for Problems 3 and 4. For Problem 3, its true frequency is 1. In the numerical study, we select two frequencies 1.2 ( Fig. 3) and 1 (Fig. 4). From Figs. 3-4, we can see the accuracies are quite different. The accuracy of 1 is much higher than that of 1.2. In this problem, we know its true frequency. But when it comes to applications, the true frequency is often unpredictable. Therefore, we need to try some different candidates. For Problems 2 and 3, we find that ISSEFMRKN2 is much more accurate and efficient than our methods considered in this paper. As is pointed out in introduction, ISSEFRKN2 is not a complete EF method. However, ISSEFMRKN2 is completely EF. The solutions of Problems 2 and 3 are all in the form of a triangular function. This is just in line with EF methods. So, ISSEFMRKN2 performs very well.

Conclusions
In this paper a two-stage symmetric, symplectic IEFMRKN integrator has been derived. Like the existing EFRKN integrators such as [34,35], the coefficients of the new method depend on the product of the dominant frequency ω and the step size h. When the parameter z approaches zero, the ISSEFMRKN2 method reduces to the classical RKN method.
The numerical experiments carried out show that the new method is more efficient than some implicit RKN methods. In every experiment, the method ISSEFMRKN2 is shown to be the most efficient one among the methods used for comparison. However, like the ISSEFRKN2 method in [35], we derive only one method, not a class of methods whose coefficients can be dependent on one or more parameters. In the future, we will consider deriving a class of ISSEFRKN methods.