Numerical methods for Chaotic ODE

TL;DR

The paper investigates backward error analysis for chaotic ODEs through residual analysis, the method of modified equations, and shadowing, arguing that numerical solutions effectively solve a persistently disturbed system $\dot{Y}=F(Y)+\varepsilon v(t)$ and that many attractor statistics remain robust despite exponential trajectory divergence. It demonstrates with the Lorenz and Henon–Heiles models how residuals and Hamiltonian structure influence backward error, and discusses the role and limits of shadowing and structure-preserving integrators in preserving correct statistics. A key message is practical: use variable steps and tolerance checks, monitor residuals and conserved quantities to detect spurious behavior, and recognize that exact trajectory convergence cannot be guaranteed in chaotic regimes, though statistics can still be trustworthy. The work also highlights open questions about the typicality of shadowing orbits in long-time simulations and the interplay between numerical perturbations and physical invariants. These insights have practical impact for scientists performing long-time chaotic simulations, guiding the choice of integrators and diagnostic checks to assess reliability of computed statistics.

Abstract

This paper explores backward error analysis for numerical solutions of ordinary differential equations, particularly focusing on chaotic systems. Three approaches are examined: residual assessment, the method of modified equations, and shadowing. We investigate how these methods explain the success of numerical simulations in capturing the behavior of chaotic systems, even when facing issues like spurious chaos introduced by numerical methods or suppression of chaos by numerical methods. Finally, we point out an open problem, namely to explain why the statistics of long orbits are usually correct, even though we do not have a theoretical guarantee why this should be so.

Figures (4)

• Figure 1: The computed solution to the Lorenz system on the interval $0 \le t \le 50$ with parameters $\sigma=10$, $\rho = 28$, and $\beta=8/3$, as computed by the Julia solver "Tsit5" with tolerances $10^{-8}$. With tolerances $10^{-9}$ the figure looks no different, to the casual eye.
• Figure 2: 2025 samples of the residual in the computed solution to the Lorenz system on the interval $0 \le t \le 50$ with parameters $\sigma=10$, $\rho = 28$, and $\beta=8/3$, as computed by the Julia solver "Tsit5" with tolerances $10^{-10}$.
• Figure 3: Solving the Henon--Heiles equations with the given initial conditions using the Størmer--Verlet method with $h=1.175$ (top left) and $h=1.18$ (top right) for $N=16,000$ steps. Both have bounded state and nearly-constant energy, though resonance effects make the plots look quite different. Intermittent spurious chaos ensues with $h$ slightly larger, say by $h=79/64 = 1.234375$; by $h=81/64=1.265625$ stability has apparently returned. In the graph, $q_2$ is plotted against $q_1$. For slightly larger $h$, all structure vanishes. The energy of the solution remains reasonably constant over this time interval for three out of the four graphs. The energy starts varying noticeably when the trajectory is chaotic, so the spurious chaos is detectable.
• Figure 4: All the possible orbits of $G(x)$ in quarter-precision IEEE floating-point.