advanced32 min studyLesson 922 in Physics

Symplectic integrators for orbital motion

Preserve Hamiltonian geometry while assessing energy, phase and close-encounter errors separately.

On this page

What you’ll understand

  • Construct a kick–drift–kick Verlet step.
  • Distinguish symplecticity from exact conservation of the original Hamiltonian.
  • Identify why fixed steps and resolved orbital time scales matter.

An orbit can look plausible for a few revolutions while its numerical phase or energy drifts steadily. Symplectic integration addresses the geometry of Hamiltonian flow, which often improves long-time behavior in conservative mechanics. It does not force the original energy to be exactly constant at every step. A useful evaluation therefore checks geometric structure, energy error, orbital phase and spatial resolution together. A tiny harmonic benchmark makes these distinctions visible before a singular gravitational potential introduces close encounters and additional modeling choices.

Build the map from exact subflows

For a separable Hamiltonian, kinetic motion changes position at fixed momentum and potential motion changes momentum at fixed position. Each subflow is canonical. Composing a half potential kick, a full kinetic drift and another half kick gives velocity Verlet in momentum form. The symmetric composition has second-order global accuracy over a fixed time interval for sufficiently smooth forces. The corresponding one-step state error is third order; those two error statements refer to different comparisons.

For a gravitational test particle, the potential is minus GMm divided by radius. The kicks then use the force at the old and new positions. Rotationally invariant pair kicks can preserve angular momentum in exact arithmetic, but energy generally oscillates around a nearby value. Arbitrarily replacing the force or applying inconsistent updates to different particles can destroy properties of the intended Hamiltonian map. One must specify the actual force calculation as carefully as the time integrator.

pn+1/2=pn−h2∇V(qn),qn+1=qn+hmpn+1/2,pn+1=pn+1/2−h2∇V(qn+1)p_{n+1/2}=p_n-\frac h2\nabla V(q_n),\quad q_{n+1}=q_n+\frac h m p_{n+1/2},\quad p_{n+1}=p_{n+1/2}-\frac h2\nabla V(q_{n+1})
For H = p²/(2m) + V(q), q is position, p momentum and h time. Every kick has momentum units. The formula assumes constant positive mass and a differentiable potential along the step.↔ Scroll wide equations horizontally.

Preserve the two-form rather than impose an energy identity

A symplectic map preserves the canonical two-form. Its derivative obeys a matrix identity involving the canonical antisymmetric matrix J. In a single position–momentum plane this implies area preservation; in higher dimensions volume preservation alone is weaker than symplecticity. Explicit Euler applied simultaneously to position and momentum fails even the planar area test for a harmonic oscillator. Small time steps reduce its error but do not make that particular map exactly symplectic.

Backward error analysis explains why a suitable fixed-step symplectic scheme approximately follows a modified Hamiltonian. Under smoothness, bounded-motion and sufficiently small-step assumptions this can yield small energy error over long intervals. The conclusion is conditional, especially near singularities, resonances or poorly resolved fast motion. It also does not imply uniformly small phase error: a slightly incorrect orbital frequency can accumulate a substantial angular displacement while the energy remains well behaved.

(DΦh)TJ(DΦh)=J,H(Φh(q,p))≠H(q,p) in general,Hh=H+O(h2) for Verlet(D\Phi_h)^{\mathsf T}J(D\Phi_h)=J,\qquad H(\Phi_h(q,p))\ne H(q,p)\ \text{in general},\qquad H_h=H+O(h^2)\ \text{for Verlet}
Φh is the fixed-step map and J the canonical symplectic matrix in the chosen coordinate ordering. The modified Hamiltonian statement is an asymptotic description under the stated regularity and step restrictions.↔ Scroll wide equations horizontally.

Resolve the orbit before trusting long-time diagnostics

An eccentric orbit has a short dynamical time near periapsis, so a step adequate near apoapsis may miss its strongest acceleration. Ordinary state-dependent step-size substitution generally sacrifices the symplectic property; special extended-phase-space constructions require their own derivation. Collisions and close encounters may need regularization or a deliberately different integrator. Softening a force avoids a singular acceleration by changing the potential, and therefore changes the physical model being integrated.

Useful diagnostics include convergence of periapsis location, phase, conserved quantities and trajectories over a stated interval. In chaotic many-body motion, pointwise trajectories can separate even when statistical observables converge, so the target observable matters. The oscillator lab supplies the exact analytic motion of a single damped or undamped oscillator. It can benchmark expected periods, but it contains no finite time step, kick–drift update or computed orbit and cannot demonstrate a numerical energy-conservation theorem.

V(r)=−GMmr,tdyn(r)∼r3GM,h≪tdyn is a resolution aim, not a universal error boundV(r)=-\frac{GMm}{r},\qquad t_{\rm dyn}(r)\sim\sqrt{\frac{r^3}{GM}},\qquad h\ll t_{\rm dyn}\ \text{is a resolution aim, not a universal error bound}
r is separation, GM has length³/time² units and tdyn is a local time scale. A quantitative step requirement depends on the orbit, method, tolerance and desired observable.↔ Scroll wide equations horizontally.

One Verlet step changes the true energy slightly

Use dimensionless harmonic units m = ω = 1, so H = (p² + q²)/2. Starting at q = 1, p = 0, take one kick–drift–kick step h = 0.20. Compute the new state and its original Hamiltonian.

  1. The first kick gives p1/2=0−0.10(1)=−0.10p_{1/2}=0-0.10(1)=-0.10. The initial energy was one half.

  2. The drift gives q1=1+0.20(−0.10)=0.98q_1=1+0.20(-0.10)=0.98. The second kick therefore uses the new position and gives p1=−0.10−0.10(0.98)=−0.198p_1=-0.10-0.10(0.98)=-0.198.

  3. The original energy becomes (0.982+0.1982)/2=0.499802(0.98^2+0.198^2)/2=0.499802. Its change is −0.000198, despite the symplectic construction. This single step is evidence of the distinction, not a long-time error bound for every orbit.

The result

The new state is (q, p) = (0.98, −0.198) and H = 0.499802, slightly below the initial 0.5.

Explore the model

Interactive lab

Interactive lab · oscillator

An analytic period benchmark

Reset and set BOTH Damping ratio ζ = 0 and Natural frequency ω₀ = 1; predict period 2π. Reset and set ζ = 0 and ω₀ = 2; predict period π. Compare with Play or Step over several cycles. The lab evaluates analytic cosine motion from x0 = 1, v0 = 0; it does not run Verlet, vary a numerical step or calculate gravitational orbits.

An analytic period benchmarkUnit mass, x(0) = 1 and v(0) = 0. The analytic solution includes underdamped, critically damped and overdamped cases; time is measured in seconds. Current readings: Displacement x(t): 1; Velocity v(t): 0; Regime: underdamped. All parameters are adjustable with the labeled controls below.02.557.510−1.15−0.5700.571.15Time · sDisplacementunderdamped
Displacement x(t)Current time
Displacement x(t)
1
Velocity v(t)
0
Regime
underdamped
t = 0 s

Unit mass, x(0) = 1 and v(0) = 0. The analytic solution includes underdamped, critically damped and overdamped cases; time is measured in seconds. Parameter changes pause playback and recompute the model at the current elapsed time. Reset restores both parameters and returns time to zero. Playback starts paused. Adjust either parameter, then choose Play or Step.

Look forAn exact reference solution supplies a benchmark, while numerical geometric properties require examining the integrator.

Check your understanding

Which statement correctly describes a sufficiently resolved fixed-step symplectic integrator?

Try it yourself

For the same dimensionless harmonic oscillator and initial state, apply simultaneous explicit Euler with h = 0.20. Find the new energy and the determinant of its position–momentum update matrix. Does it preserve planar area?

Take a moment to make a prediction or work it out on paper. Open a hint when you need one.

Give me a hint

Euler uses qnew = q + hp and pnew = p − hq, both evaluated at the old state.

Compare with the solution

Euler gives qnew = 1 and pnew = −0.20, hence Hnew = 0.52. The update matrix is (1h−h1)\left(\begin{smallmatrix}1&h\\-h&1\end{smallmatrix}\right), with determinant 1+h2=1.041+h^2=1.04. A small phase-space area therefore grows by four percent in one step in this harmonic test. It is not area preserving or symplectic. The larger energy here is a concrete defect for this method and data; a general comparison of integrators would also examine accumulated error, force evaluations and the requested orbital observable.

Take it with you

  • Verlet composes canonical kicks and drifts.
  • Symplectic geometry differs from exact conservation of H.
  • Phase accuracy, close encounters and step selection need separate checks.

Read further

Original sources and deeper explanations to follow at your own pace.

Optional progress, saved in this browser when available. Revisit any lesson at any time.

Lesson not marked complete.