Symplectic Runge–Kutta (SRK) methods form the second main family of geometric integrators after splitting methods. They act directly on the system of ODEs instead of decomposing the vector field. Their key virtues:
The classical example: the Gauss–Legendre collocation integrators, which are implicit, A-stable, symmetric, symplectic, and of order 2s for s stages.
Consider an autonomous Hamiltonian system on \(\mathbb{R}^{2m}\):
\[ \dot{y} = J^{-1} \nabla H(y), \qquad J = \begin{pmatrix} 0&-I\\ I&0 \end{pmatrix}. \]
A Runge–Kutta method with Butcher tableau:
| c_i | a_{ij} | ||
| b_1 | b_2 | b_3 | |
is symplectic if and only if:
\[ b_i a_{ij} + b_j a_{ji} - b_i b_j = 0 \qquad \text{for all } i,j. \tag{5.1.1} \]
This set of bilinear conditions characterises exactly those RK methods that preserve the canonical two-form \(\omega = dq \wedge dp\).
The Gauss–Legendre nodes c_i in [0,1] are roots of:
\[ P_s(2c-1)=0, \]
where P_s is the Legendre polynomial of degree s.
The associated collocation method has:
The Butcher coefficients are:
\[ a_{ij} = \int_0^{c_i} \ell_j(\tau)\, d\tau,\qquad b_j = \int_0^1 \ell_j(\tau)\, d\tau, \]
where \(\ell_j\) are the Lagrange basis polynomials at the Gauss nodes.
These are the highest-order (per stage) RK methods known, and are optimal for many Hamiltonian problems.
SRK methods arise naturally from a discrete variational principle. For a Lagrangian L(q,\dot{q}):
\[ S[q] = \int_0^h L(q,\dot q)\,dt. \]
Approximate the trajectory inside each timestep by a polynomial interpolant defined at stage points:
\[ q(t) \approx \sum_{j=1}^s q_j \ell_j(t/h). \]
Then extremise the discrete action:
\[ S_h(q_1,\dots,q_s) = h \sum_{i=1}^s b_i\, L(q_i, \dot{q}_i). \]
The resulting discrete Euler–Lagrange equations produce exactly the SRK update.
This variational origin guarantees symplecticity, momentum conservation, and excellent energy behaviour.
For s=2, the Gauss nodes are:
\[ c_{1,2} = \tfrac{1}{2} \mp \tfrac{\sqrt{3}}{6}. \]
The resulting method has:
This is one of the best all-purpose geometric ODE solvers available.
A hallmark of Gauss collocation methods: their superior phase accuracy for Hamiltonian oscillatory systems. For a harmonic oscillator:
\[ y' = J^{-1} \nabla \tfrac{1}{2}(q^2+p^2), \]
the Gauss-2 method yields:
Comparison of energy drift:
Gauss methods are algebraically stable and A-stable. Thus they handle:
They are often the preferred integrator when:
Runge–Kutta (RK) methods admit an algebraic representation of their action on differential equations through B-series. These are formal expansions indexed by rooted trees, which encode precisely how an RK method approximates the Taylor expansion of the exact flow.
This section introduces the B-series formalism, describes order conditions, and explains why symplectic Runge–Kutta methods naturally restrict to symplectic B-series, leading to:
For an autonomous ODE \(\dot{y}=f(y)\), the Taylor expansion of the exact flow can be written using rooted trees.
\[ y(h) = y(0) + h f + \frac{h^2}{2} f'f + \frac{h^3}{6} \big(f''(f,f) + f'(f'f)\big) + \cdots. \]
Each term corresponds to a rooted tree. For example:
For a tree \(\tau\), the associated elementary differential is denoted
\[ F(\tau)(y) \in \mathbb{R}^m. \]
The exact flow has a B-series:
\[ \phi_h(y) = y + \sum_{\tau\in\mathcal{T}} \frac{h^{|\tau|}}{\sigma(\tau)} F(\tau)(y). \tag{5.2.1} \]
Any (consistent) RK method produces a B-series:
\[ \Psi_h(y) = y + \sum_{\tau\in\mathcal{T}} h^{|\tau|}\alpha(\tau)F(\tau)(y), \tag{5.2.2} \]
where \(\alpha(\tau)\) are algebraic expressions in a_{ij},b_i,c_i.
The method has order p if and only if:
\[ \alpha(\tau) = \frac{1}{\sigma(\tau)} \qquad \text{for all } |\tau|\le p. \tag{5.2.3} \]
B-series form a group under composition: the Butcher group \(G_B\).
Key facts:
Gauss collocation methods correspond to the exponential of the truncated pre-Lie product associated with the vector field.
For Hamiltonian ODEs \(\dot{y}=J^{-1}\nabla H(y)\), the exact B-series satisfies:
\[ \phi_h^*(\omega) = \omega, \]
where \(\omega\) is the symplectic form.
An RK B-series is symplectic if and only if the Butcher coefficients satisfy:
\[ b_i a_{ij} + b_j a_{ji} - b_i b_j = 0, \tag{5.1.1 revisited} \]
These constraints force the B-series coefficients into the Hamiltonian sub-group of the Butcher group.
Consequences:
A Gauss–Legendre collocation method with s stages satisfies all order conditions up to order 2s. This follows because:
This gives the celebrated order-barrier optimality result:
\[ \text{For any s-stage RK method: } \quad p \le 2s. \]
Gauss methods achieve this bound with equality.
Rooted trees form a pre-Lie algebra:
\[ \tau_1 \triangleright \tau_2 = \text{sum over graftings of the root of } \tau_1 \text{ onto each node of } \tau_2. \]
The exact solution is the Lie exponential of the vector field in this algebra:
\[ \phi_h = \exp_{\triangleright}(h f). \]
Gauss methods are the truncated exponential of this pre-Lie product, preserving symplecticity by restriction to the Hamiltonian subalgebra.
This provides the deepest algebraic explanation for their optimal order, symmetry, and geometric behaviour.
Symmetry (time reversibility), symplecticity, and energy conservation are three of the most important geometric properties for numerical integrators of Hamiltonian systems. While no Runge–Kutta method can preserve all invariants exactly, certain RK families preserve the symplectic form, quadratic invariants, and/or time-reversal symmetry, which in turn guarantees excellent long-time behaviour.
This section develops:
A one-step method \(\Phi_h\) is symmetric (or time reversible) if:
\[ \Phi_{-h} \circ \Phi_h = \mathrm{id}. \tag{5.3.1} \]
For RK methods with tableau \((A,b,c)\), symmetry requires:
\[ c_i = 1 - c_{s+1-i}, \qquad b_i = b_{s+1-i}, \qquad a_{ij} = b_{j} - a_{s+1-i,s+1-j}. \tag{5.3.2} \]
Examples:
Key fact:
A symmetric method automatically has even order.
When an RK method is both symplectic and symmetric:
This is why Gauss methods are the gold standard:
\[ \text{Gauss–Legendre RK} \quad \text{= symmetric + symplectic + order } 2s. \]
This combination yields qualitatively correct behaviour over astronomical time spans.
A symplectic RK method applied to \(\dot{y}=J^{-1}\nabla H(y)\) exactly solves a nearby Hamiltonian system:
\[ \dot{y} = J^{-1}\nabla \tilde{H}(y), \]
where the modified Hamiltonian has a formal expansion:
\[ \tilde{H} = H + h^p H_{p+1} + h^{p+2} H_{p+3} + \cdots. \tag{5.3.3} \]
Key consequences:
No ordinary RK method (explicit or implicit) can preserve a general Hamiltonian exactly. But several specialized families achieve exact energy conservation by modifying the discretisation:
\[ \frac{y_{n+1}-y_n}{h} = \int_0^1 f\big(\theta y_{n+1} + (1-\theta)y_n\big)\, d\theta. \]
For Hamiltonian systems, AVF preserves energy exactly:
\[ H(y_{n+1}) = H(y_n). \]
Choose a discrete gradient \(\bar{\nabla}H(y_{n+1},y_n)\) satisfying:
\[ H(y_{n+1}) - H(y_n) = \bar{\nabla}H(y_{n+1},y_n)^{\!\top} (y_{n+1}-y_n). \tag{5.3.4} \]
Then define an update:
\[ \frac{y_{n+1} - y_n}{h} = J^{-1} \bar{\nabla}H(y_{n+1},y_n). \]
This preserves energy exactly for any Hamiltonian.
There is a fundamental incompatibility: an RK method cannot be both symplectic and energy-preserving unless the Hamiltonian is quadratic.
Thus one chooses:
| Method | Symplectic? | Symmetric? | Exact Energy? | Order |
|---|---|---|---|---|
| Gauss–s | Yes | Yes | No | 2s |
| Radau IIA | No | No | No | 2s-1 |
| Implicit midpoint | Yes | Yes | No | 2 |
| AVF / DG Methods | No | Depends | Yes | 2 |
Many important ODEs have a natural partitioned structure. Examples include:
For such systems, Partitioned Runge–Kutta (PRK) methods assign different RK tableaux to each component, while maintaining a coupled update.
When the system is Hamiltonian, special algebraic conditions give rise to Symplectic PRK (SPRK) methods.
Consider a partitioned system:
\[ \begin{aligned} q' &= f(q,p), \\ p' &= g(q,p). \end{aligned} \tag{5.4.1} \]
A partitioned Runge–Kutta method is defined by two tableaux:
|
|
The stage values satisfy:
\[ \begin{aligned} Q_i &= q_n + h\sum_{j=1}^s a_{ij}\,f(Q_j,P_j), \\ P_i &= p_n + h\sum_{j=1}^s \tilde{a}_{ij}\,g(Q_j,P_j), \end{aligned} \tag{5.4.2} \]
and the step update is:
\[ \begin{aligned} q_{n+1} &= q_n + h\sum_{i=1}^s b_i\,f(Q_i,P_i), \\ p_{n+1} &= p_n + h\sum_{i=1}^s \tilde{b}_i\,g(Q_i,P_i). \end{aligned} \tag{5.4.3} \]
For a Hamiltonian system in canonical form:
\[ q' = +\nabla_p H(p,q), \qquad p' = -\nabla_q H(p,q), \tag{5.4.4} \]
the exact flow is symplectic:
\[ \Phi_h^*(dq \wedge dp) = dq \wedge dp. \]
A partitioned RK method is symplectic if and only if the coefficients satisfy:
\[ b_i \tilde{a}_{ij} + \tilde{b}_j a_{ji} = b_i \tilde{b}_j, \qquad 1\le i,j\le s. \tag{5.4.5} \]
This is the symplectic PRK condition.
Special cases:
For separable systems:
\[ H(p,q) = T(p) + V(q), \tag{5.4.6} \]
the ODE is:
\[ q' = +T_p(p), \qquad p' = -V_q(q). \tag{5.4.7} \]
Since the two components depend only on one variable each, SPRK methods simplify dramatically.
Two 1-stage PRK methods:
Verlet is the most widely used SPRK of order 2.
Taking the tableaux:
\[ A = \begin{pmatrix}0\end{pmatrix},\quad b=\begin{pmatrix}1\end{pmatrix}, \qquad \tilde{A} = \begin{pmatrix}\frac12\end{pmatrix},\quad \tilde{b}=\begin{pmatrix}1\end{pmatrix}, \tag{5.4.8} \]
yields the scheme:
\[ \begin{aligned} p_{n+\frac12} &= p_n - \tfrac{h}{2} \nabla V(q_n),\\ q_{n+1} &= q_n + h\, \nabla_p T(p_{n+\frac12}),\\ p_{n+1} &= p_{n+\frac12} - \tfrac{h}{2} \nabla V(q_{n+1}). \end{aligned} \tag{5.4.9} \]
This is exactly the Stoermer–Verlet integrator, which is:
PRK order conditions are described by two-coloured rooted trees:
Exactly as in the classical B-series case:
\[ \Psi_h = y + \sum_{\tau\in \mathcal{T}_{2}} h^{|\tau|}\,\alpha(\tau)\,F(\tau), \tag{5.4.10} \]
where \(\mathcal{T}_{2}\) is the set of two-coloured rooted trees.
Using coloured trees splits the order conditions into mixed and diagonal classes, allowing high-order symmetric & symplectic PRK methods to be constructed systematically.
Examples include:
A typical high-order symmetric composition:
\[ \Psi_h^{(4)} = \Phi_{\gamma_1 h}\circ \Phi_{\gamma_2 h}\circ \Phi_{\gamma_1 h}, \tag{5.4.11} \]
with \(\gamma_1 = 1/(2-2^{1/3})\), \(\gamma_2 = -2^{1/3}\gamma_1\).
Such compositions retain symplecticity and symmetry by construction.
Many evolution equations have the semilinear structure:
\[ y' = Ly + N(y), \qquad y(0)=y_0 \tag{5.5.1} \]
where:
Typical examples include:
In such cases, classical RK methods suffer order reduction and stability problems. Exponential Runge–Kutta (ERK) and integrating factor / Lawson methods resolve this by treating the stiff linear part exactly.
The exact flow of the linear part is:
\[ \phi_L^h(y) = e^{hL}y. \]
Introduce the change of variables:
\[ v(t) = e^{-tL} y(t). \tag{5.5.2} \]
Then:
\[ v' = e^{-tL}N\!\big(e^{tL}v\big), \tag{5.5.3} \]
which is nonstiff if N is nonstiff. Ordinary RK applied to v'(t) gives the classical Lawson method.
This is the foundation of all exponential integrators.
Apply an s-stage RK method to (5.5.3):
\[ v_{n+1} = v_n + h\sum_{i=1}^s b_i\,e^{-c_i h L} N\!\big(e^{c_i h L} V_i\big), \tag{5.5.4} \]
with stages:
\[ V_i = v_n + h\sum_{j=1}^s a_{ij}\,e^{-c_j h L} N\!\big(e^{c_j hL} V_j\big). \tag{5.5.5} \]
Return to y via:
\[ y_{n+1} = e^{hL} v_{n+1}. \tag{5.5.6} \]
Advantages:
ERK methods express the solution using the variation-of-constants formula:
\[ y(t+h)=e^{hL}y(t) +\int_0^h e^{(h-\tau)L} N(y(t+\tau))\,d\tau. \tag{5.5.7} \]
This integral naturally leads to the phi functions:
\[ \varphi_0(z)=e^z,\qquad \varphi_1(z)=\frac{e^z-1}{z},\qquad \varphi_2(z)=\frac{e^z - z - 1}{z^2}, \dots \tag{5.5.8} \]
Higher-order ERK methods are linear combinations of \(\varphi_k(hL)\) acting on nonlinear terms.
\[ Y_i = e^{c_i h L} y_n + h\sum_{j=1}^s a_{ij}(hL)\,N(Y_j), \tag{5.5.9} \]
and
\[ y_{n+1} = e^{hL}y_n + h\sum_{i=1}^s b_i(hL)\,N(Y_i), \tag{5.5.10} \]
where the coefficients are matrix functions:
\[ a_{ij}(hL) = \int_0^{c_i}e^{(c_i-\tau)hL}\ell_j(\tau)\,d\tau, \tag{5.5.11} \]
and similarly for b_i(hL).
If the linear part L generates a symmetry group of the system, exponential integrators inherit invariances automatically:
\[ L^T = -L \;\Longrightarrow\; \|e^{hL}y\|=\|y\|. \]
If L is a Hamiltonian matrix:
\[ e^{hL} \in \mathrm{Sp}(2n). \tag{5.5.12} \]
Then exponential integrators preserve linear symplectic structure automatically.
Example: rotation, Schrödinger flow, Maxwell–Bloch dynamics.
Because ERK treats the stiff part exactly, its stability region includes the entire left half-plane for linear test problems y' = \lambda y with stiff \(\lambda\). This eliminates order reduction for parabolic PDEs.
ERK is closely linked with:
For stiff linear operators, ERK methods are effectively splitting methods with exact L-flow and approximated nonlinear flow.
For non-autonomous linear systems
\[ y'(t) = A(t)\, y(t), \tag{5.6.1} \]
the exact solution is given by a time-ordered exponential:
\[ y(t+h) = \mathcal{T} \exp\!\left(\int_t^{t+h}A(\tau)\,d\tau\right)y(t), \tag{5.6.2} \]
which lies in the Lie group generated by the Lie algebra containing A(t). This structure must be preserved for:
U(t) ∈ U(n))R(t) ∈ SO(3))M(t) ∈ Sp(2n))This section develops three families of integrators that preserve the Lie-group geometry:
Magnus (1954) showed that the solution of (5.6.1) can be written exactly as:
\[ y(t+h) = \exp\!\left(\Omega(t,h)\right)y(t), \tag{5.6.3} \]
where
\[ \Omega(t,h)=\int_t^{t+h}A(\tau)d\tau - \frac12 \int_t^{t+h}\!\!\int_t^{\tau_1} [A(\tau_1),A(\tau_2)]\,d\tau_2 d\tau_1 + \cdots. \tag{5.6.4} \]
Properties:
The 2nd-order Magnus method:
\[ \Omega^{[2]} = h\, A(t+\tfrac{h}{2}). \tag{5.6.5} \]
The 4th-order Magnus method:
\[ \Omega^{[4]} = h\Big(\tfrac12(A_1+A_2) - \tfrac{\sqrt{3}}{12}h [A_2, A_1]\Big), \tag{5.6.6} \]
where A_1, A_2 are samples at Gauss–Legendre nodes.
Standard Magnus integrators require evaluating matrix commutators [A_i, A_j] and nested commutators, which can be expensive.
CFM integrators avoid commutators entirely: instead of \(\Omega\) being a sum of Lie brackets, we use sums of exponentials:
\[ y_{n+1} = \exp\!\left(h\sum_k \alpha_k A(t_n + c_k h)\right) y_n, \tag{5.6.7} \]
with coefficients chosen to match the Magnus expansion up to order p.
Advantages:
\[ y_{n+1} = \exp\!\big(h(\tfrac12 A_1 + \tfrac12 A_2)\big)\; y_n, \tag{5.6.8} \]
with A_1, A_2 chosen at Gauss nodes. No commutators appear, yet the method is 4th order.
For semilinear systems:
\[ y' = A(t) y + N(y,t), \tag{5.6.9} \]
one can combine a Magnus integrator for the linear part with a Runge–Kutta method for the nonlinear part:
\[ y_{n+1} = \exp\!\left(\Omega(t_n,h)\right) \left( y_n + h\sum_{i=1}^s b_i\, \tilde{N}_i \right), \tag{5.6.10} \]
where \(\tilde{N}_i\) are nonlinear stages evaluated in a transformed frame:
\[ \tilde{N}_i = N\Big( \exp(\Omega_i)\big(y_n + \cdots\big),\, t_n + c_i h \Big). \tag{5.6.11} \]
These hybrid methods:
Blanes, Casas, Iserles, and others developed high-order methods combining:
A typical 6th-order CFM–RK method uses:
\[ y_{n+1} = \exp\!\left( h(\alpha_1 A_1 + \alpha_2 A_2 + \alpha_3 A_3) \right) \Big( y_n + h\sum b_i N(Y_i) \Big). \tag{5.6.12} \]
All coefficients arise from matching the Magnus series through order 6, but no commutators are used.
If the exact flow satisfies:
\[ y(t) \in G, \qquad G = \exp(\mathfrak{g}), \]
then Magnus and CFM–RK methods ensure:
\[ y_{n+1} = \exp(\xi_n) y_n \in G, \tag{5.6.13} \]
because \(\xi_n \in \mathfrak{g}\) for all steps.
Thus: