Strang Splitting: Theory & Applications

Updated 4 January 2026

Strang splitting is a second-order time integration method that employs symmetric operator splitting to solve nonlinear and stiff evolution equations with improved efficiency and stability.
The method decomposes operators so that sub-flows are integrated individually, preserving time symmetry and achieving global error of O(Δt²) for complex PDEs.
Advanced implementations leverage FFT acceleration and boundary correction techniques to enhance performance in multidimensional and fractional models while maintaining discrete maximum principles.

Strang splitting is a widely used second-order time integration method for systems of evolution equations with multiple nonlinear or stiff components. It provides a systematic operator-splitting framework that is symmetric and achieves second-order accuracy in time under broad conditions. The method has had extensive impact in computational mathematics, especially for parabolic, hyperbolic, and dispersive PDEs, and has been the subject of significant theoretical and algorithmic development across diverse application domains, including fractional phase field models, reaction-diffusion equations, quantum dynamics, and more.

1. Fundamental Principles and Classical Construction

Strang splitting was introduced for the numerical solution of semilinear evolution equations of the form

$u'(t) = (A + B)u(t), \qquad u(0) = u^0$

where $A$ and $B$ represent operators (possibly nonlinear) that can each be integrated separately, often more efficiently or accurately than their sum.

The classical (second-order, symmetric) Strang step of size $\Delta t$ from $t_n$ to $t_{n+1} = t_n + \Delta t$ is

$u^{n+1} = e^{(\Delta t/2) B} \, e^{\Delta t A} \, e^{(\Delta t/2) B} u^{n}.$

This scheme is time-reversible and achieves global error $O(\Delta t^2)$ provided both $A$ and $B$ generate sufficiently regular flows. For nonlinear split systems, the scheme can be interpreted using the flows of the individual vector fields, i.e., via compositions of the analytic or numerically exact evolution maps for each operator.

Key properties:

Time symmetry: The method is self-adjoint; applying the steps backward undoes the advance.
Second-order accuracy: The global error is $O(\Delta t^2)$ as shown by Baker–Campbell–Hausdorff expansions and commutator estimates.
Extensible to multiple operators: Through recursive application or composition.

2. Algorithmic Variants and Two-Level Strang Splitting

In many applications, spatial discretization of PDEs yields effective ODEs with matrix or operator structure, motivating fast or structure-preserving extensions of Strang splitting.

A prominent example is the two-level Strang splitting for multidimensional spatial fractional Allen–Cahn equations (Cai et al., 2022). After semi-discretization of the fractional Riesz derivative, the system becomes

$u'(t) = A u(t) + f(u(t)),$

where $A$ is a symmetric Toeplitz matrix and $f(u) = u - u^3$ is the Allen–Cahn nonlinearity. Direct computation with Toeplitz $A$ is expensive ( $O(M^2)$ ), but $A$ can be decomposed:

$A = C + S$ , with $C$ a circulant and $S$ a skew-circulant matrix.

The corresponding two-level Strang splitting involves:

Nonlinear (reaction) half-step: $u_1 = Q(\Delta t/2, u^n)$ —exact for the ODE $v' = v - v^3$ .
Linear (diffusion) inner Strang splitting on $C + S$ : $e^{(\Delta t/2)C} \, e^{\Delta t S} \, e^{(\Delta t/2)C}$ , each efficiently computed using the FFT, reducing cost to $O(M \log M)$ .
Another nonlinear half-step.

Pseudocode:

def TwoLevelStrang(u, Δt):
    u1 = Q(Δt/2, u)
    y = FFT(u1)
    y = exp((Δt/2) Λ_C) * y
    u2 = IFFT(y)
    z = Ψ * u2
    z_hat = FFT(z)
    z_hat = exp(Δt Λ_S) * z_hat
    z = IFFT(z_hat)
    u3 = Ψ^* z
    w = FFT(u3)
    w = exp((Δt/2) Λ_C) * w
    u4 = IFFT(w)
    u_next = Q(Δt/2, u4)
    return u_next

This approach is provably second-order accurate in both time and space, preserves a discrete maximum principle for

\|u^n\|_\infty \le 1

, and enables substantial acceleration and stability even in high dimensions (Cai et al., 2022).

3. Convergence, Stability, and Structure-Preserving Properties

Extensive theoretical work underpins the stability and accuracy of Strang splitting:

Convergence: Provided the full problem and each split flow are sufficiently regular, the global error is $O(\Delta t^2)$ , established via local error expansions and commutator estimates.
Discrete Maximum Principle: For Allen–Cahn-type equations, under suitable matrix properties, each sub-flow is a contraction in $\ell^\infty$ , and the Strang composition preserves maximum principles for all steps (Cai et al., 2022).
Energy Dissipation and Invariants: For gradient flow models (Allen–Cahn, matrix-valued analogs), a modified energy functional can be constructed under which the Strang method is (unconditionally) dissipative, guaranteeing stability without time-step restriction (Quan et al., 6 Aug 2025, Li et al., 2021). Determinant bounds and $H^1$ -stability can also be enforced in matrix-valued settings.

Furthermore, the principle generalizes to variable step sizes and to cases with inhomogeneous or mixed boundary conditions, with dedicated boundary corrections ensuring the preservation of formal order (Dang et al., 14 Apr 2025, 1711.02193, Einkemmer et al., 2016).

4. Order Reduction and Boundary Correction Techniques

A crucial practical challenge for Strang splitting is order reduction when applied to PDEs with inhomogeneous, Neumann, or mixed boundary conditions. While the method is formally second order for problems with periodic or homogeneous boundary data, in the presence of Dirichlet, Neumann, or Robin boundaries, classical splitting may degenerate to first or $1.5$th order global accuracy (Dang et al., 14 Apr 2025, 1711.02193, Nakano et al., 2019, Einkemmer et al., 2016).

Mitigating strategies include:

Direct compatibility-enforcing corrections: Adding a correction term $q(x)$ so that the sub-flows' boundary behavior matches that of the full problem.
Initial-boundary correction: Reformulating the problem with homogeneous boundary conditions at each step by subtracting a time-dependent or Taylor-expanded lift, so that splitting applies to boundary-homogeneous variables and then reconstructing the original solution (Dang et al., 14 Apr 2025).
Time-dependent boundary correction: Modifying temporal boundary prescriptions in the linear sub-flows to match Taylor expansions up to the correct order (Einkemmer et al., 2016).
Elliptic and multigrid corrections: Numerically solving for smooth compatibility functions via elliptic or multigrid methods at each step (1711.02193).

These techniques restore $O(\Delta t^2)$ global convergence for general boundaries, validated theoretically and empirically (Dang et al., 14 Apr 2025, 1711.02193, Nakano et al., 2019, Einkemmer et al., 2016).

5. Computational Efficiency and Advanced Implementations

Strang splitting enables specialized and highly efficient computational implementations:

FFT Acceleration: When split operators admit diagonalization in Fourier (circulant/skew-circulant) bases, as in the two-level Strang for Toeplitz-structured problems, each sub-flow step reduces to FFT-based multiplication, with $O(M \log M)$ complexity (Cai et al., 2022).
Lattice and Spectral Schemes: For high-dimensional Schrödinger equations, combining Strang splitting with lattice-based pseudo-spectral discretization achieves dimension-independent error constants and efficiency, conserving discrete norm and energy (Suzuki et al., 2018).
Bound/Mass Preservation: For nonlinear convection–diffusion, coupling Strang splitting with high-order compact difference spatial schemes and Karush–Kuhn–Tucker corrections enforces discrete conservation of bounds and mass without weakening temporal or spatial accuracy (Kuang et al., 2024).
Adaptivity and Non-uniform Steps: Adaptive or variable-step Strang splitting retains stability and accuracy for Allen–Cahn and related equations with Neumann, periodic, or more general nonlinear boundary conditions (Quan et al., 10 Apr 2025, Li et al., 2021).

6. Domain-Specific Applications and Generalizations

Strang splitting is pivotal across computational disciplines:

Phase-Field and Fractional Models: Widely used for Allen–Cahn, Cahn–Hilliard, and associated nonlocal or fractional PDEs to efficiently resolve complex interface dynamics under strict theoretical control (Cai et al., 2022, Quan et al., 6 Aug 2025, Li et al., 2021).
Quantum and Spectroscopic Dynamics: Central in time-dependent Schrödinger equation solvers (split-operator methods), conical intersection modeling, and high-dimensional pseudo-spectral quantum simulations (Suzuki et al., 2018, Gu, 2023).
Reaction–Hydrodynamics, Astrophysics: Used to couple fast reactions and fluid dynamics with guaranteed second-order accuracy only if energy/temperature is advanced together with species (Zingale et al., 2021).
General ODE/PDE Solvers: Foundation for high-order and composition methods, including almost-symmetric, iterated, and compositionally-boosted schemes that attain arbitrary even temporal orders (Einkemmer et al., 2013, Einkemmer et al., 2013).

7. Limitations, Extensions, and Theoretical Insights

While Strang splitting provides a robust and extensible framework, its performance is subject to:

Order reduction under certain boundary conditions, ameliorated by boundary correction schemes.
Sensitivity to inexact flows: Loss of symmetry and limitations in constructing high-order compositions if sub-flows are not solved exactly; the iterated almost-symmetric Strang splitting recovers high symmetry and enables composition to arbitrary order (Einkemmer et al., 2013, Einkemmer et al., 2013).
Optimality in Stability: Strang splitting and its repeated compositions (e.g., Störmer–Verlet, RKR/KRK integrators) maximize the stability region for certain second-order systems compared to other symmetric or palindromic splittings; analytical results establish Chebyshev-optimality for these cases (Casas et al., 2022).

The method continues to evolve, with research directions including structure preservation under more general operator algebras, treatment of stiff/nonlocal couplings, and automated adaptivity for large-scale scientific simulations.

References:

"A fast two-level Strang splitting method for multi-dimensional spatial fractional Allen-Cahn equations with discrete maximum principle" (Cai et al., 2022)
"Practical Effects of Integrating Temperature with Strang Split Reactions" (Zingale et al., 2021)
"An initial-boundary corrected splitting method for diffusion-reaction problems" (Dang et al., 14 Apr 2025)
"Unconditional energy dissipation of Strang splitting for the matrix-valued Allen-Cahn equation" (Quan et al., 6 Aug 2025)
"Strang splitting in combination with rank-$1$ and rank- $r$ lattices for the time-dependent Schrödinger equation" (Suzuki et al., 2018)
"Efficient boundary corrected Strang splitting" (1711.02193)
"Strang splitting structure-preserving high-order compact difference schemes for nonlinear convection diffusion equations" (Kuang et al., 2024)
"Stability and Convergence of Strang Splitting Method for the Allen-Cahn Equation with Homogeneous Neumann Boundary Condition" (Quan et al., 10 Apr 2025)
"Stability and convergence of Strang splitting. Part I: Scalar Allen-Cahn equation" (Li et al., 2021)
"A New Optimality Property of Strang's Splitting" (Casas et al., 2022)
"An almost symmetric Strang splitting scheme for the construction of high order composition methods" (Einkemmer et al., 2013)
"An almost symmetric Strang splitting scheme for nonlinear evolution equations" (Einkemmer et al., 2013)