A Simulation Can Conserve Energy Perfectly and Still Be Wrong

A double pendulum run for 10^6 steps can hold its energy to many digits and still give a wrong trajectory after about 5 s of simulated time. That is my estimate for a plausible case, and I derive it below. The thesis: small energy drift from a symplectic integrator is a necessary check on a chaotic trajectory, not a sufficient one.

I hold the energy claim on the other side too. A second-order symplectic method at step h should keep energy error below RK4 at step h/4 over 10^6 steps. I put that at 0.8 confidence, and the Convergence Lab run will test it. This essay says why winning that contest proves less than people think. I state the claim before the run, so the run can embarrass me.

What the energy check certifies

Stormer-Verlet is symplectic. Backward error analysis (Hairer, Lubich and Wanner [1][2]) shows that the numerical map is, up to a tiny remainder, the exact flow of a modified Hamiltonian $\tilde H = H + h^2 H_2 + h^4 H_4 + \dots$ Here h is the step size in seconds and H is the true energy in joules. The modified energy is conserved by the numerical flow almost exactly. The true energy then differs from it by $O(h^2)$ and does not drift. The same paper lists the results that follow: long-time energy conservation, linear error growth for integrable problems, and preserved invariant tori in near-integrable ones [1].

Note what this says. The numerical trajectory is an accurate trajectory of the wrong Hamiltonian $\tilde H$. The energy check confirms that $\tilde H$ is close to H. It says nothing about whether the numerical orbit stays close to the true orbit with the same initial data.

Where this fails: the theory assumes the numerical solution stays in a bounded region where the series for $\tilde H$ behaves. It also assumes h is small compared with the fastest time scale. For a pendulum flipping through large angles, the fast scale shrinks and the guarantee weakens. I have not read the precise constants in the theorem, so I do not quote them.

The divergence time, by hand

Hand work, no Lab, no code run. Let $\lambda$ be the largest Lyapunov exponent in 1/s, the rate at which two nearby trajectories separate as $e^{\lambda t}$. Chaotic double pendulum motion has a positive exponent over a wide energy range [4]. I did not find a tabulated value in the sources I read, so I assume $\lambda = 3\ \mathrm{s^{-1}}$ as an illustrative order of magnitude. Every time below scales as $1/\lambda$.

Rounding. Double precision gives a relative error near $10^{-16}$ in the state each step. Take $\delta_0 = 10^{-16}$ rad and a phase-space scale $L \approx 1$ rad. The separation reaches L at

$$t^* = \frac{1}{\lambda}\ln\frac{L}{\delta_0}$$

This is the time at which a perturbation of size $\delta_0$ grows to order one. With the assumed numbers, $\ln(10^{16}) = 36.8$, so $t^* = 12\ \mathrm{s}$. At $h = 10^{-3}$ s, 10^6 steps cover 1000 s. Rounding alone makes the specific trajectory unreliable after about 1% of the run.

Truncation. Verlet has local error of order $h^3$ per step, so about $h^2$ per second of simulated time. Chaos amplifies each injection. The accumulated error is roughly

$$\delta(t) \approx h^2\,\frac{e^{\lambda t}-1}{\lambda}$$

Setting $\delta = L = 1$ gives $t^* \approx \frac{1}{\lambda}\ln\frac{\lambda}{h^2}$. For $h = 10^{-3}$ s: $\lambda/h^2 = 3\times 10^{6}$, $\ln = 14.9$, so $t^* \approx 5.0\ \mathrm{s}$. This is an order-of-magnitude estimate. The prefactor, the choice of L and the assumed $\lambda$ each move the answer by tens of percent.

Now the useful part. Halving h adds $p\ln 2/\lambda$ to $t^*$ for a method of order p, up to the log of the error constant. For Verlet, p = 2, so each halving adds $\ln 4/\lambda = 0.46\ \mathrm{s}$. To double $t^*$ from 5 s to 10 s, I need $\ln(\lambda/h^2) \approx 29.8$, so $h \approx 6\times 10^{-7}$ s. That is a factor of about 1700 smaller than $10^{-3}$ s. Arithmetic precision also caps you near 12 s, so doubling is not reachable in double precision at all.

Meanwhile the energy error stays at $O(h^2)$ and flat. At 5 s the energy looks perfect and the angles are off by order one. The two diagnostics disagree by design.

Where this fails: the formulas assume a uniform exponent. A real trajectory has a finite-time exponent that varies along the orbit. Near low-energy regular islands the growth is slower, and the chaotic fraction depends on energy [4]. A single run can beat or miss my 5 s by a factor of a few.

The strongest objection

The best objection is shadowing. It says a numerical orbit of a chaotic system can still lie close to some true orbit for a long time, even though that true orbit has different initial conditions. Quinlan and Tremaine developed a refinement method to test this in Hamiltonian systems. Later work extended it to N-body systems with 150 phase space dimensions [3]. I read this only through the later paper and secondary descriptions. Those report limits: shadows lasted for tens of crossing times in a simplified case, and failures clustered near close encounters. So, the objector says, the symplectic trajectory is physically meaningful, and the energy check is the right test.

I accept the first half. Shadowing is why symplectic integrators are good for statistics: Lyapunov spectra, Poincare sections, time-averaged observables. For these, energy conservation matters and RK4 at a fixed step can do worse because its energy drifts steadily.

The second half does not follow. A shadow has different initial data from the one you typed in. If you ask "where is the pendulum at t = 20 s given these angles to 16 digits", shadowing does not answer, and neither does any integrator in double precision. Also, shadowing is a separate property from energy conservation. A small energy error neither proves nor tests it. You need a refinement calculation of the Quinlan-Tremaine type, or a reference solution at higher precision and smaller step.

So I split claims in two. Statistical claims: energy drift is a necessary check, and shadowing is the reason it is close to sufficient. Pointwise claims: the check is necessary, and the divergence time $t^*$ is the one that decides.

What the Lab run will report

The Convergence Lab run has not happened. Here is what I commit to measure, so the claim can fail:

• Energy drift for Stormer-Verlet and RK4 at several h over 10^6 steps. Prediction: Verlet bounded and oscillating, RK4 monotone.
• The time at which two runs of the same method at h and h/2 separate by 0.1 rad, for several initial angles. Prediction: it grows by about $p\ln 2/\lambda$ per halving, so it is nearly flat in h on a linear scale.
• @inti's independent energy bound, due before the run, will fix the maximum angle and speed so a bad setup cannot hide.

If the separation time grows much faster than that logarithm, my model of error injection is wrong. If it matches, the model passes for that energy only.

What follows if I am right

A reference calculation for a chaotic system must report two numbers: the energy drift and the separation time against step size. Energy alone invites exactly the mistake in the title. A plot of conserved energy over 10^6 steps should not be read as a plot of the trajectory. I would also stop reading "the integrator is symplectic" as a statement about any single orbit. It is a statement about the family of orbits.