Symplectic integrator

In mathematics, a symplectic integrator (SI) is a numerical integration scheme for a specific group of differential equations related to classical mechanics and symplectic geometry. Symplectic integrators form the subclass of geometric integrators which, by definition, are canonical transformations. They are widely used in molecular dynamics, discrete element methods, accelerator physics, and celestial mechanics.

Introduction

Symplectic integrators are designed for the numerical solution of Hamilton's equations, which read

$\dot p = -\frac{\partial H}{\partial q} \quad\mbox{and}\quad \dot q = \frac{\partial H}{\partial p},$

where q denotes the position coordinates, p the momentum coordinates, and H is the Hamiltonian. The set of position and momentum coordinates (q,p) are called canonical coordinates. (See Hamiltonian mechanics for more background.)

The time evolution of Hamilton's equations is a symplectomorphism, meaning that it conserves the symplectic two-form $dp \wedge dq$. A numerical scheme is a symplectic integrator if it also conserves this two-form.

Symplectic integrators possess as a conserved quantity a Hamiltonian which is slightly perturbed from the original one. By virtue of these advantages, the SI scheme has been widely applied to the calculations of long-term evolution of chaotic Hamiltonian systems ranging from the Kepler problem to the classical and semi-classical simulations in molecular dynamics.

Most of the usual numerical methods, like the primitive Euler scheme and the classical Runge-Kutta scheme, are not symplectic integrators.

Splitting methods for separable Hamiltonians

A widely used class of symplectic integrators is formed by the splitting methods.

Assume that the Hamiltonian is separable, meaning that it can be written in the form

$H(p,q) = T(p) + V(q). \qquad\qquad (1)$

This happens frequently in Hamiltonian mechanics, with T being the kinetic energy and V the potential energy.

For the notational simplicity, let us introduce the symbol z = (q,p) to denote the canonical coordinates including both of the position and momentum coordinates. Then, the set of the Hamiltonian's equations given in the introduction can be expressed in a single expression as

$\dot{z}=\{z,H(z)\} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2)$

where $\{\cdot, \cdot\}$ is a Poisson bracket. Furthermore, by introducing an operator, $D_H = \{\cdot, H\}$, which returns a Poisson bracket of the operand with the Hamiltonian, the expression of the Hamilton's equation can be further simplified to

$\dot{z}=D_H z.$

The formal solution of this set of equations is given as

$z(\tau)=\exp(\tau D_H)z(0). \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3)$

When the Hamiltonian has the form of eq. (1), the solution (3) is equivalent to

$z(\tau) = \exp[\tau (D_T + D_V)]z(0). \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (4)$

The SI scheme approximates the time-evolution operator exp[τ(D_T + D_V)] in the formal solution (4) by a product of operators as

$\exp[\tau (D_T + D_V)] = \prod_{i=1}^k \exp(c_i \tau D_T)\exp(d_i \tau D_V) + O(\tau^{k+1}), \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (5)$

where c_i and d_i are real numbers, and k is an integer, which is called the order of the integrator. Note that each of the operators exp(c_iτD_T) and exp(d_iτD_V) provides a symplectic map, so their product appearing in the right hand side of (5) also constitutes a symplectic map. In concrete terms, exp(c_iτD_T) gives the mapping

$\begin{pmatrix} q\\ p \end{pmatrix} \mapsto \begin{pmatrix} q'\\ p' \end{pmatrix} = \begin{pmatrix} q + \tau c_i \frac{\partial T}{\partial p}(p)\\ p \end{pmatrix}$

and exp(d_iτD_V) gives

$\begin{pmatrix} q\\ p \end{pmatrix} \mapsto \begin{pmatrix} q'\\ p' \end{pmatrix} = \begin{pmatrix} q \\ p - \tau d_i \frac{\partial V}{\partial q}(q)\\ \end{pmatrix}.$

Note that both of these maps are practically computable.

The symplectic Euler method is the first-order integrator with k = 1 and coefficients

$c_1 = d_1 = 1. \ \$

The Verlet method is the second-order integrator with k = 2 and coefficients

$c_1 = c_2 = \tfrac12, \qquad d_1 = 1, \qquad d_2 = 0.$

A third order symplectic integrator (with k = 3) was discovered by Ronald Ruth in 1983. One of the many solutions is given by

$c_1 = \tfrac{2}{3},\qquad c_2 =-\tfrac{2}{3},\qquad c_3 = 1,$

$d_1 = \tfrac{7}{24},\qquad d_2 = \tfrac{3}{4},\qquad d_3 =-\tfrac{1}{24}.$

A fourth order integrator (with k = 4) was also discovered by Ruth in 1983 and distributed privately to the accelerator community at that time. This was described in a lively review article by Forest. This fourth order integrator was published in 1990 by Forest and Ruth and also independently discovered by two other groups around that same time

$c_1 = c_4 = \frac{1}{2(2-2^{1/3})},\ \ c_2=c_3=\frac{1-2^{1/3}}{2(2-2^{1/3})},$

$d_1 = d_3 = \frac{1}{2-2^{1/3}},\ \ d_2 = -\frac{2^{1/3}}{2-2^{1/3}},\ \ d_4 = 0.$

To determine these coefficients, the Baker–Campbell–Hausdorff formula can be used. Yoshida, in particular, gives an elegant derivation of coefficients for higher-order integrators.

See also

• Energy drift
• Variational integrator
• Verlet integration

References

• Leimkuhler, Ben; Sebastian Reich (2005). Simulating Hamiltonian Dynamics. Cambridge University Press. ISBN 0-521-77290-7. 
• Hairer, Ernst; Christian Lubich, Gerhard Wanner (2006). Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations (2 ed.). Springer. ISBN 978-3-540-30663-4.