Skip to content
Computational Lab.

Double Pendulum

Chaos in two joints.

Intermediate / integration / chaos

Two rigid links under gravity obey equations with no closed-form solution, so the motion comes from numerical integration. Run two pendulums from initial angles that differ by a thousandth of a radian and watch them diverge completely: deterministic, predictable in principle, practically unpredictable.

Tagged: chaos, RK4, Lagrangian, sensitivity

Live simulation

The ghosted pendulum started 0.000001 radians away from the primary.

Speed

The ghosted pendulum started 0.000001 radians away from the primary.

Mathematics

The double pendulum's equations of motion come from the Lagrangian and have no closed-form solution. They are also chaotic: two initial states separated by one millionth of a radian eventually look nothing alike.

θ̈₁ = [ −g(2m₁+m₂)sinθ₁ − m₂g sin(θ₁−2θ₂) − 2sin(θ₁−θ₂) m₂ (ω₂²l₂ + ω₁²l₁cos(θ₁−θ₂)) ] / [ l₁ (2m₁+m₂−m₂cos(2θ₁−2θ₂)) ]

θ̈₂ = [ 2sin(θ₁−θ₂)(ω₁²l₁(m₁+m₂) + g(m₁+m₂)cosθ₁ + ω₂²l₂m₂cos(θ₁−θ₂)) ] / [ l₂ (2m₁+m₂−m₂cos(2θ₁−2θ₂)) ]

E = ½m₁l₁²ω₁² + ½m₂(l₁²ω₁² + l₂²ω₂² + 2l₁l₂ω₁ω₂cos(θ₁−θ₂)) − (m₁+m₂) g l₁ cosθ₁ − m₂ g l₂ cosθ₂

Energy must be conserved: a growing drift would mean the integrator, not the physics, is changing the system. The ghosted pendulum is a twin integrated from the same state plus 10⁻⁶ radians; the twin-gap metric is the divergence you have heard described but probably never watched.

How it works

Classical RK4 integrates the exact equations each tick, with the energy recomputed directly from the state so drift is measured rather than assumed. The twin uses a cheaper semi-implicit step; its own integration error is negligible next to the separation it exists to measure. The trace records the lower bob's path in a ring buffer.

Set both initial angles to 90 degrees, run for a while, then reset and add 0.001 degrees. The metrics panel shows when the twins part ways.

Implementation

src/engine/experiments/double-pendulum.ts · 264 lines · the module running above

TypeScript · fixed-timestep engine

1import { BaseExperiment } from "../core/base-experiment";
2import type { ParameterDef, ThemeColors } from "../core/types";
3
4const TRACE_LENGTH = 260;
5const PX_PER_M = 130;
6
7/**
8 * Double pendulum integrated on the exact Lagrangian equations with RK4.
9 * A "shadow" twin started 1e-6 rad away makes chaos measurable: the
10 * separation metric is the real Lyapunov-style divergence of the two.
11 */
12export class DoublePendulumExperiment extends BaseExperiment {
13 readonly id = "double-pendulum";
14
15 private m1 = 1;
16 private m2 = 1;
17 private l1 = 1;
18 private l2 = 0.8;
19 private g = 9.81;
20 private th1 = 0;
21 private th2 = 0;
22 private w1 = 0;
23 private w2 = 0;
24 private s1 = 0;
25 private s2 = 0;
26 private sw1 = 0;
27 private sw2 = 0; // shadow twin
28 private time = 0;
29 private energy0 = 0;
30 private traceX = new Float32Array(TRACE_LENGTH);
31 private traceY = new Float32Array(TRACE_LENGTH);
32 private traceHead = 0;
33
34 protected params(): ParameterDef[] {
35 return [
36 { key: "m1", label: "Upper mass", min: 0.5, max: 5, step: 0.1, defaultValue: 1, unit: "kg" },
37 { key: "m2", label: "Lower mass", min: 0.5, max: 5, step: 0.1, defaultValue: 1, unit: "kg" },
38 { key: "l1", label: "Upper arm", min: 40, max: 160, step: 5, defaultValue: 100, unit: "cm" },
39 { key: "l2", label: "Lower arm", min: 40, max: 160, step: 5, defaultValue: 80, unit: "cm" },
40 { key: "g", label: "Gravity", min: 1, max: 25, step: 0.1, defaultValue: 9.8, unit: "m/s²" },
41 { key: "a1", label: "Initial θ1", min: -180, max: 180, step: 1, defaultValue: 120, unit: "°" },
42 { key: "a2", label: "Initial θ2", min: -180, max: 180, step: 1, defaultValue: -20, unit: "°" },
43 {
44 key: "integrator",
45 label: "Integrator",
46 options: [
47 { value: "rk4", label: "RK4" },
48 { value: "semi", label: "Semi-implicit Euler" },
49 ],
50 defaultValue: "rk4",
51 },
52 ];
53 }
54
55 protected onReset(): void {
56 this.m1 = this.num("m1");
57 this.m2 = this.num("m2");
58 this.l1 = this.num("l1") / 100;
59 this.l2 = this.num("l2") / 100;
60 this.g = this.num("g");
61 this.th1 = (this.num("a1") * Math.PI) / 180;
62 this.th2 = (this.num("a2") * Math.PI) / 180;
63 this.w1 = 0;
64 this.w2 = 0;
65 this.s1 = this.th1 + 1e-6;
66 this.s2 = this.th2;
67 this.sw1 = 0;
68 this.sw2 = 0;
69 this.time = 0;
70 this.traceX.fill(0);
71 this.traceY.fill(0);
72 this.traceHead = 0;
73 this.energy0 = this.energy(this.th1, this.th2, this.w1, this.w2);
74 }
75
76 /** Exact equations of motion for the planar double pendulum. */
77 private accelerations(t1: number, t2: number, w1: number, w2: number): [number, number] {
78 const { m1, m2, l1, l2, g } = this;
79 const d = t1 - t2;
80 const denom = 2 * m1 + m2 - m2 * Math.cos(2 * d);
81 const a1 =
82 (-g * (2 * m1 + m2) * Math.sin(t1) -
83 m2 * g * Math.sin(t1 - 2 * t2) -
84 2 * Math.sin(d) * m2 * (w2 * w2 * l2 + w1 * w1 * l1 * Math.cos(d))) /
85 (l1 * denom);
86 const a2 =
87 (2 *
88 Math.sin(d) *
89 (w1 * w1 * l1 * (m1 + m2) + g * (m1 + m2) * Math.cos(t1) + w2 * w2 * l2 * m2 * Math.cos(d))) /
90 (l2 * denom);
91 return [a1, a2];
92 }
93
94 private energy(t1: number, t2: number, w1: number, w2: number): number {
95 const { m1, m2, l1, l2, g } = this;
96 const kinetic =
97 0.5 * m1 * l1 * l1 * w1 * w1 +
98 0.5 * m2 * (l1 * l1 * w1 * w1 + l2 * l2 * w2 * w2 + 2 * l1 * l2 * w1 * w2 * Math.cos(t1 - t2));
99 const potential = -(m1 + m2) * g * l1 * Math.cos(t1) - m2 * g * l2 * Math.cos(t2);
100 return kinetic + potential;
101 }
102
103 protected onUpdate(dt: number): void {
104 if (this.str("integrator") === "rk4") {
105 const [a1, a2] = this.accelerations(this.th1, this.th2, this.w1, this.w2);
106 const half = dt / 2;
107 const [b1, b2] = this.accelerations(
108 this.th1 + this.w1 * half,
109 this.th2 + this.w2 * half,
110 this.w1 + a1 * half,
111 this.w2 + a2 * half,
112 );
113 const [c1, c2] = this.accelerations(
114 this.th1 + this.w1 * half + a1 * half * half,
115 this.th2 + this.w2 * half + a2 * half * half,
116 this.w1 + b1 * half,
117 this.w2 + b2 * half,
118 );
119 const [d1, d2] = this.accelerations(
120 this.th1 + this.w1 * dt + b1 * half * dt,
121 this.th2 + this.w2 * dt + b2 * half * dt,
122 this.w1 + c1 * dt,
123 this.w2 + c2 * dt,
124 );
125 const nw1 = this.w1 + (dt / 6) * (a1 + 2 * b1 + 2 * c1 + d1);
126 const nw2 = this.w2 + (dt / 6) * (a2 + 2 * b2 + 2 * c2 + d2);
127 const nt1 = this.th1 + (dt / 6) * (this.w1 + 2 * (this.w1 + a1 * half) + 2 * (this.w1 + b1 * half) + (this.w1 + c1 * dt));
128 const nt2 = this.th2 + (dt / 6) * (this.w2 + 2 * (this.w2 + a2 * half) + 2 * (this.w2 + b2 * half) + (this.w2 + c2 * dt));
129 this.th1 = nt1;
130 this.th2 = nt2;
131 this.w1 = nw1;
132 this.w2 = nw2;
133 // Shadow twin uses the cheaper semi-implicit step; its phase error is
134 // negligible next to the 1e-6 rad separation it exists to measure.
135 const [sa1, sa2] = this.accelerations(this.s1, this.s2, this.sw1, this.sw2);
136 this.sw1 += sa1 * dt;
137 this.sw2 += sa2 * dt;
138 this.s1 += this.sw1 * dt;
139 this.s2 += this.sw2 * dt;
140 } else {
141 const [a1, a2] = this.accelerations(this.th1, this.th2, this.w1, this.w2);
142 this.w1 += a1 * dt;
143 this.w2 += a2 * dt;
144 this.th1 += this.w1 * dt;
145 this.th2 += this.w2 * dt;
146 }
147 this.time += dt;
148 const bob = this.bobPosition(this.th1, this.th2);
149 this.traceX[this.traceHead] = bob.x;
150 this.traceY[this.traceHead] = bob.y;
151 this.traceHead = (this.traceHead + 1) % TRACE_LENGTH;
152 }
153
154 private pivot(): { x: number; y: number } {
155 return { x: this.width / 2, y: this.height * 0.34 };
156 }
157
158 private bobPosition(t1: number, t2: number): { x: number; y: number } {
159 const origin = this.pivot();
160 const jointX = origin.x + Math.sin(t1) * this.l1 * PX_PER_M;
161 const jointY = origin.y + Math.cos(t1) * this.l1 * PX_PER_M;
162 return {
163 x: jointX + Math.sin(t2) * this.l2 * PX_PER_M,
164 y: jointY + Math.cos(t2) * this.l2 * PX_PER_M,
165 };
166 }
167
168 render(ctx: CanvasRenderingContext2D, theme: ThemeColors, alpha: number): void {
169 ctx.fillStyle = theme.bg;
170 ctx.fillRect(0, 0, this.width, this.height);
171 const origin = this.pivot();
172
173 // Trace of the lower bob (most recent segment strongest).
174 ctx.lineWidth = 1;
175 for (let i = 1; i < TRACE_LENGTH; i++) {
176 const a = ((this.traceHead - i + TRACE_LENGTH) % TRACE_LENGTH);
177 const b = ((this.traceHead - i + 1 + TRACE_LENGTH) % TRACE_LENGTH);
178 if (this.traceX[a] === 0 && this.traceY[a] === 0) continue;
179 ctx.strokeStyle = theme.viz[1];
180 ctx.globalAlpha = (i / TRACE_LENGTH) * 0.6;
181 ctx.beginPath();
182 ctx.moveTo(this.traceX[a], this.traceY[a]);
183 ctx.lineTo(this.traceX[b], this.traceY[b]);
184 ctx.stroke();
185 }
186 ctx.globalAlpha = 1;
187
188 const inter1 = this.th1 - this.w1 * 0.016 * (1 - alpha);
189 const inter2 = this.th2 - this.w2 * 0.016 * (1 - alpha);
190 const jointX = origin.x + Math.sin(inter1) * this.l1 * PX_PER_M;
191 const jointY = origin.y + Math.cos(inter1) * this.l1 * PX_PER_M;
192 const bobX = jointX + Math.sin(inter2) * this.l2 * PX_PER_M;
193 const bobY = jointY + Math.cos(inter2) * this.l2 * PX_PER_M;
194
195 // Shadow twin, ghosted.
196 const shadow = this.bobPosition(this.s1, this.s2);
197 const jointS = {
198 x: origin.x + Math.sin(this.s1) * this.l1 * PX_PER_M,
199 y: origin.y + Math.cos(this.s1) * this.l1 * PX_PER_M,
200 };
201 ctx.strokeStyle = theme.fgTertiary;
202 ctx.globalAlpha = 0.5;
203 ctx.beginPath();
204 ctx.moveTo(origin.x, origin.y);
205 ctx.lineTo(jointS.x, jointS.y);
206 ctx.lineTo(shadow.x, shadow.y);
207 ctx.stroke();
208 ctx.fillStyle = theme.fgTertiary;
209 ctx.beginPath();
210 ctx.arc(shadow.x, shadow.y, 4, 0, Math.PI * 2);
211 ctx.fill();
212 ctx.globalAlpha = 1;
213
214 // Main pendulum
215 ctx.strokeStyle = theme.fgSecondary;
216 ctx.lineWidth = 2;
217 ctx.beginPath();
218 ctx.moveTo(origin.x, origin.y);
219 ctx.lineTo(jointX, jointY);
220 ctx.lineTo(bobX, bobY);
221 ctx.stroke();
222 ctx.fillStyle = theme.fgTertiary;
223 ctx.beginPath();
224 ctx.arc(origin.x, origin.y, 3, 0, Math.PI * 2);
225 ctx.fill();
226 ctx.fillStyle = theme.fg;
227 const r1 = 6 + this.m1 * 1.6;
228 ctx.beginPath();
229 ctx.arc(jointX, jointY, r1, 0, Math.PI * 2);
230 ctx.fill();
231 ctx.fillStyle = theme.accent;
232 const r2 = 6 + this.m2 * 1.6;
233 ctx.beginPath();
234 ctx.arc(bobX, bobY, r2, 0, Math.PI * 2);
235 ctx.fill();
236 }
237
238 getMetrics() {
239 const energy = this.energy(this.th1, this.th2, this.w1, this.w2);
240 const drift = this.energy0 !== 0 ? ((energy - this.energy0) / Math.abs(this.energy0)) * 100 : 0;
241 let sep = Math.abs(this.s1 - this.th1);
242 sep = (sep * 180) / Math.PI;
243 return {
244 theta1: `${((((this.th1 * 180) / Math.PI) % 360) + 540) % 360 - 180 | 0}°`,
245 theta2: `${((((this.th2 * 180) / Math.PI) % 360) + 540) % 360 - 180 | 0}°`,
246 energyDrift: `${drift >= 0 ? "+" : ""}${drift.toExponential(2)} %`,
247 twinGap: `${sep < 0.01 ? sep.toExponential(1) : sep.toFixed(1)}°`,
248 time: `${this.time.toFixed(1)} s`,
249 };
250 }
251
252 describe(): string {
253 const sep = (Math.abs(this.s1 - this.th1) * 180) / Math.PI;
254 if (sep < 0.01) {
255 return "Two pendulum bobs swinging; the 0.000001-radian twin is still indistinguishable from the primary.";
256 }
257 return `Two pendulum bobs swinging; the twin started 0.000001 rad away has diverged to ${sep.toFixed(1)} degrees.`;
258 }
259
260 entities(): number {
261 return 3; // primary + shadow twin (2 bodies each) as the tracked set
262 }
263}
264

Performance

Time per tick
O(1) per tick: a handful of trig calls, four RK4 stages
Space
A 260-sample trace ring buffer
Notes
Chaos is not randomness: same seed, same trajectory, every time. The twin differs by one millionth of a radian, not by chance.

Live numbers (FPS, frame time, steps per second, microseconds per step, entity count) are in the performance overlay on the simulation above: the pulse icon in its top-right corner.