Time-symmetry, symplecticity and stability of Euler-Maclaurin and Lanczos-Dyche integration

TL;DR

This paper introduces the Lanczos-Dyche (LD) formula, a time-symmetric, two-point Taylor-based integration rule that yields high-order, multi-derivative approximations for evolving $du/dt=f(t,u)$. For linear (quadratic) Hamiltonian systems, LD can preserve both energy and symplectic structure, a property not available to general explicit methods, with LD reducing to Padé approximants for linear cases. When used in a method of lines framework for PDEs, LD-based schemes are unconditionally stable, enabling large time steps without CFL limits, while maintaining fidelity to continuum invariants; for nonlinear or dissipative systems, LD remains highly accurate and exhibits bounded, oscillatory deviations in conserved quantities, behaving as “symplectic on average.” The findings position LD as a powerful, general-purpose alternative to traditional integrators, especially for linear or quasi-linear problems and PDEs, with broad implications for long-time simulations in Hamiltonian and wave-like systems.

Abstract

Numerical evolution of time-dependent differential equations via explicit Runge-Kutta or Taylor methods typically fails to preserve symmetries of a system. It is known that there exists no numerical integration method that in general preserves both the energy and the symplectic structure of a Hamiltonian system. One is thus normally forced to make a choice. Nevertheless, a symmetric integration formula, obtained by Lanczos-Dyche via two-point Taylor expansion (or Hermite interpolation), is shown here to preserve both energy as well as symplectic structure for linear systems. This formula shares similarities with the Euler-Maclaurin formula, but is superconvergent rather than asymptotically convergent. For partial differential equations, the resulting evolution methods are unconditionally stable, i.e, not subject to a Courant-Friedrichs-Lewy limit. Although generally implicit, these methods become explicit for linear systems.

Figures (6)

• Figure 1: Illustrative comparison between remainder terms of the Lanczos-Dyche (LD), Euler-Maclaurin (EM), Taylor Expansion (T) formulas, using derivatives of the integrand of order up to $15^{\rm th}$. In this example, the test function $f(t) = e^{-t^2}$ is integrated in the interval $(t_1,t_2) = (-1/2,1/2)$ with a single, large time-step $\Delta t=t_2-t_1=1.$ The solid blue, dashed green, and dashed-dot lines show the absolute difference between the exact integral $I(f)=\int_{t_1}^{t_2} f(t) dt$ and its numerical approximation $I(\mathcal{P}_n)=\int_{t_1}^{t_2} \mathcal{P}_n (t) dt$ computed by integrating an interpolating polynomial, that is, via the LD, EM, and T formulas, given by Eqs. \ref{['eq:EM']}, \ref{['eq:LD']} and \ref{['eq:Taylorint']} respectively. This absolute difference is in agreement with the estimates of the respective remainder terms $R_n(f)$ for these formulas. As the value of $\tau$ in the remainder terms makes an insignificant impact, it is chosen rather arbitrarily to be the midpoint of the interval $\tau \in (t_{1},t_{2})$.
• Figure 2: A-stability regions for different time-stepping schemes.
• Figure 3: Phase portraits of various numeric methods used to evolve the simple harmonic oscillator. Note that the Euler and RK2 methods exhibit unphysical growth. The Verlet methods create closed trajectories, but energy is not conserved throughout a period. The LD methods are symplectic and energy conserving, so they form perfect circles, the expected trajectories. To exaggerate these effects a large step size, $\Delta t = 0.75$, was used.
• Figure 4: The relative error $\delta E /E$ in numerical energy compared with exact energy of a harmonic oscillator as a function of time. Observe that the error in the LD methods is of the order of machine precision, indicating that energy is numerically conserved. Contrast this with the RK methods, where the error grows polynomially in time, indicating that energy is not conserved.
• Figure 5: The relative error $\delta E/E$ in the numerical energy compared to the exact energy as a function of time for the pendulum equation. In all simulations $\Delta t = 0.1$. Note that for the LD methods the signal is locally a periodic function, strongly suggesting that the method is effectively symplectic, even if not exactly symplectic. Note that that the magnitude of the relative error is small. Contrast this with the RK methods which quickly display polynomial growth.