Symplectic Integrators
Summer research 2024: the three-body problem · Part 2 of 11
A weekly research log from my summer project at the University of Edinburgh, supervised by Prof. Jenni Smillie.
← Plotting and Visualising the Three-Body Problem · Ellipses and the Two-Body Problem →
This post explains what makes an integrator symplectic, and compares second-, third- and fourth-order symplectic schemes on the figure-eight orbit.
What is a symplectic integrator?
Start from Liouville’s theorem in Hamiltonian mechanics: the flow of a Hamiltonian system preserves phase-space volume. For one degree of freedom, a small region of the \((x, p)\) plane keeps the same area as it evolves, \(dA = dA'\):
So evolving Newton’s equations forward by a finite time defines an area-preserving (more generally, symplectic) map of phase space. A symplectic integrator is a numerical method whose single step is itself a symplectic map.
This matters for long simulations. A symplectic integrator does not conserve the energy exactly, but it exactly conserves a slightly modified “shadow” Hamiltonian, so the energy error stays bounded and oscillates instead of drifting. Non-symplectic methods such as the explicit Euler method and RK4 show a steady energy drift.
Checking whether a scheme is symplectic
For one degree of freedom, write a step of a linear scheme as
\[\begin{pmatrix} x_{n+1} \\ p_{n+1} \end{pmatrix}=A\begin{pmatrix} x_n \\ p_n \end{pmatrix}.\]The scheme is symplectic if and only if \(\det A = 1\). For a nonlinear scheme the same test is applied to the Jacobian of the map, and in higher dimensions the condition becomes \(M^{T}JM = J\) for the Jacobian \(M\) and the standard symplectic matrix \(J\).
Higher-order schemes
A family of explicit symplectic schemes applies, for \(i = 1, \ldots, n\) within each time step \(\Delta t\), the updates
\[\begin{aligned} \mathbf{v}_{i+1}&=\mathbf{v}_i+d_i\, \mathbf{a}\left(\mathbf{x}_i\right) \Delta t, \\ \mathbf{x}_{i+1}&=\mathbf{x}_i+c_i\, \mathbf{v}_{i+1}\, \Delta t, \end{aligned}\]with suitably chosen coefficients. A third-order scheme (Ruth) uses
\[c = \left(1, -\tfrac{2}{3}, \tfrac{2}{3}\right), \qquad d = \left(-\tfrac{1}{24}, \tfrac{3}{4}, \tfrac{7}{24}\right),\]and a fourth-order scheme (Forest–Ruth) uses
\[\begin{aligned} c &= \biggl(\frac{1}{2(2 - 2^{1/3})},\ \frac{1 - 2^{1/3}}{2(2 - 2^{1/3})}, \\ &\qquad \frac{1 - 2^{1/3}}{2(2 - 2^{1/3})},\ \frac{1}{2(2 - 2^{1/3})}\biggr), \\ d &= \left(\frac{1}{2 - 2^{1/3}},\ -\frac{2^{1/3}}{2 - 2^{1/3}},\ \frac{1}{2 - 2^{1/3}},\ 0\right). \end{aligned}\]Each sub-step is a “kick” followed by a “drift”, both symplectic maps, so their composition is symplectic too.
Accuracy of the higher-order schemes
I integrated the figure-eight orbit for \(50\,000\) steps with \(\Delta t = 0.01\). The energy deviations were:
| Scheme | Order | Energy deviation |
|---|---|---|
| Velocity Verlet | 2 | \(5.9 \times 10^{-5}\) |
| Ruth | 3 | \(4.1 \times 10^{-8}\) |
| Forest–Ruth | 4 | \(1.9 \times 10^{-8}\) |
The higher-order schemes improve the energy conservation by more than three orders of magnitude.
A remaining problem: momentum
For pairwise forces like gravity, the total linear momentum should be conserved to round-off precision by all of these schemes. In my simulation it drifted: the error was small, but it grew steadily with time.
Gallery
The animations are a highlight of this project:
And an animation of a perturbed figure-eight:
Credits
The section on the definition of symplectic integrators, including the figure, is based on these notes on the leapfrog method by Peter Young.
Changes on 25 October 2026
What was wrong. Nothing in the results.
What changed. The definition of a symplectic map is given in general, the coefficients are laid out more clearly and the energy errors are rounded to two figures.