Why Can’t Anyone Forecast the Weather a Month Ahead?

• A millionth off. The lines are one until step 1,784. Move the Step slider to 1,000 and read Gap: about 2 × 10⁻⁴, 160 times the start, still far too small to see.
• A billionth off, a thousand times closer. Surely this one stays together a thousand times longer? It parts at step 2,648, only 864 steps later.
• A trillionth off, another thousand times closer: step 3,367, only 719 steps more.
• Rounded. The second run starts from the first run’s numbers rounded to three decimal places, the way Edward Lorenz typed them into his computer in 1961 (the story is below). They part at step 1,233.

A pan of water as a tiny weather

The ground warms the air above it. Warm air rises, cools high up and sinks again, and the air goes round in great rolls: this is convection, the engine of most weather, from an afternoon thunderstorm to the monsoon. A pan of water on a stove does the same.

In 1963 Edward Lorenz, a meteorologist at the Massachusetts Institute of Technology in the United States, cut convection down to the smallest model he could find that still behaved like weather. He borrowed it from his colleague Barry Saltzman’s equations and kept only three numbers:

• x: how fast the roll of air turns. When x is positive it turns one way; when it is negative, the other way.
• y: how different the temperature is between the rising and the sinking air.
• z: how far the temperature from bottom to top differs from a steady, even change.

And three rules for how fast each one changes, with time measured in a unit of the model’s own:

$$\text{change in } x = \sigma (y - x) \qquad \text{change in } y = x(\rho - z) - y \qquad \text{change in } z = xy - \beta z$$

with σ = 10, ρ = 28 and β = 8⁄3, the values Lorenz chose (ρ is how strongly the pan is heated). That is the whole weather of this topic, with no randomness anywhere in it: start it from the same three numbers and it does exactly the same thing every time.

One step at a time

In Grade 11 you moved a thrown ball forward with Euler’s method: work out how fast things are changing now, and move that far in one small time step. Lorenz used a slightly better step, which he called a double approximation:

1. Take a trial step, with the rates of change at the start, exactly as Euler would.
2. Work out the rates of change again at the end of the trial step.
3. Take the real step from the start, with the average of the two rates.

Averaging the start and the end of the step makes up for the rates changing during it, so each step’s error is far smaller than Euler’s. Each step is a hundredth of the model’s unit of time; Lorenz’s computer took about a second for each.

Why a tiny gap grows

In Figure 13.1, two runs that start very close follow almost the same path. But near the middle, where the path decides which loop to go round next, nearby paths are pulled apart: one goes left, its neighbor right. Elsewhere they are squeezed together again; on average the stretching wins.

This program watches the gap between two runs a billionth apart. It imports the model from lorenz.py, a file this page keeps for you (you will write the same code yourself in Step 1):

from lorenz import step x, y, z = 0, 1, 0for n in range(1000):    x, y, z = step(x, y, z)a = (x, y, z)b = (x + 1e-9, y, z)                     # a billionth offfor n in range(2801):    if n % 200 == 0:        print(f'step {n:4}: gap {abs(a[0] - b[0]):.1e}')    a = step(*a)    b = step(*b)

step    0: gap 1.0e-09step  200: gap 7.2e-10step  400: gap 2.2e-09step  600: gap 3.9e-09step  800: gap 5.3e-08step 1000: gap 1.6e-07step 1200: gap 5.3e-07step 1400: gap 6.2e-06step 1600: gap 1.0e-04step 1800: gap 2.6e-04step 2000: gap 1.0e-03step 2200: gap 1.2e-02step 2400: gap 2.0e-01step 2600: gap 6.9e-01step 2800: gap 2.4e+00

1e-9 is how Python writes 1 × 10⁻⁹, a billionth, and :.1e prints a number the same way. (step(*a) hands the three numbers in a to step one by one.)

The gap doesn’t grow evenly; over the first 200 steps it even shrinks. But every few hundred steps it is about ten times bigger, and by step 2,800 it is as big as the swings of the weather itself. Growth by the same factor in the same time is called exponential, like a population doubling: tiny for a long time and then, suddenly, everything.

All in Python, with pixels, this page’s own module for showing a list of rows as a picture (on your own computer you would use a library such as Pillow).

Step 1: Lorenz’s computer, again

Here are the model and its steps, as the lists above say. To check it is Lorenz’s, compare it with the table his computer printed: x, y and z every five steps from his start (0, 1, 0), each times ten, and only the digits before the decimal point (37.009 printed as 37).

SIGMA = 10RHO = 28BETA = 8 / 3DT = 0.01                     # one step is a hundredth of the model's time def slopes(x, y, z):    # how fast x, y and z are changing now    return SIGMA * (y - x), x * (RHO - z) - y, x * y - BETA * z def step(x, y, z):    dx, dy, dz = slopes(x, y, z)    tx, ty, tz = x + dx * DT, y + dy * DT, z + dz * DT    # a trial step    ex, ey, ez = slopes(tx, ty, tz)                        # how fast things change at its end    # the real step uses the average of the two    return x + (dx + ex) / 2 * DT, y + (dy + ey) / 2 * DT, z + (dz + ez) / 2 * DT # Lorenz's Table 1 of 1963: x, y and z times ten, whole part onlyLORENZ = {0: (0, 10, 0), 5: (4, 12, 0), 10: (9, 20, 0), 15: (16, 36, 2), 20: (30, 66, 7),          25: (54, 115, 24), 30: (93, 192, 74), 35: (150, 268, 201), 40: (195, 234, 397)} x, y, z = 0, 1, 0                      # Lorenz's startprint('step      yours        Lorenz, 1963')for n in range(41):    if n in LORENZ:        yours = (int(10 * x), int(10 * y), int(10 * z))        his = LORENZ[n]        mark = '' if yours == his else '  <- differs'        print(f'{n:4} {yours[0]:5}{yours[1]:5}{yours[2]:5}   {his[0]:5}{his[1]:5}{his[2]:5}{mark}')    x, y, z = step(x, y, z)

step      yours        Lorenz, 1963   0     0   10    0       0   10    0   5     4   12    0       4   12    0  10     9   20    0       9   20    0  15    16   37    2      16   36    2  <- differs  20    30   66    7      30   66    7  25    54  116   24      54  115   24  <- differs  30    94  192   75      93  192   74  <- differs  35   150  269  201     150  268  201  <- differs  40   196  234  398     195  234  397  <- differs

The same weather, to within one in the last printed figure. In each row that differs, yours is a little over a whole number and his a little under it: your y at step 15 is 37.009 times ten; his was below 37. Your laptop or phone keeps about 16 digits of every number; Lorenz’s computer, a Royal McBee LGP-30 that he got in about 1958, worked with fewer and rounded them its own way. The two agree to about two decimal places, and it doesn’t matter which is a little off. For now.

Step 2: your computer against his

Lorenz printed a second table: the state at every step from 0 to 6,000 at which z reached a peak. At each one, x says which way the roll was turning. Here are his 83 peaks (x times ten), in a file the page keeps, lorenz1963.py:

Run your model for 6,000 steps, and at each of his steps ask: was your roll turning the same way as his?

from lorenz import stepfrom lorenz1963 import TABLE2 x, y, z = 0, 1, 0mine = []                         # your x after each step, 0 to 6,000for n in range(6001):    mine.append(x)    x, y, z = step(x, y, z) # at each step Lorenz printed: did his roll and yours turn the same way?agree = [(his > 0) == (mine[n] > 0) for n, his in TABLE2]k = agree.index(False)            # the first time they didn'tn, his = TABLE2[k]print(f'Step {n}: Lorenz {his / 10:.1f}, you {mine[n]:.1f}')print(f'Before step {n}: the same way at all {k} of his steps')rest = agree[k:]print(f'From step {n} on: the same way at {sum(rest)} of {len(rest)}')

Step 2110: Lorenz 14.6, you -8.4Before step 2110: the same way at all 31 of his stepsFrom step 2110 on: the same way at 26 of 52

For more than 2,000 steps your computer and Lorenz’s tell the same story, peak after peak, though the numbers were drifting apart: at step 1,643 his x was −13.9 and yours −13.0; at step 2,029, −13.5 and −9.2. Then, at step 2,110, his roll is turning one way and yours the other, and from then on they agree only half the time: 26 of 52, what two coins tossed side by side would do. The equations and the start are the same. The difference is how two computers, over 60 years apart, round the digits they can’t keep. Lorenz knew this: “As simple a change as transferring the experiment to a different computer with a different rounding-off procedure can cause the experiments to diverge,” he said in 1991.

Step 3: two runs, a gap apart

Now the widget’s experiment, as your own program. run gives back the list of x’s of one run; parted finds the first step at which two runs’ x’s differ by more than 1, a difference anyone would see (x itself swings between about −20 and 20).

from lorenz import step def run(x, y, z, steps):    xs = []    for n in range(steps + 1):        xs.append(x)        x, y, z = step(x, y, z)    return xs def parted(a, b):    # the first step at which the two runs' x differ by more than 1    for n in range(len(a)):        if abs(a[n] - b[n]) > 1:            return n    return None x, y, z = 0, 1, 0for n in range(1000):                  # partway through, like Lorenz's printout    x, y, z = step(x, y, z)print(f'Step 1000: x = {x}, y = {y}, z = {z}') first = run(x, y, z, 4000)second = run(round(x, 3), round(y, 3), round(z, 3), 4000)print('Rounded to 3 decimals: parted at step', parted(first, second)) last = Nonefor gap in [1e-3, 1e-6, 1e-9, 1e-12]:    n = parted(first, run(x + gap, y, z, 4000))    more = '' if last is None else f'   ({n - last} steps more)'    print(f'x off by {gap:g}: parted at step {n}{more}')    last = n

Step 1000: x = -5.923047633983275, y = -5.899584473251274, z = 24.000352340705373Rounded to 3 decimals: parted at step 1233x off by 0.001: parted at step 1080x off by 1e-06: parted at step 1784   (704 steps more)x off by 1e-09: parted at step 2648   (864 steps more)x off by 1e-12: parted at step 3367   (719 steps more)

These are the widget’s numbers: it does the same arithmetic in the same order. Rounding to three decimals moved the start by less than a thousandth (z was 24.00035, typed as 24.0), and the runs still parted after 1,233 steps.

Now draw it: a graph made of squares, one column for every 20 steps, one row for each whole number of x, the first run blue, the second orange, black where both share a square.

from lorenz import stepfrom pixels import show_grid def run(x, y, z, steps):    xs = []    for n in range(steps + 1):        xs.append(x)        x, y, z = step(x, y, z)    return xs x, y, z = 0, 1, 0for n in range(1000):    x, y, z = step(x, y, z)first = run(x, y, z, 4000)second = run(x + 1e-6, y, z, 4000)          # a millionth off # one column every 20 steps; one row for each whole number of x, 20 at the toppicture = [['.'] * 201 for row in range(41)]for col in range(201):    picture[20][col] = '-'                      # the line x = 0    for xs, mark in [(first, 'a'), (second, 'b')]:        row = 20 - max(-20, min(20, round(xs[col * 20])))        if picture[row][col] in '.-':            picture[row][col] = mark        elif picture[row][col] != mark:            picture[row][col] = 'c'             # both runs in this squareshow_grid(picture, {'.': 'white', '-': 'lightgray', 'a': 'steelblue', 'b': 'darkorange', 'c': 'black'})

Picture

Change 1e-6 to 1e-12, Run, and watch the black line last longer, but not much longer.

Step 4: a forecast from 50 runs

If nobody can measure today’s weather exactly, run the model many times from starts a tiny bit different, and see how many runs agree: an ensemble. Here are 50 runs, each starting up to a thousandth away; x > 0 might stand for “rain” and x < 0 for “dry”.

import randomfrom lorenz import stepfrom pixels import show_grid random.seed(2083)x, y, z = 0, 1, 0for n in range(1000):    x, y, z = step(x, y, z) # 50 runs, each starting up to a thousandth away in xruns = [(x + random.uniform(-0.001, 0.001), y, z) for member in range(50)]picture = [[] for member in runs]          # one row for each runfor n in range(3001):    if n % 20 == 0:        for k in range(len(runs)):            picture[k].append('+' if runs[k][0] > 0 else '-')    if n % 250 == 0:        plus = sum(1 for r in runs if r[0] > 0)        print(f'step {n:4}: x > 0 in {2 * plus:3}% of the runs')    runs = [step(*r) for r in runs]show_grid(picture, {'+': 'darkorange', '-': 'steelblue'})

step    0: x > 0 in   0% of the runsstep  250: x > 0 in   0% of the runsstep  500: x > 0 in   0% of the runsstep  750: x > 0 in   0% of the runsstep 1000: x > 0 in   0% of the runsstep 1250: x > 0 in  46% of the runsstep 1500: x > 0 in  40% of the runsstep 1750: x > 0 in  60% of the runsstep 2000: x > 0 in  44% of the runsstep 2250: x > 0 in  48% of the runsstep 2500: x > 0 in  42% of the runsstep 2750: x > 0 in  54% of the runsstep 3000: x > 0 in  40% of the runs

Picture

Each row is a run, left to right in time. For over a thousand steps the stripes line up: all 50 runs tell the same story, and the forecast is sure. Then they stop lining up, and from step 1,250 the runs split roughly half and half: past that point the forecast is a coin toss. (random.seed(2083) gives the page and your run the same random starts; take it out and the percentages change, but not the picture’s shape.)

• Your computer and Lorenz’s agreed for over 2,000 steps, then went their own ways, with the same equations and the same start: only the rounding differed.
• A thousand times closer buys only a few hundred steps: 704, 864 and 719 more. Ten doublings make a thousandfold (2¹⁰ = 1,024), so the gap doubles roughly every 76 steps.
• An ensemble is sure, and then it isn’t: 50 runs agreed for over a thousand steps, then split about half and half.

This is the answer to the question. A forecast starts from today’s weather as measured, and no measurement is exact: there are only so many weather stations, balloons and satellites, and the air between them is guessed. In a system like this, the error in the start grows exponentially until it is as big as the weather. Better instruments push the limit back a little; nothing pushes it back a lot. This sensitive dependence on the start is called chaos: not randomness, since every step follows a rule, but a rule that turns differences too small to measure into differences too big to miss.

A smaller step doesn’t help

In lorenz.py, make DT = 0.005 and double every step count in Step 3. Each thousandfold better start now buys about 1,100 to 2,300 more steps, but these steps are half as long: about the same stretch of the model’s time as before. The parting belongs to the equations, not the method.

A calm sky

Make RHO = 20 in lorenz.py and run Step 3. The runs never part: the roll settles to turning steadily one way, and every small difference dies away. The equations are chaotic only when the heating (ρ) is strong enough.

Smaller than the computer can hold

In Step 3, make the list of gaps [1e-16] and Run. The runs never part (None), because x + 1e-16 is exactly x: on a number near 6 the gap is smaller than the last of the 16 digits the computer keeps Grade 11: Grade 11’s floating point. Now try [1e-15]. That is a different number, one step in the last digit, yet the runs never part either. This shows why:

from lorenz import step x, y, z = 0, 1, 0for n in range(1000):    x, y, z = step(x, y, z)a, b = (x, y, z), (x + 1e-15, y, z)for n in range(4):    print(n, b[0] - a[0])    a, b = step(*a), step(*b)

0 8.881784197001252e-161 8.881784197001252e-162 0.03 0.0

Two steps on, the step’s own arithmetic has rounded the gap away, and the runs are one run. Down here, the computer’s rounding decides the weather.

An open question

Is there a hard limit of about two weeks, however good the measurements and computers become? The 2019 study says one probably exists, built into the atmosphere, but the atmosphere can’t be run twice to check. And slow things can be forecast further ahead: the sea’s temperature, and with it whether a whole monsoon season is likely to be wetter or drier than usual.

Three challenges, each checked on the page: write your function in the program above the checker (press Edit), then press Check my program. The examples show what your function should give; there are more tests you can’t see.

★ Where do two lists part?

Write first_apart(a, b, limit): the first place in two lists of numbers where they differ by more than limit (as a size: −3 and 3 differ by 6). If they never do, give back −1. If one list is shorter, only compare as far as it goes.

def first_apart(a, b, limit):    # Give back the first place where a and b differ by more than limit, or -1.    return -1 print(first_apart([1, 2, 3, 4], [1, 2.5, 5, 4], 1))

★★ When do they part?

Write parting_step(start, gap, most). start is a list [x, y, z]. Run Lorenz’s model from it, and a second run from the same start with gap added to x. Give back the number of steps after which the two x’s first differ by more than 1 (0 if they already do at the start), or −1 if they haven’t parted after most steps. step is given: don’t change it.

SIGMA, RHO, BETA, DT = 10, 28, 8 / 3, 0.01 def step(x, y, z):    dx, dy, dz = SIGMA * (y - x), x * (RHO - z) - y, x * y - BETA * z    tx, ty, tz = x + dx * DT, y + dy * DT, z + dz * DT    ex, ey, ez = SIGMA * (ty - tx), tx * (RHO - tz) - ty, tx * ty - BETA * tz    return x + (dx + ex) / 2 * DT, y + (dy + ey) / 2 * DT, z + (dz + ez) / 2 * DT def parting_step(start, gap, most):    # Give back the step at which the two runs' x first differ by more than 1, or -1.    return -1 print(parting_step([-5.923047633983275, -5.899584473251274, 24.000352340705373], 1e-6, 6000))

★★★ How many decimals?

Lorenz typed in his numbers rounded to three decimals. Write decimals_needed(start, steps): the fewest decimals d, from 1 to 15, such that a run started from start rounded to d decimals (each of x, y and z, with round(v, d)) has not parted from the run started from start itself (x’s differing by more than 1) within steps steps.

SIGMA, RHO, BETA, DT = 10, 28, 8 / 3, 0.01 def step(x, y, z):    dx, dy, dz = SIGMA * (y - x), x * (RHO - z) - y, x * y - BETA * z    tx, ty, tz = x + dx * DT, y + dy * DT, z + dz * DT    ex, ey, ez = SIGMA * (ty - tx), tx * (RHO - tz) - ty, tx * ty - BETA * tz    return x + (dx + ex) / 2 * DT, y + (dy + ey) / 2 * DT, z + (dz + ez) / 2 * DT def decimals_needed(start, steps):    # Give back the fewest decimals that keep the rounded run with the real one for steps steps.    return 15 print(decimals_needed([-5.923047633983275, -5.899584473251274, 24.000352340705373], 2000))

• Grade 11, Modelling and Simulation: the update rule, Euler’s method and why a step’s errors can drift, and Monte Carlo simulations, the idea behind an ensemble.
• Grade 11, Numbers with fractions: why a computer keeps only about 16 digits, and what happens to the rest.
• Grade 12, How sure? Uncertainty by resampling: giving an answer with how far to trust it, as an ensemble does.
• Where Does the Rain Come From? (Science, Going Further 6–8): how real weather works and what a “70% chance of rain” means.
• How can one formula be calm, then wild?: chaos from a single formula for a population, and the famous branching picture. (Coming soon.)
• When does a forest fire stop spreading?: another simulation where a small change makes all the difference. (Coming soon.)