Systematic Improvement of Splitting Methods for the Hamilton Equations

We show how the standard (St{"o}rmer-Verlet) splitting method for differential equations of Hamiltonian mechanics (with accuracy of order $\tau^2$ for a timestep of length $\tau$) can be improved in a systematic manner without using the composition method. We give the explicit expressions which increase the accuracy to order $\tau^8$, and demonstrate that the method work on a simple anharmonic oscillator.

💡 Research Summary

The paper addresses a long‑standing limitation of the standard Störmer‑Verlet (also known as leap‑frog or symplectic splitting) integrator for Hamiltonian systems. While the Verlet scheme is second‑order accurate (error ∼ τ²) and preserves the symplectic structure, many applications—such as long‑term celestial mechanics, molecular dynamics, and wave propagation—require much higher accuracy without sacrificing the geometric properties. The usual route to higher order is the composition method (e.g., Suzuki‑Yoshida or triple‑jump schemes), which builds a high‑order map by concatenating several low‑order steps with carefully chosen coefficients. Although effective, composition quickly becomes cumbersome: the number of sub‑steps grows exponentially with the desired order, and the resulting coefficients can be large, leading to increased round‑off error and computational cost.

The authors propose a fundamentally different strategy: instead of composing many copies of the basic Verlet step, they augment the single “half‑kick – drift – half‑kick’’ sequence with explicit high‑order correction operators that cancel the leading error terms identified by backward error analysis. The key steps are:

Modified Hamiltonian Analysis – By applying the Baker‑Campbell‑Hausdorff (BCH) expansion to the Verlet map, they derive the modified Hamiltonian (\tilde H = H + τ²H₂ + τ⁴H₄ + τ⁶H₆ + …). The extra terms (H₂, H₄, H₆) are expressed in terms of nested Poisson brackets of the kinetic and potential parts, (T(p)) and (V(q)). These terms are responsible for the τ²‑order global error and for the slow drift of conserved quantities.
Construction of Correction Operators – They introduce a correction exponential (\exp(τ³C₃ + τ⁵C₅ + τ⁷C₇)) that is inserted between the two half‑kicks. The operators (C₃, C₅, C₇) are chosen such that the BCH expansion of the full map (kick‑correction‑drift‑kick) eliminates the unwanted terms up to τ⁸. Explicitly, \