Operator Splitting

Pdes

Operator splitting is the practical answer to a common situation: you have a time-dependent PDE whose right-hand-side decomposes naturally into PIECES, each easy to integrate on its own, but the COMBINATION is awkward. The split-step strategy is to alternate between integrating each piece separately — advance with operator $A$ for a short time, then with operator $B$ , then back to $A$ , etc. — and let the alternation approximate the full coupled evolution. The leading error is proportional to the COMMUTATOR $[A, B]$ ; in the limit of small time step, splitting is exact when the operators commute and second-order accurate (with the right symmetric splitting) when they don't.

Ubiquitous in time-dependent PDEs: reaction-diffusion equations (split reaction from diffusion), advection-diffusion (split advection from diffusion), time-dependent Schrödinger (split kinetic from potential), Navier-Stokes via pressure projection (split convection-diffusion from pressure correction). The technique trades MULTIPLICATIVE complexity (each piece's solver vs the full coupled solver) for an ADDITIVE structure that decouples the difficulty.

The Trotter formula

The mathematical foundation is the Trotter product formula. For (suitably nice) operators $A$ and $B$ :

exp (T (A + B)) = n \to \infty lim (exp (T A / n) exp (TB / n))^{n} .

For finite $n$ , the approximation is $exp (T (A + B)) \approx (exp (Δ t A) exp (Δ tB))^{n}$ with $Δ t = T / n$ . This is the LIE-TROTTER splitting. The error per step is governed by the BAKER-CAMPBELL-HAUSDORFF formula:

exp (Δ t A) exp (Δ tB) = exp (Δ t (A + B) + \frac{Δ t ^{2}}{2} [A, B] + O (Δ t^{3})) .

So the per-step error is $O (Δ t^{2})$ in the EXPONENT, contributing $O (Δ t^{2})$ to the unitary one-step error. Over $n = T /Δ t$ steps that accumulates to $O (Δ t)$ — FIRST-ORDER convergence. The error coefficient is $[A, B] /2$ ; when $A$ and $B$ commute, the commutator vanishes and the splitting is exact at any time step.

The Strang splitting

Symmetrizing the splitting gives second-order accuracy at modest extra cost:

exp (Δ t (A + B)) \approx exp (Δ t A /2) exp (Δ tB) exp (Δ t A /2) .

Strang (1968) showed that this expression matches $exp (Δ t (A + B))$ through $O (Δ t^{2})$ : the leading $[A, B]$ error cancels by symmetry, and the next-order error is $O (Δ t^{3})$ per step, accumulating to $O (Δ t^{2})$ over the simulation — second-order CONVERGENCE.

Cost: ONE extra half-step of $A$ per outer time step compared to Lie-Trotter. But the two half-steps at consecutive Strang iterations can be COMBINED into a single full step except at the very first and very last — so the "amortised" cost of Strang is identical to Lie-Trotter (one A-step and one B-step per outer time step). For essentially no overhead, you get one extra order of accuracy. This is why Strang is the default splitting choice in practice.

Higher-order splittings exist (Yoshida fourth-order, Suzuki sixth-order) at the cost of more sub-steps with carefully tuned coefficients. They are useful for very-high-accuracy simulations of Hamiltonian systems (long-time molecular dynamics, planetary mechanics).

Code: the textbook test

The cleanest demonstration uses 2×2 non-commuting matrices, where the "exact" answer can be computed by scipy.linalg.expm of $(A + B) T$ and the split approximations by exponentials of $A$ and $B$ separately. The convergence orders fall out cleanly.

# Lie-Trotter and Strang splitting on a canonical non-commuting test:
#   v(t) = exp((A + B) t) v_0,
# with A and B non-commuting 2x2 matrices. The exact solution uses
# scipy.linalg.expm of the SUM; the split approximations use the
# matrix exponentials of A and B individually.
#
#   Lie-Trotter:  ( exp(A dt) exp(B dt) )^nsteps  → 1st order in dt
#   Strang:       ( exp(A dt/2) exp(B dt) exp(A dt/2) )^nsteps → 2nd order

import numpy as np
from scipy.linalg import expm

# Two non-commuting matrices (test problem)
A = np.array([[0.0,  1.0], [-1.0, 0.0]])     # rotation generator
B = np.array([[1.0,  0.0], [0.0, -1.0]])     # boost-like
v0 = np.array([1.0, 0.0])
T  = 2.0    # final time

def split_solve(T, dt, method='strang'):
    nsteps = int(T / dt)
    dt = T / nsteps
    if method == 'lie':
        step = expm(A * dt) @ expm(B * dt)
    else:
        step = expm(A * dt / 2) @ expm(B * dt) @ expm(A * dt / 2)
    v = v0.copy()
    for _ in range(nsteps):
        v = step @ v
    exact = expm((A + B) * T) @ v0
    return np.max(np.abs(v - exact))

print(f"Non-commuting 2x2 matrix splitting, T = {T}:")
print(f"{'method':>10s}  {'dt':>8s}  {'err':>14s}  {'order':>8s}")
for method in ['lie', 'strang']:
    errs = []
    for dt in [0.4, 0.2, 0.1, 0.05, 0.025]:
        e = split_solve(T, dt, method=method)
        errs.append(e)
        if len(errs) > 1:
            order = np.log2(errs[-2] / errs[-1])
            print(f"  {method:>8s}  {dt:>8.4f}  {e:>14.4e}  {order:>8.2f}")
        else:
            print(f"  {method:>8s}  {dt:>8.4f}  {e:>14.4e}  {'--':>8s}")
print("(Lie observed order ≈ 1, Strang observed order ≈ 2 — as theory says.)")

Output:

Non-commuting 2x2 matrix splitting, T = 2.0:
    method        dt             err     order
       lie    0.4000      8.0646e-01        --
       lie    0.2000      4.0525e-01      0.99
       lie    0.1000      2.0177e-01      1.01
       lie    0.0500      1.0050e-01      1.01
       lie    0.0250      5.0132e-02      1.00
    strang    0.4000      1.2056e-01        --
    strang    0.2000      3.0862e-02      1.97
    strang    0.1000      7.7621e-03      1.99
    strang    0.0500      1.9435e-03      2.00
    strang    0.0250      4.8605e-04      2.00
(Lie observed order ≈ 1, Strang observed order ≈ 2 — as theory says.)

Two things to read off. Lie-Trotter observed convergence order is 1.00 across all refinements — error halves when $Δ t$ halves. Strang observed order is 2.00 — error quarters when $Δ t$ halves. At $Δ t = 0.025$ , Strang has error $5 \times 1 0^{- 4}$ while Lie-Trotter has $5 \times 1 0^{- 2}$ — TWO orders of magnitude better for the same number of operator exponentials, just rearranged. This is the empirical content of "one extra order of accuracy for essentially no extra work."

The canonical PDE applications

Time-dependent Schrödinger: split-step Fourier

The Schrödinger equation $i \partial_{t} ψ = (\hat{T} + \hat{V}) ψ$ has the kinetic operator $\hat{T} = - \nabla^{2} /2$ (diagonal in MOMENTUM space) and potential $\hat{V} = V (x)$ (diagonal in POSITION space). Neither is sparse in the OTHER's natural representation, but each is trivial in its own. Strang splitting gives:

ψ (t + Δ t) \approx e^{- i \hat{V} Δ t /2} e^{- i \hat{T} Δ t} e^{- i \hat{V} Δ t /2} ψ (t) .

Each application: multiply by a complex phase in position space, FFT, multiply in momentum space, IFFT, multiply by a phase in position space. $O (N lo g N)$ per step, unconditionally stable, exactly unitary (norm-conserving). See the time-dependent Schrödinger page for the full implementation.

Reaction-diffusion equations

$u_{t} = D \nabla^{2} u + R (u)$ with diffusion operator $D \nabla^{2}$ and (possibly stiff, nonlinear) reaction $R (u)$ . Splitting decouples the two: a diffusion sub-step solved by an implicit solver (Crank-Nicolson or BDF) and a reaction sub-step that becomes an ODE at every grid point (often analytically integrable, or solved with a stiff ODE solver). Standard for chemical kinetics, biological pattern formation (Turing patterns), and combustion simulations.

Navier-Stokes via projection

Chorin's projection method (1968): split incompressible Navier-Stokes $\partial_{t} u + (u \cdot \nabla) u = - \nabla p + ν \nabla^{2} u$ with $\nabla \cdot u = 0$ into a CONVECTION-DIFFUSION sub-step (compute $u^{*}$ from current velocity) followed by a PRESSURE-PROJECTION sub-step (project onto the divergence-free subspace by solving a Poisson equation for $p$ ). The framework underlying most modern CFD codes.

When splitting fails

Two failure modes worth knowing.

Large step at order-1. Lie-Trotter with too-large $Δ t$ can lose monotonicity properties (e.g. the maximum principle for reaction-diffusion) that the original PDE possesses. Strang is more forgiving; some applications use even higher-order splittings to maintain positivity at large steps.
Stiff B with explicit A. If the B-step is stiff and you can only handle it implicitly, but A is fast and best done explicitly, splitting still works — but the step size is bounded by the EXPLICIT A's stability, not the unconditional B. Make sure the cheap operator drives the step size, not the expensive one.
Discontinuous coefficients. Splitting can produce spurious oscillations across discontinuities; if your A or B has a jump (interface problem), Strang's nice symmetry doesn't help and order reduction is common.

Why splitting is so widely used

The deep reason is that real physical problems are usually written as a SUM of physical effects — advection plus diffusion plus reaction, kinetic plus potential, transport plus collision. Each piece has a well-developed solver community, and the FULL coupled problem rarely does. Splitting lets you compose specialist solvers without re-inventing each one. As an organizing principle, it is to time-dependent PDEs what divide-and-conquer is to algorithms: SEPARATE the difficulty, then handle each piece with the best tool for it.

Method of lines — the alternative organizing principle: discretize space first, then integrate as one big ODE.
Time-dependent Schrödinger — the prototype application of Strang splitting in quantum dynamics.
Spectral methods — Fourier spectral representation pairs naturally with split-step methods for the Schrödinger equation.
Runge-Kutta methods — the natural choice for each individual sub-step when the operators don't have closed-form exponentials.