How it works

Take a single pendulum, a weight on a string. You can tell exactly where it will be tomorrow, in a year, any time you want. A physicist can write a formula and calculate it.

Now connect a second pendulum to the bottom of the first. Start two copies from positions that are nearly identical, a millimetre apart, say. For the first second or two they look the same. Then one dips a little earlier, the other swings wide, and within maybe ten seconds they are doing completely different things. After a minute, the starting positions are completely irrelevant.

This is not randomness. Both pendulums are following Newton's laws exactly. The issue is that the system is extremely sensitive to starting conditions. Any tiny difference grows. Double it by a small factor, wait, it grows again. That compounding amplification is what physicists call chaos, and it only shows up at large angles where the coupling between the two rods is strongly nonlinear. For small starting angles the system actually behaves predictably, like two linked springs.

Positions (direct trigonometry from pivot):

$$x_1 = L_1 \sin\theta_1 \qquad y_1 = -L_1 \cos\theta_1$$ $$x_2 = x_1 + L_2 \sin\theta_2 \qquad y_2 = y_1 - L_2 \cos\theta_2$$

Velocities, differentiating with respect to time:

$$\dot{x}_1 = \dot\theta_1 L_1 \cos\theta_1 \qquad \dot{y}_1 = \dot\theta_1 L_1 \sin\theta_1$$ $$\dot{x}_2 = \dot{x}_1 + \dot\theta_2 L_2 \cos\theta_2 \qquad \dot{y}_2 = \dot{y}_1 + \dot\theta_2 L_2 \sin\theta_2$$

Accelerations, differentiating again (product rule):

$$\ddot{x}_1 = -\dot\theta_1^2 L_1 \sin\theta_1 + \ddot\theta_1 L_1 \cos\theta_1$$ $$\ddot{y}_1 = \phantom{-}\dot\theta_1^2 L_1 \cos\theta_1 + \ddot\theta_1 L_1 \sin\theta_1$$ $$\ddot{x}_2 = \ddot{x}_1 - \dot\theta_2^2 L_2 \sin\theta_2 + \ddot\theta_2 L_2 \cos\theta_2$$ $$\ddot{y}_2 = \ddot{y}_1 + \dot\theta_2^2 L_2 \cos\theta_2 + \ddot\theta_2 L_2 \sin\theta_2$$

Applying Newton's second law to the geometry above gives the angular accelerations. Let $\Delta = \theta_1 - \theta_2$ and $D = 2m_1 + m_2 - m_2\cos(2\Delta)$:

$$\ddot\theta_1 = \frac{-g(2m_1+m_2)\sin\theta_1 - m_2 g\sin(\theta_1-2\theta_2) - 2\sin\Delta\,m_2(\dot\theta_2^2 L_2 + \dot\theta_1^2 L_1\cos\Delta)}{L_1 D}$$ $$\ddot\theta_2 = \frac{2\sin\Delta\,\bigl[\dot\theta_1^2 L_1(m_1+m_2) + g(m_1+m_2)\cos\theta_1 + \dot\theta_2^2 L_2 m_2\cos\Delta\bigr]}{L_2 D}$$

No small-angle approximations are made. These are exact for any angle.

The equations of motion give $\ddot\theta_1$ and $\ddot\theta_2$. To integrate numerically, introduce $\omega_1 = \dot\theta_1$ and $\omega_2 = \dot\theta_2$ so everything is first-order:

$$\mathbf{s} = (\theta_1,\, \omega_1,\, \theta_2,\, \omega_2) \qquad \dot{\mathbf{s}} = f(\mathbf{s}) = (\omega_1,\; \ddot\theta_1,\; \omega_2,\; \ddot\theta_2)$$

Each time step advances $\mathbf{s}$ using Runge-Kutta 4, which evaluates $f$ four times per step:

$$k_1 = f(\mathbf{s}_n)$$ $$k_2 = f\!\left(\mathbf{s}_n + \tfrac{h}{2}k_1\right)$$ $$k_3 = f\!\left(\mathbf{s}_n + \tfrac{h}{2}k_2\right)$$ $$k_4 = f(\mathbf{s}_n + h\,k_3)$$ $$\mathbf{s}_{n+1} = \mathbf{s}_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$

With $h = \Delta t = 0.005\,\text{s}$:

In practice, total energy drifts by less than $10^{-8}$ J over hundreds of simulated seconds.

This simulator runs physics in two separate places. The Pendulum page runs on the CPU in JavaScript, which uses IEEE 754 64-bit floats (f64) by default. The Phase Map runs on the GPU in a WGSL compute shader, which uses 32-bit floats (f32), the native type for GPU arithmetic. f64 on GPU hardware is rare and typically 16–32× slower, so f32 is the practical choice.

CPU (f64): Machine $\varepsilon \approx 10^{-16}$ is four orders of magnitude below the RK4 per-step truncation error $\approx 10^{-12}$. Rounding is negligible; integration error dominates.

GPU (f32): Machine $\varepsilon \approx 10^{-7}$ is five orders of magnitude above the f64 truncation error. Each arithmetic operation rounds at the 7th decimal digit. The Phase Map shader runs roughly 50 floating-point operations per RK4 step, giving round-off noise per step of $\approx 50 \times 10^{-7} = 5 \times 10^{-6}$. Accumulation over time:

At 60,000 steps individual pixel trajectories can be off by $\sim 0.3\,\text{rad}$, no longer physically exact. This is acceptable because the Phase Map visualises relative behaviour: which starting conditions flip quickly versus slowly. The fractal boundary survives f32 precision; absolute flip-time values should not be read as precise beyond ~2 significant figures for long runs.

The Phase Map simulates an $N \times N$ grid of pendulums in parallel. Each pixel is one independent starting condition $(\theta_1, \theta_2)$ with $\omega_1 = \omega_2 = 0$.

All state lives in a single GPU storage buffer. Each pendulum occupies 8 × f32 (32 bytes):

One WGSL compute thread per pendulum, workgroup size 256. At $800 \times 800$ (640,000 pendulums) that is 2,500 workgroups per dispatch, all running in parallel on the GPU. Each thread runs its own RK4 loop for $N$ steps and writes results back in-place. No inter-thread communication is needed.

Flip detection runs inside the per-step loop. A "flip" occurs when the lower bob passes over the top, detected when $\sin(\theta_2)$ changes sign while $\cos(\theta_2) < 0$:

$$\sin(\theta_2^{\text{prev}}) \cdot \sin(\theta_2^{\text{new}}) < 0 \;\text{ and }\; \cos(\theta_2^{\text{new}}) < 0$$

firstFlipTime is recorded once and never overwritten. A separate reinit shader resets the buffer on resolution or region change, avoiding a 20 MB CPU→GPU upload.

Coloring modes:

Time to first flip: each pixel's color encodes firstFlipTime. The boundary between fast-flip and slow-flip regions is fractal: a one-pixel shift in initial angle can double the flip time.

Live $\theta_2$: pixels colored by the current lower angle through a palette. Colors evolve continuously as the simulation advances, revealing the flow structure of phase space.