Deep
Patterns

Hidden Patterns asked whether you could spot the algorithm. This asks whether you can build it, understand why it scales the way it does, and find it in the real world.

This will take several sessions. That's the point.

The Starry Night — Van Gogh

The Starry Night

Van Gogh · 1889

Flow fields → O(n) advection ↓
The Great Wave — Hokusai

The Great Wave

Hokusai · c.1831

Fluid dynamics → $1M open ↓
The Scream — Munch

The Scream

Munch · 1893

Reaction–diffusion → Turing ↓
01 / III paintings Van Gogh · 1889

The Starry Night

B — The algorithm

Particles advected through a curl-noise flow field

Van Gogh's sky is not random. The brushstrokes follow coherent streamlines that swirl, split and rejoin — the visual signature of a vector field. We reproduce it by seeding thousands of massless particles and pushing each one along the local flow, frame after frame.

The field itself is the curl of a smooth scalar noise function. Taking the curl guarantees the flow is divergence-free — particles circulate rather than pile up — which is exactly why the result reads as turbulence rather than noise.

pt+1  =  pt  +  (∇ × ψ)(pt) · Δt

Eq. 1 — explicit Euler advection step

# one frame — O(n), n particles, no interaction
for p in particles:
    vx, vy = curl(psi, p.x, p.y)
    p.x += vx * dt
    p.y += vy * dt
    if offscreen(p):
        reseed(p)   # keep density constant

C — The complexity

Linear — and that is the whole point

Each particle reads the field once and moves. No particle looks at any other. That independence makes one frame O(n) — double the particles, double the work, no worse. The table below contrasts that with the O(n²) cost you would pay if every particle had to consider every other, as a naïve fluid solver does.

particles n O(n) O(n log n) O(n²) naïve
10010066410,000
1,0001,0009,9661,000,000
10,00010,000132,877100,000,000
1,000,0001,000,00019,931,5691,000,000,000,000
rows land · then numbers count · easeOutCubic

D — In the real world

Where linear flow fields earn their keep

VFX & games Smoke, fire and magic in real-time engines are millions of particles advected through exactly this kind of field — cheap because the cost stays linear, expensive if you needed exact inter-particle forces.
Meteorology Wind-map visualisations trace tracer particles over forecast velocity grids — the same advect-and-reseed loop, scaled to the globe. ECMWF runs particle ensembles over 137 vertical atmospheric levels.
Medical imaging Diffusion-tensor tractography follows particles along neural fibre directions to reconstruct white-matter pathways in the brain — critical for surgical planning around tumours.
Your simulation runs a few thousand particles at 60 fps in a browser. A weather model runs billions on a supercomputer. The loop is identical — so what, exactly, did the supercomputer buy?
02 / III paintings Hokusai · c.1831

The Great Wave

B — The algorithm

A grid of heights evolving under the wave equation

Hokusai's wave captures the moment before a crest breaks — the peak of a travelling disturbance through water. The 2D wave equation describes how height disturbances propagate through a fluid surface at speed c. Each point in the grid depends only on its state at the two previous time steps and its four nearest neighbours.

Discretised with finite differences, this becomes an update rule that requires no matrix operations, no global solves. Each cell is independent; the wave emerges from the aggregate of local interactions.

∂²η/∂t²  =  c²∇²η

Eq. 1 — 2D wave equation (height η, speed c)

ηt+1  =  2ηt − ηt−1  +  c²Δt²∇2hηt

Eq. 2 — explicit finite-difference (Verlet integration)

# one time step — O(n), n = grid cells
for y in range(1, H-1):
    for x in range(1, W-1):
        lap = (h[y,x-1] + h[y,x+1]
             + h[y-1,x] + h[y+1,x] - 4*h[y,x])
        h_next[y,x] = (2*h[y,x] - h_prev[y,x]
                     + c**2 * lap) * damping
h_prev, h, h_next = h, h_next, h_prev

C — The complexity

Linear updates — but an unsolved problem hides inside

The update rule above is O(n) per time step, where n is the number of grid cells. Double the grid, double the work. No surprises there.

But notice what we quietly skipped. The simulation uses a simplified wave equation — not the Navier-Stokes equations that govern real fluid dynamics. Those are different in kind: non-linear, fully coupled, and despite over 170 years of effort, not fully understood.

grid cells n O(n) O(n log n) O(n²) direct
10010066410,000
1,0001,0009,9661,000,000
10,00010,000132,877100,000,000
1,000,0001,000,00019,931,5691,000,000,000,000
rows land · then numbers count

Clay Millennium Prize · $1,000,000

Navier–Stokes existence and smoothness

Given a smooth initial velocity field for a 3D incompressible fluid, do the Navier–Stokes equations have a smooth solution for all future time — or can they develop a singularity (a point where velocity becomes infinite) in finite time? The Clay Mathematics Institute designated this one of seven Millennium Prize Problems in the year 2000. It remains unsolved. Your simulation sidesteps the question by using a linear approximation. Real ocean-wave models do the same, for the same reason.

D — In the real world

Where wave equations earn their keep

Tsunami prediction NOAA's DART buoy network detects pressure disturbances in the deep ocean; wave propagation models identical in structure to your simulation project arrival times and wave heights. The 2004 Indian Ocean tsunami was the first major event partially predicted by such systems in real time.
Aeroacoustics Aircraft noise regulations require predicting how sound waves — governed by the wave equation — propagate from engines. The O(n) scaling per step makes full-aircraft simulation tractable on modern supercomputers. Engine placement on modern aircraft is partly an acoustic optimisation problem.
Seismology Earthquake ground motion is modelled by wave equations across 3D geological grids covering continents, with tens of billions of cells. The finite-difference scheme is the same; the engineering is in stable boundary conditions and parallelising across thousands of CPUs.
Your simulation shows a 2D surface. Real ocean waves are three-dimensional, non-linear, and break. Which of those three complications do you think introduces the most error into a tsunami prediction model — and how would you measure it?
03 / III paintings Munch · 1893

The Scream

B — The algorithm

Two chemicals, one unstable equilibrium, and Turing patterns

In 1952, Alan Turing — better known for his work on computation and codebreaking — published a paper arguing that certain chemical systems could spontaneously produce spatial patterns from a nearly uniform initial state. The turbulent sky of The Scream, with its swirling organisation, has something of this quality: a system driven from equilibrium, trying to arrange itself.

The Gray-Scott model captures the essence with two interacting chemicals U and V. U is fed in from outside at rate F. When U and V meet, they produce more V (autocatalytic reaction). V also decays at rate k. Crucially, the two diffuse at different rates — U fast, V slow. This asymmetry is the mechanism: V concentrates where it forms, while U replenishes the gaps, producing self-sustaining spatial structure.

∂u/∂t = Du∇²u − uv² + F(1−u)

∂v/∂t = Dv∇²v + uv² − (F+k)v

Eq. 1 — Gray-Scott reaction–diffusion model

# one Gray-Scott step — O(n), n = grid cells
for y in range(H):
    for x in range(W):
        lu = laplacian(u, x, y)   # 4-neighbour sum
        lv = laplacian(v, x, y)
        uvv = u[y,x] * v[y,x]**2  # autocatalytic term
        u_next[y,x] = u[y,x] + Du*lu - uvv + F*(1 - u[y,x])
        v_next[y,x] = v[y,x] + Dv*lv + uvv - (F+k)*v[y,x]
u, v = u_next, v_next

C — The complexity

Linear in space — but how long until the pattern?

The update rule is O(n) per step: each cell reads its four neighbours and applies a fixed formula. The table looks identical to the others. The interesting question is different: you cannot predict, from the initial conditions, how many steps will pass before a coherent pattern emerges.

The Gray-Scott system has a transient phase — a period of apparent disorder before the pattern crystallises. The length of that transient depends on F, k, and the initial perturbations in ways that resist analytical bounds. In biological contexts, this corresponds to the actual developmental timeline: a zebrafish's stripe pattern takes as long as it takes, regardless of how many cells are computing it.

grid cells n O(n) per step O(n log n) O(n²) naïve
10010066410,000
1,0001,0009,9661,000,000
10,00010,000132,877100,000,000
1,000,0001,000,00019,931,5691,000,000,000,000
rows land · then numbers count

D — In the real world

Where reaction–diffusion shows up, uninvited

Animal markings Turing's 1952 paper predicted the mathematical existence of stripes and spots before biologists had evidence for the mechanism. In 2012, researchers confirmed that the chemical signals governing mouse digit formation follow a Turing-type reaction-diffusion process — 60 years after the prediction.
Cardiac arrhythmia Spiral waves in the cardiac muscle — rotating patterns of electrical excitation — are a form of reaction-diffusion instability. They underlie ventricular fibrillation. The FitzHugh-Nagumo model used in cardiology is structurally similar to Gray-Scott, sharing the same autocatalytic mechanism.
Materials science Regular nanostructures in certain chemical vapour deposition processes form by reaction-diffusion dynamics. Engineers exploit this to produce self-assembling surface patterns at scales too small for conventional lithography, without directing each feature individually.
The patterns depend sensitively on F and k. What experiment could you design to test whether a real biological pattern — say, an angelfish's stripes — is genuinely produced by a Turing mechanism, and not by a developmental clock or some other process?

Downloads

Run it on your own machine

starry_night.py ~80 lines · numpy · matplotlib ↓ Download
great_wave.py ~90 lines · numpy · matplotlib ↓ Download
the_scream.py ~100 lines · numpy · matplotlib ↓ Download

Open on a desktop to download and run the scripts.

# one-time setup
pip install numpy matplotlib

# run any script
python starry_night.py
python great_wave.py
python the_scream.py

Challenge progression · four tiers per painting

The Starry Night

01 OBSERVERun it. Change one parameter. Predict what the output will look like before you press enter.
02 EXTENDAdd a second noise octave to the field. Measure whether the cost per frame is still O(n) — time it for n = 1000 and n = 10000.
03 REBUILDReplace the synthetic noise field with one sampled from a real wind dataset. ECMWF provides forecast data free via CDS API.
04 INVESTIGATEMake particles interact — each pushes its neighbours slightly. Watch O(n) become O(n²). Find the n at which your machine drops below 10 fps.

The Great Wave

01 OBSERVERun it. Add a wave source by clicking. Watch how two sources interfere — where do you see constructive vs destructive interference?
02 EXTENDImplement a reflective wall along one edge. How does changing the angle of incidence affect the reflected wave?
03 REBUILDFind the Courant limit experimentally: increase c until the simulation produces numerical artefacts, then prove the stability condition c·Δt/Δx ≤ 1/√2.
04 INVESTIGATEThe wave equation is linear — waves pass through each other unchanged. Add a non-linear term (η³ dispersion). What new phenomena appear?

The Scream

01 OBSERVERun it. Change F or k by 0.002. Does the resulting pattern change gradually or appear all at once? Record the parameter pairs that produce stripes vs spots vs nothing.
02 EXTENDMap the (F, k) parameter space: run a grid of 10×10 combinations and save a thumbnail for each. Which region of the space corresponds to which pattern type?
03 REBUILDImplement the FitzHugh-Nagumo model — a related two-variable reaction-diffusion system used in cardiac modelling. Do Turing patterns emerge with the standard cardiac parameters?
04 INVESTIGATECan you find (F, k) values where the system never reaches a stable pattern? Provide either a proof or a practical test for detecting non-convergence.

About the scripts

Each script is a self-contained implementation of the algorithm from the browser page, written for readability not performance. The code matches the pseudocode in section B — you should recognise every line.

matplotlib's FuncAnimation handles the update loop. The simulation logic is separate from the display — you can swap in a different renderer without touching the physics.

Requires Python 3.9+ · numpy ≥ 1.22 · matplotlib ≥ 3.5