Double Pendulum Simulation

Deterministic, yet unpredictable

Nothing random happens in this simulation. The equations are exact Newtonian mechanics, the numerical method is fixed, and the same parameters always produce the same run, in the browser and in the Python port. What makes the double pendulum hard to predict is sensitive dependence on initial conditions: a tiny difference in the starting state grows roughly exponentially, like e^λt, until it is as large as the motion itself.

• The separation plot is on a logarithmic scale, so exponential growth is a straight line. From δ = 10⁻⁹ rad it climbs at a nearly steady slope, with wiggles, for many seconds, then saturates: two angles cannot be more than about 4.4 rad apart (π√2).
• The angle plot shows why an exponential is so unforgiving. The curves are indistinguishable for most of the run; the difference only becomes visible in the last few doublings.
• Precision buys little time. With the default settings, making δ a thousand times smaller (10⁻¹²) delays the divergence only from about 18.7 s to 25.8 s. Knowing the initial state better adds time to the forecast only in proportion to the logarithm of the improvement.

Not every starting point is chaotic. Release both arms at 10° and the motion is a regular combination of two normal modes: the separation stays near 10⁻⁹ rad for the whole run and the fitted growth rate is essentially zero. Chaos needs enough energy for the arms to exchange it irregularly; at 120° and −20° (the default) the lower arm flips over the top repeatedly.

The equations of motion

Two point masses m₁, m₂ hang on massless rigid rods of lengths l₁, l₂. Measuring both angles from the downward vertical and writing Δ = θ₁ − θ₂, the Lagrangian gives two coupled equations that are linear in the angular accelerations:

The code solves this 2 × 2 system for θ̈₁ and θ̈₂ at every evaluation; the determinant l₁l₂(m₁ + m₂ sin²Δ) is always positive. The total energy

is conserved exactly by these equations, which makes it the accuracy check.

Numerical method and the growth rate

Both pendulums are integrated with the classical fourth-order Runge–Kutta method (RK4) in 1 ms steps. When the arms whip around fast (a light upper mass yanked by a heavy lower one, say) a step is split into as many as 50 equal sub-steps, chosen from the current angular speeds and accelerations; both pendulums always take the same sub-steps, so they differ only through δ. The reported energy drift is the largest |E − E₀| over the run, divided by the energy scale (m₁ + m₂)gl₁ + m₂gl₂. A chaotic system amplifies integration error just as it amplifies δ, so no integrator can follow the true trajectory forever; a small drift shows that the integrator respects the physics even where it cannot follow one exact trajectory.

The separation is the distance between the two pendulums in angle space, √(Δθ₁² + Δθ₂²), with each difference wrapped into [−π, π). The growth rate λ is the least-squares slope of ln(separation) against time while the separation is still below 10⁻³ rad, sampled every 10 ms; it needs at least a second of such data. It is a finite-time estimate in the spirit of the largest Lyapunov exponent, which is defined in the limit of infinitely long times and infinitesimal separations, so expect it to vary with the starting angles and δ.

Bit-for-bit reproducibility needs one more precaution. Chaos amplifies even a difference in the last bit of sin(x), and JavaScript and Python compute sine with different library code. The model therefore uses its own sine and cosine, built only from addition, multiplication, division and rounding down, so both languages produce identical numbers.

Example result

Released from rest at θ₁ = 120° and θ₂ = −20° with a shadow pendulum offset by δ = 10⁻⁹ rad, the two pendulums stay together until the separation exceeds 1 rad at t = 18.68 s. While it was small, the separation grew at λ ≈ 0.99 per second, doubling every 0.70 s. Over the 30 s run the lower arm flips over the top 11 times and the energy drifts by only 8.2e-10 of the energy scale; the final separation is 1.68e+0 rad.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator. The excerpts show the core of each.

// State [θ₁, θ₂, ω₁, ω₂], angles from the downward vertical. The Lagrangian equations of

// motion are linear in the angular accelerations, a 2 × 2 system solved by Cramer's rule:

//   (m₁+m₂)·l₁·α₁ + m₂·l₂·cos Δ·α₂ = −m₂·l₂·ω₂²·sin Δ − (m₁+m₂)·g·sin θ₁

//   l₁·cos Δ·α₁   + l₂·α₂          =  l₁·ω₁²·sin Δ − g·sin θ₂,      Δ = θ₁ − θ₂

export function derivatives(s, p) {
  const [t1, t2, w1, w2] = s;
  const d = t1 - t2, sd = sin(d), cd = cos(d), M = p.m1 + p.m2;
  const a11 = M * p.l1, a12 = p.m2 * p.l2 * cd, a21 = p.l1 * cd, a22 = p.l2;
  const b1 = -p.m2 * p.l2 * w2 * w2 * sd - M * p.g * sin(t1);
  const b2 = p.l1 * w1 * w1 * sd - p.g * sin(t2);
  const det = a11 * a22 - a12 * a21;
  return [w1, w2, (b1 * a22 - a12 * b2) / det, (a11 * b2 - a21 * b1) / det];
}

export function rk4Step(s, h, p) {
  const k1 = derivatives(s, p);
  const k2 = derivatives(s.map((y, i) => y + (h / 2) * k1[i]), p);
  const k3 = derivatives(s.map((y, i) => y + (h / 2) * k2[i]), p);
  const k4 = derivatives(s.map((y, i) => y + h * k3[i]), p);
  return s.map((y, i) => y + (h / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]));
}

export function energy(s, p) {
  const [t1, t2, w1, w2] = s, M = p.m1 + p.m2;
  const kinetic = 0.5 * M * p.l1 * p.l1 * w1 * w1 + 0.5 * p.m2 * p.l2 * p.l2 * w2 * w2 + p.m2 * p.l1 * p.l2 * w1 * w2 * cos(t1 - t2);
  return kinetic - M * p.g * p.l1 * cos(t1) - p.m2 * p.g * p.l2 * cos(t2);
}

def derivatives(s, p):
    """State (theta1, theta2, omega1, omega2), angles from the downward vertical.
    The Lagrangian equations are linear in the angular accelerations: solve the 2x2 system."""
    t1, t2, w1, w2 = s
    d = t1 - t2
    sd, cd, M = sin(d), cos(d), p["m1"] + p["m2"]
    a11, a12, a21, a22 = M * p["l1"], p["m2"] * p["l2"] * cd, p["l1"] * cd, p["l2"]
    b1 = -p["m2"] * p["l2"] * w2 * w2 * sd - M * p["g"] * sin(t1)
    b2 = p["l1"] * w1 * w1 * sd - p["g"] * sin(t2)
    det = a11 * a22 - a12 * a21
    return [w1, w2, (b1 * a22 - a12 * b2) / det, (a11 * b2 - a21 * b1) / det]


def rk4_step(s, h, p):
    k1 = derivatives(s, p)
    k2 = derivatives([y + h / 2 * k1[i] for i, y in enumerate(s)], p)
    k3 = derivatives([y + h / 2 * k2[i] for i, y in enumerate(s)], p)
    k4 = derivatives([y + h * k3[i] for i, y in enumerate(s)], p)
    return [y + h / 6 * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) for i, y in enumerate(s)]


def energy(s, p):
    t1, t2, w1, w2 = s
    M = p["m1"] + p["m2"]
    kinetic = 0.5 * M * p["l1"] * p["l1"] * w1 * w1 + 0.5 * p["m2"] * p["l2"] * p["l2"] * w2 * w2 + p["m2"] * p["l1"] * p["l2"] * w1 * w2 * cos(t1 - t2)
    return kinetic - M * p["g"] * p["l1"] * cos(t1) - p["m2"] * p["g"] * p["l2"] * cos(t2)

Run the full model from a terminal:

node simulation.js --theta1 120 --theta2 -20
python simulation.py --theta1 120 --theta2 -20

Download the model

Assumptions & limitations

• Ideal pendulum. Point masses on rigid, massless rods with frictionless pivots and no air resistance, so energy is conserved. A real double pendulum loses energy and eventually settles.
• Planar motion. Both arms swing in one vertical plane.
• Finite precision. The integrator's own error is amplified exactly like δ. Beyond the divergence time the run shows typical chaotic motion, not the exact trajectory of these initial conditions.
• A finite-time growth rate. λ is a least-squares fit over a finite window, not the true largest Lyapunov exponent.

Related simulations

• OS-035Harmonic oscillatorA mass on a spring with friction and a periodic push: damping regimes, phase portraits and resonance.
• OS-034Projectile motionThe vacuum parabola against a trajectory with quadratic air drag, and why 45° stops being the best angle.
• OS-030Lotka–VolterraTwo coupled equations whose predator and prey populations circle forever on closed orbits.
• OS-012Random walkCoin-flip steps on a line or grid: distance grows like √n, and in 1-D and 2-D every walker comes back.

Background: Euler vs Runge–Kutta: Time Stepping in Simulations

References

1. Shinbrot, T., Grebogi, C., Wisdom, J. & Yorke, J. A. (1992). Chaos in a double pendulum. American Journal of Physics, 60(6), 491–499. doi.org/10.1119/1.16860
2. Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos, 2nd ed. Westview Press.
3. Taylor, J. R. (2005). Classical Mechanics. University Science Books.