Which numerical results can you trust? The Lorenz system as a test case
A computational physics course talk, From the convection equations to a trustworthy numerical solution. The slides run in the browser: soundadam.github.io/lorenz63. Use the arrow keys to move through them and N for the speaker notes. The slides and notes are in Chinese; the plots, equations and live demos read the same in any language.
A paradox
Integrate the Lorenz equations (σ=10, ρ=28, β=8/3) with DOP853 from
scipy.solve_ivp, same initial condition, once with rtol=1e-6 and once with
rtol=1e-12. For the first ten time units the two trajectories agree closely.
Then they come apart until they are as far apart as the attractor is wide.
Which one is right? The intuitive answer is the tighter tolerance. That holds for t<10 and means nothing for t>40.
Precision buys only logarithmic time
Sweep the tolerance from 1e-4 to 1e-12 and record when each run first moves more than 1 away from a reference solution (rtol=1e-13):
| rtol | Divergence time |
|---|---|
| 1e-4 | 5.69 |
| 1e-6 | 10.62 |
| 1e-8 | 16.11 |
| 1e-10 | 20.71 |
| 1e-12 | 26.88 |
Divergence time against ln(1/rtol) is a straight line, and its slope is one over the largest Lyapunov exponent (λ₁ ≈ 0.904). Per digit: each extra digit of precision buys only ln 10 / λ₁ ≈ 2.55 time units. Double precision bottoms out at 10⁻¹⁶, so a single trajectory cannot be predicted much past t ≈ 41, however tight the tolerance.
But the statistics do not move
Take three runs at rtol 1e-6, 1e-9 and 1e-13. As trajectories, the first two are lost by t=10.5 and t=18.3. Yet their Lorenz maps (each local maximum of z against the next) coincide to within the line width, with a thickness of 0.0086% of the range of z. The invariant densities of z differ by at most 1% in total variation, and most of that is finite sampling. The Lyapunov spectrum is just as insensitive to the tolerance: the three exponents sum to the phase space divergence −(σ+1+β) to within about 10⁻⁹.
The same data are worthless as trajectories and intact as statistics.
The test depends on the kind of system
| System | How to check it | One long trajectory |
|---|---|---|
| Conservative (Hamiltonian) | Invariants stay bounded; use a symplectic scheme | Trustworthy, error grows polynomially |
| Dissipative, not chaotic | Converges to the right steady state or limit cycle | Trustworthy |
| Dissipative, chaotic | Statistics insensitive to method and tolerance | Not trustworthy beyond order 1/λ₁ |
Conservative systems are checked against their invariants, and Lorenz has no conserved quantity to check. In a chaotic system a single trajectory cannot be trusted, but its statistics can. The talk works through the first row on the Kepler problem: the energy error of Störmer–Verlet stays flat over 1000 orbits while that of RK4 drifts linearly, and for both the time bought per unit of cost is a power law. For Lorenz it is a logarithm.
What to try in the slides
- The cover replays two trajectories whose starting points differ by 10⁻⁵; |Δ| reaches 1 at t ≈ 13.6
- The opening paradox: pick an rtol and see at which t it leaves the reference solution
- Drag ρ and watch the convection roll turn
- Shrink the initial separation and watch the divergence time recede only logarithmically
The demos integrate in the browser with DOP853, rewritten in JavaScript from the scipy coefficient table. The two printouts used at the lectern are on the same site: a mind map (A4 landscape) and a board cheat sheet (A4 portrait).
Every number here was computed by a Python package, and the key claims are tests in it, for instance that the three Lyapunov exponents sum to the divergence and that divergence time is linear in ln(1/rtol).