Chaos lab

1The double pendulum

Two rods of 1 m and two bobs of 1 kg, with no friction. Lagrange's method turns the energy of the swing into two equations for the angular accelerations. They are completely deterministic: give the same θ₁, θ₂, ω₁ and ω₂ and you get the same future every time.

The catch is the coupling. The second rod is flung around by the first and tugs back on it, and the sin(θ₁ − θ₂) terms make every small difference feed on itself. At small angles the pendulum is nearly two linked simple pendulums and stays calm; start it above about 90° and it tumbles.

Energy has to stay fixed, so it doubles as a check on the integrator below.

θ₁″ = [−3g sin θ₁ − g sin(θ₁ − 2θ₂) − 2 sin Δ (ω₂² + ω₁² cos Δ)] / (3 − cos 2Δ)Δ = θ₁ − θ₂, g = 9.81 m/s², equal masses and lengthsθ₁ = …, ω₁ = … → θ₁″ = …

θ₂″ = 2 sin Δ (2ω₁² + 2g cos θ₁ + ω₂² cos Δ) / (3 − cos 2Δ)the second bobθ₂ = …, ω₂ = … → θ₂″ = …

E = ω₁² + ½ω₂² + ω₁ω₂ cos Δ − 2g cos θ₁ − g cos θ₂kinetic plus potential energy, in joulesE = …, drift since release …

2Stepping forward: Runge-Kutta 4

There is no formula for where a double pendulum will be in ten seconds, so the page walks there in tiny steps. The simplest walk, Euler's method, uses the slope at the start of each step and overshoots; for a swinging pendulum it pumps in energy until the motion is nonsense.

The classic Runge-Kutta method samples the slope four times per step (at the start, twice in the middle and at the end) and blends them, so its error per step shrinks like h⁵. Every pendulum and every Lorenz rider here uses it, with … steps per frame.

Even RK4 cannot beat chaos: rounding errors of 10⁻¹⁶ grow just like the offsets you set. What it can do is keep the energy right, which is why the curves still look like a real pendulum long after the exact path is lost.

k₁ = f(yₙ), k₂ = f(yₙ + h k₁ / 2), k₃ = f(yₙ + h k₂ / 2), k₄ = f(yₙ + h k₃)four slopes per step of size hh = …, steps so far …

yₙ₊₁ = yₙ + h (k₁ + 2k₂ + 2k₃ + k₄) / 6local error O(h⁵), global error O(h⁴)

Euler: yₙ₊₁ = yₙ + h f(yₙ)the same pendulum, the same h, for 10 senergy drift: RK4 …, Euler …

3The Lorenz system

In 1963 Edward Lorenz boiled a model of convection, warm air rising and cool air sinking, down to three variables: x is the speed of the rolling air, y and z are temperature differences. With σ = 10, β = 8/3 and ρ = 28 a single point traces a butterfly that never repeats and never crosses itself.

The two wings circle the fixed points C±, where the air would roll steadily one way or the other. Above ρ ≈ 24.74 they are unstable, so the point loops round one, gets flung across to the other, and the number of loops before each switch is unpredictable. Below it, riders spiral into a fixed point and the chaos dies.

The same equations drive the Lorenz scene in the visualiser, where the music pushes ρ about. The attractor is a fractal too, with dimension about 2.06: more than a surface, less than a solid.

dx/dt = σ(y − x), dy/dt = x(ρ − z) − y, dz/dt = xy − βzσ = 10, β = 8/3ρ = …; (x, y, z) = …→ (…)

C± = (±√(β(ρ − 1)), ±√(β(ρ − 1)), ρ − 1)the two fixed points, the eyes of the wingsC± = …

ρ_H = σ(σ + β + 3) / (σ − β − 1) ≈ 24.74above this, C± are unstable and the motion is chaotic…

4Three bodies

Three equal stars pulling on each other by Newton's gravity, borrowed whole from N-body gravity: the same starting orbits, the same force sum and the same Euler, Verlet and RK4 steps, so the two pages always agree. Each copy moves star 1 a little further along x, and the page measures how fast the copies part.

The figure-eight is one of the rare stable choreographies: copies drift apart only slowly and λ stays near zero. The Lagrange triangle with equal masses is unstable, so its copies peel apart exponentially. Raise the nudge to knock the eight off its loop, or switch to Euler and watch numerical error look like chaos. Live sky's planets are the calm case: one sun so heavy that each orbit is nearly a two-body ellipse.

aᵢ = Σⱼ G mⱼ (rⱼ − rᵢ) / (|rⱼ − rᵢ|² + ε²)^3/2G = m = 1, three stars…

|δ| = √Σ (Δx² + Δy² + Δvₓ² + Δvᵧ²)the gap between two copies in phase space, all three starsλ = …

5The Lyapunov exponent

Measure the gap |δ(t)| between two nearby starts. In a chaotic system it grows exponentially, |δ(t)| ≈ |δ₀| e^(λt), so its logarithm climbs in a straight line: the first chart. The slope of that line is the largest Lyapunov exponent λ, and the page fits it by least squares from just after release until the gap saturates at the size of the whole system.

The fit only has a short straight stretch to work with, so the second chart uses Benettin's trick: a shadow start is kept a fixed tiny distance from the first, and after every frame the page notes how much the gap stretched, then pulls the shadow back. Averaging the logs of those stretches converges on λ without ever saturating. For the Lorenz butterfly at ρ = 28 the answer is about 0.906.

λ above zero means chaos, zero means a regular orbit that drifts slowly, below zero means everything settles to the same place. Many bodies at once, as in N-body, are chaotic the same way, and its three-body mode runs this same renormalised estimate live.

λ = lim_t→∞ lim_δ₀→0 (1/t) ln(|δ(t)| / |δ₀|)the largest Lyapunov exponentln |δ₀| = …, ln |δ(t)| = …

ln |δ(t)| ≈ ln |δ₀| + λtleast-squares line over the straight stretchfit over …: λ = …, R² = …

λ ≈ (1 / Nτ) Σ ln(dₖ / d₀)Benettin: stretch over each frame τ, then renormalise to d₀ = 10⁻⁸after … frames: λ = …

6The prediction horizon

If the gap grows like e^(λt), the time until a start error δ₀ grows to a tolerance Δ is only logarithmic in how good the measurement was. Measuring ten times better buys a fixed ln 10 / λ of extra time, no matter how good you already were. That is why weather forecasts stop at about two weeks: better instruments add hours, not months.

Try it: slide how much better the start is measured and watch the horizon barely move.

a three-line rule and a hinged rod already hold more future than any computer can reach.

T ≈ (1/λ) ln(Δ / δ₀)time until an error δ₀ reaches the tolerance ΔΔ = …, δ₀ = … → T = …measured … better → T = … (…)

t₂ = ln 2 / λdoubling time of a small errort₂ = …

The maths is the same kind of thing that makes maths art's Clifford attractor: a simple rule, repeated, that never settles.