Symmetric and symplectic ERKN methods for oscillatory Hamiltonian systems

https://doi.org/10.1016/j.cpc.2011.09.002Get rights and content

Abstract

The ERKN methods proposed by H. Yang et al. [Comput. Phys. Comm. 180 (2009) 1777] are an important improvement of J.M. Francoʼs ARKN methods for perturbed oscillators [J.M. Franco, Comput. Phys. Comm. 147 (2002) 770]. This paper focuses on the symmetry and symplecticity conditions for ERKN methods solving oscillatory Hamiltonian systems. Two examples of symmetric and symplectic ERKN (SSERKN) methods of orders two and four respectively are constructed. The phase and stability properties of the new methods are analyzed. The results of numerical experiments show the robustness and competence of the SSERKN methods compared with some well-known methods in the literature.

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{q+ω2q=f(t,q),t[t0,T],q(t0)=q0,q(t0)=q0, where q:RRd is a vector-valued function, f(t,q) is a small perturbing force. This kind of problems are known to have oscillatory solutions. Here, ω>0, the principal frequency of the problem, is assumed to be known or can be accurately estimated in advance. In applications the function f(t,q) often satisfy f(t,q)ω2q=MVq(t,q) for some positive definite constant d×d matrix M and some smooth function V(t,q). Then the problem (1) can be considered as the IVP of a Hamiltonian system in the form{q=Mp,p=Vq(t,q),q(t0)=q0,p(t0)=p0. In mechanics, q:RRd and p:RRd are generalized positions and generalized momenta, respectively, M is a symmetric and positive definite matrix. The Hamiltonian of the system (2) is H(t,p,q)=12pTMp+V(t,q). 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 q+ω2q=0. 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 q=g(t,q) (particularly for the system (2) or (1) g(t,q)=f(t,q)ω2q=MV(t,q)) by introducing suitable functions of frequencies into the coefficients and requiring that the resulting method{Qi=β¯i(ν)qn+γ¯i(ν)cihqn+h2j=1sa¯ij(ν)g(tn+cjh,Qj),i=1,,s,qn+1=β¯s+1(ν)qn+γ¯s+1(ν)hqn+h2i=1sb¯i(ν)g(tn+cih,Qi),qn+1=γs+1(ν)qn+hi=1sbi(ν)g(tn+cih,Qi) integrate exactly the differential equations whose solutions are linear combinations of functions exp(λt) and exp(λt), or equivalently, cosωt and sinωt when λ=iω, i2=1. 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 f(q)=εg(q), ε1, Franco following Bettisʼ line [11], proposes to modify the updates by introducing ϕ-functions ϕ0(ν),ϕ1(ν) and obtained his ARKN methods (Runge–Kutta–Nyström methods adapted to the numerical integration of perturbed oscillators), where ν=ωh, ϕ0(ν)=cos(ν) and ϕ1(ν)=sin(ν)/ν (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 ω2q 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 q reads{Qi=qn+cihqn+h2j=1sa¯ij(f(tn+cjh,Qj)ω2Qj),i=1,,s,qn+1=qn+hqn+h2i=1sb¯i(f(tn+cih,Qi)ω2Qi),qn+1=qn+hi=1sbi(f(tn+cih,Qi)ω2Qi),n=0,1,,N1, where h is the step size, N is the number of computation steps, and a¯ij,b¯i,bi,i,j=1,,s, 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{b¯T(ν)Φ(τ)ρ(τ)!γ(τ)ϕρ(τ)+1(ν)=O(hpρ(τ)),ρ(τ)=1,,p1,bT(ν)Φ(τ)ρ(τ)!γ(τ)ϕρ(τ)(ν)=O(hp+1ρ(τ)),ρ(τ)=1,,p, where τ is the Nyström tree corresponding to an elementary differential F(q,q) of the function f˜(q):=f(q)ω2q at (qn,qn), 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 equationy(t)=λ2y(t),λ>0. Here, for our SSERKN methods, we follow the approach of Franco [14] and use the following test equationy(t)+ω2y(t)=ε2y(t), where ε>0 represents the error of the estimate ω of the principal frequency λ. Applying an SSERKN method to Eq. (23) yields{Y=ϕ0(cν)yn+cϕ1(cν)hynz2A¯(ν)Y,yn+1=ϕ0(ν)yn+ϕ1(ν)hyn

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)

  • M. Calvo et al.

    J. Comput. Appl. Math.

    (2008)
  • A. Tocino et al.

    Math. Comput. Modell.

    (2005)
  • J.M. Franco

    Comput. Phys. Comm.

    (2002)
  • J.M. Franco

    J. Comput. Appl. Math.

    (2003)
  • J.M. Franco

    Appl. Numer. Math.

    (2006)
  • H. Yang et al.

    Comput. Phys. Comm.

    (2009)
  • M.P. Calvo et al.

    J. Comput. Appl. Math.

    (2009)
  • X. Wu et al.

    Comput. Phys. Comm.

    (2009)
  • H. Van de Vyver

    Comput. Phys. Comm.

    (2005)
  • H. Van de Vyver

    New Astron.

    (2005)
  • J.M. Franco

    Comput. Phys. Comm.

    (2007)
  • H. Van de Vyver

    Phys. Lett. A

    (2007)
  • Cited by (0)

    View full text