4.2.5 Symplectic Algorithms

Hamiltonian dynamics is an important field of application for ODE algorithms. One simple performance test is the conservation of total mechanical energy.

Another quantity that should be conserved in Hamiltonian dynamics is the symplectic form. A number of algorithms have been devised which by construction conserve this quantity.

In a classical system with $M$ degrees of freedom, let the complete set of (generalized) coordinates be denoted by $\mbox{$\bf q$}$, and the conjugate momenta by $\mbox{$\bf p$}$. Then

$$\frac{d \mbox{$\bf q$}}{dt} = \nabla_{p} \, H(\mbox{$\bf q$}... ...\bf p$}}{dt} = -\nabla_{q} \, H(\mbox{$\bf q$},\mbox{$\bf p$})$$

where $H(\mbox{$\bf q$},\mbox{$\bf p$})$ is the Hamiltonian.

Rewrite this by defining the phase space vector $\mbox{$\bf z$}$ $\equiv \{ \mbox{$\bf q$},\mbox{$\bf p$}\}$:

$$\frac{d \mbox{$\bf z$}}{dt} = \mbox{${\bf J}$} \cdot \nabla_{z} \, H(\mbox{$\bf z$})$$

where

$$\mbox{${\bf J}$} \equiv \mbox{$\left( \begin{array}{cc}\mbox... ...-9pt}\\ -\mbox{${\bf I}$}&\mbox{${\bf0}$}\end{array} \right)$}$$

is the ``symplectic matrix''. (``symplectic'' means ``intertwined'').

The canonical transformation

$$\mbox{$\bf z$}(t_{0}) \; \Longrightarrow \; \mbox{$\bf z$}(t)$$

conserves not only the energy (= numerical value of the Hamiltonian), but also the symplectic form

$$s(\mbox{$\bf z$}_{1}, \mbox{$\bf z$}_{2}) \equiv \mbox{$\bf z$}_{1}^{T} \cdot \mbox{${\bf J}$} \cdot \mbox{$\bf z$}_{2}$$

EXAMPLE: Harmonic oscillator:

$$H(\mbox{$\bf z$}) \equiv H(q,p) = \frac{k}{2} q^{2} + \frac{p^{2}}{2m}$$

Canonical transformation:

$$\mbox{$\bf z$}(t) = \mbox{$\left( \begin{array}{r} q \\ \vsp... ...y} \right)$} \equiv \mbox{${\bf A}$} \cdot \mbox{$\bf z$}(0)$$

(with $\omega^{2} = k/m$.)

Energy:

$$\frac{k}{2} q^{2} + \frac{p^{2}}{2m} = \frac{k}{2} q^{2}(0) + \frac{p^{2}(0)}{2m}$$

Symplectic structure:
Let $\{q_{1}(0), p_{1}(0) \}$ and $\{q_{2}(0), p_{2}(0) \}$ be two different initial conditions; then

$$s(q_{1}(0) \dots p_{2}(0)) \equiv (q_{1}(0), p_{1}(0)) \cdot \mbox{$\left( \begin{array}{cc}0&1\\ ... ... \begin{array}{r} q_{2}(0) \\ \vspace{-9 pt}\\ p_{2}(0) \end{array} \right)$}$$

$$= q_{1}(0) p_{2}(0) - p_{1}(0) q_{2}(0)$$

Geometric interpretation: $s$ is the area of a parallelogram defined by the two initial state vectors $\mbox{$\bf z$}_{1,2}$.
$s$ is indeed constant:

$$s(\mbox{$\bf z$}_{1}(t),\mbox{$\bf z$}_{2}(t)) = \mbox{$\bf z$}_{1}^{T}(t) \cdot \mbox{${\bf J}$} \cdot \mbox{$\bf z$}_{2}(t)$$

$$= \mbox{$\bf z$}_{1}^{T}(0) \cdot \mbox{${\bf A}$}^{T} \cdot \mbox{${\bf J}$} \cdot \mbox{${\bf A}$} \cdot \mbox{$\bf z$}_{2}(0)$$

$$= \mbox{$\bf z$}_{1}^{T}(0) \cdot \mbox{${\bf J}$} \cdot \mbox{$\bf z$}_{2}(0)$$

This means that $\mbox{${\bf A}$}^{T} \cdot \mbox{${\bf J}$} \cdot \mbox{${\bf A}$} = \mbox{${\bf J}$}$.

Let us check the performance of the Euler-Cauchy scheme:

$$\mbox{$\bf z$}(t) = \mbox{$\left( \begin{array}{r} q \\ \vsp... ...y} \right)$} \equiv \mbox{${\bf E}$} \cdot \mbox{$\bf z$}(0)$$

The EC algorithm enhances the symplectic form by a factor $1+ (\omega \Delta t)^{2}$ at each time step. It is therefore not a symplectic algorithm.

Symplectic integrator of fourth order (Neri and Candy-Rozmus):

Third-order scheme (Ruth): same structure as Candy's algorithm, but with coefficients

$$(a_{1},a_{2},a_{3}) = (2/3,-2/3,1)$$

$$(b_{1},b_{2},b_{3}) = (7/24, 3/4, -1/24)$$

Of the algorithms discussed previously only the Størmer-Verlet formula is symplectic.

Some theory:
For non-integrable Hamiltonians there can be no algorithm that conserves both energy and symplectic structure. But symplectic integrators conserve a Hamiltonian function that is different from, but close to, the given Hamiltonian (YOSHIDA). Therefore they display no long-time energy trend. (Contrast: e.g. Runge-Kutta integrators have good short-time accuracy but show a regularly increasing deviation in energy.

Example: the first-order symplectic algorithm (Euler-Cromer)

$$\mbox{$\bf p$}_{n+1}=\mbox{$\bf p$}_{n}+\mbox{$\bf F$}(\mbox... ...mbox{$\bf q$}_{n}+\mbox{$\bf P$}(\mbox{$\bf p$}_{n+1})\Delta t$$

exactly conserves

$$\tilde{H} \equiv H + H_{1}\Delta t + H_{2}(\Delta t)^{2}+ H_{3}(\Delta t)^{3}+ \dots$$

where

$$H_{1}=\frac{1}{2}H_{p}H_{q}\, , \;\;\; H_{2}=\frac{1}{12}(H_... ...)\, , \;\;\; H_{3}=\frac{1}{12}H_{pp}H_{qq}H_{p}H_{q}\,\;\dots$$

( $H_{q} =\nabla_{q}H$ etc.) To take the harmonic oscillator, the perturbed Hamiltonian

$$\tilde{H}=H_{ho}+\frac{\textstyle \omega^{2} \Delta t}{\textstyle 2}pq$$

is conserved exactly.

EXERCISE: Apply the (non-symplectic) RK method and the (symplectic) Størmer-Verlet algorithm (or the Candy procedure) to the one-body Kepler problem with elliptic orbit. Perform long runs to assess the long-time performance of the integrators. (For RK the orbit should eventually spiral down towards the central mass, while the symplectic procedures should only give rise to a gradual precession of the perihelion.)