Many differential equations encountered in physics, chemistry, astronomy, and PDE simulation can be written as the sum of simpler vector fields:
\[ \dot{z} = (A + B)(z) \]
where the flows of \(A\) and \(B\) can be computed exactly or at low computational cost.
For example, for separable Hamiltonian systems with
\[ H(q,p) = T(p) + V(q), \]
the Hamiltonian vector field splits into:
\[ A : \dot{q} = \nabla_p T(p),\quad \dot{p} = 0, \qquad B : \dot{q} = 0,\quad \dot{p} = -\nabla_q V(q). \]
Each part is trivially solvable; the full system is not.
Splitting methods use the exact flows of these pieces to build highly accurate and structure-preserving integrators.
The flow of \( \dot{z} = A(z) \) over time \(h\) is denoted:
\[ \varphi_h^A = e^{hA}. \]
Similarly for \(B\). These exponentials satisfy:
\[ \varphi_h^{A+B} = e^{h(A+B)}. \]
In general:
Splitting provides accurate approximations to the combined flow using compositions of the partial flows.
The simplest splitting is the Lie–Trotter product formula:
\[ e^{h(A+B)} = \lim_{n\to\infty} \left( e^{hA/n} e^{hB/n} \right)^n. \]
For numerical integration with step size \(h\), we approximate:
\[ \Phi_h^{\mathrm{LT}} = e^{hA} e^{hB}. \]
This is a first-order integrator.
Error analysis uses the Baker–Campbell–Hausdorff (BCH) formula:
\[ e^{hA} e^{hB} = e^{h(A+B) + \frac{h^2}{2}[A,B] + O(h^3)}. \]
Thus the local error is:
\[ \mathcal{O}(h^2). \]
Lie–Trotter splitting is widely used in:
The classic Strang splitting is:
\[ \Phi_h^{\mathrm{S}} = e^{\frac{h}{2}A}\, e^{hB}\, e^{\frac{h}{2}A}. \]
BCH gives:
\[ e^{\frac{h}{2}A}\, e^{hB}\, e^{\frac{h}{2}A} = e^{h(A+B) + \frac{h^3}{12}[A,[A,B]] - \frac{h^3}{24}[B,[A,B]] + \cdots}. \]
Therefore Strang splitting is:
It is the basis of the Verlet method in molecular dynamics.
Splitting methods are canonical transformations when applied to Hamiltonian systems where each partial flow is symplectic.
In coordinates:
\[ \varphi_h^{T(p)} : \begin{cases} q \mapsto q + h\, \nabla_p T(p),\\[4pt] p \mapsto p, \end{cases} \qquad \varphi_h^{V(q)} : \begin{cases} q \mapsto q,\\ p \mapsto p - h\, \nabla_q V(q). \end{cases} \]
The geometric picture:
In effect, Strang splitting alternates between geodesics of A and B.
Splitting methods are ubiquitous.
\[ i\partial_t \psi = (T + V)\psi. \]
Lie–Trotter and Strang splitting:
Splitting methods naturally apply to:
Below we compare first-order Lie–Trotter vs second-order Strang for the simple separable harmonic oscillator.
Having introduced Lie–Trotter (first order) and Strang (second order) splittings, we now explain how to construct arbitrarily high-order geometric integrators by composition.
The fundamental idea is:
A symmetric second-order method can be composed with itself (with carefully chosen coefficients) to produce a higher-order symmetric method.
This is a cornerstone of modern geometric integration and of Trotter–Suzuki decompositions in quantum simulation.
Start with a symmetric 2nd-order method \(S_h\), such as Strang splitting.
A symmetric composition takes the form:
\[ \Psi_h = S_{\alpha_1 h} S_{\alpha_2 h} \cdots S_{\alpha_s h} S_{\alpha_s h} \cdots S_{\alpha_2 h} S_{\alpha_1 h}. \]
Symmetry → automatic cancellation of all even-order error terms. Therefore:
\[ \Psi_h \text{ is of order } p \; \Longleftrightarrow\; \text{all odd-order error terms up to } p-1 \text{ vanish}. \]
Order conditions reduce to a small algebraic system in the coefficients \(\alpha_i\).
Let the local error of the 2nd-order symmetric method be:
\[ S_h = \exp\big( h(A+B) + h^3 E_3 + h^5 E_5 + \cdots \big). \]
For the composition:
\[ \Psi_h = \prod_{i=1}^s S_{\alpha_i h} \prod_{i=s}^1 S_{\alpha_i h}, \]
the BCH formula implies the 3rd-order terms vanish automatically due to symmetry.
Imposing fourth order reduces to:
\[ \sum_{i=1}^s \alpha_i^3 = \frac{1}{24}. \]
And the consistency condition:
\[ 2\sum_{i=1}^s \alpha_i = 1. \]
Similar polynomial conditions determine 6th, 8th, 10th order, etc.
Take:
\[ \Psi_h^{(4)} = S_{\gamma h} S_{\delta h} S_{\gamma h}, \]
with:
\[ \gamma = \frac{1}{2 - 2^{1/3}}, \qquad \delta = -\,\frac{2^{1/3}}{2 - 2^{1/3}}. \]
Properties:
Suzuki discovered a recursive formula for producing order \(2k+2\) methods from order \(2k\) ones:
\[ S^{(2k+2)}_h = S^{(2k)}_{\alpha h} S^{(2k)}_{\beta h} S^{(2k)}_{\alpha h}, \]
with coefficients chosen to cancel higher-order errors. Example for 4th order:
\[ \alpha = \frac{1}{2 - 2^{1/3}},\quad \beta = -\frac{2^{1/3}}{2 - 2^{1/3}}. \]
The order hierarchy continues indefinitely:
But coefficients rapidly become:
Crucial theorem:
If all coefficients αᵢ are real and non-negative, no splitting method of order > 2 exists for general operators A+B.
Thus:
Modern developments use complex coefficients:
These are particularly useful in:
One common 4th-order complex scheme uses:
\[ \alpha = 0.5 + \frac{i}{2\sqrt{3}},\qquad \bar{\alpha} = 0.5 - \frac{i}{2\sqrt{3}}. \]
These satisfy necessary order conditions with purely nonnegative real parts.
Below we compare:
The example is a perturbed harmonic oscillator; error in phase is displayed.
Higher-order splitting methods obtained by composition (Yoshida–Suzuki) grow in cost and often suffer from large error constants, negative coefficients, and stability limitations. Optimised splitting methods aim to reduce the leading error terms by solving constrained nonlinear optimisation problems on the splitting coefficients.
This allows one to build practical high-order schemes with:
These methods appear in molecular dynamics, quantum simulation, semiclassical analysis, and geometric PDE solvers.
Let the exact flow be:
\[ \varphi_h = e^{h(A+B)}. \]
A splitting with substeps \(\alpha_i\) approximates:
\[ \Phi_h = \prod_{i=1}^s e^{\alpha_i h A} e^{\beta_i h B}. \]
Using BCH:
\[ \Phi_h = \exp\!\Big( h(A+B) + h^3 C_3 + h^5 C_5 + \cdots \Big). \]
A method of order \(p\) cancels all error terms \(C_3,\dots,C_{p-1}\). Optimised methods attempt to minimise the norm of the first nonzero term \(C_{p+1}\).
Given a target order (often 4, 6, or 8), one seeks coefficients that:
The resulting optimal compositions are much more efficient than the straightforward Yoshida–Suzuki iteration.
Strang splitting (order 2) has error \(\mathcal{O}(h^3)\). A 4th-order symmetric composition uses the sequence:
\[ \Phi_h = e^{a_1 h A} e^{b_1 h B} e^{a_2 h A} e^{b_2 h B} e^{a_2 h A} e^{b_1 h B} e^{a_1 h A}. \]
Conditions for order 4 include:
\[ a_1 + a_2 = \frac{1}{2},\qquad b_1 + b_2 = 1, \]
and cancellation of:
\[ [A,[A,B]],\quad [B,[A,B]]. \]
Optimised values (Blanes–Casas, 2002):
\[ a_1 = 0.5153528374311229364,\quad a_2 = -0.085782019412973646,\quad b_1 = 0.4415830236164665242,\quad b_2 = 0.1184634424268122861. \]
These give significantly smaller error constants than Yoshida’s triple jump.
For 6th-order splittings, order conditions involve nested commutators:
\[ [A,[A,[A,B]]],\quad [B,[A,[A,B]]],\quad [B,[B,[A,B]]],\dots \]
The number of polynomial constraints grows rapidly:
Blanes–Casas–Murua and Kahan–Li constructed optimised solutions with far superior performance compared to Suzuki compositions.
Example: an optimised 6th-order, 7-stage symmetric splitting:
\[ \Phi_h^{(6)} = S_{\alpha_1 h} S_{\alpha_2 h} \dots S_{\alpha_7 h} S_{\alpha_7 h} \dots S_{\alpha_2 h} S_{\alpha_1 h}, \]
with the \(\alpha_i\) chosen by solving a constrained optimisation problem.
For Hamiltonian systems, the modified Hamiltonian expands as:
\[ \tilde{H} = H + h^p H_{p+1} + h^{p+2} H_{p+3} + \cdots. \]
Optimised methods minimise:
\[ \|H_{p+1}\| \quad \text{or} \quad \sum_i |\alpha_i|. \]
The result is far better long-time energy behaviour than naïve compositions.
For PDEs such as:
the flow of \(A\) (typically linear) is unbounded; large negative coefficients make compositions highly unstable.
Thus optimised schemes seek:
Special splittings exploit system-specific structure:
In these cases “optimal” means minimising error relative to a problem-specific asymptotic regime.
Below is a visual comparison of:
We display energy error for the harmonic oscillator.
In many applications — quantum control, laser-driven Schrödinger systems, NMR, accelerator physics, molecular dynamics with thermostats, celestial mechanics with perturbations, and time-dependent PDEs — the governing equation is non-autonomous:
\[ \dot{u}(t) = \big(A(t) + B(t)\big)\,u(t),\qquad u(0)=u_0. \]
The exact solution involves the time-ordered exponential:
\[ u(t) = \mathcal{T}\exp\!\left(\int_0^t (A(\tau)+B(\tau))\,d\tau\right)u_0. \]
Because \(A(t)\) and \(B(t)\) fail to commute at different times, standard autonomous splittings do not apply directly.
This section establishes the modern theory of splitting for non-autonomous systems, combining:
The time-ordered exponential can be rewritten using the Magnus expansion:
\[ u(t) = \exp\big(\Omega(t)\big)\,u(0), \]
where:
\[ \Omega(t) = \int_0^t A(\tau)\,d\tau + \tfrac{1}{2}\!\int_0^t\!\!\int_0^{\tau_1}[A(\tau_1),A(\tau_2)]\,d\tau_2\,d\tau_1 + \cdots. \]
This provides a universal autonomous lifting of non-autonomous systems: a time-dependent problem is equivalent to an autonomous flow generated by a (possibly infinite) Lie series.
In practice, truncations provide high-order approximations preserving unitarity or symplecticity exactly.
For autonomous systems:
\[ e^{h(A+B)} \approx e^{hA} e^{hB}. \]
For non-autonomous systems:
\[ e^{\int A(t)\,dt} e^{\int B(t)\,dt} \neq e^{\int (A+B)(t)\,dt}. \]
Commutators like [A(t_1),A(t_2)] and [A(t),B(s)] prevent naive splitting.
Modern solutions involve:
Consider:
\[ \dot{u} = A(t)u + B(t)u. \]
A robust second-order method evaluates operators at midpoints:
\[ \Phi_h^{(2)} = \exp\!\big(\tfrac{h}{2}A(t+\tfrac{h}{2})\big) \exp\!\big(h B(t+\tfrac{h}{2})\big) \exp\!\big(\tfrac{h}{2}A(t+\tfrac{h}{2})\big). \]
This coincides with a truncated Magnus expansion and produces a time-reversible method for symmetric systems.
Standard Magnus methods require evaluating commutators. CFMIs avoid this by expressing solutions as products of exponentials:
\[ \Phi_h = \exp\!\big(a_1 h\, A(t_1)\big)\, \exp\!\big(a_2 h\, A(t_2)\big)\cdots \]
The idea is to choose quadrature nodes \(t_i\) and coefficients \(a_i\) so that the method matches the Magnus expansion up to order \(p\) without computing commutators.
These methods are central in quantum computing, laser physics, and time-dependent PDEs.
For systems decomposing as:
\[ \dot{u} = A(t)u + B(t)u, \]
where each subflow is easily computable, hybrid methods combine Magnus approximations for time dependence with splitting for operator structure.
Example (4th order):
\[ \Phi_h^{(4)} = \exp\!\left( a_1 h A(t_1) \right) \exp\!\left( b_1 h B(t_1) \right) \exp\!\left( a_2 h A(t_2) \right) \exp\!\left( b_1 h B(t_1) \right) \exp\!\left( a_1 h A(t_1) \right), \]
where the nodes \(t_1,t_2\) approximate the Gauss–Legendre nodes.
These methods:
A beautiful geometric trick introduces an extended phase-space:
\[ (q,p,t,p_t), \]
with extended Hamiltonian:
\[ K(q,p,t,p_t) = H(q,p,t) + p_t. \]
The resulting system is autonomous:
\[ \dot{t} = 1,\qquad \dot{p_t} = -\frac{\partial H}{\partial t},\qquad \dot{q},\dot{p}\text{ as usual}. \]
Any autonomous splitting method — symplectic, symmetric, or optimised — now applies.
For the Schrödinger equation with time-dependent Hamiltonian:
\[ i\frac{d\psi}{dt} = H(t)\psi, \]
splitting + Magnus is the dominant method for:
CFMIs give unitarity exactly, with no commutator cost.
Celestial mechanics and plasma physics contain Hamiltonians of type:
\[ H(q,p,t) = T(p) + V(q) + \epsilon\,W(q,t), \]
where \(W\) is a periodic or quasiperiodic driving term. Magnus–splitting hybrids offer:
We solve:
\[ \ddot{q} + q = 0.1\sin(3 t). \]
In first-order form:
\[ \dot{q}=p,\qquad \dot{p}=-q + 0.1\sin(3t). \]
After establishing splitting methods (Sections 4.1–4.4), we now move to the frontier: designing high-order, low-error, structure-preserving integrators for non-autonomous systems.
These schemes combine:
Their purpose: achieve high accuracy at minimal computational cost while preserving geometric structure (symplecticity, reversibility, unitarity, invariants).
In the non-autonomous setting, the Baker–Campbell–Hausdorff (BCH) expansion acquires time-dependent commutators. For operators A(t), B(t):
\[ e^{h A(t_1)} e^{h B(t_2)} = \exp\!\Big( h(A(t_1)+B(t_2)) + \tfrac{h^2}{2} [A(t_1),B(t_2)] + \cdots \Big). \]
High-order accuracy requires matching not only the autonomous BCH conditions but the generalised non-autonomous Magnus expansion.
This leads to order conditions expressed in terms of time-integrated commutators, e.g.:
\[ \int_0^h \!\! \int_0^{\tau_1} [A(\tau_1), A(\tau_2)] d\tau_2 d\tau_1. \]
Direct enforcement is intractable → hence the need for:
CFQI methods approximate the Magnus operator \(\Omega(h)\) via a product of exponentials:
\[ \Phi_h = \exp\!\left( h \sum_j a_{1j} A(t_j) \right) \exp\!\left( h \sum_j a_{2j} A(t_j) \right) \cdots \]
with quadrature nodes t_j chosen by Gauss–Legendre, Gauss–Lobatto, or optimised nodes.
The advantages:
A GL-Magnus integrator uses Gauss–Legendre quadrature to approximate:
\[
\Omega(h) \approx h \sum_j w_j A(t_j)
+ \frac{h^3}{2} \sum_{i
In practice:
GL-Magnus methods provide excellent error-per-cost behaviour for:
A non-autonomous symmetric composition:
\[
\Phi_h^{(p)}
= \prod_{i=1}^s
\exp\!\left( a_i h A(t+c_i h) \right)
\exp\!\left( b_i h B(t+c_i h) \right)
\]
satisfies:
Yoshida’s method generalises by choosing c_i, a_i, b_i to satisfy non-autonomous order conditions.
Modern optimised values (Blanes–Casas–Murua 2022+) minimise:
For stiff PDEs, negative real coefficients make schemes unstable.
A breakthrough: use complex coefficients with positive real part.
Let:
\[
a_i = \alpha_i + i\beta_i,\qquad
\operatorname{Re}(a_i) > 0.
\]
Then each flow remains stable for operators like:
These schemes maintain symmetry and reach orders 6–12 with no negative coefficients at all.
In highly oscillatory regimes (e.g. laser pulses, RF fields), time dependence can vary rapidly.
Adaptive schemes choose:
\[
\{t_j\},\quad \{w_j\}
\]
dynamically based on:
These are the current state of the art for quantum control and Hamiltonian microlocal analysis.
Another modern improvement: perform a low-order splitting step, then use:
\[
u^{(k+1)} = u^{(k)} + \text{SDC correction}
\]
using a spectrally accurate quadrature rule.
This yields arbitrarily high order while preserving the preferred splittings
in each correction stage.
This is powerful for:
The main error sources:
Modern methods choose coefficients to minimise a weighted combination:
\[
\varepsilon
= \|C_{p+1}\| + \lambda \|Q_{\text{err}}\| + \mu \sum_i |a_i|.
\]
This produces exceptionally efficient non-autonomous integrators for large-scale simulations.
4.5.4 High-Order Non-Autonomous Splitting via Symmetric Composition
4.5.5 Complex Time Steps and Forward-Time Constraints
4.5.6 Adaptive Quadrature Magnus–Splitting Hybrids
4.5.7 Splitting with Spectral Deferred Correction (SDC)
4.5.8 Error Budget and Practical Optimisation
4.5.9 Summary of the State of Practice
4.5.10 References for Section 4.5