loading...
running...
How to use
Pick a preset, or set masses, lengths, gravity and the starting state. The orange, red and blue trails belong to the three bobs, the cyan ghost is a twin started 10^-n rad away. Watch the twin-separation graph rise exponentially, then saturate. Space pauses, R resets. Drag to orbit, wheel to zoom, camera buttons for presets.
Models and baked PBR textures made in Blender; Three.js r160 animates every part from the shared physics module.
The science
Three point masses m1, m2, m3 hang from rigid, massless rods of length l1, l2, l3 in a vertical plane. With angles measured from the downward vertical, the kinetic and potential energies are
T = ½ Σ m_i |v_i|², v_1 = l1 w1 (cos th1, sin th1), v_2 = v_1 + l2 w2 (cos th2, sin th2), v_3 = v_2 + l3 w3 (cos th3, sin th3)
V = g (m1 y1 + m2 y2 + m3 y3), y_1 = -l1 cos th1, y_2 = y_1 - l2 cos th2, y_3 = y_2 - l3 cos th3
The Euler-Lagrange equations take the compact form A(th) α = f(th, w) with the symmetric mass matrix (M = m1+m2+m3, M23 = m2+m3)
A = [ M l1² , M23 l1 l2 cos(th1-th2) , m3 l1 l3 cos(th1-th3) ;
M23 l1 l2 cos(th1-th2) , M23 l2² , m3 l2 l3 cos(th2-th3) ;
m3 l1 l3 cos(th1-th3) , m3 l2 l3 cos(th2-th3) , m3 l3² ]
f_1 = -M23 l1 l2 sin(th1-th2) w2² - m3 l1 l3 sin(th1-th3) w3² - M g l1 sin th1
f_2 = M23 l1 l2 sin(th1-th2) w1² - m3 l2 l3 sin(th2-th3) w3² - M23 g l2 sin th2
f_3 = m3 l1 l3 sin(th1-th3) w1² + m3 l2 l3 sin(th2-th3) w2² - m3 g l3 sin th3
A 3x3 linear solve (Gaussian elimination with pivoting) gives the three angular accelerations at every integrator stage; an optional viscous joint term -b w is added to each. The state (th1, th2, th3, w1, w2, w3) is advanced with classical RK4 (or RK2 / Euler for comparison).
Small angles. Linearising about the hanging state gives the generalised eigenproblem det(K - w² M) = 0 with K = g diag(M l1, M23 l2, m3 l3). Three normal modes exist; for equal masses and lengths their shapes are about (0.43, 0.56, 0.71), (0.38, 0.13, -0.91) and (0.48, -0.79, 0.37) with periods 3.11 s, 1.32 s and 0.80 s at l = 1 m. The first check below solves that determinant independently and compares it with periods measured from the full nonlinear integrator.
Chaos. Away from the linear regime the system is non-integrable and nearby states separate roughly as d(t) = d0 e^(lambda t). The largest Lyapunov exponent lambda is estimated with the Benettin method (a shadow trajectory is renormalised to d0 every 0.05 s and the log growth is averaged). The cyan twin starts 10^-n rad away in th1; its separation saturates once it reaches the size of the accessible phase space, so the exponential law only holds on the rising part of the graph. Typical values for strongly chaotic starts are a few per second and depend on the energy; small-angle starts give lambda near 0 and a quasi-periodic Poincare section drawn as closed curves. A third link adds a dimension to phase space (6 instead of 4) and usually a larger lambda than the double pendulum, but there is no universal factor.
Numerics. RK4 with a 1/480 s step conserves total energy to better than 1e-5 relative over 20 s even on a chaotic orbit; explicit Euler gains energy steadily and RK2 drifts visibly, switch the integrator to see it. The time-reversal check integrates forward, flips the velocities and integrates back, recovering the start to about 1e-7; this proves integrator accuracy over that horizon, it is not a symmetry of the chaos (rounding errors would eventually be amplified).
Simplifications. Point masses, rigid massless rods, frictionless pivots unless damping is set, planar motion, no air drag, no rod mass or inertia. A real compound pendulum would need distributed-mass terms.