Symmetric and symplectic ERKN methods for oscillatory Hamiltonian systems
Introduction
Hamiltonian systems constitute an important category of differential equations which arise in a variety of applied fields such as classical and quantum mechanics, astronomy, chemistry and biology. These systems have some pronounced structures like the invariance of Hamiltonian and symplecticity of flow, which require special features of integrators. Until today among the most popular methods designed deliberately for Hamiltonian systems are symplectic integrators of Runge–Kutta type or Runge–Kutta–Nyström type which can preserve the symplecticity of the exact solutions (see [1], [2], [3]). On the other hand, in a multitude of applications like the mathematical pendulum and the two-body problem, the Hamiltonian systems of interest have oscillatory solutions. One is then motivated to construct numerical integrators which are symplectic and can preserve the oscillatory properties of the exact solutions simultaneously.
This paper is concerned with the time integration of initial value problems (IVPs) of systems of second order differential equations in the form where is a vector-valued function, is a small perturbing force. This kind of problems are known to have oscillatory solutions. Here, , the principal frequency of the problem, is assumed to be known or can be accurately estimated in advance. In applications the function often satisfy for some positive definite constant matrix M and some smooth function . Then the problem (1) can be considered as the IVP of a Hamiltonian system in the form In mechanics, and are generalized positions and generalized momenta, respectively, M is a symmetric and positive definite matrix. The Hamiltonian of the system (2) is . For time integration of the system (1), Runge–Kutta–Nyström (RKN) methods, which are universally adopted for second order differential equations (see [1], [2], [4]), are of course applicable directly. Unfortunately, the general-purpose RKN methods fail to take into account the oscillatory character possessed by many Hamiltonian systems like the problem (1) or (2).
One fact that motivates one to consider symplectic methods is that symplectic methods are zero-dissipative which is a prominent property for solving oscillatory problems. Simos and Vigo-Aguiar [5] introduce a symplectic modified RKN method of order two with frequency-dependent coefficients which integrates exactly the harmonic oscillators . The method is proved to be much more efficient than the classical symplectic RKN method of order two given by Calvo and Sanz-Serna [6].
Exponential fitting is a recently developed technique for dealing with oscillatory problems that are shown more effective and more efficient than traditional approaches. A systematic exposition with an extensive bibliography on this subject can be found in Ixaru and Vanden Berghe [7]. Calvo et al. [8] give a complete characterization of EFRK methods that preserve linear or quadratic invariants and derive a family of symplectic EFRK two-stage methods of order four. The basic idea behind the exponentially fitted methods is modifying traditional RKN methods for the system (particularly for the system (2) or (1) ) by introducing suitable functions of frequencies into the coefficients and requiring that the resulting method integrate exactly the differential equations whose solutions are linear combinations of functions and , or equivalently, and when , . In the pioneering work [9], Tocino and Vigo-Aguiar give the symplectic conditions for a general class of exponentially fitted Runge–Kutta–Nyström methods (EFRKN) solving second order differential equations that are equivalent to a special form of Hamiltonian systems, but they fail to give a practical symplectic EFRKN integrator. Van de Vyver [10] derives sufficient conditions on symplecticity of EFRK methods and constructs an effective fourth order symplectic exponentially fitted integrator.
We note that the approach of exponential fitting does not take into account the particular structure of the Eq. (1). Considering (1) as a perturbed oscillatory system with the perturbing force of the form , , Franco following Bettisʼ line [11], proposes to modify the updates by introducing ϕ-functions and obtained his ARKN methods (Runge–Kutta–Nyström methods adapted to the numerical integration of perturbed oscillators), where , and (see [12], [13], [14]). Based on Francoʼs work, Yang et al. [15] make full use of the structure of the system (1) brought by the term to improve ARKN methods and obtained more effective ERKN (extended RKN) methods, which are naturally of exponentially-fitted type.
On the other hand, conservative mechanical systems are time reversible. Accordingly symmetric integrators are then designed to retain this significant property in the long-term numerical integration. For instance, Calvo, Franco, Montijano and Rández [16] construct three-stage sixth-order symmetric and symplectic EFRK methods of Gauss type which can preserve linear and quadratic invariants.
The purpose of this paper is to investigate the symmetry and symplecticity conditions for ERKN methods when they are applied to Hamiltonian systems in the form (2) or (1). The remainder of this paper is organized as follows. In Section 2 we briefly recall the ideas in the formulation of the ERKN methods solving the system (1). Section 3 derives the symmetry and symplecticity conditions for ERKN methods. In Section 4 we construct two symmetric and symplectic ERKN methods of orders two and four. Section 5 examines the stability and phase properties of the new methods. In Section 6 we carry out some numerical experiments to check the effectiveness and efficiency of the new methods comparing with some highly efficient codes in the literature. Section 7 is devoted to conclusions.
Section snippets
Formulation of ERKN methods
An s-stage RKN method solving the perturbed oscillators (1) with the perturbing force f independent of the derivative reads where h is the step size, N is the number of computation steps, and , are constants. Considering the oscillatory character of the solutions of the problem (1), Franco [12] incorporates the special
Order conditions
As explained in Yang et al. [15], in order to derive the order conditions of ARKN or ERKN methods, one need only to consider the autonomous equations. Theorem 3.1 An ARKN method for the system (1) is of order p if and only if where τ is the Nyström tree corresponding to an elementary differential of the function at , the functions (the density of τ) and (the
Construction of SSERKN methods
For convenience, we restrict ourselves to explicit methods.
Stability and phase properties of the new SSERKN methods
Stability analysis for a new integrator is an important step before it is put into use. To this end, one usually applies the method to the second order homogeneous linear equation Here, for our SSERKN methods, we follow the approach of Franco [14] and use the following test equation where represents the error of the estimate ω of the principal frequency λ. Applying an SSERKN method to Eq. (23) yields
Numerical experiments
In this section, in order to test the effectiveness of the new SSERKN methods derived in Section 4, we apply these methods as well as some other codes from the literature to six practical model problems. The numerical integrators we choose for comparison are as follows:
- •
SSERKN2: the two-stage SSERKN method of order two derived in Section 4 of this paper;
- •
SSERKN4: the three-stage SSERKN method of order four derived in Section 4 of this paper;
- •
SEFRKN2: the two-stage symplectic exponentially fitted
Conclusions
The ERKN methods proposed by Yang et al. [15] form an essential improvement of the ARKN methods of Franco. We have investigated the symmetry and symplecticity of ERKN. A critical feature of the ERKN methods is that they are exact for the unperturbed harmonic oscillator in both internal stages and the updates. This makes the ERKN methods promising when applied to perturbed oscillators. Conditions for ERKN methods to be symmetric and symplectic are formulated. Symplectic ERKN methods are shown to
Acknowledgements
The authors are grateful to the anonymous referees for their careful reading of the manuscript and for their many invaluable comments and suggestions, which largely help to improve this paper. We thank Professor Xinyuan Wu and Professor Jinxi Zhao for their beneficial discussions and continual encouragement. The work was supported by Youth Innovation Fund of Naning Agricultural University (grant No. KJ09023), the Fundamental Research Fund for the Central Universities (grant No. Y0201100265) and
References (24)
- et al.
J. Comput. Appl. Math.
(2008) - et al.
Math. Comput. Modell.
(2005) Comput. Phys. Comm.
(2002)J. Comput. Appl. Math.
(2003)Appl. Numer. Math.
(2006)- et al.
Comput. Phys. Comm.
(2009) - et al.
J. Comput. Appl. Math.
(2009) - et al.
Comput. Phys. Comm.
(2009) Comput. Phys. Comm.
(2005)New Astron.
(2005)