Chapter 6Do it yourself

Every picture on this site comes out of a few lines of plain Python with nothing imported but the standard library. Here is the heart of it. Copy it, run it, change the numbers, break it.

Download chaos.py   The full source on GitHub

#!/usr/bin/env python3
"""chaos.py — the whole site's mathematics, no libraries, runs in a second.

    python3 chaos.py

Prints: a chaotic run you can't forecast, the average of a chaotic pile that never moves,
the weather model's Lyapunov exponent, and one feedback loop taming a run of shocks.
"""
import math, random


def logistic(r, x, n):
    "The rabbit rule: where you are, times the growth knob, times the room left."
    for _ in range(n):
        x = r * x * (1 - x)
        yield x


def lorenz_step(s, dt=0.006, sig=10, rho=28, beta=8/3):
    x, y, z = s
    def d(st):
        x, y, z = st
        return (sig*(y-x), x*(rho-z)-y, x*y-beta*z)
    k1 = d(s)
    k2 = d((x+dt/2*k1[0], y+dt/2*k1[1], z+dt/2*k1[2]))
    k3 = d((x+dt/2*k2[0], y+dt/2*k2[1], z+dt/2*k2[2]))
    k4 = d((x+dt*k3[0],   y+dt*k3[1],   z+dt*k3[2]))
    return tuple(s[i] + dt/6*(k1[i]+2*k2[i]+2*k3[i]+k4[i]) for i in range(3))


def lyapunov(dt=0.006, n=60000, d0=1e-9):
    "Benettin's method: run a twin a hair away, measure the growth, pull it back, repeat."
    a = (1.0, 1.0, 1.0)
    for _ in range(4000):            # settle onto the attractor
        a = lorenz_step(a, dt)
    b = (a[0]+d0, a[1], a[2])
    total = 0.0
    for _ in range(n):
        a, b = lorenz_step(a, dt), lorenz_step(b, dt)
        gap = math.dist(a, b)
        total += math.log(gap/d0)
        f = d0/gap                   # rescale the twin back to distance d0
        b = (a[0]+(b[0]-a[0])*f, a[1]+(b[1]-a[1])*f, a[2]+(b[2]-a[2])*f)
    return total/(n*dt)


def die_average(n, seed=70118):
    rng = random.Random(seed)
    s = 0
    for i in range(1, n+1):
        s += rng.randint(1, 6)
    return s/n


def governor(k=0.28, n=240, seed=4669):
    "One feedback loop against a run of shocks: subtract k of the error each step."
    rng = random.Random(seed)
    drift = held = 0.0
    for _ in range(n):
        shock = rng.gauss(0, 1) * 0.5
        drift += shock                       # no correction
        held += shock; held -= k * held      # pull back k of the error
    return drift, held


if __name__ == "__main__":
    run = list(logistic(4.0, 0.4, 12))
    print("chaotic run (unforecastable):", " ".join(f"{v:.3f}" for v in run))

    pile = list(logistic(4.0, 0.31415926, 300000))
    print("average of 300,000 chaotic values:", round(sum(pile)/len(pile), 4), "(theory 0.5)")

    print("Lyapunov exponent of the weather model:", round(lyapunov(), 3), "(known 0.906)")

    print("die average over 3,000 rolls:", round(die_average(3000), 3), "(true 3.5)")

    drift, held = governor()
    print(f"shock left alone drifted to {drift:.2f}; one feedback loop held it to {held:.2f}")

What it prints

Run python3 chaos.py and you get, in order: a dozen steps of the fully chaotic rabbit map (you will not be able to forecast the thirteenth); the average of three hundred thousand of those same chaotic values, which sits on 0.5 every time; the weather model's Lyapunov exponent, measured the way chapter three measured it, landing near the textbook 0.906; the average of three thousand die rolls, near 3.5; and a single feedback loop pulling a run of shocks back to centre.

Six functions, one file, no dependencies. The chaos, the two tools, and the proof that the pile holds still while the grains jump — all of it fits on a page.