Conventional numerical ODE methods (Runge–Kutta, multistep, etc.) approximate solutions to differential equations with high accuracy over short time intervals. But for Hamiltonian systems, long-time behaviour—not short-time accuracy—is critical.
Symplectic integrators are designed explicitly to preserve:
This section explains the fundamental reasons why symplectic integrators dramatically outperform conventional schemes for Hamiltonian dynamics.
Consider a Hamiltonian system
\[ \dot{z} = X_H(z), \qquad H(z) = \text{constant along solutions}. \]
Standard integrators rarely preserve this exactly. Even if a method has high order, the computed energy \(H(z_n)\) typically satisfies:
\[ H(z_n) = H(z_0) + C h t + \mathcal{O}(h), \]
showing a linear drift in time.
\[ z_{n+1} = z_n + h X_H(z_n) \]
Explicit Euler is not symplectic and produces spiralling trajectories even for simple systems like the harmonic oscillator.
A numerical map \(\Phi_h: M\to M\) is symplectic if:
\[ \Phi_h^* \omega = \omega. \]
In canonical coordinates, this reads:
\[ (D\Phi_h)^T J (D\Phi_h) = J, \qquad J=\begin{pmatrix}0 & I \\ -I & 0\end{pmatrix}. \]
Symplectic methods conserve:
But the most striking property involves backward error analysis.
For symplectic integrators, one can show that:
The numerical solution exactly follows the flow of a modified Hamiltonian \(\tilde{H}\), where \[ \tilde{H} = H + h^p H_{p+1} + h^{p+2} H_{p+2} + \cdots. \]
This implies:
\[ \tilde{H}(z_n) = \tilde{H}(z_0) \quad \text{exactly}. \]
Therefore, the true energy satisfies:
\[ |H(z_n) - H(z_0)| \le C h^p \]
for times that grow like:
\[ t = \mathcal{O}(e^{c/h}). \]
This is an exponentially long interval of near-conservation — a result that no non-symplectic method can match.
For nearly integrable systems, the motion lies on invariant tori (Kolmogorov–Arnold–Moser theory).
A symplectic integrator maps invariant tori to slightly deformed tori, preserving quasiperiodic structure.
A non-symplectic method destroys tori, causing:
Below is a JavaScript comparison for the harmonic oscillator:
This section motivates the use of symplectic integrators by explaining:
These ideas form the conceptual foundation for the concrete algorithms presented in the remainder of Chapter 3.
The simplest symplectic integrators arise from the splitting of a separable Hamiltonian
\[ H(q,p) = T(p) + V(q), \]
where \(T\) is kinetic energy and \(V\) potential energy. This structure allows the Hamiltonian vector field to be decomposed into two exactly integrable flows: the drift and the kick.
From these flows one constructs the Symplectic Euler (two variants) and the Störmer–Verlet (or Leapfrog) integrator.
These are the foundation for higher-order symplectic composition methods introduced later.
For \(H = T(p) + V(q)\) the Hamiltonian vector field decomposes:
\[ X_H = X_T + X_V. \]
Their flows are exactly solvable:
The exact time-\(h\) flows are therefore:
\[ \Phi_T^h(q,p) = (q + h\,\nabla_p T(p),\, p), \] \[ \Phi_V^h(q,p) = (q,\; p - h\,\nabla_q V(q)). \]
Symplectic integrators are obtained by composing these exact flows.
There are two first-order symplectic Euler methods:
\[ \begin{aligned} p_{n+1} &= p_n - h\,\nabla_q V(q_n), \\ q_{n+1} &= q_n + h\,\nabla_p T(p_{n+1}). \end{aligned} \]
This corresponds to:
\[ \Phi_h^{\mathrm{A}} = \Phi_T^h \circ \Phi_V^h. \]
\[ \begin{aligned} q_{n+1} &= q_n + h\,\nabla_p T(p_n), \\ p_{n+1} &= p_n - h\,\nabla_q V(q_{n+1}). \end{aligned} \]
Corresponding to:
\[ \Phi_h^{\mathrm{B}} = \Phi_V^h \circ \Phi_T^h. \]
Although only first order, these methods:
The second-order Störmer–Verlet method is obtained by symmetric composition:
\[ \Phi^{\mathrm{SV}}_h = \Phi_T^{h/2} \circ \Phi_V^h \circ \Phi_T^{h/2}. \]
The update rule is:
\[ \begin{aligned} p_{n + 1/2} &= p_n - \frac{h}{2}\nabla_q V(q_n), \\ q_{n+1} &= q_n + h\,\nabla_p T(p_{n+1/2}), \\ p_{n+1} &= p_{n+1/2} - \frac{h}{2}\nabla_q V(q_{n+1}). \end{aligned} \]
Properties:
Since each subflow is symplectic, the composed method is symplectic.
Below, we integrate the harmonic oscillator with:
This section developed the simplest—and most fundamental—symplectic integrators:
These schemes constitute the backbone of geometric integration, molecular dynamics, celestial mechanics, HMC (Hamiltonian Monte Carlo), and quantum simulation splitting methods.
Symplectic Euler and Störmer–Verlet can be obtained in many ways:
The variational viewpoint is particularly elegant: one discretises the action integral rather than the differential equations. The resulting numerical method is automatically:
Let \(Q\) be the configuration manifold and \(L: TQ \to \mathbb{R}\) the Lagrangian, typically
\[ L(q,\dot{q}) = T(q,\dot{q}) - V(q). \]
The action of a curve \(q:[t_0,t_1]\to Q\) is
\[ \mathcal{S}[q] = \int_{t_0}^{t_1} L\big(q(t),\dot{q}(t)\big)\,dt. \]
Hamilton’s principle states that the actual trajectory satisfies
\[ \delta \mathcal{S}[q] = 0 \]
for variations \(\delta q\) vanishing at the endpoints. This yields the Euler–Lagrange equations:
\[ \frac{d}{dt}\frac{\partial L}{\partial \dot{q}} - \frac{\partial L}{\partial q} = 0. \]
For a mechanical Lagrangian \(L = \tfrac{1}{2}\dot{q}^{\mathsf T}M\dot{q}-V(q)\), this is equivalent to Hamilton’s equations after a Legendre transform.
We now discretise time:
\[ t_k = t_0 + k h, \qquad k = 0,\dots,N. \]
A discrete Lagrangian is a function
\[ L_d : Q\times Q \to \mathbb{R}, \]
intended as an approximation to the exact action along the true trajectory joining \(q_k\) and \(q_{k+1}\):
\[ L_d(q_k,q_{k+1}) \approx \int_{t_k}^{t_{k+1}} L\big(q(t),\dot{q}(t)\big)\,dt. \]
The discrete action is defined as the sum over all segments:
\[ \mathcal{S}_d[q_0,\dots,q_N] = \sum_{k=0}^{N-1} L_d(q_k,q_{k+1}). \]
Discrete Hamilton’s principle:
A discrete trajectory \((q_k)_{k=0}^N\) is a solution if \(\mathcal{S}_d\) is stationary with respect to variations of \(q_1,\dots,q_{N-1}\) with fixed endpoints \(q_0,q_N\).
Vary \(\mathcal{S}_d\) with respect to \(q_k\) for \(1 \le k \le N-1\). Each \(q_k\) appears only in \(L_d(q_{k-1},q_k)\) and \(L_d(q_k,q_{k+1})\):
\[ \delta\mathcal{S}_d = \sum_{k=1}^{N-1} \Big( D_2 L_d(q_{k-1},q_k) + D_1 L_d(q_k,q_{k+1}) \Big)\cdot \delta q_k, \]
where \(D_1\) and \(D_2\) denote partial derivatives with respect to the first and second argument.
Stationarity for all \(\delta q_k\) yields the discrete Euler–Lagrange equations:
\[ D_2 L_d(q_{k-1},q_k) + D_1 L_d(q_k,q_{k+1}) = 0, \qquad k=1,\dots,N-1. \]
If the mixed derivative \(D_{12}^2 L_d\) is nondegenerate, these equations implicitly define an update map \((q_{k-1},q_k)\mapsto(q_k,q_{k+1})\).
Define discrete momenta via the discrete Legendre transforms:
\[ p_k^- := - D_1 L_d(q_k,q_{k+1}), \qquad p_k^+ := D_2 L_d(q_{k-1},q_k). \]
The discrete Euler–Lagrange equation is exactly
\[ p_k^- = p_k^+. \]
Hence we simply write \(p_k := p_k^- = p_k^+\) and obtain an update rule \((q_k,p_k) \mapsto (q_{k+1},p_{k+1})\).
A key theorem:
Theorem 3.1 (Symplecticity of Variational Integrators). The map \((q_k,p_k)\mapsto(q_{k+1},p_{k+1})\) defined by the discrete Euler–Lagrange equations is symplectic with respect to the canonical symplectic form on \(T^*Q\).
Proof (outline): view the discrete action as a generating function of type I, and show that the pullback of the canonical 1-form by the update map differs by an exact form, implying preservation of the symplectic 2-form.
Consider a mechanical Lagrangian (mass matrix \(M\) constant, for simplicity):
\[ L(q,\dot{q}) = \frac{1}{2}\dot{q}^{\mathsf T}M\dot{q} - V(q). \]
A natural second-order accurate discrete Lagrangian:
\[ L_d(q_k,q_{k+1}) = h\,L\left( \frac{q_k+q_{k+1}}{2}, \frac{q_{k+1}-q_k}{h} \right). \]
Compute:
\[ D_1 L_d(q_k,q_{k+1}) = -M\,\frac{q_{k+1}-q_k}{h} -\frac{h}{2}\nabla V\!\left(\frac{q_k+q_{k+1}}{2}\right), \] \[ D_2 L_d(q_k,q_{k+1}) = M\,\frac{q_{k+1}-q_k}{h} -\frac{h}{2}\nabla V\!\left(\frac{q_k+q_{k+1}}{2}\right). \]
Discrete Euler–Lagrange:
\[ D_2 L_d(q_{k-1},q_k) + D_1 L_d(q_k,q_{k+1}) = 0 \]
yields (after rearrangement) the familiar Störmer–Verlet scheme. In the simplest case \(M=I\) and using equivalent algebra, one recovers:
\[ q_{k+1} - 2q_k + q_{k-1} = - h^2 \nabla V(q_k), \]
which is the position form of the Verlet method.
Suppose a Lie group \(G\) acts on \(Q\) and leaves the discrete Lagrangian invariant:
\[ L_d(g\cdot q_k, g\cdot q_{k+1}) = L_d(q_k,q_{k+1}), \qquad \forall g\in G. \]
Then there exists a discrete momentum map \(J_d\) such that the discrete flow preserves \(J_d\):
\[ J_d(q_{k},q_{k+1}) = J_d(q_{k-1},q_k). \]
This is the discrete analogue of Noether’s theorem. In particular:
For the harmonic oscillator with
\[ L(q,\dot{q}) = \tfrac{1}{2}\dot{q}^2 - \tfrac{1}{2}\omega^2 q^2, \]
the Verlet position update is
\[ q_{k+1} - 2q_k + q_{k-1} = -h^2\omega^2 q_k. \]
Below we iterate this scheme and plot phase space \((q_k,\dot{q}_k)\) approximated by \(\dot{q}_k \approx (q_{k+1}-q_{k-1})/(2h)\).
Variational integrators arise from a discrete action principle and inherit:
Störmer–Verlet is a prime example, providing a unifying viewpoint that connects splitting, variational, and Hamiltonian perspectives.
Splitting methods are among the most fundamental geometric integrators, applicable to Hamiltonian, Poisson, Lie–Poisson, and PDE systems. They are constructed by decomposing the vector field or operator into subproblems with exactly computable flows.
This section explains:
These techniques generalise Störmer–Verlet to arbitrary Hamiltonian or non-Hamiltonian systems.
Consider an ODE:
\[ \dot{z} = (A + B)z, \]
where \(A\) and \(B\) are vector fields (or linear operators). If the flows
\[ \Phi_A^t = e^{tA},\qquad \Phi_B^t = e^{tB} \]
are easily computed, then the flow of the sum is approximated by composing \(\Phi_A^t\) and \(\Phi_B^t\).
The Lie–Trotter splitting is:
\[ \Phi_{A+B}^h \approx \Phi_A^h \circ \Phi_B^h, \]
which is first-order accurate.
The error of splitting is governed by the BCH formula. For noncommuting operators:
\[ e^{hA}e^{hB} = \exp\left( h(A+B) + \frac{h^2}{2}[A,B] + \frac{h^3}{12}[A,[A,B]] - \frac{h^3}{12}[B,[A,B]] + \cdots \right). \]
Key consequences:
A symmetric composition gives second order:
\[ \Phi_{A+B}^h \approx \Phi_A^{h/2}\circ\Phi_B^h\circ\Phi_A^{h/2}. \]
Because the composition is symmetric:
\[ S_h^{-1} = S_{-h}, \]
all odd BCH terms cancel, yielding second-order accuracy.
The Störmer–Verlet method is exactly the Strang splitting for a separable Hamiltonian \(H=T(p)+V(q)\).
To obtain order \(p>2\), one composes several Strang splittings:
\[ \Psi_h = S_{\gamma_1 h}\circ S_{\gamma_2 h}\circ\cdots\circ S_{\gamma_s h}, \]
where the coefficients \(\gamma_i\) must satisfy nonlinear algebraic order conditions.
Well-known families:
\[ \Psi_h^{(4)} = S_{\alpha h}\circ S_{\beta h}\circ S_{\alpha h}, \]
where:
\[ \alpha = \frac{1}{2 - 2^{1/3}},\qquad \beta = 1 - 2\alpha. \]
This cancels BCH error terms up to order \(h^4\).
Higher-order splitting methods require some coefficients \(\gamma_i\) to be:
Real composition methods of order > 2 must contain negative coefficients (Sheng–Suzuki no-go theorem). Complex-coefficient methods can achieve:
If the basic method \(S_h\) is:
This follows because the properties are preserved under composition and inversion.
We compare Strang splitting (2nd order) versus Yoshida 4th order on the harmonic oscillator. Strang produces a deformed ellipse; Yoshida produces a nearly perfect ellipse.
Splitting and composition techniques provide a powerful framework for constructing symplectic, volume-preserving, and structure-preserving integrators of arbitrary order.
These methods underpin modern geometric integration for:
Traditional error analysis estimates the difference between a numerical solution and the exact flow of the original differential equation. This yields useful short-time bounds but fails to explain the remarkable long-time stability of structure-preserving (especially symplectic) integrators.
Backward error analysis (BEA) takes the opposite perspective:
Instead of asking: “How far is the numerical solution from the true solution of the original ODE?”, BEA asks: “Is the numerical solution the exact solution of a nearby differential equation?”
This leads to the concept of a modified (or shadow) differential equation. For symplectic methods applied to Hamiltonian systems, the modified equation is itself Hamiltonian with a modified Hamiltonian.
This yields exponentially accurate energy conservation for exponentially long times, and explains the preservation of quasiperiodic motion, tori, and other qualitative features.
Consider an ODE
\[ \dot{z} = f(z), \]
and a one-step numerical integrator \(\Phi_h\). If \(\Phi_h\) is smooth in \(h\), one may formally write
\[ \Phi_h = \exp(h F_h), \]
where \(F_h\) is a formal power series of differential operators:
\[ F_h = f + h f_1 + h^2 f_2 + h^3 f_3 + \cdots. \]
This defines the modified vector field.
Truncating at order \(h^N\) gives a modified ODE that the method integrates exactly up to errors of order \(h^{N+1}\).
For splitting methods, the BCH formula gives explicit modified equations. Example for Strang splitting of \(\dot{z}=(A+B)z\):
\[ \Phi_h = e^{\frac{h}{2}A} e^{hB} e^{\frac{h}{2}A} = \exp\left( h(A+B) + \frac{h^3}{24}([A,[A,B]] - [B,[A,B]]) + \mathcal{O}(h^5) \right). \]
Thus the modified vector field is:
\[ F_h = A + B + \frac{h^2}{24}\big([A,[A,B]] - [B,[A,B]]\big) + \mathcal{O}(h^4). \]
Higher-order terms involve nested commutators; symbolic or automatic computation is standard in modern implementations.
For a symplectic integrator applied to a Hamiltonian system
\[ \dot{z} = X_H(z), \]
the modified vector field is also Hamiltonian:
\[ F_h = X_{\tilde{H}}, \]
with the modified Hamiltonian
\[ \tilde{H} = H + h^p H_{p+1} + h^{p+2} H_{p+2} + \cdots. \]
This follows from:
Since symplectic maps satisfy
\[ (\Phi_h)^*\omega = \omega, \]
\(\Phi_h\) lies in the Lie group of canonical transformations whose Lie algebra consists of Hamiltonian vector fields.
Since \(\tilde{H}\) is exactly preserved:
\[ \tilde{H}(z_n) = \tilde{H}(z_0), \]
and the difference between \(H\) and \(\tilde{H}\) is small:
\[ |\tilde{H}(z) - H(z)| \le C h^p, \]
the true energy exhibits only small, bounded oscillations:
\[ |H(z_n) - H(z_0)| \le C h^p \]
for times up to:
\[ t = \mathcal{O}(e^{c / h}), \]
an exponentially long interval.
No non-symplectic method exhibits this property.
For the harmonic oscillator:
\[ H = \frac{1}{2}(p^2 + \omega^2 q^2), \]
symplectic Euler and Verlet modify \(\omega\) to a nearby \(\tilde{\omega}\), explaining the slight deformation of the ellipse in phase space.
Example for Verlet:
\[ \tilde{\omega} = \frac{1}{h}\arccos\left(1 - \frac{h^2\omega^2}{2}\right) = \omega + \frac{\omega^3 h^2}{24} + \mathcal{O}(h^4). \]
So Verlet is exactly periodic with period \(2\pi/\tilde{\omega}\).
For PDEs (e.g. Schrödinger, diffusion–advection, Vlasov–Poisson), splitting methods correspond to Trotter–Suzuki decompositions of linear or nonlinear operators.
The modified PDE typically contains commutators of differential operators:
\[ [\mathcal{A},\mathcal{B}] = \mathcal{A}\mathcal{B} - \mathcal{B}\mathcal{A}. \]
The same cancellation logic applies:
Below we simulate energy evolution of:
Backward error analysis explains the extraordinary long-time fidelity of structure-preserving methods:
This provides the theoretical core for Hamiltonian Monte Carlo, molecular dynamics, symplectic PDE solvers, and high-order geometric integrators.
Symmetry plays an essential role in geometric integration. A numerical method is symmetric (or time-reversible) if its inverse is obtained by flipping the sign of the time step:
\[ \Phi_h^{-1} = \Phi_{-h}. \]
Symmetric integrators have error expansions containing only odd powers of \(h\), giving them “free” advantages:
Moreover, processing (conjugation by a near-identity map) allows a low-order method to behave like a high-order one, without additional cost per step. This is widely used in molecular dynamics and celestial mechanics.
A numerical integrator \(\Phi_h\) is symmetric if:
\[ \Phi_h = \Phi_{-h}^{-1}. \]
Equivalently:
\[ \Phi_{-h} = \Phi_h^{-1}. \]
Example: the Strang splitting for \(A+B\)
\[ S_h = e^{\frac{h}{2}A} e^{hB} e^{\frac{h}{2}A} \]
satisfies:
\[ S_h^{-1} = S_{-h}. \]
This symmetry forces the modified vector field to contain only odd powers:
\[ F_h = A+B + h^2 E_3 + h^4 E_5 + \cdots. \]
Thus a symmetric method of order \(p\) actually behaves like an order \(p+1\) method over long times (when \(p\) is even), thanks to the absence of even-order error terms.
Compose a method with its adjoint:
\[ \Phi_h^* := (\Phi_h)^{-1}|_{h\mapsto -h}. \]
A simple composition:
\[ \Psi_h = \Phi_{a h} \circ \Phi_{b h}^*, \]
is symmetric if \(a=b\).
In general, if the sequence of coefficients \(\{a_i\}\) is palindromic:
\[ (a_1,a_2,\dots,a_s,a_s,\dots,a_2,a_1), \]
then the resulting method is symmetric.
This is used in:
The central idea:
A low-order integrator can often be conjugated by a near-identity map to produce an integrator of much higher effective order.
Let \(\Phi_h\) be any integrator and let \(P_h\) be a near-identity map (the processor).
\[ \Psi_h := P_h^{-1} \circ \Phi_h \circ P_h. \]
Properties:
The modified vector fields transform as:
\[ e^{hF_h^\Psi} = e^{-G_h} e^{hF_h} e^{G_h}, \]
with \(G_h\) corresponding to the generator of \(P_h\).
Suitable choice of \(G_h\) cancels many low-order error terms.
Let \(\Phi_h^{\mathrm{SV}}\) be the Verlet method. One can design a processor:
\[ P_h = \exp(h^2 G), \]
such that the processed method:
\[ \Psi_h = P_h^{-1}\circ \Phi_h^{\mathrm{SV}} \circ P_h \]
has a modified Hamiltonian:
\[ \tilde{H}_\Psi = H + h^4 K + \mathcal{O}(h^6), \]
even though the unprocessed Verlet method has error terms at order \(h^2\).
That is, processed Verlet becomes a fourth-order method in terms of phase accuracy for very little cost.
The modified flow is:
\[ \dot{z} = F_h(z). \]
One may integrate a truncated version:
\[ \dot{z} = f(z) + h^p f_{p+1}(z), \]
where \(f_{p+1}\) is computed symbolically or automatically.
This can dramatically reduce global error (similar to Richardson or extrapolation), but is more compatible with geometric integration because the added term preserves (in Hamiltonian case) the symplectic structure.
Below we illustrate phase-space trajectories of:
In this section we established:
Symmetry + processing are two of the most powerful engineering ideas in geometric integration.
Runge–Kutta (RK) methods form a very large and flexible family of numerical integrators for ODEs. For Hamiltonian and symplectic problems, certain RK methods possess the remarkable property of being symplectic, meaning they exactly preserve the canonical symplectic form.
The most important examples are:
These methods play a fundamental role in modern geometric integration, especially in achieving arbitrarily high order while preserving symplectic structure and long-time fidelity.
For the ODE \(\dot{z} = f(z)\), an \(s\)-stage RK method has the form:
\[ \begin{aligned} Z_i &= z_n + h\sum_{j=1}^{s}a_{ij} f(Z_j), \qquad i=1,\dots,s, \\ z_{n+1} &= z_n + h\sum_{i=1}^{s}b_i f(Z_i). \end{aligned} \]
The coefficients \((a_{ij}, b_i, c_i)\) define the method and are often expressed in a Butcher tableau:
c₁ | a₁₁ a₁₂ ... a₁ₛ
c₂ | a₂₁ a₂₂ ... a₂ₛ
⋮ | ⋮ ⋮ ⋱ ⋮
cₛ | aₛ₁ aₛ₂ ... aₛₛ
---|---------------------
| b₁ b₂ ... bₛ
A Runge–Kutta method is symplectic if and only if its coefficients satisfy:
\[ b_i a_{ij} + b_j a_{ji} = b_i b_j, \qquad \forall i,j. \]
This is a deep condition derived from preservation of the canonical 2-form:
\[ \omega = dq\wedge dp. \]
Consequences include:
The constraint essentially says: A symplectic RK method must be implicit.
Consider the collocation method where the stage values satisfy:
\[ \dot{z}(t_n + c_i h) = f(Z_i), \]
and where the collocation points \(c_i\) are the Gauss–Legendre points on \([0,1]\). These produce:
The 1-stage Gauss method is the implicit midpoint rule:
\[ z_{n+1} = z_n + h f\!\left(\frac{z_n + z_{n+1}}{2}\right). \]
It is symplectic, second order, symmetric, and fundamental.
Gauss–Legendre methods can be interpreted as:
Their modified Hamiltonian \(\tilde{H}\) has unusually small error constants, leading to:
For separable Hamiltonian systems:
\[ H(q,p) = T(p) + V(q), \]
PRK methods apply different Butcher tableaux to the \(q\) and \(p\) components:
\[ \begin{aligned} Q_i &= q_n + h\sum_j a_{ij}^{(q)} p_j,\\ P_i &= p_n - h\sum_j a_{ij}^{(p)} \nabla V(Q_j),\\ q_{n+1} &= q_n + h\sum_i b_i^{(q)} P_i,\\ p_{n+1} &= p_n - h\sum_i b_i^{(p)} \nabla V(Q_i). \end{aligned} \]
Symplecticity requires:
\[ b_i^{(q)} a_{ij}^{(p)} + b_j^{(p)} a_{ji}^{(q)} = b_i^{(q)} b_j^{(p)}. \]
Many well-known integrators (e.g. ‘velocity Verlet’) are PRK methods.
Symplectic RK and PRK methods are widely used in:
The Gauss–Legendre methods are especially effective for:
Below we compare:
The harmonic oscillator shows energy drift for RK4 and near-perfect bounded energy for midpoint.
This section established:
Gauss–Legendre schemes represent the “Rolls-Royce” of high-order symplectic ODE integrators.
While symplectic integrators preserve the symplectic 2-form and thus induce near-conservation of Hamiltonian energy, they do not generally conserve the Hamiltonian exactly.
Some systems, however, require exact energy conservation:
In these cases, one turns to energy-preserving integrators such as:
Another important structure is volume-preservation, which appears in:
We now survey these ideas and their geometric foundations.
Let \(I(z)\) be a conserved quantity of a system:
\[ \frac{d}{dt} I(z(t)) = \nabla I(z)^\top f(z) = 0. \]
A numerical method preserves this invariant if:
\[ I(z_{n+1}) = I(z_n) \quad \text{for all }n. \]
In general this is impossible for arbitrary high-order RK methods.
But it is possible when the method is built using:
Let \(H(z)\) be the Hamiltonian. A discrete gradient of \(H\) is a function
\[ \bar{\nabla}H(z_n,z_{n+1}) \]
satisfying:
\[ H(z_{n+1}) - H(z_n) = \bar{\nabla}H(z_n,z_{n+1})^\top (z_{n+1} - z_n). \]
This ensures that any update of the form:
\[ \frac{z_{n+1} - z_n}{h} = S(z_n,z_{n+1})\, \bar{\nabla} H(z_n,z_{n+1}) \]
preserves energy exactly whenever the matrix \(S\) is skew-symmetric.
Typical choices:
For \(\dot{z} = f(z)\), the AVF method is:
\[ z_{n+1} = z_n + h \int_0^1 f\!\left( (1-\xi) z_n + \xi z_{n+1} \right)\, d\xi. \]
For Hamiltonian systems with \(\dot{z} = S \nabla H(z)\), AVF preserves energy exactly:
\[ H(z_{n+1}) = H(z_n). \]
It is implicit but only requires one integral over the line segment between \(z_n\) and \(z_{n+1}\).
Advantages:
A simple idea:
\[ z_{n+1} = z^* - \lambda \nabla H(z^*) \]
where \(\lambda\) solves the constraint:
\[ H(z_{n+1}) = H(z_n). \]
This technique is general, easy to implement, but:
Many systems satisfy:
\[ \nabla \cdot f(z) = 0. \]
By Liouville’s theorem, their flow preserves volume:
\[ \det D\varphi_t = 1. \]
A numerical method is volume-preserving if:
\[ \det D\Phi_h = 1. \]
\[ f = f^{(1)} + \cdots + f^{(k)}, \quad \nabla \cdot f^{(i)} = 0. \]
Symplectic ↔ Energy-preserving methods form two distinct families:
In general:
Only in rare cases does a method preserve both (AVF on certain quadratic Hamiltonians).
Below we compare energy error for the harmonic oscillator: