Fact-checked by Grok 2 weeks ago

Strang splitting

Strang splitting is a numerical integration technique used to solve time-dependent differential equations of the form \dot{x} = f(x) = f_1(x) + f_2(x), where the vector field f is decomposed into two subproblems whose exact flows can be computed efficiently. Introduced as a symmetric composition method, it approximates the exact solution operator by applying half-steps of the first flow, a full step of the second flow, and another half-step of the first flow:

x_{n+1} = \varphi_{h/2}^{{{grok:render&&&type=render_inline_citation&&&citation_id=1&&&citation_type=wikipedia}}} \circ \varphi_h^{{{grok:render&&&type=render_inline_citation&&&citation_id=2&&&citation_type=wikipedia}}} \circ \varphi_{h/2}^{{{grok:render&&&type=render_inline_citation&&&citation_id=1&&&citation_type=wikipedia}}}(x_n)

, achieving second-order accuracy O(h^2) while preserving key structural properties of the original system.^[1] The method was proposed independently by Gilbert Strang and Gury I. Marchuk in 1968 as a device for constructing accurate finite difference schemes, particularly for nonlinear partial differential equations in multiple dimensions, building on earlier operator splitting ideas like the Lie-Trotter product formula from 1959.^[1]^[2] In its mathematical foundation, Strang splitting relies on the Baker-Campbell-Hausdorff formula to analyze the local truncation error, which is given by \mathrm{e}^{h(F_1 + F_2)} - \mathrm{e}^{h/2 F_1} \mathrm{e}^{h F_2} \mathrm{e}^{h/2 F_1} = C h^3 + O(h^4), where C involves nested commutators of the operators F_1 and F_2. This symmetric arrangement ensures time-reversibility, making it a foundational example of a geometric integrator. Key advantages of Strang splitting include its computational efficiency for problems where subflows are solvable explicitly or via fast algorithms, such as in stiff systems or high-dimensional settings, and its ability to maintain long-term qualitative behavior like symplecticity in Hamiltonian dynamics. It is second-order accurate without requiring higher-order derivatives, outperforming first-order Lie splitting in precision while remaining simple to implement. Extensions to more than two operators or non-autonomous equations further enhance its versatility, though order reduction can occur in certain semilinear cases unless modified. Strang splitting finds wide applications in scientific computing, including the simulation of Hamiltonian systems like the pendulum or N-body problems, where it preserves energy over long times; quantum mechanics for the time-dependent Schrödinger equation via split-step Fourier methods; and reaction-diffusion partial differential equations such as the Allen-Cahn or Burgers' equations in chemistry and fluid dynamics. Its use in these fields has been pivotal for developing higher-order variants and adaptive schemes, influencing modern numerical libraries for partial differential equations.

Background on Operator Splitting

Lie-Trotter Product Formula

The Lie-Trotter product formula approximates the solution operator for the sum of two non-commuting operators A and B in the context of differential equations. For a time step t > 0, it states that

e^{t(A + B)} = \lim_{n \to \infty} \left( e^{t A / n} e^{t B / n} \right)^n,

where the finite-n approximation (e^{t A / n} e^{t B / n})^n achieves first-order accuracy in t, meaning the local error is O(t^2).^[3] This formula originates from Sophus Lie's 1875 work on continuous transformation groups, where it was established for finite-dimensional matrices in the study of Lie algebras and their exponentials.^[4] In 1959, H. F. Trotter extended the result to strongly continuous semigroups of bounded operators on Banach spaces, providing a rigorous convergence proof for the infinite product limit and enabling applications to infinite-dimensional evolution equations.^[5]^[3] The derivation relies on the Baker-Campbell-Hausdorff (BCH) formula, which expresses the logarithm of a product of exponentials as a series involving nested commutators: \log(e^X e^Y) = X + Y + \frac{1}{2}[X, Y] + \ higher\ order\ terms. Substituting X = tA/n and Y = tB/n yields \log(e^{tA/n} e^{tB/n}) = \frac{t}{n}(A + B) + \frac{t^2}{2n^2}[A, B] + O(1/n^3), so raising to the nth power gives

\left( e^{tA/n} e^{tB/n} \right)^n = \exp\left( t(A + B) + \frac{t^2}{2n}[A, B] + O\left(\frac{t^3}{n^2}\right) \right).

The leading error term \frac{t^2}{2n}[A, B] vanishes as n \to \infty, confirming convergence, while for fixed n=1, the O(t^2) error arises from the non-zero commutator [A, B] when A and B do not commute.^[3] As an illustrative example, consider the linear ordinary differential equation (ODE) \frac{dy}{dt} = (A + B)y on a vector space, where A and B are linear operators and y(0) is the initial condition. The exact solution is y(t) = e^{t(A + B)} y(0). The Lie-Trotter formula approximates this via the product (e^{tA/n} e^{tB/n})^n y(0), which alternates applications of the semigroups generated by A and B; for n=1, it simplifies to e^{tA} e^{tB} y(0), capturing the first-order dynamics but introducing an O(t^2) discrepancy due to non-commutativity.

First-Order Splitting Methods

First-order splitting methods, derived from the Lie-Trotter product formula, provide practical numerical approximations for solving systems of ordinary differential equations (ODEs) of the form \frac{dy}{dt} = L_1(y) + L_2(y), where L_1 and L_2 are operators that may not commute. These methods, also known as fractional step or operator splitting techniques, decompose the problem into sequential subproblems, each involving only one operator, to leverage computational efficiencies in cases where individual operators are easier to handle. The core approximation in the first-order Lie-Trotter splitting is given by

y_{n+1} = e^{\Delta t L_2} \left( e^{\Delta t L_1} y_n \right),

which advances the solution from y_n at time t_n to y_{n+1} at t_{n+1} = t_n + \Delta t. This formulation arises from the exponential map in the Lie-Trotter product, where the non-commutativity of L_1 and L_2 introduces a local truncation error of \mathcal{O}(\Delta t^2).^[6] Algorithmically, the method proceeds in two steps: first, solve the subproblem \frac{dy_1}{dt} = L_1(y_1) with initial condition y_1(t_n) = y_n over the interval [t_n, t_{n+1}] to obtain y_{n+1/2} = y_1(t_{n+1}); second, solve \frac{dy_2}{dt} = L_2(y_2) with initial condition y_2(t_n) = y_{n+1/2} over the same interval to yield y_{n+1} = y_2(t_{n+1}). In partial differential equations (PDEs), this corresponds to sequentially applying operators, such as solving the diffusion subproblem followed by the advection subproblem.^[7] These methods offer significant computational advantages when the sub-operators can be solved independently, as each step can employ tailored solvers— for instance, implicit schemes for stiff terms like diffusion while using explicit methods for non-stiff advection.^[8] This modularity reduces overall complexity and enables parallelization in multi-dimensional problems.^[9]^[6] However, the global accuracy remains first-order, resulting in an accumulated error of \mathcal{O}(\Delta t) over many time steps, which can limit their use in long-time integrations without refinement. A representative example is the one-dimensional advection-diffusion equation \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial x^2}, where a > 0 is the advection velocity and \nu > 0 is the diffusion coefficient.^[7] Here, the advection operator L_1(u) = -a \frac{\partial u}{\partial x} and diffusion operator L_2(u) = \nu \frac{\partial^2 u}{\partial x^2} are split sequentially: first, solve the advection equation \frac{\partial u_1}{\partial t} + a \frac{\partial u_1}{\partial x} = 0 explicitly (e.g., via upwind differencing), then solve the diffusion equation \frac{\partial u_2}{\partial t} = \nu \frac{\partial^2 u_2}{\partial x^2} implicitly using a Crank-Nicolson scheme, yielding the updated solution u_{n+1}.^[7] This approach demonstrates first-order convergence while efficiently handling the disparate physical scales.^[7]

Formulation of Strang Splitting

Symmetric Splitting for Two Operators

The Strang splitting method, proposed by Gilbert Strang in 1968 for the numerical solution of partial differential equations, provides a second-order accurate approximation for evolution equations of the form \frac{dy}{dt} = L_1(y) + L_2(y), where L_1 and L_2 are typically non-commuting operators. Over a time step \Delta t, the method advances the solution y_n to y_{n+1} through a symmetric composition of exponentials:

\hat{y} = e^{(\Delta t/2) L_1} y_n, \quad \hat{z} = e^{\Delta t L_2} \hat{y}, \quad y_{n+1} = e^{(\Delta t/2) L_1} \hat{z}.

This formulation improves upon first-order splitting by centering the full step of L_2 between two half-steps of L_1. The symmetric structure of Strang splitting ensures time-reversibility, a desirable property for preserving the qualitative behavior of reversible systems such as Hamiltonian dynamics. Applying the method with a negative time step -\Delta t to y_{n+1} recovers y_n exactly, due to the adjoint nature of the exponential operators and the balanced half-steps. This symmetry enhances long-term stability in simulations of oscillatory or periodic phenomena.^[10] In practice, implementing Strang splitting requires solving two half-step evolutions under L_1 and one full-step evolution under L_2 per time step, often leveraging analytical solutions or efficient numerical solvers for each subproblem when the operators are linear or separable.^[11] The computational cost is thus comparable to two first-order splittings, but with doubled accuracy. A representative example is the simple harmonic oscillator, modeled by the Hamiltonian H(q,p) = \frac{1}{2}(p^2 + q^2) with operators L_1(q,p) = p \frac{\partial}{\partial q} (kinetic) and L_2(q,p) = -q \frac{\partial}{\partial p} (potential). Applying Strang splitting over multiple periods with \Delta t = 0.1 yields a relative phase error of approximately O(\Delta t^3), significantly lower than the O(\Delta t^2) phase lag of first-order Lie-Trotter splitting, as demonstrated in phase-space trajectories that remain close to the exact circular orbits.

Extension to Multiple Operators

The Strang splitting method, originally formulated for two operators, can be generalized to problems governed by the evolution equation \dot{u} = (L_1 + \cdots + L_N) u with N > 2 operators by constructing a symmetric composition of the corresponding flow maps in a sequenced manner that preserves the second-order accuracy of the original scheme. One standard generalization alternates half time steps for the odd-indexed operators and full time steps for the even-indexed ones, or employs a cyclic permutation of the operators to ensure overall symmetry in the splitting sequence.^[12] A representative example for N=3 is the symmetric splitting

\exp(\Delta t (L_1 + L_2 + L_3)) \approx \exp\left( \frac{\Delta t}{2} L_1 \right) \exp\left( \frac{\Delta t}{2} L_2 \right) \exp(\Delta t L_3) \exp\left( \frac{\Delta t}{2} L_2 \right) \exp\left( \frac{\Delta t}{2} L_1 \right),

where the central operator L_3 receives a full step flanked by symmetric half-steps of the preceding operators. This formulation achieves second-order accuracy in \Delta t.^[13]^[12] As N increases, the computational cost scales linearly with the number of subproblems solved per time step, since each exponential \exp(\Delta t L_i) typically requires integrating a simplified equation. Moreover, preserving the requisite symmetry for second-order convergence becomes more intricate, and deviations in sequencing can degrade the order of accuracy if the operator commutators do not align favorably.^[13]^[12] In multi-dimensional partial differential equations, this extension is particularly valuable for decomposing multidimensional operators into unidirectional components, enabling efficient one-dimensional solvers for each axis. For the two-dimensional heat equation \partial_t u = \partial_{xx} u + \partial_{yy} u, the diffusion operator is split as L_x = \partial_{xx} and L_y = \partial_{yy}, with the Strang sequence advancing the solution via alternating one-dimensional heat equations along the x- and y-directions to capture cross-dimensional interactions accurately while reducing dimensionality in each substep. This dimensional splitting strategy extends seamlessly to higher dimensions by cycling through the axes.^[12]

Theoretical Properties

Order of Accuracy and Error Analysis

Strang splitting achieves second-order accuracy, meaning its local truncation error is of order O(\Delta t^3) and the global error over a fixed time interval T = n \Delta t is O(\Delta t^2), assuming stability conditions hold.^[14] In contrast to the first-order Lie-Trotter splitting, which incurs a local error of O(\Delta t^2) due to the uncanceled commutator term [L_1, L_2], the symmetric structure of Strang splitting eliminates this leading error. The local truncation error can be derived using Taylor expansion of the exact and approximate solutions. Consider the evolution equation \frac{du}{dt} = (L_1 + L_2) u, where the exact solution over a time step \Delta t is u(t + \Delta t) = e^{\Delta t (L_1 + L_2)} u(t). The Strang approximation is S(\Delta t) u(t) = e^{(\Delta t/2) L_1} e^{\Delta t L_2} e^{(\Delta t/2) L_1} u(t). Expanding both in Taylor series around t, the exact solution satisfies

u(t + \Delta t) = u + \Delta t (L_1 + L_2) u + \frac{(\Delta t)^2}{2} (L_1 + L_2)^2 u + \frac{(\Delta t)^3}{6} (L_1 + L_2)^3 u + O((\Delta t)^4).

The Strang operator S(\Delta t) matches this expansion up to the quadratic term, as

S(\Delta t) u = u + \Delta t (L_1 + L_2) u + \frac{(\Delta t)^2}{2} (L_1 + L_2)^2 u + O((\Delta t)^3),

with the cubic term differing by the leading error involving double commutators. Specifically, the local error is

\tau(\Delta t) u = \frac{(\Delta t)^3}{24} \left( [[L_1, L_2], L_2] - [[L_2, L_1], L_1] \right) u + O((\Delta t)^4),

where [L_1, L_2] = L_1 L_2 - L_2 L_1. This confirms the O(\Delta t^3) local truncation error.^[14] An elegant proof of the second-order accuracy uses the Baker-Campbell-Hausdorff (BCH) formula, which expresses the logarithm of the product of exponentials as a single exponential with correction terms involving nested commutators. For the Strang splitting,

\log \left( e^{(\Delta t/2) L_1} e^{\Delta t L_2} e^{(\Delta t/2) L_1} \right) = \Delta t (L_1 + L_2) + \frac{(\Delta t)^3}{24} \left( [L_1, [L_1, L_2]] + [L_2, [L_2, L_1]] \right) + O((\Delta t)^5).

The BCH expansion cancels the first-order commutator term \frac{1}{2} [L_1, L_2] that appears in the Lie-Trotter product, leaving the next term at O((\Delta t)^3). Alternative proofs rely on rooted tree analysis, which systematically enumerates the order conditions for splitting methods by associating terms with labeled trees representing nested Lie brackets, or on formal power series expansions of the operators. These approaches generalize the BCH result to nonlinear settings and higher-order schemes. For the global error, repeated application of the Strang splitting over n = T / \Delta t steps yields an accumulation of local errors bounded by O(\Delta t^2), provided the method is stable (e.g., the semigroup generated by L_1 + L_2 is bounded). This follows from standard error recursion arguments in numerical ODE theory. The second-order accuracy assumes exact solutions to the subproblems e^{\Delta t L_i} u; if approximate solvers are used with errors larger than O((\Delta t)^3), the overall order reduces accordingly.

Stability Conditions

The Strang splitting method is consistent provided that the numerical solvers for the individual subproblems are consistent, with the overall order of consistency aligning with the second-order nature of the splitting scheme itself.^[14] For linear problems, stability can be assessed via von Neumann analysis, which reveals the need for CFL-like conditions restricting the time step \Delta t based on the eigenvalues of the split operators to prevent amplification of high-frequency modes.^[15] Due to its symmetric structure, Strang splitting preserves the stability properties of underlying first-order splitting methods, often resulting in larger stability regions compared to non-symmetric alternatives.^[16] Under assumptions of local Lipschitz continuity for the operators and satisfaction of the stability conditions, the Strang splitting method converges with a global error of order O(\Delta t^2) over a finite time interval, as established by semigroup theory in Banach spaces.^[17] In special cases, such as semilinear parabolic partial differential equations solved with implicit subproblem integrators, Strang splitting exhibits unconditional stability, allowing arbitrary time steps without instability.^[18] However, in non-stiff regimes where operator commutators do not dominate, the inherent splitting errors can become relatively more prominent, potentially degrading the observed accuracy. In Hamiltonian systems, the method preserves symplecticity when applied to separable potentials, maintaining long-term energy behavior.

Applications

Hamiltonian and Quantum Systems

Strang splitting finds extensive application in quantum mechanics for approximating the time evolution of wave functions governed by the time-dependent Schrödinger equation, i \hbar \frac{\partial \psi}{\partial t} = \hat{H} \psi, where the Hamiltonian operator is decomposed as \hat{H} = \hat{T} + \hat{V}, with \hat{T} = -\frac{\hbar^2}{2m} \nabla^2 representing the kinetic energy and \hat{V}(\mathbf{r}) the potential energy. The Strang splitting scheme approximates the unitary evolution operator \exp(-i \hat{H} \Delta t / \hbar) over a time step \Delta t by symmetrically composing the propagators for the kinetic and potential parts: \exp(-i \hat{T} \Delta t / 2\hbar) \exp(-i \hat{V} \Delta t / \hbar) \exp(-i \hat{T} \Delta t / 2\hbar). This decomposition ensures a second-order accurate approximation that preserves unitarity up to the local truncation error, making it suitable for long-time simulations of quantum dynamics without introducing artificial dissipation. In classical Hamiltonian systems, defined by the canonical equations \dot{q} = \frac{\partial H}{\partial p} and \dot{p} = -\frac{\partial H}{\partial q} for a separable Hamiltonian H(q, p) = T(p) + V(q), Strang splitting generates a symplectic integrator. The method alternates exact flows of the kinetic and potential sub-Hamiltonians, preserving the symplectic structure of phase space and thus maintaining volume preservation over extended simulations. This property is crucial for avoiding artificial energy drift in conservative systems, such as molecular dynamics or celestial mechanics. A representative example is the propagation of Gaussian wave packets in one-dimensional quantum systems, such as a particle in a potential well. The kinetic evolution step corresponds to free propagation, efficiently computed via fast Fourier transform (FFT) in momentum space, while the potential step involves pointwise multiplication in position space. This reduces the full operator exponential to simple, non-stiff substeps, enabling accurate tracking of wave packet dynamics over many time steps. Compared to explicit Runge-Kutta methods, Strang splitting exhibits superior long-term energy conservation in Hamiltonian systems due to its geometric properties, with bounded energy errors over exponentially long times rather than linear growth. This advantage is particularly evident in simulations requiring fine time steps for accuracy, where symplectic integrators like Strang maintain near-constant total energy fluctuations. Split-operator methods based on Strang splitting have been integral to time-dependent density functional theory (TDDFT) since the 1980s, facilitating real-time simulations of electronic excitations in molecules and solids by propagating Kohn-Sham orbitals under time-dependent potentials.^[19]

Reaction-Diffusion and Fluid Dynamics

Strang splitting is applied to reaction-diffusion equations of the form \partial_t u = D \Delta u + R(u), where D > 0 is the diffusion coefficient, \Delta is the Laplacian, and R(u) represents nonlinear reaction terms. The method decomposes the evolution operator into diffusion and reaction substeps, advancing the solution over a time step \tau via half a diffusion step, a full reaction step, and another half diffusion step: u^{n+1} = e^{\tau D \Delta / 2} \circ e^{\tau \mathcal{L}_R} \circ e^{\tau D \Delta / 2} (u^n), where \mathcal{L}_R solves \partial_t u = R(u). This symmetric splitting achieves second-order accuracy in time while allowing the stiff diffusion term to be treated implicitly and the reaction explicitly or implicitly as needed.^[20] A representative example is the Allen-Cahn equation \partial_t u = \varepsilon^2 \Delta u - (u^3 - u), modeling phase separation in materials, where the splitting separates the linear diffusion \varepsilon^2 \Delta u from the nonlinear reaction -(u^3 - u). Numerical implementations confirm second-order convergence, with error bounds O(\tau^2) in appropriate Sobolev norms, and stability under variable time steps for initial data in H^k(\Omega). The approach is particularly beneficial for multi-scale problems involving stiff reactions, such as fast chemical kinetics versus slow diffusion, as it enables efficient treatment of disparate timescales without compromising global accuracy.^[20] In fluid dynamics, Strang splitting facilitates the solution of the Navier-Stokes equations by separating advection, viscous diffusion, and pressure projection into distinct operators. For incompressible flows, the velocity update often involves a half-advection step, full diffusion and pressure correction, and another half-advection step, ensuring the divergence-free condition via projection onto a solenoidal space. This decomposition, \mathbf{u}^{n+1} = \Phi_{\text{adv}}(\tau/2) \circ \Phi_{\text{diff-pressure}}(\tau) \circ \Phi_{\text{adv}}(\tau/2) (\mathbf{u}^n), yields second-order accuracy and allows specialized solvers for each subproblem, such as spectral methods for diffusion.^[21] An illustrative case is the viscous Burgers' equation \partial_t u + u \partial_x u = \nu \partial_{xx} u, a simplified model for fluid shock formation, where splitting isolates the nonlinear advection -u \partial_x u from diffusion \nu \partial_{xx} u, demonstrating convergence at O(\tau^2) for smooth initial data. In reactive turbulent flows governed by compressible Navier-Stokes with chemical source terms, the method splits advection-diffusion from stiff reactions, treating the latter implicitly to handle fast chemical timescales efficiently. Such applications benefit from managing multi-scale stiffness in combustion simulations, where reaction rates vastly exceed flow velocities.^[22]^[23] Numerical considerations arise from non-commuting operators, such as advection and reaction, where the symmetric ordering in Strang splitting mitigates but does not eliminate local truncation errors influenced by commutator terms like [ \mathbf{u} \cdot \nabla, R(u) ]. While the method maintains second-order global accuracy, operator sequencing can affect practical stability and dissipation in non-linear regimes, often requiring careful choice of implicitness for the stiffest components. For complex systems with multiple operators, extensions to Strang splitting briefly incorporate additional terms like pressure without altering the core symmetric structure.^[22]

Extensions

Higher-Order Splitting Schemes

Higher-order splitting schemes extend the second-order accuracy of Strang splitting by composing multiple instances of the basic Strang operator or related exponential maps, solving for coefficients that satisfy higher-order conditions derived from the Baker-Campbell-Hausdorff formula. These methods are particularly useful for preserving symplectic structure in Hamiltonian systems while achieving improved long-term accuracy. Seminal contributions include Yoshida's recursive composition technique, which builds integrators of arbitrary even order using symmetric compositions, and Suzuki's general framework for both symmetric and nonsymmetric decompositions of exponential operators.^[24] Yoshida's method constructs higher-order integrators by composing basic second-order Strang splittings, denoted as S_h = e^{A h/2} e^{B h} e^{A h/2}, where A and B are non-commuting operators. For fourth-order, Yoshida uses a triple composition

S^{{{grok:render&&&type=render_inline_citation&&&citation_id=4&&&citation_type=wikipedia}}}_h = S_{\gamma_1 h} \circ S_{\gamma_2 h} \circ S_{\gamma_1 h}

, with coefficients \gamma_1 = \gamma_3 = \frac{1}{2(2 - 2^{1/3})} \approx 0.6756, \gamma_2 = 1 - 2\gamma_1 \approx -0.3513, ensuring all error terms up to O(h^4) vanish.^[24]^[12] In general, higher-order schemes rely on Suzuki-Trotter compositions, where the exponential e^{(A+B)h} is approximated by products of the form \prod_{i=1}^s e^{a_i A h} e^{b_i B h}, with coefficients a_i, b_i solved from BCH expansion conditions to achieve order p. This "forest of exponentials" approach allows for orders beyond two by increasing the number of stages s, often using recursive constructions like Suzuki's triple-jump formula \mathcal{S}_{2k+2}(h) = \mathcal{S}_k(w_k h) \circ \mathcal{S}_k((1-2w_k)h) \circ \mathcal{S}_k(w_k h), where w_k are determined iteratively (e.g., w_1 = 1/(2-2^{1/3}) for order 4). These methods are versatile for multi-operator problems and maintain symplecticity when the base Strang splitting is symmetric.^[12] The computational cost of these schemes scales with the number of sub-steps, as each Strang composition requires three exponential evaluations, but adjacent flows of the same operator can be merged, reducing the total. For instance, a fourth-order scheme using the triple composition demands 7 evaluations per time step, compared to 3 for the base second-order Strang method, while sixth-order recursive constructions require up to 27 stages. This trade-off favors higher-order methods for long-time simulations where fewer large steps reduce global error accumulation, but limits practicality beyond order 4-6 due to the rapid growth in stages and potential ill-conditioning of coefficient equations.^[24]^[12] As an example, the fourth-order Yoshida scheme has been applied to separable Hamiltonian integration, such as the Kepler problem H = p^2/2 + (-1)/|q|, where it reduces phase errors by factors of $10^3 to $10^6 over long integrations compared to second-order Strang, enabling accurate orbital simulations with step sizes up to 0.5 (versus 0.1 for second-order) while preserving energy within $10^{-8}. Limitations include the exponential increase in sub-steps, which can make orders above 6 inefficient even on modern hardware, and the introduction of negative coefficients that may amplify oscillations in non-symplectic or dissipative systems. Optimal performance is typically achieved up to order 4-6, beyond which alternative techniques like optimized Runge-Kutta compositions are preferred.^[24]^[12]

Non-Uniform Time Stepping Variants

Non-uniform time stepping variants of Strang splitting extend the classical symmetric scheme to handle systems with varying temporal scales, stiffness, or time-dependent coefficients, allowing for adaptive adjustment of step sizes to optimize accuracy and efficiency while preserving second-order convergence where possible. These methods address limitations of fixed-step approaches in multi-scale problems, such as those arising in stiff evolutionary partial differential equations (PDEs), by incorporating error estimation and control mechanisms inspired by adaptive Runge-Kutta techniques.^[25] Adaptive Strang splitting adjusts the time step Δt dynamically based on local error estimates, typically by embedding a higher-order variant or companion scheme to predict and control truncation errors. For instance, in solving stiff multi-scale evolutionary PDEs, an adaptive time splitting technique monitors the splitting error through a predictor-corrector framework, where a preliminary step with a modified Strang composition estimates the local error, enabling step size reduction or enlargement to maintain a prescribed tolerance. This approach guarantees effective error control, achieving global second-order accuracy O(Δt²) with overhead comparable to fixed-step methods, as demonstrated in numerical experiments on reaction-diffusion systems.^[26]^[27] Multi-rate splitting variants apply different step sizes to distinct operators within the Strang framework, particularly beneficial for systems where components evolve at disparate rates, such as fast stiff reactions coupled with slower diffusion processes. In these schemes, the fast operator (e.g., reaction term) is integrated with finer sub-steps within each coarse step of the slow operator (e.g., diffusion), effectively resolving stiffness without uniformly refining the global time grid. A modified Strang splitting with variable step sizes for degenerate reaction-diffusion equations exemplifies this, where the Lie-Trotter substeps for the degenerate operator use adaptive refinement to avoid singularity issues, retaining favorable stability and second-order accuracy for the overall scheme.^[28] For time-dependent operators, variable coefficient Strang splitting symmetrizes the integration around the midpoint of the interval, adapting the classical form to account for evolving coefficients. In the context of the variable-coefficient Burgers' equation, ut = a(t)uxx + b(x,t)uux, the scheme splits into diffusion and nonlinear advection subproblems, applying half-steps of the diffusion operator at the interval endpoints and a full step of the advection at the midpoint, with coefficients evaluated at appropriate intermediate times to preserve symmetry and second-order accuracy. This adaptation ensures robustness for non-autonomous systems without introducing additional splitting errors beyond the standard commutator terms.^[29] An illustrative application appears in molecular dynamics simulations with varying potentials, where adaptive Strang splitting—equivalent to variable-step velocity Verlet—adjusts steps to mitigate resonance errors in oscillatory systems, such as diatomic molecules under time-varying forces. By estimating local errors from embedded higher-order compositions, the method dynamically varies Δt to balance energy conservation and computational cost, demonstrating improved long-term stability over fixed-step alternatives in benchmarks involving stiff potentials.^[30] Implementation of these variants often employs embedded schemes for error control, where a pair of Strang-like integrators (one second-order and one higher-order) computes the step size adjustment with minimal extra cost, ensuring the global O(Δt²) error remains bounded. This is particularly effective in operator splitting for stiff chemistry in combustion models, where dynamic adaptation of the Strang substeps for transport and reaction fractions optimizes efficiency without compromising accuracy.^[31]

References

[1]
On the Construction and Comparison of Difference Schemes
Strang-Type Difference Schemes for Multidimensional Problems · David Gottlieb. SIAM Journal on Numerical Analysis, Vol. 9, No. 4 | 14 July 2006. AbstractPDF ...
[2]
[PDF] 1 Trotter's product formula - Niels Benedikter
Jun 8, 2017 · Trotter's product formula and the BCH formula. Emil Frølich. June ... Now consider vectors lying in D(A) n D(B). We have that lims→0. Ss ...
[3]
[PDF] Numerical Solutions of a System of ODEs Based on Lie-Trotter and ...
Abstract In the present paper, we consider a first-order exponential splitting method (or exponential Lie-Trotter splitting) and second order exponential ...
[4]
Numerical Solution of Advection-Diffusion Equation Using Operator ...
In this study, effects of operator splitting methods to the solution of advection-diffusion equation are examined. Within the context of this work two ...
[5]
A time-reversible integrator for the time-dependent Schrödinger ...
Dec 16, 2019 · Recovery of time reversibility by a combination of the splitting method and Ehrenfest theorem. The Ehrenfest theorem52 states that the time ...
[6]
Balanced Splitting and Rebalanced Splitting - SIAM.org
In the situation described the Strang splitting procedure is not more accurate than Lie splitting. In this paper, we propose a modified Lie/Strang splitting ...
[7]
https://dergipark.org.tr/en/pub/ijeas/issue/32676/357237
[8]
https://www.sciencedirect.com/science/article/abs/pii/S0169809504001012
[9]
https://www.researchgate.net/publication/251552236_An_operator-splitting_finite_element_method_for_the_efficient_parallel_solution_of_multidimensional_population_balance_systems
[10]
[PDF] OPERATOR SPLITTING METHODS FOR DIFFERENTIAL EQUATIONS
May 13, 2010 · The classical splitting methods are the Lie-Trotter splitting, the Strang splitting. (Dimov et al., 2001), (Strang, 1968), (Faragó & Havasi, ...
[11]
https://epubs.siam.org/doi/10.1137/120878641
[12]
https://doi.org/10.1017/S0962492923000077
[13]
[PDF] Stability and convergence analysis of Strang splitting method ... - arXiv
Apr 10, 2025 · In. [32], Jahnke and Lubich derive the error bounds of using Strang splitting to solve the linear initial problem ∂tu = (A + B)u with A ...
[14]
Propagators for the Time-Dependent Kohn–Sham Equations
Apr 19, 2018 · An old-time favorite in condensed matter physics is the split-operator formula. (35) It belongs to the wide class of splitting techniques, ...3. Results · 3.4. Runge--Kutta Schemes · 3.4. 2. Exponential Rk...
[15]
Stability and Convergence of Strang Splitting Method for the Allen ...
Apr 10, 2025 · The Strang splitting method has been widely used to solve nonlinear reaction-diffusion equations, with most theoretical convergence analysis ...
[16]
Operator splitting for two-dimensional incompressible fluid equations
Feb 5, 2011 · We analyze splitting algorithms for a class of two-dimensional fluid equations, which includes the incompressible Navier-Stokes equations and ...
[17]
Operator splitting for partial differential equations with Burgers ...
Feb 21, 2011 · Particular examples include the viscous Burgers' equation, the ... We show that the Strang splitting method converges with the expected ...
[18]
[PDF] Efficient time stepping for reactive turbulent simulations with stiff ...
Dec 5, 2017 · The Strang splitting method is strongly stable and second-order accurate in time. By splitting the reaction operator from the advection-.
[19]
https://arxiv.org/abs/1803.02113
[20]
[1104.3697] Adaptive time splitting method for multi-scale ... - arXiv
Apr 19, 2011 · This paper introduces an adaptive time splitting technique for the solution of stiff evolutionary PDEs that guarantees an effective error control of the ...
[21]
https://arxiv.org/abs/1102.1094
[22]
[PDF] Adaptive time splitting method for multi-scale evolutionary ... - HAL
Apr 19, 2011 · This paper introduces an adaptive time splitting technique for the solution of stiff evo- lutionary PDEs that guarantees an effective error ...
[23]
Solving Degenerate Reaction-Diffusion Equations via Variable Step ...
In this paper, we propose a modified Lie/Strang splitting procedure that, while retaining all the favorable properties of the original method, does not suffer ...
[24]
[PDF] Strang Splitting for the Variable-coefficient Burgers' Equation
Jul 31, 2018 · The convergence of Strang splitting allows us to numerically solve the variable- coefficient Burgers' equation by solving its two constituent ...
[25]
A dynamic adaptive method for hybrid integration of stiff chemistry
In particular, the Strang splitting scheme [26] is among the most widely used operator splitting schemes for practical combustion simulations. It features ...
[26]
Dynamic adaptive chemistry with operator splitting schemes for ...
Strang-based splitting schemes are used to separate the governing equations into transport fractional substeps and chemical reaction fractional substeps. The ...<|control11|><|separator|>