The Lorenz System and Chaos

Learning Objectives

By the end of this section you will be able to:

1. Write down the three Lorenz equations and explain what each term models physically (Prandtl convection in a two-dimensional atmospheric slab).
2. Find the three fixed points $P_0, C_+, C_-$ and classify their stability as $\rho$ is varied.
3. Demonstrate sensitive dependence on initial conditions and read the predictability horizon directly off a log-separation plot.
4. Integrate the system with explicit fourth-order Runge–Kutta in Python, and reproduce the same trajectory with PyTorch while differentiating through the integrator.
5. Estimate the largest Lyapunov exponent $\lambda \approx 0.906$ from a single simulation and connect it to a predictability time scale $T \sim 1/\lambda$.

The Question Lorenz Was Asking

"Predictability: Does the flap of a butterfly's wings in Brazil set off a tornado in Texas?" — Edward Lorenz, AAAS 1972.

Edward Lorenz was a meteorologist at MIT in 1961 running a simple numerical weather model on an LGP-30. To save time he restarted a simulation from the middle, typing the printed state $0.506$ instead of the full internal precision $0.506127$. The new run agreed with the old one for a short while, then diverged completely. The difference of three parts per thousand had erased the future.

This was not a numerical bug. Lorenz had stumbled into a feature of deterministic non-linear systems that nobody had named yet: sensitive dependence on initial conditions. Two years later he distilled the issue into the three coupled ODEs we are about to study — a system simple enough to fit on a napkin yet rich enough to break every classical notion of long-range predictability.

The Three Equations

Lorenz simplified a Boussinesq fluid-convection model down to three Fourier modes. Calling them $x, y, z$ (no longer physical coordinates — they are mode amplitudes) gives the Lorenz system:

The three parameters are dimensionless ratios pinned down by the physics:

Three things to notice in the equations themselves:

1. The system is autonomous. Time $t$ appears only through the derivatives — the right-hand side depends only on the current state $(x, y, z)$.
2. It is non-linear. The terms $x z$ in $\dot y$ and $x y$ in $\dot z$ are products of state variables — there is no way to decouple them with a linear change of basis.
3. It is dissipative. The divergence of the velocity field is $\nabla\!\cdot\!\mathbf{F} = -\sigma - 1 - \beta = -41/3 < 0$. Volumes in phase space shrink at a constant rate $e^{-41 t / 3}$ — yet the trajectory is bounded, so the limit set has zero volume but non-zero fractal dimension (≈ 2.06).

Fixed Points and Their Stability

The fixed points are the simultaneous zeros of the three right-hand sides. From $\dot x = 0$ we get $y = x$. Substituting into $\dot y = 0$ gives $x(\rho - 1 - z) = 0$, and $\dot z = 0$ then yields $z = x^2/\beta$. Two branches emerge:

• Origin $P_0 = (0, 0, 0)$, valid for every $\rho$.
• Convection pair $C_\pm = \bigl(\pm\sqrt{\beta(\rho-1)},\ \pm\sqrt{\beta(\rho-1)},\ \rho - 1\bigr)$, valid only when $\rho > 1$.

For Lorenz's canonical values $\sigma = 10,\ \rho = 28,\ \beta = 8/3$, the Jacobian

evaluated at each fixed point has the following eigenvalue structure:

Every fixed point is unstable. The trajectory cannot settle anywhere. Instead it spirals out from $C_+$, gets flung past the saddle near the origin, captured by $C_-$, spirals out there, gets flung back — forever. The two lobes of the butterfly are the slow-manifold signatures of the two complex-conjugate eigenvalue pairs.

Sensitive Dependence: The Butterfly Effect

Take two initial conditions that differ by a tiny vector $\delta_0$. For a generic chaotic system the separation grows like

where $\lambda > 0$ is the largest Lyapunov exponent. For the standard Lorenz parameters numerical experiments give $\lambda \approx 0.906$. So even a microscopic perturbation of size $10^{-6}$ grows to order 1 in

Past $T$ the two trajectories are effectively independent samples from the same attractor: a forecast is useless. This is the precise mathematical content of the butterfly metaphor.

Interactive 3D Attractor

Two trails are integrated live in your browser with the SAME RK4 algorithm the rest of this section uses. They start $10^{-5}$ apart in the $x$ coordinate. Drag the sliders to leave the chaotic regime and come back; drag the canvas to rotate.

A few experiments worth running before you read on:

• Set $\rho = 0.5$. The origin is now globally attracting — both trails collapse to a single point and stay there. The system is stable. Boring. Predictable.
• Set $\rho = 15$. The trails wind into one of the two convection points $C_\pm$ and stop. Two basins of attraction; deterministic but no chaos.
• Crank $\rho$ back to 28. The trails now weave the butterfly pattern, visibly diverge after a few seconds, and end up on different lobes.
• Push $\rho = 100$. The chaotic regime breaks down into long periodic windows interlaced with bursts of chaos — a phenomenon called intermittency.

Worked Example — One RK4 Step by Hand

The simulator above runs thousands of RK4 steps per second. To make absolutely sure you can see how the system advances even by a single step, work the following micro-example by hand.

Lyapunov Exponent and Predictability Horizon

The Lyapunov exponent is the average exponential growth rate of an infinitesimal perturbation:

In practice we estimate it by integrating two nearby trajectories and fitting a straight line to $\ln|\delta(t)|$ over the window where the separation is small enough that the linearised dynamics still apply.

Given $\lambda$ and an initial measurement precision $\varepsilon_0$, the time until the perturbation grows to a tolerable bound $\varepsilon_{\text{tol}}$ is

Improving the initial measurement by a factor of $10$ only buys you $\ln(10)/\lambda \approx 2.5$ extra seconds of predictability. Improving it by a factor of $10^{6}$ only buys you $\approx 15$ seconds. This is the iron law of chaotic forecasting.

Sensitivity Visualizer

Drag the slider to shrink the initial offset all the way down to $10^{-12}$. Watch the log-separation plot on the right: the curve always grows linearly with the same slope — only its starting height shifts. That slope IS the Lyapunov exponent.

The left panel shows the actual $x(t)$ signals of the two trajectories. They are visually indistinguishable for a while, then disagree on which lobe of the butterfly they are on, then independently wander. The red dashed line marks the predictability horizon where the separation first exceeds 1 — the point past which an averaged forecast is the best you can do.

Python: Integrating the Lorenz System with RK4

Below is the entire engine that drove the 3D viewer above — pure NumPy, around 30 lines. Click any line to see what is in memory at that point.

1import numpy as np
2
3def lorenz_rhs(state, sigma=10.0, rho=28.0, beta=8.0/3.0):
4    """Right-hand side of the Lorenz system: returns [dx/dt, dy/dt, dz/dt]."""
5    x, y, z = state
6    return np.array([
7        sigma * (y - x),
8        x * (rho - z) - y,
9        x * y - beta * z,
10    ])
11
12def rk4_step(f, state, h):
13    """One classical 4th-order Runge-Kutta step for an autonomous ODE."""
14    k1 = f(state)
15    k2 = f(state + 0.5 * h * k1)
16    k3 = f(state + 0.5 * h * k2)
17    k4 = f(state + h * k3)
18    return state + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
19
20def integrate(state0, t_final=40.0, h=0.005):
21    """Integrate the Lorenz system and return the full trajectory."""
22    n_steps = int(t_final / h) + 1
23    traj = np.zeros((n_steps, 3))
24    traj[0] = state0
25    state = np.array(state0, dtype=float)
26    for k in range(1, n_steps):
27        state = rk4_step(lorenz_rhs, state, h)
28        traj[k] = state
29    return traj
30
31# Two trajectories from nearly identical starts
32A = integrate([1.0,        1.0, 1.0])
33B = integrate([1.0 + 1e-5, 1.0, 1.0])
34
35# When do they decisively diverge?
36sep = np.linalg.norm(A - B, axis=1)
37horizon = np.argmax(sep > 1.0) * 0.005
38print(f"Predictability horizon ≈ {horizon:.2f} s")  # ≈ 8 s
39
40# Estimate the largest Lyapunov exponent from the linear-growth phase
41mask = (sep > 1e-12) & (sep < 1e-2)
42t = np.arange(len(sep)) * 0.005
43slope, _ = np.polyfit(t[mask], np.log(sep[mask]), 1)
44print(f"λ ≈ {slope:.3f}   (textbook value ≈ 0.906)")

The single line that produced the headline number $\lambda \approx 0.906$ is slope, _ = np.polyfit(t[mask], np.log(sep[mask]), 1). Everything else is bookkeeping. The lesson is that Lyapunov exponents are not exotic theoretical objects; they are simple linear regressions performed on a log-transformed diagnostic from a single forward simulation.

PyTorch: Autodiff Through Chaos

The same RK4 routine rewritten in PyTorch buys us something NumPy cannot: gradients of the trajectory with respect to the system parameters, computed by automatic differentiation through every single step.

Why would you want this? Because the inverse problem — given an observed time series, what are the underlying parameters? — is exactly a gradient-descent problem on a loss like $\sum_t \|x_t^{\text{obs}} - x_t(\theta)\|^2$. Backpropagating through the integrator is how modern data assimilation systems learn parameters of climate, neural, and epidemiological models.

1import torch
2
3def lorenz_rhs(state, params):
4    """state: (3,) tensor. params: (3,) tensor [sigma, rho, beta]."""
5    x, y, z = state[0], state[1], state[2]
6    sigma, rho, beta = params[0], params[1], params[2]
7    return torch.stack([
8        sigma * (y - x),
9        x * (rho - z) - y,
10        x * y - beta * z,
11    ])
12
13def rk4_step(state, params, h):
14    k1 = lorenz_rhs(state,                    params)
15    k2 = lorenz_rhs(state + 0.5 * h * k1,     params)
16    k3 = lorenz_rhs(state + 0.5 * h * k2,     params)
17    k4 = lorenz_rhs(state + h * k3,           params)
18    return state + (h / 6.0) * (k1 + 2*k2 + 2*k3 + k4)
19
20# Make σ, ρ, β trainable so we can ask "how sensitive is the trajectory to them?"
21params = torch.tensor([10.0, 28.0, 8.0/3.0], requires_grad=True)
22state  = torch.tensor([1.0, 1.0, 1.0])
23
24# Roll the system forward for 200 steps (t = 1.0 s, well inside the predictable window)
25for _ in range(200):
26    state = rk4_step(state, params, h=0.005)
27
28# A scalar "observable" — say the height z at t = 1.0
29observable = state[2]
30observable.backward()
31
32print("z(1.0)            =", observable.item())
33print("∂z/∂σ, ∂z/∂ρ, ∂z/∂β =", params.grad.tolist())

Where Chaos Shows Up in the Real World

Once you can recognise the Lorenz signature — bounded, non-periodic, sensitive to initial conditions — you start to see it everywhere:

The deepest take-away is statistical predictability replacing pointwise prediction. We cannot say where the trajectory will be at $t = 20$, but we can compute the density of where many trajectories land — the SRB measure of the attractor. Most modern climate science lives in this statistical regime.

Common Pitfalls

• Trusting pointwise values past the horizon. Anything you print at $t = 50$ from a single Lorenz simulation is essentially a random sample from the attractor, not a deterministic answer.
• Using forward Euler. Euler is $O(h)$ and energy-non-conserving on oscillatory pieces. It will produce a trajectory that drifts off the attractor within a few seconds even at $h = 10^{-4}$. Always use RK4 or higher.
• Mistaking transients for the attractor. The first 5 seconds of any simulation depend strongly on the initial state. Discard the first $\sim 100$ time units before computing statistics.
• Fitting $\lambda$ over the saturation window. Once $|\delta|$ reaches the attractor diameter, exponential growth STOPS — the slope of $\ln|\delta|$ goes to zero. Only fit the linear-growth window.
• Comparing trajectories in single precision. float32 round-off introduces perturbations near $10^{-7}$ at every step. For chaos studies always use float64.

Summary

Lorenz's three equations look harmless. Their behaviour is anything but.

1. Deterministic chaos exists. A continuous, autonomous, three-dimensional ODE can produce trajectories that never repeat and are exponentially sensitive to initial conditions.
2. The Lyapunov exponent $\lambda \approx 0.906$ sets a hard limit on predictability. Time horizons grow only logarithmically in initial-condition precision.
3. The attractor is a fractal of dimension ≈ 2.06. Zero volume, non-zero structure, infinitely intricate.
4. RK4 is the right tool. It gives you trajectories that are statistically faithful — exactly the resolution chaos permits.
5. Modern computation re-opens the inverse problem. Autodiff through the integrator turns parameter inference into a gradient-descent exercise on a chaotic flow.

Final picture. The Lorenz attractor is the bridge between two cultures: the classical, equation-driven world of Newton and the modern, simulation-driven world of climate models and differentiable physics. Calculus gave us the equations. Computation gave us the courage to take them seriously.