Visualizing Double Pendulum Chaos — RK4 Integration and Sensitivity to Initial Conditions

#Why a Double Pendulum

A simple pendulum is covered in high school physics. Its period is determined solely by length and gravity, and slightly changing the initial conditions barely changes the behavior. Predictable, and boringly obedient.

But the moment you attach another pendulum, the world changes. The double pendulum is a textbook example of deterministic chaos — the equations of motion are completely deterministic (zero randomness), yet shifting the initial angle by just 0.001 radians produces an entirely different trajectory within seconds.

Implementing the DoublePendulum component (362 lines) meant repeatedly confronting the core of chaotic dynamics. This article covers, from an implementer's perspective, Lagrangian mechanics, RK4 integration, NaN prevention, and visualization of sensitivity to initial conditions.

#From Lagrangian Mechanics to Equations of Motion

#Why Lagrangian Instead of Newtonian Mechanics

Newtonian mechanics (F = ma) works well for free particles, but becomes cumbersome for constrained systems. A pendulum has the constraint "the rod length is constant." Solving this with Newtonian mechanics introduces tension as an additional unknown, complicating the equations.

Lagrangian mechanics takes the approach of "deriving forces from energy." Define the Lagrangian L = T - V (the difference between kinetic energy T and potential energy V), substitute into the Euler-Lagrange equations, and constraints are handled automatically.

#State Vector

The state of a double pendulum is completely described by four variables:

• θ₁ (a1): angle of the upper pendulum
• θ₂ (a2): angle of the lower pendulum
• ω₁ (w1): angular velocity of the upper pendulum
• ω₂ (w2): angular velocity of the lower pendulum

interface State {
  a1: number  // θ₁
  a2: number  // θ₂
  w1: number  // ω₁ = dθ₁/dt
  w2: number  // ω₂ = dθ₂/dt
}

Position (x, y) is uniquely determined from angles and rod lengths, so it's not included in the state. This is the advantage of generalized coordinates.

#Deriving Angular Acceleration

The angular acceleration equations derived from the Lagrangian, shown in the form used in the actual code. Let Δa = θ₁ - θ₂, m = M₁ + M₂:

dω₁/dt:

dω₁/dt = [-g·m·sin(θ₁) - M₂·g·sin(θ₁ - 2θ₂)·0.5
          - sin(Δa)·M₂·(ω₂²·L₂ + ω₁²·L₁·cos(Δa))]
         / [L₁·(m - M₂·cos²(Δa))]

dω₂/dt:

dω₂/dt = [sin(Δa)·(ω₁²·L₁·m + g·m·cos(θ₁)
          + ω₂²·L₂·M₂·cos(Δa))]
         / [(L₂/L₁)·L₁·(m - M₂·cos²(Δa))]

The equations are long, but what they do is the rotational version of Newton's second law. The numerator corresponds to "total torque acting on the system," and the denominator to "effective moment of inertia."

In the implementation, these are packaged into a derivatives function:

function derivatives(s: State, g: number): State {
  const { a1, a2, w1, w2 } = s
  const da = a1 - a2
  const sinDa = Math.sin(da)
  const cosDa = Math.cos(da)
  const m = M1 + M2

  const dw1 =
    (-g * m * Math.sin(a1) -
      M2 * g * Math.sin(a1 - 2 * a2) * 0.5 -
      sinDa * M2 * (w2 * w2 * L2_RATIO + w1 * w1 * L1_RATIO * cosDa)) /
    (L1_RATIO * (m - M2 * cosDa * cosDa))

  const dw2 =
    (sinDa * (w1 * w1 * L1_RATIO * m + g * m * Math.cos(a1)
      + w2 * w2 * L2_RATIO * M2 * cosDa)) /
    ((L2_RATIO / L1_RATIO) * L1_RATIO * (m - M2 * cosDa * cosDa))

  return { a1: w1, a2: w2, w1: dw1, w2: dw2 }
}

Note the return value a1: w1, a2: w2. The time derivative of angle is angular velocity itself. This function returns the "time derivative of the state," which serves as input to the integrator.

#Implementing RK4 Integration

#Why Euler's Method Isn't Enough

The simplest numerical integration method is Euler's method:

x(t + dt) = x(t) + f(x, t) * dt

Looks reasonable, but Euler's method has only first-order accuracy. Halving the time step dt only halves the error. In chaotic systems, small errors are exponentially amplified, so Euler's method quickly leads to energy divergence. The pendulum keeps spinning or suddenly stops — physically impossible behavior.

#Runge-Kutta 4th Order (RK4)

RK4 achieves "fourth-order accuracy with four function evaluations." Halving dt reduces the error by a factor of 16. Excellent cost-effectiveness.

The idea is to evaluate the slope at the beginning, midpoint, and end of the time step, then take a weighted average:

1. k₁: Slope at the current state
2. k₂: Slope at the state advanced half a step by k₁
3. k₃: Slope at the state advanced half a step by k₂
4. k₄: Slope at the state advanced a full step by k₃
5. Final update: (k₁ + 2k₂ + 2k₃ + k₄) / 6

function rk4Step(s: State, g: number, dt: number): State {
  const k1 = derivatives(s, g)
  const s2: State = {
    a1: s.a1 + k1.a1 * dt * 0.5,
    a2: s.a2 + k1.a2 * dt * 0.5,
    w1: s.w1 + k1.w1 * dt * 0.5,
    w2: s.w2 + k1.w2 * dt * 0.5,
  }
  const k2 = derivatives(s2, g)
  const s3: State = {
    a1: s.a1 + k2.a1 * dt * 0.5,
    a2: s.a2 + k2.a2 * dt * 0.5,
    w1: s.w1 + k2.w1 * dt * 0.5,
    w2: s.w2 + k2.w2 * dt * 0.5,
  }
  const k3 = derivatives(s3, g)
  const s4: State = {
    a1: s.a1 + k3.a1 * dt,
    a2: s.a2 + k3.a2 * dt,
    w1: s.w1 + k3.w1 * dt,
    w2: s.w2 + k3.w2 * dt,
  }
  const k4 = derivatives(s4, g)

  return {
    a1: s.a1 + (k1.a1 + 2*k2.a1 + 2*k3.a1 + k4.a1) * (dt / 6),
    a2: s.a2 + (k1.a2 + 2*k2.a2 + 2*k3.a2 + k4.a2) * (dt / 6),
    w1: s.w1 + (k1.w1 + 2*k2.w1 + 2*k3.w1 + k4.w1) * (dt / 6),
    w2: s.w2 + (k1.w2 + 2*k2.w2 + 2*k3.w2 + k4.w2) * (dt / 6),
  }
}

#Substeps for Stability

requestAnimationFrame fires at roughly 60fps. At about 16.7ms per frame, this time step is too large for a physics simulation.

The implementation uses DT = 0.02 and STEPS_PER_FRAME = 8. Running 8 RK4 steps per frame keeps the effective time step sufficiently small:

const DT = 0.02
const STEPS_PER_FRAME = 8

// Inside the animation loop
for (let i = 0; i < STEPS_PER_FRAME; i++) {
  stateRef.current = rk4Step(stateRef.current, g, DT)
}

The choice of DT = 0.02 was empirical. At 0.05, energy drift is visually noticeable. At 0.01, it's stable but requires more STEPS_PER_FRAME, increasing CPU load. 0.02 x 8 = 0.16 effective step creates a "slightly faster flow of time" relative to the 60fps frame interval (~0.0167s), which also produces a visually satisfying pace.

#Fighting NaN

#The Problem

Floating-point arithmetic has finite precision. IEEE 754 64-bit floats have about 15 significant digits, but in chaotic systems, errors accumulate exponentially.

Specifically, when the denominator of angular acceleration m - M₂·cos²(Δa) approaches 0, the division result becomes enormous. If this continues for several steps, Infinity is born. Calculations involving Infinity produce NaN. NaN propagates through all operations, and Canvas rendering stops completely.

#The Fix

The implementation uses a simple Number.isFinite() guard:

const { a1, a2, w1, w2 } = stateRef.current

if (
  !Number.isFinite(a1) ||
  !Number.isFinite(a2) ||
  !Number.isFinite(w1) ||
  !Number.isFinite(w2)
) {
  const angles = randomAngles()
  stateRef.current = { ...angles, w1: 0, w2: 0 }
  trailRef.current = []
  return
}

When NaN/Infinity is detected, the state is reset with random initial angles and the trail buffer is cleared. To the user, it looks like "the pendulum restarted with new initial conditions."

Chaotic systems are numerically unstable too. Without this guard, leaving the simulation running for a few minutes with high gravity is enough to cause a crash. This kind of defensive programming is essential for physics simulations.

Number.isFinite() catches both NaN and ±Infinity. !Number.isNaN() alone would miss Infinity, making it insufficient.

#Visualizing Sensitivity to Initial Conditions

The defining property of chaos is sensitivity to initial conditions (also known as sensitive dependence on initial conditions). A positive Lyapunov exponent means nearby trajectories diverge exponentially.

#Speed-to-Color Mapping

The trail color corresponds to the second mass point's speed (magnitude of angular velocity):

function speedToHue(speed: number): number {
  const t = Math.min(1, speed / 15)
  return 240 - t * 240
}

• Low speed (speed ≈ 0) → hue = 240 (blue)
• High speed (speed ≥ 15) → hue = 0 (red)

The blue-to-red gradient intuitively maps to "cold → hot" and "slow → fast," making energy distribution visible just by looking at the trail.

#Trail Buffer

The trail is implemented as a circular buffer of TrailPoint[]:

interface TrailPoint {
  x: number
  y: number
  speed: number
}

const MAX_TRAIL = 2000

// Every frame
trailRef.current.push({ x: x2, y: y2, speed })
while (trailRef.current.length > MAX_TRAIL) {
  trailRef.current.shift()
}

2000 points represents about 33 seconds of history at 60fps. Making it longer increases memory and rendering cost, but making it shorter hides the beautiful chaotic patterns. 2000 is a balance point found experimentally.

During rendering, older points are made more transparent so the trail naturally fades out:

for (let i = 1; i < trail.length; i++) {
  const age = i / trail.length
  const alpha = age * 0.7 + 0.02
  const hue = speedToHue(trail[i].speed)

  ctx.beginPath()
  ctx.moveTo(trail[i - 1].x, trail[i - 1].y)
  ctx.lineTo(trail[i].x, trail[i].y)
  ctx.strokeStyle = `hsla(\${hue}, 80%, 60%, \${alpha})`
  ctx.stroke()
}

age is a normalized value from 0 (oldest) to 1 (newest), making newer points more opaque. Even the oldest points have alpha = 0.02, leaving a faint trace to prevent the trail from abruptly vanishing.

#"A Difference of 0.001 Radians"

In practice, shifting the initial angle by 0.001 radians (about 0.057 degrees) produces a completely different trajectory after 5-10 seconds. To the human eye, the initial conditions look identical, but the results are dramatically different. This is chaos.

In this implementation, pressing the Reset button generates initial angles with Math.random(), so getting the same trajectory twice is practically impossible. Every trajectory is unique — a once-in-a-lifetime encounter.

#Parameters and Feel

#Constant Design

const M1 = 10       // Mass of upper pendulum
const M2 = 10       // Mass of lower pendulum
const L1_RATIO = 0.28  // Rod length as viewport ratio
const L2_RATIO = 0.28
const MAX_TRAIL = 2000
const DT = 0.02
const STEPS_PER_FRAME = 8

Equal masses (M1 = M2 = 10) were chosen because chaos is most pronounced in this configuration. When the upper mass is extremely large, the lower pendulum becomes a mere appendage and chaos diminishes. Equal masses let both pendulums influence each other equally, producing the most complex motion.

L1_RATIO = L2_RATIO = 0.28 is the ratio relative to Math.min(width, height) of the viewport. 0.28 x 2 = 0.56, so even when both rods are fully extended, they fit within 56% of the screen. This was set as the maximum value that keeps the pendulum tip from going off-canvas.

#Gravity Slider

Gravity is adjustable from 0.5 to 3.0:

• g = 0.5: Slow, graceful motion. The trail draws large curves
• g = 1.5 (default): Scaled standard value of Earth's gravity
• g = 3.0: Intense motion. Angular velocity increases, and NaN risk rises

Changing gravity changes the energy scale of the system, which also changes how chaos "looks." Low gravity produces leisurely spirals; high gravity produces violent whip-like snaps.

#Pivot Position

const pivotX = w / 2
const pivotY = h * 0.3

Placing the pivot at 30% from the top maximizes the pendulum's range of motion. Centering it leaves only the bottom half usable; placing it too high makes the rods too short.

#What Implementation Taught Me

#Practical Value of Lagrangian Mechanics

Lagrangian mechanics is the approach of "deriving forces from energy." It automatically handles constraint conditions (constant rod length), eliminating the need to separately calculate tension as in Newtonian mechanics. Describing the system in generalized coordinates (angles) makes equation derivation a mechanical process.

Conversely, for systems with many degrees of freedom (like N-body problems), Newtonian mechanics can be simpler. Choose the tool to fit the problem.

#Debugging Chaotic Systems

Chaotic systems have no "correct trajectory." Two implementations with the same initial conditions producing different results doesn't tell you which is right (because minuscule rounding error differences are amplified).

Instead, verify using qualitative properties:

1. Energy conservation: Without friction, total energy (kinetic + potential) should be constant. RK4 still drifts slightly, but far less than Euler
2. Symmetry: Left-right mirrored initial conditions should produce left-right mirrored trajectories
3. Limiting cases: At small angles (θ ≈ 0), the motion should approximate a simple pendulum's periodic motion

The implementation omitted numerical energy monitoring, but during development, T + V values were tracked via console.log to confirm drift was within acceptable bounds.

#Cost-Effectiveness of RK4

RK4 achieves "fourth-order accuracy with four function evaluations." Symplectic integrators (like Verlet) excel at energy conservation, but are difficult to apply directly to non-separable Hamiltonians like the double pendulum. RK4 is general-purpose, easy to implement, and sufficiently accurate — the best balance for physics simulations in the browser.

The "4" in fourth-order means it accurately reproduces terms up to the fourth order of the Taylor expansion. The error is O(dt⁵), so at dt = 0.02, the error per step is approximately 3.2 x 10⁻⁹. Even accumulated over 8 steps per frame at 60fps, energy is practically conserved for several minutes.

#Canvas API Drawing Techniques

For trail rendering, beginPath() → stroke() is called per segment. This looks inefficient, but since each segment has a different color (hue) and transparency (alpha), they cannot be combined into a single Path.

If performance becomes an issue, the trail could be replaced with WebGL line strips or drawn to an offscreen Canvas and cached, but for ~2000 segments, Canvas 2D is fast enough.

#Summary: Implementing Double Pendulum Chaos

The double pendulum lets you experience deterministic chaos in just 362 lines. Derive equations of motion with Lagrangian mechanics, integrate with RK4, ensure numerical stability with NaN guards, and visualize sensitivity to initial conditions with speed-to-color mapping.

Chaos is terrifying in its "unpredictability," but behind it lies a completely deterministic equation. The equations are beautiful, and the trajectories are beautiful. You just can't see that beauty the same way twice.