Original Article
Splitting schemes and energy preservation for separable Hamiltonian systems

https://doi.org/10.1016/j.matcom.2013.11.002Get rights and content

Abstract

It is known that symplectic algorithms do not necessarily conserve energy even for the harmonic oscillator. However, for separable Hamiltonian systems, splitting and composition schemes have the advantage to be explicit and can be constructed to preserve energy. In this paper we describe and test an integrator built on a one-parameter family of symplectic symmetric splitting methods, where the parameter is chosen at each time step so as to minimize the energy error. For second-degree polynomial Hamiltonian functions as the one describing the linear oscillator, we build up second and fourth order symmetric methods which are symplectic, energy-preserving and explicit. For non-linear examples, it is possible to construct schemes with minimum error on energy conservation. The methods are semi-explicit in the sense that they require, as additional computational effort, the search for a zero of a scalar function with respect to a scalar variable. Therefore, our approach may represent an effective alternative to energy-preserving implicit methods whenever multi-dimensional problems are dealt with as is the case of many applications of interest.

Introduction

We consider conservative Hamiltonian differential equations x˙=J1H(x), where J=OInIn0 is a skew-symmetric constant matrix, In is the n-dimensional identity matrix, the vector x = (p, q)T provides the positions and momenta variables q,pn. Symplectic methods for the numerical treatment of Hamiltonian systems have been extensively studied in literature (see, for example [32], [10], [28], [13] and all the references therein contained). When the Hamiltonian function H(p, q) = T(p) + V(q) is separable in two parts T(p) and V(q), one depending on p and the other on q, the vector field of the differential system can be decomposed in two parts, each of them exactly integrable. In this case, splitting and composition methods [2], [25] represent symplectic methods which follow the solution providing, in the long run, a nearly-preservation of the energy.

Since the exact preservation of the energy is a desirable property when considering a numerical integrator of an Hamiltonian system, it would be optimal to have symplectic energy-preserving methods. However, symplectic algorithms do not necessarily retain energy even for the harmonic oscillator (see, for example [12]). In the general case of non linear Hamiltonian, the backward error analysis proves that it is not possible to have a symplectic method which exactly preserve the energy [17], [14] at least for the fixed time steps [15]. To attain this goal, most use variable time steps and a discrete total variation calculus developed in [23], [24], [21], [7]. Their basic idea was to construct a discrete action functional with variable time steps and then apply a discrete total variation calculus. In this way, two-step symplectic integrators and their associated energy conservation laws were derived. However, for fixed time steps, the resulting integrators are equivalent to the symplectic integrators derived directly from the Hamiltonian systems in some special cases and do not preserve the energy anymore.

Several energy-preserving with fixed time steps (hence, not symplectic) methods have been proposed in literature: the discrete gradient method [30], [26], its modified versions methods for one-dimensional Hamiltonian [9], [8], energy-preserving B-series methods [29], [6], extended collocation methods (see [19], [20]) as well an energy-preserving variant of collocation methods [18], [4]. However, since these scheme are fully implicit, they could be less competitive with respect to semi-explicit integrators if applied to multi-dimensional Hamiltonian systems, including high oscillatory systems as Fermi–Pasta–Ulam and Toda lattice, N-body problems arising in molecular dynamics and in outer solar systems, as well as multi species population models.

Splitting and composition schemes have the advantage to be explicit schemes for separable Hamiltonian systems; however, lying in the class of symplectic integrators, they are affected by an error in the preservation of the energy. We briefly recall the idea on which splitting and composition schemes are built. Let f[A](x) = (Tp(p), 0) and f[B](x) = (0, − Vq(q)), where Tp and Vq are the gradients of T and V with respect to p and q, respectively. Then, the evolution of the variable x can be written as x˙=f[A](x)+f[B](x). Let φt denote the exact flow, i.e. x(t) = φt(x0), and φt[A], φt[B] represent the exact flows associated to x˙=f[A](x) and x˙=f[B](x). We are going to focus attention on splitting and composition methods which are defined as discrete mappings xn = ψh(xn−1) whereψh=φbs+1h[B]φash[A]φbsh[B]φas1h[A]φb2h[B]φa1h[A]φb1h[B]and the coefficients ai,bi ensure the order p of accuracy on the approximation i.e. xn=ψh(xn1)+O(hp+1) (see e.g. [1]).

In this paper we propose to compose symplectic maps of the same order each defined by coefficients that depends on a parameter ρn i.e. ai = ai(ρn), bi = bi(ρn) and to formulate methods in the form xn=ψh(ρn)(xn1) whereψh(ρn)=φbs+1(ρn)h[B]φas(ρn)h[A]φbs(ρn)h[B]φas1(ρn)h[A]φb2(ρn)h[B]φa1(ρn)h[A]φb1(ρn)h[B]The aim is to find, in a one-parameter family of methods of a given order, the method defined by the parameter ρn which minimizes, at each timestep tn, the absolute error on the preservation of the energy |En|, where En = H(xn)  H(x0). The idea of considering parametric generalizations of existing symplectic formulae is not new in the literature. Among the first attempts are [31], [21], while more recent approaches have been investigated in [3], [22].

The advantage of our approach is the possibility to define schemes based on the composition of (different) symplectic maps featured by minimum error on energy. These schemes are semi-explicit since require, at each step, the search for a zero of a scalar function |En| with respect to the scalar variable ρn, independently of the dimension of the Hamiltonian systems. We will show by numerical tests that such schemes share good qualitative properties and turns out to be of great interest whenever high dimensional systems have to be solved and extremely small errors in energy preservation (smaller than the ones provided by symplectic schemes) are required. Indeed, even if the introduction of a parameter in a symplectic formula in order to attain energy conservation, destroys the original symplecticity property of the underlying method, there could be properties which are a consequence of symplecticity that are not lost in the parametric variant of the method. Among that, the most important is the conservation of any quadratic invariant. So, for example, the conservation property of the total angular momentum of a mechanical isolated system, is surely satisfied by the new formulae, since for any fixed choice of the parameter, the underlying method is symplectic.

In the following we restrict our search to the class of symmetric methods. In this case, the 2s + 1 free parameters reduce to s + 1 parameters, since by satisfying odd orders of accuracy p, the conditions for the even order p + 1 are automatically satisfied.

Section snippets

Second order methods

For a second order method, two conditions have to be imposed

i=1s+1bi=1,i=1sai=1

It is enough to take s = 1 and solve a1 = 1 and 2b1 = 1 with b1 = b2 for defining the well-known Störmer–Verlet scheme. In order to preserve one degree of freedom, we take s = 2 and we solve the system2b1+b2=1,2a1=1.with respect to b2 and a1, considering ρ = b1 as free parameter. We obtain the class of second order methods defined as xn+1=ψh(ρ)(xn) whereψh(ρ)=φρh[B]φh/2[A]φ(12ρ)h[B]φh/2[A]φρh[B]In [25] the author

Splitting methods with minimum error in energy-preservation

In order to preserve the energy during the evolution, instead of minimizing the truncation error, we search for the value of ρ which minimizes the error in energy at each timestep tn. In more details, starting from x0 = (p0, q0) = (p(t0), q(t0)), for n = 1, …, N, we define ρn asρn=argmin|H(xn)H(x0)|at tn = t0 + n h, where xn=ψh(ρn)(xn1) is given in (2). Whenever it is possible to find a zero as a minimum for the error in energy, at each tn, the method provides an energy-preserving integrator. For second

Second order methods

The harmonic oscillator is defined by a quadratic Hamiltonian as followsx˙=J1H(x),H(x)=xTx2,x2.Let us evaluate the numerical solution xn+1=ψh(ρ)(xn) where ψh(b) is the second order method defined in (4) x1=ψh(ρ)(x0)=Kh(ρ)x0 withKh(ρ)=18ph(ρ)eh(ρ)+qh(ρ)eh(ρ)qh(ρ)ph(ρ)where ph(ρ),eh(ρ) and qh(ρ) are the polynomials with respect to the parameter ρph(ρ)=84h2+2h4ρ4h4ρ2eh(ρ)=h3(2h2ρ3(4+h2)ρ2+6ρ1)qh(ρ)=h(8h2)+2h3ρh3(4+h2)ρ2+2h5ρ3.Let us evaluate the error in energy asE=12x0T(Kh(ρ)TKh(ρ)I)x

Non linear oscillator

Consider as a first non linear example the simple oscillator, which is described by the Hamiltonian H(p,q)=12p2cos(q). In order to find ρn = argmin|En|, where En = H(pn, qn)  H(p0, q0), we use the Matlab code lsqnonlin, providing the derivative of the objective function with respect to the parameter ρn as described in (7). With abuse of notation we will denote with S2 and S4 the resulting methods of second and fourth order, respectively. As starting guesses we set ρ0 = ρML2 for S2 and ρ0 = ρopt4 for S4

Multi-dimensional Hamiltonian system: FUP model

As third example we choose a multi-dimensional Hamiltonian problem provided by the Fermi Ulam Pasta β-lattice as described in [16]. The motion of N interacting particles evolves according to the Hamilton's equations with Hamiltonian given by:H=i=1Npi22+i=0N(qi+1qi)22+β(qi+1qi)44The boundary values are given by q0 = qN+1 = 0 and the parameter β is set at 1.8. The initial conditions are set at qi = 0.1 and pi = 0, for i = 1, …, N.

In the first test we set N = 4 and we integrate the system in the interval

Conclusion

We propose second and fourth order simmetyric splitting schemes for separable Hamiltonians, obtained by composing, at each step, symplectic maps which depend on a parameter. The parameter is chosen in order to have methods which minimize the error on the Hamiltonian function of the system. In the case of the harmonic oscillator we prove that they are explicit, symplectic and energy-preserving integrators.

In the case of simple non-linear pendulum the schemes preserve the energy when small

Acknowledgements

We thank the anonymous reviewers for their helpful comments that contributed to improve our paper.

References (32)

  • N. Del Buono et al.

    Symplectic methods based on the matrix variational equation for Hamiltonian system

    Lecture Notes in Computer Science

    (2002)
  • F. Diele et al.

    Explicit symplectic partitioned Runge–Kutta–Nystrom methods for non-autonomous dynamics

    Applied Numerical Mathematics

    (2010)
  • K. Feng et al.

    The symplectic methods for the computation of Hamiltonian equations

    Lecture Notes in Mathematics

    (1987)
  • K. Feng et al.

    Symplectic Geometric Algorithms for Hamiltonian Systems

    (2010)
  • K. Feng

    Formal power series and numerical methods for differential equations

  • Z. Ge et al.

    Lie–Poisson–Hamilton–Jacobi theory and Lie–Poisson integrators

    Physics Letters A

    (1988)
  • Cited by (4)

    This work has been partially supported by National Research Council grant CNR-RSTL id.332/2008.

    View full text