Strang Splitting Method

• Strang Splitting Method is a numerical technique that decomposes evolution equations into subproblems for efficient and accurate integration.
• It achieves second-order accuracy through a symmetric composition that yields O(τ^3) local error and O(τ^2) global error.
• The method preserves key structures and invariants, making it effective for applications in PDEs, quantum dynamics, and stiff matrix problems.

The Strang splitting method is a fundamental operator-splitting scheme for numerically integrating evolution equations where the right-hand side is decomposable into two (or more) components, each of which may be efficiently and/or exactly solvable. Initially introduced for partial differential equations (PDEs) and later widely adopted in applied mathematics, physics, and computation, the Strang splitting method has become central in the numerical solution of reaction-diffusion systems, quantum dynamics, matrix ODEs, and other classes of evolution equations. Its particular strength is a flexible second-order accurate composition that exploits problem-specific structure, boundary conditions, and invariants, with substantial literature dedicated to its stability, error properties, and high-dimensional applicability.

1. Formal Definition and Algorithmic Formulation

Consider an evolution equation in either ODE or PDE form: $\frac{du}{dt} = A(u) + B(u), \quad u(0) = u_0,$ where $A$ and $B$ are operators (possibly differential, nonlinear, or matrix-valued).

The Strang splitting operator $S_\tau$ over a time step $\tau$ is defined as

$$S_\tau = e^{(\tau/2)A} \, e^{\tau B} \, e^{(\tau/2)A}$$

in the linear (or linearizable) case, or as the composition of exact (or high-accuracy) flows for the A- and B-subproblems. Given $u^n \approx u(n\tau)$, the method computes the next step as:

1. A-step: $u^* = \varphi_A(\tfrac{\tau}{2}, u^n)$
2. B-step: $u^{**} = \varphi_B(\tau, u^*)$
3. A-step: $u^{n+1} = \varphi_A(\tfrac{\tau}{2}, u^{**})$ where $\varphi_A$ and $\varphi_B$ are the flows generated by $A$ and $B$, respectively.

The symmetric composition produces a local error of order $O(\tau^3)$, yielding a global error of order $O(\tau^2)$ under mild regularity, commutator, and stability assumptions (Quan et al., 6 Aug 2025, Li et al., 2021).

2. Mathematical Properties: Consistency, Symmetry, and Convergence

The Strang splitting method is formally second-order accurate: for sufficiently regular solutions,

$u^n - u(n\tau) = O(\tau^2)$

uniformly as $\tau \to 0$, provided the compositions $e^{\tau A}$, $e^{\tau B}$ generate smooth flows. The method is symmetric (self-adjoint): applying the sequence in reverse order yields the adjoint method, which, combined with even-order composition, facilitates the construction of higher-order schemes (Einkemmer et al., 2013, Einkemmer et al., 2013).

Local error analysis via the Baker–Campbell–Hausdorff expansion yields

$$e^{(\tau/2)A} e^{\tau B} e^{(\tau/2)A} = e^{\tau(A+B) + (\tau^3/24)[B,[B,A]] - (\tau^3/48)[A,[A,B]] + O(\tau^5)}$$

which underpins the $O(\tau^3)$ local (one-step) defect (1711.02193). This sets the basis for global convergence analysis via discrete Grönwall inequalities and separation of truncation and stability estimates.

3. Stability, Energy Dissipation, and Invariance Preservation

Strang splitting schemes for gradient-flow systems frequently preserve qualitative geometric or analytic structures at the discrete level. For the matrix-valued Allen–Cahn equation, unconditional energy dissipation for the modified energy functional holds for all $\tau > 0$ without step-size restriction: $\widetilde E(U^{n+1}) \leq \widetilde E(U^n)$ (Quan et al., 6 Aug 2025). Similar strict dissipation results are available for scalar Allen–Cahn equations (polynomial/logarithmic), with unconditional monotonicity enforced in the modified discrete energy (Li et al., 2021). Maximum principles and mass bounds persist in bounded ODE propagators or through technical corrections for strong nonlinearities or mean constraints.

For matrix flows and quantum models, Strang schemes can be tailored to preserve Hermiticity, positivity, and norm/tracial invariants (e.g., $\operatorname{tr}\rho=1$ in the $N$-level Bloch equations) (Songolo et al., 2019). Boundedness of the determinant, e.g., $|\det U^n(x)|\le 1$ for matrix Allen–Cahn, can be rigorously enforced via singular-value analysis and algebraic inequalities (Quan et al., 6 Aug 2025).

4. Boundary Conditions, Constraints, and Order Reduction

While Strang splitting yields unconditional second-order convergence in the periodic or full-space setting, nontrivial boundary conditions typically induce order reduction unless properly handled. For Dirichlet data, unmodified splitting degrades to $O(\tau)$ accuracy, and for Neumann/Robin, to $O(\tau^{1.5})$, unless discrete corrections or compatibility-preserving strategies are adopted (1711.02193, Dang et al., 14 Apr 2025).

Several families of boundary correction methods are developed:

• Homogenization and boundary-compatible splitting via change of variables and problem reformulation (Nakano et al., 2019).
• Auxiliary correction terms (direct or elliptic), e.g., $q_n(x)$ so that $B q_n = B f(u)$ on the boundary (1711.02193, Einkemmer et al., 2016).
• Initial-boundary corrected splitting (IBC-Strang) that subtracts $z_n(t)$ so the split flows avoid commutator-induced leading errors entirely (Dang et al., 14 Apr 2025).
• Flow-based correctors or projection-based adjustment, constructing corrections via the flow of the nonlinearity; this achieves genuine $O(\tau^2)$ accuracy under general inhomogeneous boundary settings (Bertoli et al., 2019).

Further, constrained PDEs (e.g., diffusion–reaction systems with mass or moment conservation) encounter order reduction in naïve Strang schemes, resolvable by substep-wise corrections ensuring constraint consistency at transitions (Altmann et al., 2016).

5. Extensions: Nonlinear, Stiff, High-Dimensional, and Matrix Problems

Strang splitting naturally generalizes to a variety of advanced and high-dimensional settings:

• Matrix-valued and low-rank equations: For Sylvester-type stiff matrix differential equations, a combination of partial matrix exponentials (linear flow) and projection-based dynamical low-rank integrators maintains second-order convergence and robustness to stiffness (Scalone et al., 7 Feb 2026).
• Quantum dynamics and Schrödinger equations: Splitting-in-time with pseudo-spectral collocation (FFT or rank-1 lattice in space) enables efficient, high-accuracy simulation of Schrödinger and Bloch models up to high dimensions, with unitarity and mass/norm conservation preserved throughout (Suzuki et al., 2018, Songolo et al., 2019, Gu, 2023).
• Fractional and nonlocal operators: Strang schemes combined with circulant/skew-circulant splittings admit fast implementation via FFTs for problems with Toeplitz structure, unconditionally preserving maximum principles (Cai et al., 2022).
• Nonlinear and quasilinear equations: Modified or filtered Strang splittings relax spatial regularity requirements and adapt to loss of derivatives in quasilinear examples, with error rates characterized by the degree of regularity and frequency-separation methods (Wu, 2022, Lu et al., 2013).

6. High-Order, Adaptive, and Structure-Preserving Modifications

While classical Strang splitting is second-order, several approaches extend it to higher order or enhance structure preservation:

• Almost-symmetric and iterated Strang splittings: For cases where one partial flow cannot be computed exactly, iterated Picard or Newton schemes produce higher-order symmetric Strang-type integrators, which may then be composed (triple-jump and beyond) to attain arbitrary even order (Einkemmer et al., 2013, Einkemmer et al., 2013).
• High-order compact difference and ADI-based spatial discretizations: Temporal Strang splitting combines with fourth-order compact stencils and alternating direction–implicit (ADI) strategies for efficient, high-accuracy multiscale systems (Kuang et al., 2024).
• Bound-preserving and mass-conservative corrections: The introduction of Lagrange multipliers, cut-off, or constraint-enforcement mechanisms in the splitting framework yields strong invariance and physical constraint preservation without loss of formal accuracy (Kuang et al., 2024).

7. Practical Implementation and Numerical Evidence

Strang splitting is readily implemented in combination with Fourier (for periodic problems), cosine/sine spectral (for Neumann/Dirichlet), high-order finite difference, or finite element spatial discretizations. The flow maps for A and B may be realized via exponentials (matrix, operator, or FFT), Runge–Kutta schemes, or specialized solvers for each split component.

Numerical experiments across diverse settings (e.g., Allen–Cahn, Burgers, reaction–diffusion, Schrödinger, Bloch equations) robustly confirm

• second-order convergence in $L^2$, $H^1$, and energy metrics (when corrections are implemented as required);
• unconditional stability and invariance under large time steps;
• preservation of maximum principles, determinant/mass bounds, and monotonic energy dissipation for properly-constructed schemes (Quan et al., 6 Aug 2025, Li et al., 2021, Quan et al., 10 Apr 2025, Cai et al., 2022).

Corrected and extended Strang approaches consistently outperform unmodified versions, especially for boundary- and constraint-dominated regimes (Dang et al., 14 Apr 2025, 1711.02193, Bertoli et al., 2019). For high-order requirements or low-regularity data, tailor-fitted schemes and adaptive step sizes preserve both accuracy and computational efficiency (Einkemmer et al., 2013, Kuang et al., 2024).

---

Key References:

• Unconditional energy dissipation for matrix-valued Allen–Cahn: (Quan et al., 6 Aug 2025)
• Scalar Allen–Cahn stability and convergence: (Li et al., 2021)
• Boundary and constraint correction: (Dang et al., 14 Apr 2025, 1711.02193, Bertoli et al., 2019, Altmann et al., 2016)
• Low-rank splitting for stiff matrix ODEs: (Scalone et al., 7 Feb 2026)
• Structure-preserving and high-order compact methods: (Kuang et al., 2024)
• Strang splitting for fractional and high-dimensional PDEs: (Cai et al., 2022, Suzuki et al., 2018)
• Advanced symmetric and high-order compositions: (Einkemmer et al., 2013, Einkemmer et al., 2013)
• Nonlinear/constrained modifications and instability: (Wu, 2022, Lu et al., 2013)