Original ArticleSplitting schemes and energy preservation for separable Hamiltonian systems☆
Introduction
We consider conservative Hamiltonian differential equations , where 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 . 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 . Let φt denote the exact flow, i.e. x(t) = φt(x0), and , represent the exact flows associated to and . We are going to focus attention on splitting and composition methods which are defined as discrete mappings xn = ψh(xn−1) whereand the coefficients ensure the order p of accuracy on the approximation i.e. (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 whereThe 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
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 systemwith respect to b2 and a1, considering ρ = b1 as free parameter. We obtain the class of second order methods defined as whereIn [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 asat tn = t0 + n h, where 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 followsLet us evaluate the numerical solution where is the second order method defined in (4) withwhere ph(ρ),eh(ρ) and qh(ρ) are the polynomials with respect to the parameter ρLet us evaluate the error in energy as
Non linear oscillator
Consider as a first non linear example the simple oscillator, which is described by the Hamiltonian . 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:The 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)
- et al.
Splitting and composition methods for explicit time dependence in separable dynamical systems
Journal of Computational and Applied Mathematics
(2010) - et al.
Practical symplectic partitioned Runge–Kutta and Runge–Kutta–Nyström methods
Journal of Computational and Applied Mathematics
(2002) - et al.
A note on the efficient implementation of Hamiltonian BVMs
Journal of Computational and Applied Mathematics
(2011) - et al.
Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics
Computer Methods in Applied Mechanics and Engineering
(1992) - et al.
Energy- and quadratic invariants-preserving integrators based upon Gauss collocation formulae
SIAM Journal on Numerical Analysis
(2012) - et al.
Hamiltonian boundary value methods (energy preserving discrete line integral methods)
Journal of Numerical Analysis Industrial and Applied Mathematics
(2010) - et al.
Energy-preserving integrators and the structure of B-series
Foundations of Computational Mathematics
(2010) - et al.
Total variation in Hamiltonian formalism and symplectic-energy integrators
Journal of Mathematical Physics
(2003) - et al.
Long-time behaviour of discretizations of the simple pendulum equation
Journal of Physics A: Mathematical and Theoretical
(2009) - et al.
Improving the accuracy of the discrete gradient method in the one-dimensional case
Physical Review E
(2010)
Symplectic methods based on the matrix variational equation for Hamiltonian system
Lecture Notes in Computer Science
Explicit symplectic partitioned Runge–Kutta–Nystrom methods for non-autonomous dynamics
Applied Numerical Mathematics
The symplectic methods for the computation of Hamiltonian equations
Lecture Notes in Mathematics
Symplectic Geometric Algorithms for Hamiltonian Systems
Formal power series and numerical methods for differential equations
Lie–Poisson–Hamilton–Jacobi theory and Lie–Poisson integrators
Physics Letters A
Cited by (4)
A dissipation-preserving scheme for damped oscillatory Hamiltonian systems based on splitting
2021, Applied Numerical MathematicsIntegrating dissipative particle dynamics with energy conservation
2020, Physical Review EGeometric numerical integration in ecological modelling
2020, Mathematics
- ☆
This work has been partially supported by National Research Council grant CNR-RSTL id.332/2008.