In the previous chapters we developed geometric structure, symplecticity, Lie–algebraic flows, exponential integrators, and stability theory. We now introduce one of the most powerful and widely used classes of geometric integrators: splitting and composition methods.
These methods exploit the fact that many differential equations can be written as the sum of simpler vector fields whose individual flows are exactly integrable or preserve essential geometric properties.
Consider an evolution equation
\[ y' = (B + C)(y), \]
where the flows of the subfields B and C are individually computable:
\[ \Phi_B^t = e^{tB}, \qquad \Phi_C^t = e^{tC}. \]
Splitting methods approximate the combined flow by compositions of these simpler flows, e.g.
\[ e^{h(B+C)} \approx e^{hB} e^{hC}, \]
or symmetric/higher-order variants. Since each exponential preserves geometric structure (e.g. symplectic, orthogonal, or unitary), the composition inherits these invariants automatically.
Advantages of splitting approaches:
Starting from a basic (possibly low-order) splitting operator
\[ \mathcal{S}(h) = e^{hB}e^{hC}, \]
one constructs higher-order methods by nesting compositions:
\[ \mathcal{S}(a_1 h)\,\mathcal{S}(a_2 h)\cdots \mathcal{S}(a_m h). \]
Carefully chosen coefficients eliminate lower-order error terms in the Baker–Campbell–Hausdorff expansion, yielding high-order integrators with minimal additional computational cost.
Examples include:
Splitting methods reside at the intersection of:
Thus splitting and composition methods form a unifying framework for geometric time integration across physics, chemistry, applied mathematics, and computational engineering.
This chapter develops the foundational splitting methods:
Splitting is thus the backbone of a huge class of geometric integrators.
In this section we study the two fundamental splitting schemes:
We will:
Let \(y(t)\) solve an autonomous ODE on a manifold \(M\):
\[ \dot{y}(t) = f(y(t)), \qquad y(0) = y_0 \in M. \]
Suppose the vector field \(f\) admits a decomposition into simpler parts
\[ f = f^{[1]} + f^{[2]}, \]
such that the flows of the subsystems
\[ \dot{y} = f^{[1]}(y), \qquad \dot{y} = f^{[2]}(y) \]
can be computed (almost) exactly or very efficiently. Denote their flows by
\[ \varphi_h^{[1]} = \exp(h F^{[1]}), \qquad \varphi_h^{[2]} = \exp(h F^{[2]}), \]
where \(F^{[j]}\) are the corresponding differential operators \(F^{[j]} = f^{[j]}\!\cdot\nabla\). The full flow is
\[ \varphi_h = \exp\!\big(h(F^{[1]} + F^{[2]})\big). \]
A splitting method approximates \(\varphi_h\) by a composition of subflows, e.g.
\[ \Phi_h = \varphi_{a_1 h}^{[1]} \circ \varphi_{b_1 h}^{[2]} \circ \cdots \circ \varphi_{a_s h}^{[1]} \circ \varphi_{b_s h}^{[2]}, \]
with suitable coefficients \(a_j, b_j\). The simplest choices give the classical Lie–Trotter and Strang splittings.
Task marker: later sections will generalize this to multiple splits \(f = f^{[1]} + \cdots + f^{[m]}\) and to composition schemes of arbitrarily high order.
The Lie–Trotter splitting (or first-order splitting) uses the simple product of exponentials
\[ \Phi_h^{\mathrm{LT}} := \varphi_h^{[1]} \circ \varphi_h^{[2]} = \exp(h F^{[1]}) \exp(h F^{[2]}). \]
The corresponding one-step method is
\[ y_{n+1} = \Phi_h^{\mathrm{LT}}(y_n) = \varphi_h^{[1]}\big( \varphi_h^{[2]}(y_n)\big). \]
The Baker–Campbell–Hausdorff (BCH) formula gives
\[ \exp(h F^{[1]}) \exp(h F^{[2]}) = \exp\Big( h(F^{[1]} + F^{[2]}) + \tfrac{h^2}{2}[F^{[1]},F^{[2]}] + \mathcal{O}(h^3) \Big). \]
Thus the Lie–Trotter step is the exact flow of a modified vector field
\[ \tilde{F} = F^{[1]} + F^{[2]} + \tfrac{h}{2}[F^{[1]},F^{[2]}] + \mathcal{O}(h^2), \]
so the local truncation error is \(\mathcal{O}(h^2)\) and the method has global order 1.
Task marker: in BEA we will explicitly compute the modified Hamiltonian for Hamiltonian splittings and relate the commutator term to Poisson brackets.
The Strang splitting (Strang 1968) is the symmetric composition
\[ \Phi_h^{\mathrm{S}} := \varphi_{h/2}^{[1]} \circ \varphi_{h}^{[2]} \circ \varphi_{h/2}^{[1]} = \exp\!\left(\tfrac{h}{2}F^{[1]}\right) \exp\!\left(h F^{[2]}\right) \exp\!\left(\tfrac{h}{2}F^{[1]}\right). \]
The corresponding scheme is:
Using BCH and symmetry, one can show that the odd-order error terms cancel. More concretely,
\[ \Phi_h^{\mathrm{S}} = \exp\!\Big( h(F^{[1]}+F^{[2]}) + \mathcal{O}(h^3) \Big), \]
so the local error is \(\mathcal{O}(h^3)\) and the global error is \(\mathcal{O}(h^2)\): Strang splitting has order 2.
A key general principle is:
Self-adjoint (time-symmetric) one-step methods of order ≥ 1 are always of even order.
Strang splitting is self-adjoint, hence order 2 rather than 1.
Let \(\Phi_h^\star\) denote the adjoint method defined by
\[ \Phi_h^\star := \Phi_{-h}^{-1}. \]
A one-step method is self-adjoint (or time-symmetric) if \(\Phi_h^\star = \Phi_h\). For Strang splitting,
\[ \big(\Phi_h^{\mathrm{S}}\big)^{-1} = \varphi_{-h/2}^{[1]} \circ \varphi_{-h}^{[2]} \circ \varphi_{-h/2}^{[1]} = \Phi_{-h}^{\mathrm{S}}, \]
hence \(\Phi_h^{\mathrm{S}}\) is self-adjoint. This is crucial for good long-time behaviour, especially in reversible and Hamiltonian systems.
Let us summarize the main geometric features when the splitting is applied to a structured system, e.g. a Hamiltonian system with \(H = H^{[1]} + H^{[2]}\).
Task marker: in later sections we derive explicit modified Hamiltonians for separable systems \(H(q,p) = T(p) + V(q)\) under Strang splitting, and compare with standard symplectic Euler and Verlet.
A prototypical application is a separable Hamiltonian system
\[ H(q,p) = T(p) + V(q), \qquad (q,p) \in \mathbb{R}^m \times \mathbb{R}^m. \]
The Hamiltonian vector field splits as
\[ \dot{q} = \nabla_p T(p), \qquad \dot{p} = -\nabla_q V(q), \]
which we interpret as the sum of two subfields:
Their flows are explicitly available:
\[ \varphi_h^{[1]}(q,p) = \big(q + h\nabla_p T(p),\; p\big), \qquad \varphi_h^{[2]}(q,p) = \big(q,\; p - h\nabla_q V(q)\big). \]
Applying Strang splitting gives the well-known velocity Verlet / leapfrog scheme:
This is the cornerstone integrator in molecular dynamics, celestial mechanics, and Hamiltonian Monte Carlo.
In this section we:
These schemes are the building blocks for higher-order composition methods and for many modern applications (Hamiltonian Monte Carlo, quantum simulation, PDE splitting, etc.). Subsequent sections generalize these ideas to multi-part splittings and higher-order factorizations.
Strang splitting is the simplest nontrivial symmetric composition:
\[ \Phi_h^{\mathrm{S}} = \varphi^{[1]}_{h/2} \circ \varphi^{[2]}_{h} \circ \varphi^{[1]}_{h/2}, \]
which is second order because symmetry causes cancellation of all odd-order error terms in the BCH expansion. This section develops the general theory: how symmetric compositions lead to arbitrarily high-order geometric integrators.
Let \(\Psi_h\) be a basic method of order \(p\), and assume it is self-adjoint (i.e. reversible, symmetric):
\[ \Psi^{\star}_{h} = \Psi_h. \]
Consider a symmetric composition with real coefficients
\[ \Psi_h^{(m)} = \Psi_{a_1 h} \circ \Psi_{a_2 h} \circ \cdots \circ \Psi_{a_s h} \circ \Psi_{a_{s} h} \circ \cdots \Psi_{a_2 h} \circ \Psi_{a_1 h}. \tag{6.2.1} \]
Symmetry ensures the overall method is self-adjoint:
\[ \big(\Psi_h^{(m)}\big)^\star = \Psi_h^{(m)}. \]
As a consequence, all odd-order terms in the BCH expansion cancel.
Let the modified vector field of \(\Psi_h\) be
\[ \tilde{F}(h) = F + h^{p} E_{p} + h^{p+1} E_{p+1} + h^{p+2} E_{p+2} + \cdots. \tag{6.2.2} \]
Inserting scaled steps \(a_j h\) yields modified fields
\[ \tilde{F}(a_j h) = F + a_j^{p} h^{p} E_p + a_j^{p+1} h^{p+1} E_{p+1} + \cdots. \]
The BCH formula for the symmetric composition (6.2.1) has only even powers:
\[ \log\big(\Psi_h^{(m)}\big) = h F + h^{p+2} C_{p+2}(a_1,\ldots,a_s) + h^{p+4} C_{p+4}(a_1,\ldots,a_s) + \cdots. \tag{6.2.3} \]
The order conditions are the algebraic constraints on \(a_1,\ldots,a_s\) ensuring that
\[ C_{p+2} = C_{p+4} = \cdots = C_{P-1} = 0. \]
This yields an overall method of order \(P\).
Task marker: later sections will list explicit order conditions for splitting and composition up to 8th order using the Lyndon word basis for the free Lie algebra.
Let \(\Psi_h\) be a symmetric method of order \(2m\). Yoshida (1990) showed that the composition
\[ \Psi_h^{(m+1)} = \Psi_{\alpha h} \circ \Psi_{\beta h} \circ \Psi_{\alpha h} \tag{6.2.4} \]
is a symmetric method of order \(2m + 2\) provided \(\alpha, \beta\) satisfy
\[ 2\alpha + \beta = 1, \qquad 2\alpha^{2m+1} + \beta^{2m+1} = 0. \tag{6.2.5} \]
These equations have the explicit solution
\[ \alpha = \frac{1}{2 - 2^{1/(2m+1)}}, \qquad \beta = -\,\frac{2^{1/(2m+1)}}{2 - 2^{1/(2m+1)}}. \tag{6.2.6} \]
For example:
| Base Order | Resulting Order | \(\alpha\) | \(\beta\) |
|---|---|---|---|
| 2 | 4 | \(\displaystyle \frac{1}{2-2^{1/3}} \approx 0.675603\) | \(1 - 2\alpha \approx -0.351207\) |
| 4 | 6 | \(\displaystyle \frac{1}{2-2^{1/5}} \approx 0.587\) | \(1 - 2\alpha \approx -0.174\) |
Starting from Strang (\(2^{nd}\) order), repeated Yoshida triple-jumps yield methods of orders 4, 6, 8, …
Suzuki (1991) extended Yoshida’s idea to general “fractal” compositions. A common example is:
\[ \Psi_h^{(p+2)} = \Psi_{\gamma h}^{(p)} \circ \Psi_{\delta h}^{(p)} \circ \Psi_{\gamma h}^{(p)} \tag{6.2.7} \]
with coefficients chosen to annihilate additional BCH terms. These schemes form a recursive family generating arbitrarily high order from a base symmetric integrator.
When the base integrator \(\Psi_h\) is itself a splitting, such as Strang:
\[ \Psi_h = \varphi_{h/2}^{[1]} \circ \varphi_h^{[2]} \circ \varphi_{h/2}^{[1]}, \tag{6.2.8} \]
the Yoshida and Suzuki compositions give high-order splitting methods. For example, the 4th-order symmetric splitting obtained from Strang is:
\[ \Phi_h^{(4)} = \Phi_{\alpha h}^{\mathrm{S}} \circ \Phi_{\beta h}^{\mathrm{S}} \circ \Phi_{\alpha h}^{\mathrm{S}}, \tag{6.2.9} \]
which expands (after substituting the Strang components) into a 7-factor composition of subflows \(\varphi^{[1]}\) and \(\varphi^{[2]}\).
These high-order splittings are widely used in:
Symmetric compositions retain:
High-order symplectic compositions are often competitive with (or superior to) high-order Runge–Kutta methods because the dominant error terms in Hamiltonian problems are often commutators, which compositions suppress efficiently.
In this section we developed the general framework of symmetric composition methods. Starting from a symmetric base integrator (typically Strang), one can construct high-order geometric integrators using Yoshida triple-jumps, Suzuki fractal compositions, or more general optimized compositions.
Key takeaways:
Classical symmetric compositions (Strang, Yoshida, Suzuki) generate high-order methods with universal algebraic rules, but they are rarely optimal for a specific class of problems. In this section we present the modern theory of optimised splitting schemes, including:
These methods represent the contemporary state-of-the-art and often outperform all traditional splitting schemes for demanding applications.
Consider symmetric splitting of a two-part system:
\[ \Phi_h = \prod_{j=1}^{s} \exp(a_j h F^{[1]}) \exp(b_j h F^{[2]}), \]
with the symmetry constraint
\[ a_j = a_{s+1-j}, \qquad b_j = b_{s+1-j}. \tag{6.3.1} \]
The BCH expansion of \(\log \Phi_h\) contains commutators \[ [F^{[1]},F^{[2]}], \quad [F^{[1]},[F^{[1]},F^{[2]}]], \quad \dots \] and the order conditions require vanishing of specific linear combinations of these.
Key barrier (Sheng 1989, Suzuki 1991):
No symmetric splitting of order higher than 2 exists with all real coefficients \(a_j, b_j \ge 0\).
In particular:
Thus:
High-order real-coefficient splittings are constructed by solving the BCH order conditions while minimising either:
A minimal 4th-order symmetric scheme exists with three exponentials of each operator:
\[ \Phi_h^{(4)} = e^{a_1 h F^{[1]}} e^{b_1 h F^{[2]}} e^{a_2 h F^{[1]}} e^{b_1 h F^{[2]}} e^{a_1 h F^{[1]}}. \tag{6.3.2} \]
with coefficients
\[ a_1 = \frac{1}{2 - 2^{1/3}}, \qquad a_2 = -\frac{2^{1/3}}{2 - 2^{1/3}}, \] \[ b_1 = \frac{1}{2}. \tag{6.3.3} \]
These are Yoshida's coefficients but optimised variants exist where the error constant is minimised (e.g., McLachlan’s 1995 scheme).
Task marker: Section 6.4 derives optimal real splittings for Hamiltonian dynamics using backward error analysis and modified Hamiltonians.
Allowing complex coefficients dramatically changes the landscape:
\[ \Phi_h = \prod_{j=1}^s \exp(a_j h F^{[1]}) \exp(b_j h F^{[2]}), \tag{6.3.4} \]
where now \[ a_j, b_j \in \mathbb{C}, \qquad \Re(a_j), \Re(b_j) \ge 0. \]
These schemes are stable for:
The Blanes–Casas “C4” scheme:
\[ a_1 = 0.5 + 0.28867513459 i, \qquad a_2 = 0.5 - 0.28867513459 i, \tag{6.3.5} \]
yields a 4th-order splitting with only two exponentials of each operator, outperforming all real-coefficient schemes in quantum applications.
For diffusion or parabolic systems:
\[ \dot{u} = A u + B u, \]
the operator \( e^{a_j h A} \) must satisfy:
\[ \Re(a_j) \ge 0, \qquad \Re(b_j) \ge 0, \tag{6.3.6} \]
otherwise the method is unstable. Therefore:
They constructed 4th, 6th, and 8th-order splitting schemes with all coefficients in the right half-plane. These dominate the field in smoothing-dominated PDEs.
For Hamiltonian systems with \(H = T(p) + V(q)\), the goal is to minimise:
These yield:
\[ \Phi_h^{\mathrm{OMF}} = e^{\lambda h B} e^{\frac{1}{2} h A} e^{(1-2\lambda) h B} e^{\frac{1}{2} h A} e^{\lambda h B}, \tag{6.3.7} \]
with \(\lambda \approx 0.1931833275\), which dramatically reduces energy drift.
This section presented the full modern landscape of optimised splittings:
Backward error analysis (BEA) provides the fundamental theoretical explanation for the excellent long-time behaviour of geometric integrators. Rather than analysing the discrete numerical flow directly, BEA interprets a numerical scheme as the exact flow of a modified differential equation:
\[ \Phi_h = \exp\!\big(h \tilde{F}(h)\big), \]
where the modified vector field admits a formal series expansion \[ \tilde{F}(h) = F + h F_2 + h^2 F_3 + h^3 F_4 + \cdots. \] For geometric methods—especially symplectic mappings—the modified system inherits the geometric structure (e.g.~Hamiltonian geometry) to all orders in \(h\).
Consider a splitting integrator for an autonomous system
\[ y' = (A + B)(y). \tag{6.4.1} \]
A general splitting has the form
\[ \Phi_h = \prod_{j=1}^s \exp(a_j h A)\, \exp(b_j h B). \tag{6.4.2} \]
The Baker–Campbell–Hausdorff expansion provides
\[ \log(\Phi_h) = h(A+B) + h^2 \alpha_1[A,B] + h^3\big(\alpha_2[A,[A,B]] + \alpha_3[B,[A,B]]\big) + h^4 \cdots. \tag{6.4.3} \]
Let the modified vector field be defined by
\[ h\tilde{F}(h) = \log(\Phi_h). \]
Then
\[ \tilde{F}(h) = A + B + h \alpha_1[A,B] + h^2\big(\alpha_2[A,[A,B]] + \alpha_3[B,[A,B]]\big) + h^3(\cdots). \tag{6.4.4} \]
Each coefficient \(\alpha_k\) is a known polynomial in the splitting coefficients \(a_j, b_j\). Setting these equal to zero yields the order conditions. Optimised splittings (Section 6.3) choose them to minimise the leading error terms in the modified system.
If the original system is Hamiltonian with \[ A = X_T, \qquad B = X_V, \] then all nested commutators of \(A\) and \(B\) are Hamiltonian vector fields. Thus the modified vector field is Hamiltonian:
\[ \tilde{F}(h) = X_{\tilde{H}}, \]
where the modified Hamiltonian admits the expansion
\[ \tilde{H} = H + h^2 H_2 + h^4 H_4 + h^6 H_6 + \cdots. \tag{6.4.5} \]
For symmetric (time-reversible) methods, all odd-index terms vanish:
\[ H_{2k+1} = 0. \]
Thus Strang and all symmetric compositions preserve (to very high accuracy) an even-order approximation of the true Hamiltonian, which explains their exceptional long-time energy conservation.
For symplectic integrators, BEA yields the “shadow energy” theorem:
If \( \Phi_h \) is symplectic, then it exactly preserves the modified Hamiltonian \( \tilde{H} \) for exponentially long times \( T = \mathcal{O}(e^{c/h}) \), provided the solution remains in a compact set.
In practice:
This explains why symplectic integrators dominate long-time molecular dynamics, celestial mechanics, and semiclassical quantum simulations.
The leading term in the modified Hamiltonian for a symmetric splitting is
\[ H_2 = \alpha_{21}\{T,\{T,V\}\} + \alpha_{22}\{V,\{T,V\}\}. \tag{6.4.6} \]
The constants \( \alpha_{2j} \) encode the leading error of the scheme. Modern optimised splittings (Eqs.~(6.3.7) and others) choose coefficients to:
For example, the Omelyan–Mryglod–Folk (OMF) parameter \( \lambda \approx 0.1931833275 \) is chosen to minimise the modified energy error in lattice QCD simulations.
For complex splitting coefficients \(a_j, b_j \in \mathbb{C}\) with \( \Re(a_j), \Re(b_j) \ge 0 \), the modified vector field remains well-defined and admits the same BCH expansion (6.4.4).
If the original system is Hamiltonian but the complex step destroys exact symplecticity of each subflow (e.g. for real canonical variables), the method is usually interpreted in a complexified phase space. For Schrödinger-type PDEs, complex coefficients preserve unitarity exactly, and BEA yields a modified Hamiltonian which is anti-Hermitian:
\[ \tilde{H} = H + h^2 H_2 + h^4 H_4 + \cdots, \qquad \tilde{H}^\dagger = -\tilde{H}. \tag{6.4.7} \]
This ensures exact preservation of \( \|\psi(t)\|_{L^2} \) for highly accurate quantum time propagation.
Although BEA is asymptotically valid, it breaks down when:
High-order optimized splittings typically shift the dominant resonance tongues to higher order in \(h\), improving practical stability.
Key facts established in this section:
Splitting schemes extend naturally from finite-dimensional ODEs to evolution PDEs. Many PDEs can be written in the semilinear form
\[ u_t = A u + B(u), \tag{6.5.1} \]
where A is a (typically linear, stiff) operator generating a strongly continuous semigroup, and B is a (possibly nonlinear) locally Lipschitz operator. If the flows of \(A\) and \(B\) are explicitly computable or efficiently approximable, splitting schemes yield structure-preserving, high-accuracy time integrators.
This section develops the rigorous theory and SOTA implementations for splitting methods applied to:
The underlying functional-analytic framework uses semigroup theory, Fourier spectral discretization, and geometric properties of flows in infinite dimensions.
Let \(A: D(A)\subset X \to X\) be the generator of a \(C^0\)-semigroup \(e^{tA}\) on a Banach space \(X\). Let \(B: X \to X\) be locally Lipschitz so that the ODE
\[ u_t = B(u) \tag{6.5.2} \]
has a unique local flow \( \varphi_B^t\). Then the splitting schemes
are well-defined on \(X\).
Convergence follows from the classical theory: if \(A+B\) generates a (mild or strong) solution semigroup, then the Lie and Strang splittings converge with order 1 and 2 respectively under mild regularity assumptions.
\[ \|\Phi_h^n u_0 - e^{nh(A+B)}u_0\| \le C h^p, \qquad p = 1~\text{(Lie)},\ 2~\text{(Strang)}. \tag{6.5.3} \]
Many PDEs have linear parts diagonalizable by Fourier transform. Let \( \mathcal{F} \) denote the Fourier transform, then:
\[ A = \mathcal{F}^{-1} \Lambda(k)\, \mathcal{F}, \tag{6.5.4} \]
so that its flow is exactly computable:
\[ e^{tA} u = \mathcal{F}^{-1} \big( e^{t\Lambda(k)} \hat{u}(k) \big). \tag{6.5.5} \]
This leads to highly efficient splitting integrators:
This structure underlies spectral splittings for Schrödinger, NLS, KdV, and dispersive PDEs.
Consider the linear Schrödinger equation
\[ i\psi_t = -\frac{1}{2}\Delta \psi + V(x)\psi. \tag{6.5.6} \]
Define
Both flows are exactly computable:
\[ e^{tA}\psi = \mathcal{F}^{-1}\left( e^{-i t |k|^2/2} \hat{\psi}(k)\right), \tag{6.5.7} \]
\[ e^{tB}\psi = e^{-i t V(x)} \psi(x). \tag{6.5.8} \]
Thus Strang splitting gives a unitary, second-order, time-reversible method. High-order complex-coefficient splittings (Section 6.3) are widely used in quantum chemistry and semiclassical dynamics.
For
\[ i \psi_t = -\frac{1}{2}\Delta \psi + g|\psi|^2\psi, \tag{6.5.9} \]
the nonlinear flow is
\[ \varphi_B^t(\psi) = e^{-i g t |\psi(x)|^2} \psi(x), \tag{6.5.10} \]
which preserves mass exactly. Strang and higher-order splittings preserve:
These are state-of-the-art and standard in Gross–Pitaevskii solvers and Bose–Einstein condensate simulations.
For parabolic systems:
\[ u_t = \Delta u + B(u), \tag{6.5.11} \]
the semigroup \(e^{t\Delta}\) is contractive only for \(t\ge 0\). Therefore splitting coefficients must satisfy:
\[ \Re(a_j) \ge 0, \qquad \Re(b_j) \ge 0. \tag{6.5.12} \]
Real-coefficient high-order splittings require negative coefficients (Sheng–Suzuki barrier), which are unstable for diffusion.
Thus modern high-order solvers for parabolic PDEs use:
For KdV:
\[ u_t + 6 uu_x + u_{xxx} = 0, \tag{6.5.13} \]
split into:
The combination via Strang splitting yields a highly effective scheme used in soliton dynamics, dispersive shock wave, and Whitham modulation simulations.
For the Klein–Gordon equation:
\[ u_{tt} = \Delta u - m^2 u - f(u), \tag{6.5.14} \]
rewrite as a first-order Hamiltonian system in Fourier space:
\[ \frac{d}{dt} \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} v \\ \Delta u - m^2 u \end{pmatrix} + \begin{pmatrix} 0 \\ -f(u) \end{pmatrix}. \tag{6.5.15} \]
Linear part diagonalises in Fourier space; nonlinear part is pointwise. High-order splittings with FFTs are state-of-the-art for nonlinear wave propagation, spectral radiation studies, and weak turbulence.
Maxwell’s equations (in vacuum):
\[ \partial_t E = \nabla\times H, \qquad \partial_t H = -\nabla\times E. \tag{6.5.16} \]
preserve:
Splitting into:
yields symplectic, structure-preserving methods suitable for:
Key takeaways:
Splitting and composition methods, Magnus integrators, and exponential integrators all arise from the same operator-theoretic viewpoint: approximate the evolution operator of a differential equation by structured exponentials of simpler pieces. This section explains their common algebraic backbone and clarifies when each class is preferable.
We focus on:
Consider a linear (possibly time-dependent) evolution equation on a Banach space \(X\):
\[ u'(t) = A(t) u(t), \qquad u(t) \in X, \]
where \(A(t)\) is a (possibly unbounded) linear operator generating an evolution family \(U(t,s)\) such that
\[ u(t) = U(t,s) u(s), \qquad U(t,s) U(s,r) = U(t,r),\quad U(s,s)=I. \]
Magnus viewpoint. For sufficiently regular \(A(t)\), we may write
\[ U(t,s) = \exp\big(\Omega(t,s)\big), \]
where \(\Omega(t,s)\) is given by the Magnus expansion, an infinite series of time integrals of nested commutators of \(A(\tau)\). Truncating \(\Omega\) yields Magnus integrators.
Splitting viewpoint. When \(A\) decomposes into simpler pieces,
\[ A(t) = A_1(t) + A_2(t) + \cdots + A_m(t), \]
whose flows (or exponentials) can be computed efficiently, we approximate \(U(t,s)\) by compositions of exponentials of the sub-operators, e.g. Lie–Trotter or Strang splittings over a time step \(h\).
Exponential integrator viewpoint. For semilinear systems
\[ u'(t) = A u(t) + g(t, u(t)), \]
exponential integrators treat \(e^{hA}\) exactly (or nearly so) and approximate the contribution of the nonlinearity \(g\) via \(\varphi\)-functions and quadrature. These methods can be interpreted as Magnus-/splitting-type approximations in the variation-of-constants formula.
In all three cases, the core object is the exponential of an operator (or sum of operators) and its algebra under composition, governed by the Baker–Campbell–Hausdorff (BCH) and Magnus series.
For simplicity, consider an autonomous decomposition
\[ A = A_1 + A_2 \]
and a splitting step over a time increment \(h\):
\[ S(h) = e^{a_1 h A_1} e^{b_1 h A_2} e^{a_2 h A_1} e^{b_2 h A_2} \cdots e^{a_s h A_1} e^{b_s h A_2}. \tag{6.6.1} \]
The BCH formula ensures that there exists an operator \(Z(h)\) such that
\[ S(h) = \exp(Z(h)). \]
Expanding \(Z(h)\) as a formal power series in \(h\) gives
\[ Z(h) = h(A_1 + A_2) + h^2 \alpha_1 [A_1, A_2] + h^3 \big( \alpha_2 [A_1,[A_1,A_2]] + \alpha_3 [A_2,[A_1,A_2]] \big) + \cdots, \tag{6.6.2} \]
where the coefficients \(\alpha_j\) are polynomials in the composition coefficients \(\{a_k, b_k\}\). The order conditions for the splitting require that
Formally, \(Z(h)\) plays the same role as the Magnus generator \(\Omega(h)\) for the exact flow. For an autonomous problem, the exact Magnus generator is simply
\[ \Omega(h) = h(A_1 + A_2), \]
while for time-dependent problems \(A(t)\) the Magnus series contains time integrals of commutators. In this sense:
Splitting methods can be viewed as particular approximations of the Magnus generator by the BCH logarithm of a product of simpler exponentials.
The modern theory of commutator-free Magnus integrators makes this connection explicit by designing products of exponentials whose BCH logarithm reproduces a chosen truncated Magnus series.
For linear non-autonomous systems
\[ u'(t) = A(t) u(t), \]
the Magnus expansion over a step \([t_n,t_{n+1}] = [t_n, t_n+h]\) is
\[ \Omega_n(h) = \int_{t_n}^{t_{n+1}} A(\tau)\, d\tau - \frac{1}{2} \int_{t_n}^{t_{n+1}}\!\!\int_{t_n}^{\tau_1} [A(\tau_1), A(\tau_2)]\, d\tau_2\, d\tau_1 + \cdots. \tag{6.6.3} \]
Classical Magnus integrators compute or approximate commutators like \([A(\tau_1), A(\tau_2)]\), which may be expensive for large matrices or operators.
Commutator-free Magnus (CFM) methods avoid explicit commutators by using products of exponentials of linear combinations of evaluations of \(A\) at quadrature nodes:
\[ U_{n+1} = \prod_{j=1}^s \exp\!\left( h \sum_{\ell=1}^r \alpha_{j\ell} A(t_n + c_\ell h) \right) U_n. \tag{6.6.4} \]
Here, \(c_\ell\) are quadrature nodes (e.g. Gauss–Legendre abscissas), and the coefficients \(\alpha_{j\ell}\) are chosen so that the BCH logarithm of this product matches the truncated Magnus series up to order \(p\).
From the perspective of Chapter 4, each factor \(\exp(h \sum_\ell \alpha_{j\ell} A(t_n+c_\ell h))\) is simply another sub-flow whose exponential can be applied efficiently (e.g. via Krylov or FFT-based methods). The overall method is a high-order composition method for the time-dependent generator.
Key properties:
Thus CFM integrators can be seen as a SOTA synthesis of Magnus ideas (logarithm of the evolution) and composition/splitting ideas (products of exponentials of simpler pieces).
For semilinear problems
\[ u'(t) = A u(t) + g(t,u(t)), \tag{6.6.5} \]
the variation-of-constants formula gives
\[ u(t_{n+1}) = e^{hA} u(t_n) + \int_0^h e^{(h-\sigma)A} g(t_n + \sigma, u(t_n+\sigma))\, d\sigma. \tag{6.6.6} \]
Exponential Runge–Kutta and related exponential integrators approximate the integral term via quadrature and replace \(u(t_n+\sigma)\) by internal stage values. Typical schemes have the form
\[ U_{n+1} = e^{hA} U_n + h \sum_{j=1}^s \varphi_j(hA) G_j(U_n), \tag{6.6.7} \]
where the \(\varphi_j\) functions arise from integral representations involving \(e^{(h-\sigma)A}\).
Relation to splitting. A simple Lie–Trotter splitting for (6.6.5) reads
\[ U_{n+1} = \varphi_g^h\!\big(e^{hA}U_n\big), \]
where \(\varphi_g^h\) is the nonlinear flow generated by \(u' = g(t,u)\) (with appropriate time freezing or time-labelling). Higher-order Strang-type and composition splittings can be interpreted as particular quadrature approximations of (6.6.6), just as exponential integrators do, but organized as compositions of sub-flows rather than as direct quadrature of the integral.
Relation to Magnus. If \(g\) is linear in \(u\), then the full generator \(A(t)\) of (6.6.5) is linear, and Magnus/CFM schemes apply directly. For general nonlinear \(g\), Magnus ideas still influence the design of exponential integrators (e.g. time symmetry, commutator structure of linearized operators, etc.), but require more care.
In practice:
Modern SOTA codes often blend splitting, Magnus/CFM, and exponential ideas to balance structure preservation, efficiency, and robustness.
From the geometric standpoint, all these methods are variations on a single theme: construct a discrete flow as a product of exponentials whose BCH logarithm approximates the exact generator in a way that respects invariants, symmetries, and qualitative structure.