Splitting and composition methods in the numerical integration of differential equations

1 Introduction by examples

• Splitting composes solutions of simpler component equations; Hamiltonian kinetic/potential splitting yields first-order symplectic Euler schemes and a second-order symplectic Strang or leapfrog scheme.The component equations should be simpler to integrate than the original system.Source: §1 · pp. 2, 4
• Splitting methods preserve geometric features such as symplecticity, volume, time symmetry, and first integrals, supporting qualitative fidelity and long-time integration across diverse applications.Applications include molecular dynamics, celestial mechanics, particle accelerators, quantum physics, and spatially discretized PDEs.Source: §1 · pp. 2, 3, 7, 9
• For the harmonic oscillator, symplectic Euler preserves area unlike standard Euler and exactly advances a slightly perturbed Hamiltonian; over 100 steps its long-time behavior is visibly superior.Both methods are first-order approximations to the exact trajectory.Source: §1 · p. 5
• In the Kepler example, symplectic Euler preserves angular momentum, but its perturbed dynamics cause precession and nonclosed trajectories because the Laplace–Runge–Lenz vector is not preserved.The perturbed problem no longer has the original SO(4) symmetry.Source: §1 · p. 6
• For the Kepler problem (e=1/5) over 500 periods, all methods use 1500 force evaluations per period; at equal cost, Störmer–Verlet is more accurate than Heun, while split symplectic methods’ average energy error does not grow and position error grows linearly.Source: §1 · p. 7, Figure 3 · p. 8
• Composing component flows that belong to the same diffeomorphism group yields an approximation in that group, inheriting the exact flow’s geometric properties.The same reasoning applies, with modifications, when the exact flow forms a semigroup or symmetric space.Source: §1 · p. 10
• Splittings may lack a general construction procedure, offer more efficient alternatives, and preserve different geometric properties that cannot always all be retained.The survey focuses on combining flows for a given splitting to obtain efficient higher-order approximations.Source: §1 · p. 10

2.1 Increasing the order of an integrator by composition

• Recursive composition builds higher-order integrators from a lower-order method, and symmetric second-order bases yield schemes of order 2k.The recursive schemes are time-symmetric; a symmetric S[2] produces S[2k] of order 2k.Source: §2.1 · pp. 10–11, §2.2 · p. 13
• A general composition-based construction can require many evaluations and usually has large truncation errors.Source: §2.1 · p. 11
• For a two-part split, the splitting method has order r if and only if its corresponding composition has order r for arbitrary consistent integrators.The equivalence applies to any split into two parts.Source: §2.1 · p. 12

2.2 Integrators and series of differential operators

• A method is time-symmetric exactly when its modified series contains only odd powers of the step size, so time-symmetric methods have even order.Source: §2.2 · p. 13
• The recursive symmetric schemes S[2k] have order 2k when the starting S[2] is a symmetric second-order integrator.Source: §2.2 · p. 13
• An integrator’s order can be checked by comparing its differential-operator series with the series associated with the exact flow.Order r requires agreement through the terms relevant to order r.Source: §2.2 · pp. 12–14

3.1 Order conditions via BCH formula

• BCH order conditions come from repeatedly expanding log(Ψ(h)) and setting its terms beyond the consistent first-order term to zero; the resulting polynomial expressions can be cumbersome.For splitting schemes, the expansion uses a Hall basis of the free Lie algebra, and its coefficients are polynomials in the splitting parameters.Source: §3.1 · pp. 15–16

3.2 A set of independent order conditions

• Lyndon multi-index polynomials give independent order conditions for composition methods: consistency requires the coefficients to sum to 1, and the conditions characterize order r for arbitrary ODEs and consistent integrators.The conditions are independent when arbitrary coefficient sequences are considered.Source: §3.2 · pp. 17–18

3.3 Order conditions of compositions methods with symmetry

• Time symmetry leaves only odd-order Lyndon conditions independent; imposing the additional odd-index restriction for symmetric compositions of symmetric schemes reduces the conditions further.The paper counts the remaining conditions using n_k for general compositions and m_k for Lyndon indices with all entries odd.Source: §3.3 · pp. 18–19

3.4 Relation among different sets of order conditions of composition methods

• The independent conditions from Lyndon multi-indices and labelled rooted trees are equivalent because both generate the same graded algebra of coefficient functions.BCH-derived coefficient polynomials belong to this algebra as well.Source: §3.4 · pp. 19–20

3.5 Negative time steps

• For splitting methods of order r ≥ 3, satisfying u3 = 0 necessarily requires at least one negative a_i and one negative b_i coefficient, so backward time steps are unavoidable.Such negative steps can be problematic for semigroup equations and PDEs with unbounded operators.Source: §3.5 · p. 21

3.6 Near-integrable systems

• For near-integrable systems, splitting coefficients can be chosen to prioritize eliminating low powers of the perturbation parameter ε; each order εh^k has only one corresponding error term to eliminate.The method’s consistency order in the limit ε → 0 is s1.Source: §3.6 · p. 22

3.7 Runge–Kutta–Nyström methods

• Splitting a second-order ODE into f[a](y,v) = (v,0) and f[b](y,v) = (0,g(y)) yields schemes that are particular Runge–Kutta–Nyström methods.The component flows can be computed exactly, enabling efficient implementation.Source: §3.7 · pp. 22–23

3.7 Runge–Kutta–Nystr¨om methods

• RKN splittings require 2 + l2 + ··· + lr independent order conditions, versus 2 + n2 + ··· + nr generally; the conditions coincide through order three.For r ≥ 4, the general characterization no longer applies to the RKN case, although results for negative time steps still apply because the conditions match through order three.Source: §3.7 · p. 23
• RKN methods can be more efficient than general splittings for second-order systems, and extend to systems whose split component is linear in velocity.They are useful when the reduced equation y′′ = My′ + g1(y) is easily solvable; for the Hamiltonian generalization, exact solutions are known only in particular cases, such as the Kepler problem or harmonic oscillator.Source: §3.7 · pp. 23–24
• For generalized harmonic oscillators, only one condition raises the order from 2k − 1 to 2k, and two raise it from 2k to 2k + 1.These special features can be used to design adapted splitting methods.Source: §3.7 · p. 24

4.1 Methods with modified potentials

• Scheme (72), constructed in, is of order four, while the RKN case (68), with f[a] and f[b] of form (69), has an h-flow for F[abb].That flow maps (y, v) to (y, v + h^3g[3](y)), where g[3](y) = 2g′(y)g(y).Source: §4.1 · p. 24
• Composing flows of T with modified-potential flows adds free parameters without much added computational cost, enabling more efficient integrators when the needed derivatives are affordable to evaluate.The modified-potential maps use g and derivatives such as g[3], g[5], g[7,1], and g[7,2].Source: §4.1 · p. 25
• Modified-potential schemes need not satisfy all order conditions at orders r ≥ 3, so positive-coefficient schemes can exist, although negative coefficients appear in methods of order six.The usual negative-time-step results do not apply in this case.Source: §4.1 · p. 26

4.2 Methods with processing

• Processing can reduce the conditions the kernel must satisfy: a method has effective order r when a post-processor makes its processed method attain conventional order r.The post-processor can satisfy many of the processed method's order conditions.Source: §4.2 · p. 26
• Processed schemes can use fewer kernel conditions, but are not recommended when intermediate results are needed at every step because their total compositions per step are typically higher.A post-processor approximation derived from the kernel's intermediate stages can be virtually cost-free and introduces only local error.Source: §4.2 · p. 27
• For symmetric kernels, Table 3 compares conventional and effective order-condition counts, including for RKN splitting integrators.The reported condition counts distinguish conventional order from effective order.Source: Table 3 · p. 28, §4.2 · p. 28

5 A collection of splitting methods

• The review classifies published splitting methods by family, order, and stage count, including processed methods and schemes for near-integrable systems.Its tables collect representative methods, with near-integrable schemes organized by their orders in step size and perturbation size.Source: §5 · pp. 28–30, Table 5 · p. 32, Table 7 · p. 33
• Symmetric compositions automatically satisfy even-order conditions when odd-order conditions hold, reducing the independent conditions and usually the stage count relative to non-symmetric counterparts.The review gives the minimum stage count as 2Mr − 1 for order r.Source: §5 · p. 29
• Finding all parameter solutions is challenging, and practitioners select among obtained solutions using an objective such as minimizing the sum of absolute composition coefficients.With extra stages, free parameters can generate solution families; the described continuation procedure does not guarantee finding all solutions.Source: §5 · pp. 29–30
• Methods optimized for two-part ODE splits are not necessarily the best choice for composing a method with its adjoint, despite closely related performance.Source: §5 · pp. 30–31
• For compositions of a method and its adjoint, schemes up to order six are reported to be more efficient in most problems than compositions of the compared class at the same order.The review reports finding no order-eight methods and leaves open whether the associated polynomial system admits more efficient solutions.Source: §5 · p. 31

Runge–Kutta–Nystr¨om methods.

• The 4th-order 6-stage RKN method has global error about 0.00175 times that of classical 4th-order Runge–Kutta and computational cost about 0.31 for a given error.On the Hénon–Heiles Hamiltonian, its truncation error is contrasted with that of a leapfrog-based composition, which is about 10 times larger than the classical Runge–Kutta scheme.Source: Runge–Kutta–Nystr¨om methods · p. 33
• Modified-potential RKN schemes can reduce stage counts by using extra coefficients to satisfy order conditions, when the added evaluation cost is moderate.The review notes this setting in classical and quantum mechanics.Source: Runge–Kutta–Nystr¨om methods · p. 33

• Backward error analysis treats the numerical solution as the exact solution of a modified differential equation, enabling long-time behavior to be studied through a perturbation of the original system.For the harmonic oscillator example, the numerical solution is the exact solution of a harmonic oscillator with modified frequency approximately 1.Source: §6 · pp. 35–36
• The modified equation preserves key structure: it has only even powers of h for symmetric methods, is divergence-free for volume-preserving methods on divergence-free systems, and is locally Hamiltonian for symplectic methods on Hamiltonian systems.Source: §6 · p. 36
• The formal modified-equation series generally does not converge outside the linear case, so rigorous results require coefficient bounds, optimal truncation, and estimates of numerical-to-modified-solution error.Source: §6 · p. 37
• A symplectic method of order r with constant step size h preserves the Hamiltonian to O(h^r) over exponentially long time intervals.Source: §6 · p. 37
• Changing the step size changes the modified equation, which the authors identify as a reason for poor long-time behavior with standard variable-step strategies for symplectic schemes.Source: §6 · p. 37

7.1 Splitting methods for linear systems

• For spatially discretized Schrödinger equations, splitting H = T + V uses a diagonal potential action and FFT-based kinetic action, avoiding direct exponentiation of a potentially costly full matrix.The FFT requires O(N log N) operations.Source: §7.1 · p. 38
• The linear-system schemes also suit discretized Maxwell equations in isotropic, lossless, source-free media because those systems have a similar structure.Source: §7.1 · p. 38
• For linear systems, a method family uses s exponentials to reach orders 4, 6, 8, 10, and 12, compared with 15 and 31 stages for general compositions at orders 8 and 10.Source: §7.1 · p. 39

7.2 Splitting methods for non-autonomous systems

• Converting a non-autonomous equation to an enlarged autonomous system can alter a useful frozen-time structure, forcing more general and less efficient integrators.The conversion can also be costly when time-dependent split functions must be evaluated at every stage.Source: §7.2 · p. 40
• An alternative composes time-frozen subflows using quadrature nodes; the reviewed order-four and order-six methods are generally more efficient than standard splitting on the enlarged system.With Gauss–Legendre nodes and k evaluations of each split component, order 2k can be built when enough stages are used.Source: §7.2 · pp. 40–41

8 Numerical examples with selected methods

• The examples compare splitting methods with standard methods and examine higher-order schemes and methods selected for each problem’s structure.For the Kepler tests, the authors measure average energy error against computational cost; for Schrödinger, they measure solution error against FFT calls.Source: §8 · p. 41, §8.1 · p. 42, §8.2 · p. 45

8.1 The perturbed Kepler problem

• The Kepler system’s separable and solvable structure supports several method families, with Kepler-specific methods expected to be superior as the perturbation tends to zero.Processing and modified potentials can further improve performance; evaluating the modified perturbation costs about 10–20% more than evaluating the unmodified part.Source: §8.1 · pp. 42, 45
• For the Kepler problem, the five-stage fourth-order method is more accurate than the three-stage method at every tested computational cost, and higher-order methods perform better when high accuracy is desired.The 24-stage eighth-order method is an exception to the reported higher-order pattern; the 21-stage methods perform better.Source: §8.1 · p. 43, Figure 4 · p. 43
• Symplectic integrators’ average energy error does not increase secularly, whereas the standard fourth-order Runge-Kutta method’s energy error grows linearly and position error quadratically.The authors note that energy-error results are largely independent of final time for this comparison.Source: §8.1 · p. 44, Figure 5 · p. 45

8.2 The Schr¨odinger equation

• For the Schrödinger test, Taylor methods outperform general symmetric-symmetric splitting and RKN methods, while structure-tailored GM1212 and P382 are clearly superior.Taylor methods’ relative performance deteriorates for longer integrations.Source: §8.2 · pp. 46–47, Figure 6 · p. 46

9 Conclusions and outlook

• Processing can reduce evaluations per time step, making it perhaps the most economical path to high order despite the rapidly growing number of order conditions.The paper presents order conditions based on Lyndon words and specializes them for near-integrable systems and second-order equations y'' = g(y).Source: §9 · p. 47
• The survey collects representative splitting integrators and classifies them by equation structure, stage count, and consistency order, with references to their original proposals.Source: §9 · p. 47
• Backward error analysis explains favorable long-time error propagation and preservation of invariants by relating splitting solutions to exact solutions of modified equations with the original system’s geometric properties.The paper also notes that this interpretation underlies rigorous estimates on the numerical solution.Source: §9 · pp. 47–48
• For the illustrated problem, specially adapted fourth-order splitting schemes are up to 6 orders of magnitude more accurate than Runge–Kutta at the same computational cost.They also outperform standard higher-order composition integrators over a wide range of step sizes h.Source: §9 · p. 48
• For the generalized harmonic oscillator, the system’s structure supports an optimized processed second-order method with many stages that is more efficient than other integrators.The oscillator arises from spatial discretizations of several PDEs in quantum mechanics, optics, and electrodynamics.Source: §9 · p. 48
• No general rule is offered for splitting an arbitrary vector field; some cases permit systematic splitting, while others require case-by-case choices that may preserve different geometric properties.The paper also leaves open which property should take priority when a physical problem conserves several simultaneously.Source: §9 · p. 48
• The survey treats constant-step integration, although adaptive steps may be mandatory for close encounters; changing h changes the modified equation, so geometric preservation is not guaranteed.Purpose-designed transformations can enable variable-step splitting by making the step constant in new variables.Source: §9 · pp. 48–49
• Higher-than-second-order splitting requires negative coefficients, while large stability intervals and efficiency beyond leading-error minimization remain design challenges.Negative-time progression may be unacceptable for ill-posed PDEs; general high-order schemes with large linear stability intervals remain an open problem, and later error terms can determine efficiency.Source: §9 · pp. 49–50