Drag to rotate · scroll or pinch to zoom. Changing σ, ρ or β acts on the running trajectories immediately; Restart begins again from the start point. With many trajectories each trail is shortened automatically so the page stays smooth.

x(t), y(t) and z(t) for the same two runs (blue: from the start point, orange: the twin). Dashed lines mark the fixed points C±: the sign of x tells you which wing the orbit is circling. The thin grey line follows the animation.

What the Lorenz attractor is

In 1963 Edward Lorenz cut a model of convection down to three variables and three equations:

dx/dt = σ(y − x), dy/dt = x(ρ − z) − y, dz/dt = xy − βz

Roughly, x is the strength of the rolling motion, y the temperature difference between rising and falling fluid, and z how far the vertical temperature profile is from a straight line. σ is the Prandtl number, ρ measures how hard the layer is heated, and β is a geometric factor. With his values — σ = 10, ρ = 28, β = 8/3 — the solution never settles and never repeats: it loops around one of two centres, switches to the other at irregular moments, and stays inside a bounded region forever. That set is the attractor, and its two lobes are the butterfly people print on posters.

There is no randomness anywhere in the equations. Everything this page draws is a deterministic calculation, which is exactly what makes the butterfly effect surprising.

How to use the simulator

σ, ρ, β — drag or type. The running trajectories feel the change at once: drag ρ from 28 to 14 and the orbit spirals into a resting point.
Trajectories and Start apart by δ — trajectory i starts at x₀ + δ·i/(N − 1): 2 gives a pair exactly δ apart, 1,000 a tiny line of starts that the flow stretches over the attractor.
Start point — “On the attractor” first runs (1, 1, 20) forward for 40 time units; New start point runs it further by a random whole number of time units. “Exactly” uses the numbers you type.
Time step dt — the fixed step of the classical fourth-order Runge–Kutta integrator. The default 0.002 is the step the article used.
Trail, Colour by, Palette — only change the picture. Save PNG redraws the view at up to 4096 × 4096, Save SVG writes the trails as polylines, and Copy link to this view stores every setting, camera included, in the address.

The butterfly effect, measured

The Japanese write-up this tool comes from started two runs whose x differed by 10⁻⁹ and recorded the distance between them. With the classic parameters the article measured 0.00000026 at t = 5, 0.000035 at t = 10, 0.0025 at t = 15 and 0.027 at t = 20, and the distance first passed 1.0 at t = 23.8. On a log scale that is close to a straight line: exponential growth. For most of that time the two runs are drawn on top of each other and look like one line; then, quite suddenly, they are on different wings.

Press Two runs 10⁻⁹ apart to repeat it. The preset uses the article’s published setup — dt = 0.002, a start obtained by running (1, 1, 20) for 40 time units, and the offset added to x — and the integrator does the arithmetic in the same order as the article’s NumPy code, so the result is not just similar: the twin runs pass 1.0 at t = 23.84 (23.8 in the article) and the distances at t = 5, 10, 15, 20 round to the article’s table. Other start points give other times — four start points in this page’s check script gave 23.84, 20.13, 25.97 and 21.58 — so press New start point a few times to see how much it varies.

The Lyapunov exponent: one run is not enough

The growth rate of that distance is the largest Lyapunov exponent λ. The commonly quoted value for σ = 10, ρ = 28, β = 8/3 is about 0.9056; the article cites 0.906. Fitting a straight line to the single run above gave the article only 0.760, about 16% low. The stretching rate depends on where on the attractor the pair happens to be, and one run of 24 time units samples too few places.

The fix is Benettin’s method: follow a second trajectory a tiny distance away, and every so often measure how much the gap grew, shrink it back to its original size along the same direction, and keep a running average of the logarithm of the growth. The λ chart below the stats does exactly that, live, from the current start point. With the default settings it renews the gap every 0.5 time units starting from 10⁻⁸, and its running estimate passes through the article’s figures — 0.889 after 100 intervals, 0.917 after 500, 0.909 after 2000 and 0.9057 after 4000 — before continuing to 10,000 intervals (0.9042 in the check script). In the article the two readings differed by 0.15; the stats box shows both for whatever parameters you set.

Ten times the precision buys a fixed amount of time

If an error of size δ grows like δ·e^λt, it reaches a tolerance Δ after T = (1/λ) ln(Δ/δ). The prediction horizon depends on the logarithm of the initial precision, so each factor of ten adds the same ln 10 / λ, about 2.54 time units for λ = 0.9057. The article measured it directly for seven offsets from 10⁻³ to 10⁻¹⁵: 3.6, 8.7, 18.5, 23.8, 28.5, 33.9 and 34.0 time units until the runs were 1.0 apart — about 2.71 per factor of ten on a straight-line fit. Making the start a trillion times more precise did not make the forecast even ten times longer.

The table under the charts reruns that experiment for whatever parameters and start point you have set, puts its fitted slope next to ln 10 / λ and, when the settings match the article’s, lists the article’s own numbers alongside. Expect the slope and ln 10 / λ to disagree by a few tenths: the article got 2.71 against 2.54, and finite-time growth rates along one particular orbit are not the long-run average. Two more caveats from the article: t is dimensionless model time, not hours, and this is the toy convection model, not a weather forecast.

Fixed points, and what ρ does

The equations have a fixed point at the origin for every ρ, and for ρ > 1 two more, C± = (±√(β(ρ − 1)), ±√(β(ρ − 1)), ρ − 1) — the centres of the two wings, at (±8.485, ±8.485, 27) for the classic values. They are marked on the 3D view and as dashed lines in the x, y, z panel. C± lose stability at ρ_H = σ(σ + β + 3)/(σ − β − 1), which is 24.74 for σ = 10, β = 8/3. The presets walk through what that means; each description is what this page’s own solver shows, checked by its verification script:

ρ = 0.5 — the origin is the only fixed point and it attracts everything; the largest exponent is negative (−0.475 in closed form).
ρ = 14 — C± are stable. A fan of 100 starts spread over 10 units of x splits between the two resting points (35 and 65 of them in the check script), and λ settles at about −0.39, the real part of the eigenvalues at C±.
ρ = 21 — C± are still stable, but there is transient chaos: the cloud wanders over a butterfly-shaped region and then its members drop into C+ or C− one by one. In the check script all 300 starts of the preset settled; the median one took 49 time units and switched wings 8 times first.
ρ = 99.96 and ρ = 160 — periodic windows. After a short transient the orbit closes up (λ ≈ 0, the value for a limit cycle): at 99.96 the loop goes around the wings in the order R L L R L L before repeating, at 160 it is a symmetric R L L R loop.

Making an image of it

For wallpaper or prints, use one trajectory, a long trail and “Time along the trail” colouring, rotate to the angle you like and save a PNG at 3840 × 2160. Trails are drawn additively, so places the orbit visits often glow. A cloud of 1,000 starts with the Spectrum palette gives a different picture: the colours stay sorted along a thin filament, then fold into each other — the butterfly effect made visible.

Limits

Fixed-step RK4 is fourth order: halving dt divides the error over two time units by a factor that approaches 16 — 13.8 going from dt = 0.0025 to 0.00125 and 15.6 at the smallest steps in the check script. On a chaotic orbit even that error grows like any other, so after a few dozen time units the trajectory on screen is not “the” solution from that start — it is a different, equally valid orbit on the same attractor.
A start far from the attractor with a large step blows the integrator up — (200, 200, 200) with dt = 0.02 does it within a few steps. The page stops and says so; choose a smaller dt.
The Lyapunov estimate is a finite-time average from one start point. For a periodic orbit it approaches 0 only as 1/T, so small positive or negative values near 0 mean “periodic”.
Initial precision has a floor that depends on x₀, as the article also notes: x₀ + δ is rounded to the nearest double. At the article’s start x₀ = 3.3956 the 10⁻¹⁵ offset actually stored is 8.9 × 10⁻¹⁶, and once |x₀| reaches 16 an offset of 10⁻¹⁵ is rounded away entirely, leaving two identical runs. The table marks both cases.

Lorenz Attractor Simulator