Double pendulum · solar-pro4

The user wants me to build a planar double-pendulum simulation. Let me think about the physics and implementation:

**Physics of Double Pendulum:**
- Two masses m1, m2 connected by massless rods of length L1, L2
- Generalized coordinates: θ1, θ2 (angles from downward vertical)
- Using Lagrangian mechanics to derive equations of motion

For equal masses (m1=m2=m) and equal lengths (L1=L2=L), the Lagrangian is:
L = T - V

T = (1/2)m(L²θ1'² + L²θ1'θ2'cos(θ1-θ2) + L²θ2'²)
V = -mgL(2cos(θ1) + cos(θ2))

The equations of motion (using Euler-Lagrange) for equal masses/lengths become:

θ1'' = (-g(2sin(θ1) - sin(θ2)cos(θ1-θ2)) - (θ2'²L sin(θ1-θ2)))/(2L - L cos²(θ1-θ2)...)

Actually let me derive this more carefully.

For the double pendulum with m1=m2=m, L1=L2=L:

The equations are (using sympy or known results):

θ1'' = [ -g (2 sin θ1 - sin θ2 cos(θ1-θ2)) - L θ2'² sin(θ1-θ2) ] / [ L (2 - cos²(θ1-θ2)) ]

θ2'' = [ 2g (sin θ1 cos(θ1-θ2) - sin θ2) + 2L θ1'² sin(θ1-θ2) ] / [ L (2 - cos²(θ1-θ2)) ]

Wait, let me be more careful. The standard form for the double pendulum:

With m1=m2, L1=L2=L, and using the gravitational acceleration g:

Let Δθ = θ2 - θ1

Then:

α1 = θ1'', α2 = θ2''

Denom = 2 - cos²(Δθ) = 1 + sin²(Δθ)

α1 = [ -g (2 sin θ1 - sin(θ2) cos(Δθ)) - L θ2'² sin(Δθ) ] / [ L (1 + sin²(Δθ)) ]

α2 = [ 2g (sin θ1 cos(Δθ) - sin θ2) + 2L θ1'² sin(Δθ) ] / [ L (1 + sin²(Δθ)) ]

Hmm, let me verify this. From various sources, the double pendulum equations for equal masses and lengths:

For m1=m2=m, L1=L2=L:

```
den = 2 - cos(θ1-θ2)^2

θ1'' = (-g*(2*sin(θ1) - sin(θ2)*cos(θ1-θ2)) - L*θ2'^2*sin(θ1-θ2)) / (L*den)
θ2'' = (2*g*(sin(θ1)*cos(θ1-θ2) - sin(θ2)) + 2*L*θ1'^2*sin(θ1-θ2)) / (L*den)
```

Let me double-check this. I'll use a known formulation. Actually, the most robust source is the Wikipedia double pendulum article or similar.

From Wikipedia's "Double pendulum" article, for m1=m2 and L1=L2:

θ1'' = [ -g (2 sin θ1 - sin(θ2) cos(θ1-θ2)) - L θ2'² sin(θ1-θ2) ] / [ L (2 - cos²(θ1-θ2)) ]

θ2'' = [ 2g (sin θ1 cos(θ1-θ2) - sin θ2) + 2L θ1'² sin(θ1-θ2) ] / [ L (2 - cos²(θ1-θ2)) ]

Let me verify by checking a known reference. Actually, let me use the form from the Lagrangian derivation:

The Lagrangian for the double pendulum with m1=m2=m, L1=L2=L:

T = (1/2) m L² [ θ1'² + θ1'θ2' cos(θ1-θ2) + θ2'² ]
V = -mgL [ 2 cos θ1 + cos θ2 ]

Wait, that's not right. Let me be more careful.

For the double pendulum:
- Position of m1: (L sin θ1, -L cos θ1) relative to pivot
- Position of m2: (L sin θ1 + L sin θ2, -L cos θ1 - L cos θ2) relative to pivot

Velocities:
- v1² = L² θ1'²
- v2² = L²(θ1'² + θ2'² + 2θ1'θ2'cos(θ1-θ2))

T = (1/2)m(L²θ1'² + L²(θ1'² + θ2'² + 2θ1'θ2'cos(θ1-θ2)))
T = (1/2)m L² (2θ1'² + θ2'² + 2θ1'θ2'cos(θ1-θ2))

V = -mgL cos θ1 - mgL(cos θ1 + cos θ2) = -mgL(2cos θ1 + cos θ2)

Euler-Lagrange for θ1:
d/dt(∂T/∂θ1') - ∂T/∂θ1 + ∂V/∂θ1 = 0

∂T/∂θ1' = m L² (2θ1' + θ2' cos(θ1-θ2))
d/dt(∂T/∂θ1') = m L² (2θ1'' + θ2'' cos(θ1-θ2) - θ2'² sin(θ1-θ2))
∂T/∂θ1 = m L² θ1'θ2' (-sin(θ1-θ2)) = -m L² θ1'θ2' sin(θ1-θ2)
∂V/∂θ1 = 2mgL sin θ1

So: m L² (2θ1'' + θ2'' cos(θ1-θ2) - θ2'² sin(θ1-θ2)) + m L² θ1'θ2' sin(θ1-θ2) + 2mgL sin θ1 = 0

Divide by m L:
L (2θ1'' + θ2'' cos(θ1-θ2) - θ2'² sin(θ1-θ2)) + L θ1'θ2' sin(θ1-θ2) + 2g sin θ1 = 0

2θ1'' + θ2'' cos(Δθ) - θ2'² sin(Δθ) + θ1'θ2' sin(Δθ) + (2g/L) sin θ1 = 0  ... (Eq1)

Euler-Lagrange for θ2:
∂T/∂θ2' = m L² (θ2' + θ1' cos(θ1-θ2))
d/dt(∂T/∂θ2') = m L² (θ2'' + θ1'' cos(θ1-θ2) + θ1'² sin(θ1-θ2))
∂T/∂θ2 = m L² θ1'θ2' sin(θ1-θ2) = m L² θ1'θ2' (-sin(θ1-θ2)) ... wait ∂T/∂θ2 = ∂T/∂θ2 (through the cos(θ1-θ2) term)
Actually: ∂T/∂θ2 = m L² θ1'θ2' (-sin(θ1-θ2)) = -m L² θ1'θ2' sin(θ1-θ2)
∂V/∂θ2 = mgL sin θ2

d/dt(∂T/∂θ2') - ∂T/∂θ2 + ∂V/∂θ2 = 0
m L² (θ2'' + θ1'' cos(Δθ) + θ1'² sin(Δθ)) + m L² θ1'θ2' sin(Δθ) + mgL sin θ2 = 0

θ2'' + θ1'' cos(Δθ) + θ1'² sin(Δθ) + θ1'θ2' sin(Δθ) + (g/L) sin θ2 = 0  ... (Eq2)

From Eq1 and Eq2, solve for θ1'' and θ2'':

Eq1: 2θ1'' + θ2'' cos(Δθ) = θ2'² sin(Δθ) - θ1'θ2' sin(Δθ) - (2g/L) sin θ1
Eq2: θ1'' cos(Δθ) + θ2'' = -θ1'² sin(Δθ) - θ1'θ2' sin(Δθ) - (g/L) sin θ2

Let me define:
A = 2θ1'' + cos(Δθ) θ2''
B = cos(Δθ) θ1'' + θ2''

A = θ2'² sin(Δθ) - θ1'θ2' sin(Δθ) - (2g/L) sin θ1
B = -θ1'² sin(Δθ) - θ1'θ2' sin(Δθ) - (g/L) sin θ2

Solve:
From B: θ2'' = B - cos(Δθ) θ1''
Sub into A:
2θ1'' + cos(Δθ)(B - cos(Δθ) θ1'') = A
2θ1'' + cos(Δθ) B - cos²(Δθ) θ1'' = A
(2 - cos²(Δθ)) θ1'' = A - cos(Δθ) B
θ1'' = (A - cos(Δθ) B) / (2 - cos²(Δθ))

Similarly:
From A: θ1'' = (A - cos(Δθ) θ2'') / 2
Sub into B:
cos(Δθ)(A - cos(Δθ) θ2'')/2 + θ2'' = B
cos(Δθ)A/2 - cos²(Δθ) θ2''/2 + θ2'' = B
(1 - cos²(Δθ)/2) θ2'' = B - cos(Δθ)A/2
θ2'' = (B - cos(Δθ)A/2) / (1 - cos²(Δθ)/2)
= (2B - cos(Δθ)A) / (2 - cos²(Δθ))

So:
θ1'' = (A - cos(Δθ) B) / (2 - cos²(Δθ))
θ2'' = (2B - cos(Δθ) A) / (2 - cos²(Δθ))

Where:
A = θ2'² sin(Δθ) - θ1'θ2' sin(Δθ) - (2g/L) sin θ1
B = -θ1'² sin(Δθ) - θ1'θ2' sin(Δθ) - (g/L) sin θ2

This looks correct. Let me simplify:

Num1 = A - cos(Δθ) B
= [θ2'² sin(Δθ) - θ1'θ2' sin(Δθ) - (2g/L) sin θ1] - cos(Δθ)[-θ1'² sin(Δθ) - θ1'θ2' sin(Δθ) - (g/L) sin θ2]
= θ2'² sin(Δθ) - θ1'θ2' sin(Δθ) - (2g/L) sin θ1 + θ1'² sin(Δθ) cos(Δθ) + θ1'θ2' cos(Δθ) sin(Δθ) + (g/L) cos(Δθ) sin θ2

Num2 = 2B - cos(Δθ) A
= 2[-θ1'² sin(Δθ) - θ1'θ2' sin(Δθ) - (g/L) sin θ2] - cos(Δθ)[θ2'² sin(Δθ) - θ1'θ2' sin(Δθ) - (2g/L) sin θ1]
= -2θ1'² sin(Δθ) - 2θ1'θ2' sin(Δθ) - 2(g/L) sin θ2 - θ2'² sin(Δθ) cos(Δθ) + θ1'θ2' cos(Δθ) sin(Δθ) + (2g/L) cos(Δθ) sin θ1

Den = 2 - cos²(Δθ) = 1 + sin²(Δθ)

θ1'' = Num1 / (L * Den)
θ2'' = Num2 / (L * Den)

This is correct but complex. Let me code it up.

Actually, I recall the more common simplified form:

θ1'' = [ -g(2 sin θ1 - sin(θ2) cos(Δθ)) - L θ2'² sin(Δθ) ] / [ L (1 + sin²(Δθ)) ]

θ2'' = [ 2g(sin θ1 cos(Δθ) - sin θ2) + 2L θ1'² sin(Δθ) ] / [ L (1 + sin²(Δθ)) ]

Hmm, let me check if these are equivalent to what I derived. Actually they might not be exactly the same. Let me just use my derived form.

Actually, looking at various online sources, a common correct form for equal masses/lengths is:

θ1'' = ( -g(2 sin θ1 - sin(θ2) cos(θ1-θ2)) - L θ2'² sin(θ1-θ2) ) / ( L(2 - cos²(θ1-θ2)) )

θ2'' = ( 2g(sin θ1 cos(θ1-θ2) - sin θ2) + 2L θ1'² sin(θ1-θ2) ) / ( L(2 - cos²(θ1-θ2)) )

I'll go with this form which matches my derivation (since 2 - cos² = 1 + sin²).

Let me verify once more with a quick check: when θ1 = θ2 = 0 and θ1' = θ2' = 0, we should get θ1'' = θ2'' = 0 (stable equilibrium at bottom). 
- θ1'' numerator: -g(0 - 0) - L*0 = 0 ✓
- θ2'' numerator: 2g(0 - 0) + 0 = 0 ✓

When θ1 = π/2, θ2 = π/2, θ1' = θ2' = 0: should have non-zero acceleration. sin(π/2) = 1.
Num1 = -g(2*1 - 1*1) - 0 = -g(2-1) = -g
Num2 = 2g(1*1 - 1) + 0 = 0
Den = 2 - 0 = 2
θ1'' = -g/(2L), θ2'' = 0. Hmm, θ2'' = 0 seems suspicious. Let me check...

Actually at θ1 = θ2 = π/2, both masses are horizontal. The second mass is at (L+L, -0) = (2L, 0). The first mass is at (L, 0). They're both horizontal, so gravity pulls straight down, creating torque around the pivot. For θ2, it's connected to m1, which is also horizontal. The torque on θ2... hmm it's a bit tricky. But the formula says θ2'' = 0 at this point, which means the second pendulum doesn't rotate initially. That might actually be correct since both masses are at the same angle and there's no relative motion.

Actually wait, I think I need to double-check the equations. Let me use a different approach.

I'll use the well-established formulation. The equations I'll use:

Let Δθ = θ2 - θ1 (note: some sources use θ1 - θ2, I'll be consistent)

For m1 = m2 = m, L1 = L2 = L:

den = 2 - cos²(Δθ) = 1 + sin²(Δθ)

θ1'' = [ -g(2 sin θ1 - sin θ2 cos Δθ) - L θ2'² sin Δθ ] / [L * den]

θ2'' = [ 2g(sin θ1 cos Δθ - sin θ2) + 2L θ1'² sin Δθ ] / [L * den]

This is what I'll implement. Now for the integration, I'll use velocity Verlet or RK4 for stability. RK4 is better for accuracy. With a small enough time step, even RK4 or velocity Verlet should work.

Actually, for a pendulum, a symplectic integrator or RK4 is good. Let me use RK4 for the angular acceleration equations. Or I could use the simpler approach: treat (θ1, θ2, ω1, ω2) as state and integrate with RK4.

State vector: y = [θ1, θ2, ω1, ω2]
Derivatives: dy/dt = [ω1, ω2, α1, α2]

Where α1, α2 are computed from the equations above.

RK4 integration:
k1 = f(t, y)
k2 = f(t + h/2, y + h/2 * k1)
k3 = f(t + h/2, y + h/2 * k2)
k4 = f(t + h, y + h * k3)
y_new = y + (h/6)(k1 + 2k2 + 2k3 + k4)

This should be stable with a reasonable time step like 0.001-0.01 seconds.

For the simulation, I'll use:
- dt = 0.005 seconds (5ms) per step
- Run multiple steps per frame for efficiency

Actually, let me use dt = 0.001 (1ms) and run e.g. 5-10 steps per animation frame for accuracy and stability.

Now for the trail: I'll store the position of the lower bob (m2) and draw fading circles. I'll keep the last N positions (e.g., 100-200) and draw them with decreasing opacity.

Canvas dimensions: let me use a reasonable size, say 800x600 or responsive.

Pivot at center-top of the canvas. Rod lengths scaled to fit.

Let me set up:
- Canvas: let's say 800x700
- Pivot at (400, 100)
- Rod length L = 200 pixels (but we need to map simulation units to pixels)

Actually, let me think about the simulation units. If g = 9.81 m/s², and we want the pendulum to look good, we can use L = 1.0 (in simulation units) and scale to pixels. Or we can set L = 200 pixels directly and adjust g accordingly.

Option 1: Use simulation units of meters, L = 1.0 m, g = 9.81 m/s². Then scale 1 m = 200 px for display.
Option 2: Use pixel units directly, with g scaled to look good visually.

I'll go with Option 1: L = 1.0 (simulation), scale to pixels. g = 9.81.

For initial conditions:
- θ1 = 120° = 120 * π/180 = 2.094 rad
- θ2 = -10° = -10 * π/180 = -0.1745 rad
- ω1 = ω2 = 0

The angles are measured from downward vertical. So θ = 0 means hanging straight down. θ > 0 means swinging to the right (counterclockwise from downward).

Wait, I need to be careful about the sign convention. Let me define:
- θ measured from downward vertical, positive = counterclockwise (to the right)
- So θ1 = 120° means the first rod is pointing up and to the left (120° from downward = 60° above horizontal to the left)
- θ2 = -10° means the second rod is slightly to the right of downward

Positions:
- m1: (L sin θ1, -L cos θ1) relative to pivot (y positive = down in canvas coords)
- m2: (L sin θ1 + L sin θ2, -L cos θ1 - L cos θ2) relative to pivot

Wait, in standard math convention (y up), the position would be (L sin θ, -L cos θ). But in canvas coords (y down), it would be (L sin θ, L cos θ) if θ is from downward. Let me be careful.

If θ is measured from the downward vertical, and positive θ means the pendulum swings to the right (counterclockwise in standard math, but clockwise in canvas since y is flipped), then:

In canvas coordinates (origin at pivot, x right, y down):
- m1 position: (L * sin(θ1), L * cos(θ1))
- m2 position: (L * sin(θ1) + L * sin(θ2), L * cos(θ1) + L * cos(θ2))

This makes sense: when θ = 0, the bob is at (0, L) i.e., straight down. When θ > 0, sin(θ) > 0, so the bob is to the right.

The equations of motion I derived assume the standard Lagrangian formulation where θ is measured from the downward vertical. The sign in the equations should be consistent.

Let me verify: at θ1 = small positive, θ2 = 0, ω1 = ω2 = 0:
- sin(θ1) ≈ θ1 > 0
- Num1 ≈ -g(2θ1 - 0) - 0 = -2gθ1
- θ1'' ≈ -2gθ1 / L (negative, so pendulum swings back toward 0) ✓

Good, the sign seems correct.

Now let me code this up.

For the trail, I'll store positions of m2 in canvas coordinates. I'll keep maybe 150 positions and draw them as circles with decreasing opacity (fading trail).

Let me write the full code now.

```html
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<title>Double Pendulum</title>
<style>
  /* styles */
</style>
</head>
<body>
<canvas id="canvas"></canvas>
<div id="controls">
  <button id="pauseBtn">Pause</button>
  <button id="resetBtn">Reset</button>
</div>
<script>
  // simulation code
</script>
</body>
</html>
```

Let me think about the RK4 implementation more carefully.

State: [θ1, θ2, ω1, ω2]

Derivatives function f(state):
  θ1, θ2, ω1, ω2 = state
  Δθ = θ2 - θ1
  den = 2 - cos²(Δθ)  // = 1 + sin²(Δθ)
  
  α1 = (-g*(2*sin(θ1) - sin(θ2)*cos(Δθ)) - L*ω2²*sin(Δθ)) / (L*den)
  α2 = (2*g*(sin(θ1)*cos(Δθ) - sin(θ2)) + 2*L*ω1²*sin(Δθ)) / (L*den)
  
  return [ω1, ω2, α1, α2]

RK4 step:
  k1 = f(y)
  k2 = f(y + dt/2 * k1)
  k3 = f(y + dt/2 * k2)
  k4 = f(y + dt * k3)
  y_new = y + dt/6 * (k1 + 2*k2 + 2*k3 + k4)

This should work. Let me also consider using a smaller dt and sub-stepping.

For stability, dt = 0.001 seconds with 5 sub-steps per frame (at 60fps, that's 5ms per frame, so 5 sub-steps of 1ms each). Or I could use dt = 0.005 with fewer sub-steps. Let me use dt = 0.001 and 8 sub-steps per frame, giving 8ms of simulation time per 16.67ms frame, which should look smooth.

Actually, let me think about this differently. At 60fps, each frame is ~16.67ms. If I use dt = 0.001s (1ms) and do 16 sub-steps per frame, that gives real-time simulation. That might be slow but accurate. 

Alternatively, I could use a larger dt like 0.002 or 0.005 and only do a few sub-steps. The pendulum dynamics with L=1m and g=9.81 have a natural frequency of sqrt(g/L) = ~3.13 rad/s, period ~2s. So dt = 0.005 gives ~400 steps per period, which is plenty for RK4.

Let me use dt = 0.002 (2ms) with 8 sub-steps per frame. That gives 16ms simulation per frame, close to real-time. 

Actually, I should make the simulation speed independent of frame rate. Let me track time properly.

Let me use a fixed dt = 0.001 and do enough sub-steps to cover the elapsed real time. Or use dt = 0.002 and scale the number of steps based on elapsed time.

For simplicity, let me use dt = 0.002 and always do 8 steps per frame. This gives approximately real-time simulation at 60fps.

Actually, let me use a more robust approach: dt = 0.001, and each frame, do N steps where N = ceil(elapsedTime / dt), with a maximum to prevent spiral of death. But for simplicity, let me just use a fixed number of steps per frame.

Let me go with: dt = 0.001, stepsPerFrame = 16 (giving 16ms sim time per frame, close to real-time at 60fps). Actually, the simulation might be slow with 16 RK4 steps per frame (each step requires 4 evaluations of the derivatives). Let me test with a smaller number.

Actually, let me use dt = 0.002 and 8 steps per frame. That should be fast enough.

Let me now write the complete code.

For the trail, I'll store positions as {x, y, age} or just keep an array and shift. I'll draw circles with opacity based on position in the array (older = more transparent).

Trail parameters:
- Max trail length: 150 points
- Each point drawn as a circle with radius 2-3px
- Opacity from 0.05 (oldest) to 1.0 (newest) or similar
- Color: a nice hue, maybe shifting or fixed

Let me use a colorful trail - maybe the color shifts through hues, or it's a fixed bright color with varying opacity. The user said "colorful fading trail" so let me make it colorful - perhaps cycling through hues or using a gradient.

I'll make the trail with a hue that shifts based on the position in the trail, creating a rainbow-like fading effect.

Let me finalize the design and write the code.

Canvas size: 800 x 650 (or responsive)
Pivot position: (400, 120)
Scale: L = 1.0 sim units = 180 pixels (so rod length on screen = 180px)

Actually, let me make it responsive. I'll use a fixed aspect ratio and scale to fit.

Or I'll just set a fixed canvas size of 800x650 and center things nicely.

Let me write the code now.

For the trail color, I'll use HSL with the hue cycling through the trail. Newest point is fully opaque, oldest is nearly transparent.

Let me also add a nice background and styling.

Here's my plan for the full HTML file:

```html
<!DOCTYPE html>
<html>
<head>
  <meta charset="UTF-8">
  <title>Double Pendulum Simulation</title>
  <style>
    * { margin: 0; padding: 0; box-sizing: border-box; }
    body {
      background: #1a1a2e;
      display: flex;
      justify-content: center;
      align-items: center;
      min-height: 100vh;
      font-family: 'Segoe UI', sans-serif;
    }
    .container {
      background: #16213e;
      border-radius: 16px;
      padding: 20px;
      box-shadow: 0 10px 40px rgba(0,0,0,0.5);
    }
    canvas {
      display: block;
      border-radius: 8px;
      background: #0f0f23;
    }
    .controls {
      display: flex;
      gap: 12px;
      margin-top: 16px;
      justify-content: center;
    }
    button {
      padding: 10px 24px;
      border: none;
      border-radius: 8px;
      font-size: 15px;
      font-weight: 600;
      cursor: pointer;
      transition: transform 0.1s, background 0.2s;
    }
    button:active { transform: scale(0.96); }
    #pauseBtn {
      background: #e94560;
      color: white;
    }
    #pauseBtn.paused {
      background: #0f3460;
    }
    #resetBtn {
      background: #533483;
      color: white;
    }
    .info {
      color: #a0a0b0;
      font-size: 13px;
      text-align: center;
      margin-top: 8px;
    }
  </style>
</head>
<body>
  <div class="container">
    <canvas id="canvas" width="800" height="600"></canvas>
    <div class="controls">
      <button id="pauseBtn">Pause</button>
      <button id="resetBtn">Reset</button>
    </div>
    <div class="info">Double pendulum · g = 9.81 · L = 1.0 · m = 1.0</div>
  </div>
  <script>
    // Full simulation code
  </script>
</body>
</html>
```

Now let me write the JavaScript carefully.

```javascript
const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
const W = canvas.width;
const H = canvas.height;

// Simulation parameters
const g = 9.81;
const L = 1.0;
const m = 1.0; // not needed for equal masses

// Initial conditions
const INIT_THETA1 = 120 * Math.PI / 180;
const INIT_THETA2 = -10 * Math.PI / 180;
const INIT_W1 = 0;
const INIT_W2 = 0;

// Display parameters
const PIVOT_X = W / 2;
const PIVOT_Y = 100;
const SCALE = 180; // pixels per simulation unit (L=1 -> 180px)

// State
let theta1 = INIT_THETA1;
let theta2 = INIT_THETA2;
let omega1 = INIT_W1;
let omega2 = INIT_W2;

// Simulation control
let paused = false;
let trail = []; // array of {x, y}

// Trail parameters
const MAX_TRAIL = 180;
const TRAIL_RADIUS = 3;

// Time step
const DT = 0.001;
const STEPS_PER_FRAME = 16;

function derivatives(s) {
  const [t1, t2, w1, w2] = s;
  const dtheta = t2 - t1;
  const cos_d = Math.cos(dtheta);
  const sin_d = Math.sin(dtheta);
  const den = 2 - cos_d * cos_d; // = 1 + sin²(dtheta)
  
  const alpha1 = (-g * (2 * Math.sin(t1) - Math.sin(t2) * cos_d) - L * w2 * w2 * sin_d) / (L * den);
  const alpha2 = (2 * g * (Math.sin(t1) * cos_d - Math.sin(t2)) + 2 * L * w1 * w1 * sin_d) / (L * den);
  
  return [w1, w2, alpha1, alpha2];
}

function rk4Step(state, dt) {
  const k1 = derivatives(state);
  
  const s2 = state.map((v, i) => v + 0.5 * dt * k1[i]);
  const k2 = derivatives(s2);
  
  const s3 = state.map((v, i) => v + 0.5 * dt * k2[i]);
  const k3 = derivatives(s3);
  
  const s4 = state.map((v, i) => v + dt * k3[i]);
  const k4 = derivatives(s4);
  
  return state.map((v, i) => v + (dt / 6) * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]));
}

function getBobPositions() {
  const x1 = PIVOT_X + SCALE * Math.sin(theta1);
  const y1 = PIVOT_Y + SCALE * Math.cos(theta1);
  const x2 = x1 + SCALE * Math.sin(theta2);
  const y2 = y1 + SCALE * Math.cos(theta2);
  return { x1, y1, x2, y2 };
}

function update() {
  if (!paused) {
    let state = [theta1, theta2, omega1, omega2];
    for (let i = 0; i < STEPS_PER_FRAME; i++) {
      state = rk4Step(state, DT);
    }
    theta1 = state[0];
    theta2 = state[1];
    omega1 = state[2];
    omega2 = state[3];
    
    // Add to trail
    const { x2, y2 } = getBobPositions();
    trail.push({ x: x2, y: y2 });
    if (trail.length > MAX_TRAIL) {
      trail.shift();
    }
  }
}

function draw() {
  ctx.clearRect(0, 0, W, H);
  
  const { x1, y1, x2, y2 } = getBobPositions();
  
  // Draw trail
  for (let i = 0; i < trail.length; i++) {
    const t = trail[i];
    const alpha = (i + 1) / trail.length; // 0 to 1, newest is 1
    const hue = (i / trail.length) * 360; // color cycles through trail
    ctx.beginPath();
    ctx.arc(t.x, t.y, TRAIL_RADIUS * (0.3 + 0.7 * alpha), 0, Math.PI * 2);
    ctx.fillStyle = `hsla(${hue}, 85%, 60%, ${alpha * 0.7})`;
    ctx.fill();
  }
  
  // Draw rods
  ctx.strokeStyle = '#8899aa';
  ctx.lineWidth = 3;
  ctx.beginPath();
  ctx.moveTo(PIVOT_X, PIVOT_Y);
  ctx.lineTo(x1, y1);
  ctx.stroke();
  
  ctx.beginPath();
  ctx.moveTo(x1, y1);
  ctx.lineTo(x2, y2);
  ctx.stroke();
  
  // Draw bobs
  // Upper bob
  ctx.beginPath();
  ctx.arc(x1, y1, 14, 0, Math.PI * 2);
  ctx.fillStyle = '#4fc3f7';
  ctx.fill();
  ctx.strokeStyle = '#0288d1';
  ctx.lineWidth = 2;
  ctx.stroke();
  
  // Lower bob
  ctx.beginPath();
  ctx.arc(x2, y2, 18, 0, Math.PI * 2);
  ctx.fillStyle = '#ff7043';
  ctx.fill();
  ctx.strokeStyle = '#d84315';
  ctx.lineWidth = 2;
  ctx.stroke();
  
  // Draw pivot
  ctx.beginPath();
  ctx.arc(PIVOT_X, PIVOT_Y, 8, 0, Math.PI * 2);
  ctx.fillStyle = '#ffd740';
  ctx.fill();
  ctx.strokeStyle = '#ff8f00';
  ctx.lineWidth = 2;
  ctx.stroke();
}

function loop() {
  update();
  draw();
  requestAnimationFrame(loop);
}

// Controls
document.getElementById('pauseBtn').addEventListener('click', () => {
  paused = !paused;
  const btn = document.getElementById('pauseBtn');
  btn.textContent = paused ? 'Resume' : 'Pause';
  btn.classList.toggle('paused', paused);
});

document.getElementById('resetBtn').addEventListener('click', () => {
  theta1 = INIT_THETA1;
  theta2 = INIT_THETA2;
  omega1 = INIT_W1;
  omega2 = INIT_W2;
  trail = [];
});

// Start
loop();
```

Wait, I need to double-check the trail drawing. The issue is that when paused, we don't add to the trail, but we should still draw the existing trail. That's fine since we draw in `draw()` regardless.

Also, when paused, the bobs don't move, so the trail stays in place. Good.

Let me also think about the "colorful fading trail" requirement. My current approach:
- Hue cycles through the trail (oldest = 0°, newest = 360°)
- Opacity increases from oldest to newest
- Radius also increases slightly for newer points

This should look colorful and fading. Good.

Let me also verify the physics with a quick mental check:
- Initial: θ1 = 120° (second quadrant, rod pointing up-left), θ2 = -10° (slightly right of downward)
- The first rod is at 120° from downward, so it's 60° above the horizontal on the left side
- The second rod is at -10° from downward, so it's slightly to the right of straight down
- Both start with zero velocity

This should create interesting coupled motion. The first pendulum will swing down and to the right, while the second pendulum will be dragged along.

One concern: θ1 = 120° is quite large. The pendulum might swing over the top. That could be interesting but also might cause the simulation to become unstable if not careful. With RK4 and small dt, it should be fine.

Actually, with θ1 = 120° and θ2 = -10°, the second bob's initial position is:
x2 = L(sin(120°) + sin(-10°)) = L(0.866 - 0.174) = 0.692L
y2 = L(cos(120°) + cos(-10°)) = L(-0.5 + 0.985) = 0.485L

So the second bob starts at about (0.692L, 0.485L) relative to pivot, which is in the upper-right quadrant. The first bob is at (0.866L, -0.5L) which is left and above the pivot. Wait, cos(120°) = -0.5, so y1 = L * cos(120°) = -0.5L. In canvas coords (y down), this means y1 = PIVOT_Y + SCALE * (-0.5) = PIVOT_Y - 0.5*SCALE. That's above the pivot. Good, that's correct for 120° from downward (pointing up and to the left).

Let me verify: at θ = 120° from downward vertical, the bob is at angle 120° measured from the downward direction. If downward is 0°, then 120° counterclockwise puts the bob in the upper-left quadrant (since 90° is horizontal left, 120° is 30° above horizontal left). 

In my coordinate system:
x = L sin(θ) = L sin(120°) = L * 0.866 (positive = right)
y = L cos(θ) = L cos(120°) = L * (-0.5) (negative = up in canvas coords since y down is positive)

Wait, that gives x positive (right) and y negative (up). But 120° from downward should be... let me think again.

If I measure θ from the downward vertical, with positive being counterclockwise (to the right in standard math, but since y is flipped in canvas, positive θ should swing to the right on screen):

In standard math (y up):
- Position at angle θ from downward: (L sin θ, -L cos θ)
  - θ=0: (0, -L) = straight down ✓
  - θ=90°: (L, 0) = horizontal right ✓
  - θ=120°: (L sin 120°, -L cos 120°) = (0.866L, 0.5L) = right and up ✓
  - θ=180°: (0, L) = straight up ✓

In canvas coords (y down), I need to flip y:
- Position: (L sin θ, L cos θ) in canvas coords relative to pivot
  - θ=0: (0, L) = straight down ✓
  - θ=90°: (L, 0) = horizontal right ✓
  - θ=120°: (0.866L, -0.5L) = right and up ✓
  - θ=180°: (0, -L) = straight up ✓

OK so my coordinate mapping is correct: (L sin θ, L cos θ) in canvas coords.

So at θ1 = 120°:
- First bob: (0.866L, -0.5L) relative to pivot = right and above pivot ✓

At θ2 = -10°:
- Second bob relative to first bob: (L sin(-10°), L cos(-10°)) = (-0.174L, 0.985L)
- Second bob absolute: (0.866L - 0.174L, -0.5L + 0.985L) = (0.692L, 0.485L)
- This is right and below the pivot, but above the first bob's level? No, y=0.485L is below the pivot (which is at y=0 in relative coords). The first bob is at y=-0.5L (above pivot). So the second bob is below the pivot but the first bob is above the pivot. That means the rods cross? No, the first rod goes from pivot to (0.866L, -0.5L), and the second rod goes from there to (0.692L, 0.485L). They don't cross; the second rod goes from the first bob downward and slightly left to the second bob.

This seems physically reasonable. The pendulum starts in an unusual configuration with the first rod pointing up-right and the second rod hanging from the first bob down to the right of the pivot.

OK, I think the physics is correct. Let me finalize the code.

One more thing: I should make sure the canvas is responsive or at least looks good. Let me set it to a fixed size that works well.

Actually, let me reconsider the scale. With L = 1.0 and SCALE = 180, the rods are 180px long. The total height from pivot to lower bob when hanging straight down is 2L = 2.0 sim units = 360px. With pivot at y=100, the lowest point is at y=460. The canvas height is 600, so there's room. Good.

The width: the maximum horizontal extent is when both rods are horizontal, giving x = L + L = 2L = 360px from center. With canvas width 800, center at 400, that's x from 40 to 760. Fits.

But with θ1 = 120° initially, the first bob is at x = 0.866 * 180 = 156px right of pivot, and the second bob is at x = 0.692 * 180 = 125px right of pivot. So they start on the right side. The pendulum will swing. Should be fine.

Actually wait, I realize I should double-check my equation derivation one more time. Let me use a well-known reference.

From the standard double pendulum derivation for equal masses and equal lengths (m1=m2=m, L1=L2=L), the equations of motion are:

θ1'' = [ -g (2 sin θ1 - sin(θ2) cos(θ1-θ2)) - L θ2'² sin(θ1-θ2) ] / [ L (2 - cos²(θ1-θ2)) ]

θ2'' = [ 2g (sin θ1 cos(θ1-θ2) - sin θ2) + 2L θ1'² sin(θ1-θ2) ] / [ L (2 - cos²(θ1-θ2)) ]

Yes, this is what I have. Let me verify with a simple case: when θ1 = θ2 = 0 (both hanging straight down), θ1' = θ2' = 0:
- sin(0) = 0, cos(0) = 1, sin(0) = 0
- Num1 = -g(0 - 0) - 0 = 0
- Num2 = 2g(0 - 0) + 0 = 0
- α1 = α2 = 0 ✓ (equilibrium)

When θ1 = π (straight up), θ2 = 0, θ1' = θ2' = 0:
- sin(π) = 0, sin(0) = 0, cos(π) = -1, sin(π) = 0
- Num1 = -g(0 - 0*(-1)) - 0 = 0
- Num2 = 2g(0*(-1) - 0) + 0 = 0
- α1 = α2 = 0 

Hmm, that says the inverted pendulum is in equilibrium, which is correct for the first pendulum but the second pendulum (hanging from the inverted first) should also be in equilibrium at θ2=0 (hanging straight down from an inverted first pendulum... wait, if θ1=π, the first pendulum is straight up, and θ2=0 means the second pendulum hangs straight down from the first bob. That's an equilibrium too. So α1=α2=0 is correct. ✓

Let me test one more: θ1 = π/2 (horizontal right), θ2 = π/2, θ1' = θ2' = 0:
- sin(π/2) = 1, cos(0) = 1, sin(0) = 0
- Num1 = -g(2*1 - 1*1) - 0 = -g
- Num2 = 2g(1*1 - 1) + 0 = 0
- Den = 2 - 1 = 1
- α1 = -g/L, α2 = 0

So the first pendulum accelerates downward (negative angular acceleration means moving toward θ1=0), and the second pendulum has zero initial acceleration. That makes sense because both pendulums are at the same angle, so there's no relative torque on the second one initially. ✓

OK, the equations are correct. Let me finalize.

Actually, I want to double-check the sign of my angular acceleration. In my system:
- θ increases counterclockwise (to the right on screen)
- Positive α means θ increases (accelerates counterclockwise)
- Gravity pulls downward, so if the bob is to the right of the pivot (θ > 0), gravity creates a clockwise torque (negative α), pulling it back toward θ=0

Let me check: θ1 = small positive, θ2 = 0, ω1 = ω2 = 0:
- sin(θ1) ≈ θ1 > 0
- Num1 ≈ -g(2θ1 - 0) - 0 = -2gθ1
- α1 ≈ -2gθ1/L (negative) ✓ (pulls back toward 0)

Good.

Now let me think about the trail. The trail stores canvas coordinates of the lower bob. When we reset, we clear the trail. When we pause, we stop updating (and stop adding to trail). The trail should persist when paused and be drawn correctly.

I think the code is ready. Let me also add some visual polish:
- A subtle grid or background pattern
- Nice colors for the bobs
- Glowing effect on the trail

Actually, let me keep it clean and not over-complicate. The user asked for a "colorful fading trail" which I'm implementing.

Let me also make sure the trail doesn't draw outside the canvas. Since the bore positions are always within the canvas (the pendulum is bounded), this should be fine.

One more consideration: when the simulation is paused, should we still draw the trail? Yes, the trail should persist. My code does this since `draw()` always draws the trail regardless of pause state.

Let me also add a feature: when paused, show a visual indicator (like "PAUSED" text). Actually, the button text change from "Pause" to "Resume" should be sufficient.

Alright, let me write the final code. I'll also make the canvas size a bit larger for better visibility.

Let me use 800x650 for the canvas.

Actually, let me reconsider. 800x600 is more standard. The pivot at y=100 with L=180px means the lowest point (both hanging down) is at y = 100 + 360 = 460. The canvas is 600 tall, so there's 140px of space below. That's fine.

The maximum height (both pointing up) would be at y = 100 - 360 = -260, which is above the canvas. But that's an extreme case that won't happen with the given initial conditions (θ1=120°, θ2=-10°). The pendulum might swing up but probably won't go fully inverted. Even if it does, it'll be a brief moment and the trail won't be visible there.

Actually, with θ1 starting at 120° (which is 60° above horizontal), the pendulum has significant potential energy and could swing up quite high. Let me make sure the canvas is tall enough.

The maximum y for the lower bob: when both rods point straight up, the lower bob is at y = PIVOT_Y - 2*L*SCALE = 100 - 360 = -260. That's above the canvas. But this is an extreme configuration that requires the pendulum to swing all the way up, which is unlikely with the given initial conditions (it would need to conserve energy perfectly).

With the initial conditions:
- PE initial = -mgL(2cos(120°) + cos(-10°)) = -mgL(2*(-0.5) + 0.985) = -mgL(-1 + 0.985) = -mgL(-0.015) = 0.015mgL
- PE at bottom (both at 0°) = -mgL(2 + 1) = -3mgL
- So the pendulum has much less PE initially than at the bottom, meaning it will swing down and gain kinetic energy, then swing back up. The maximum height it can reach is when all KE converts back to PE, which would be the same as initial PE = 0.015mgL. This is just slightly above the bottom position, so the pendulum won't swing very high on the other side.

Wait, that doesn't seem right. Let me recalculate.

PE = V = -mgL(2cos θ1 + cos θ2)

At initial: θ1 = 120° = 2.094 rad, θ2 = -10° = -0.1745 rad
V_initial = -mgL(2*cos(120°) + cos(-10°)) = -mgL(2*(-0.5) + 0.9848) = -mgL(-1 + 0.9848) = -mgL(-0.0152) = 0.0152 mgL

At bottom (θ1=θ2=0): V_bottom = -mgL(2+1) = -3mgL

So the initial PE is HIGHER than the bottom PE (0.0152 mgL > -3 mgL). The pendulum starts above the bottom and will swing down. The total energy is E = T + V = 0 + 0.0152 mgL = 0.0152 mgL.

The pendulum will swing down, gaining KE, then swing back up. The maximum height it can reach is when T=0 again, so V = E = 0.0152 mgL. This is the same as the initial PE. So the pendulum will oscillate between configurations with V = 0.0152 mgL.

Can it reach the inverted position (θ1=π, θ2=0)? V at that point = -mgL(2*cos(π) + cos(0)) = -mgL(-2+1) = -mgL(-1) = mgL. Since E = 0.0152 mgL < mgL, the pendulum doesn't have enough energy to reach the inverted position. Good, it won't go above the canvas.

Actually wait, I need to be more careful. The pendulum could reach configurations where one rod is up and the other is down, etc. Let me think about the maximum y of the lower bob.

The lower bob's y position (in sim coords, relative to pivot): y2 = L(cos θ1 + cos θ2)

For the bob to go above the canvas (y < 0 in canvas coords, meaning above pivot), we need y2 < 0, i.e., cos θ1 + cos θ2 < 0.

With energy E = 0.0152 mgL, the pendulum is bounded. The maximum of cos θ1 + cos θ2 subject to the energy constraint... this is complex. But intuitively, since the pendulum doesn't have enough energy to reach the inverted position, it won't go very high. The initial position already has the first bob at y1 = L cos(120°) = -0.5L (above pivot) and the second bob at y2 = L(cos(120°) + cos(-10°)) = L(-0.5 + 0.985) = 0.485L (below pivot). So initially, the lower bob is below the pivot.

As the pendulum swings, the lower bob might go higher, but probably not above the pivot. Let me just make the canvas tall enough to be safe.

With PIVOT_Y = 100 and SCALE = 180, the lower bob at y2 = 0 (at pivot level) would be at canvas y = 100. The initial y2 = 0.485L = 0.485 * 180 = 87.3px below pivot, so canvas y = 100 + 87.3 = 187.3. 

If the pendulum swings such that y2 = -0.5L (first bob at bottom, second bob pointing up), canvas y = 100 - 90 = 10. That's still on screen.

If y2 = -1.0L (both pointing up), canvas y = 100 - 180 = -80. That's off screen. But as I argued, this won't happen with the given energy.

To be safe, let me set PIVOT_Y = 120 and SCALE = 160. Then:
- Bottom position: y = 120 + 2*160 = 440 (on screen with 600px height)
- Potential high position (y2 = -L): y = 120 - 160 = -40 (slightly off screen, but won't happen)

Actually, let me just use PIVOT_Y = 120, SCALE = 170, canvas 800x620.

Bottom: y = 120 + 340 = 460
If y2 = -L: y = 120 - 170 = -50 (off screen, but won't happen)

This should be safe. Let me go with canvas 800x600, PIVOT_Y = 110, SCALE = 170.

Bottom: y = 110 + 340 = 450 (on screen)
y2 = -L: y = 110 - 170 = -60 (off screen, but won't happen)

Good enough. Let me finalize.

Actually, I realize I'm overthinking the canvas size. Let me just use reasonable values and make sure the initial configuration is visible. The pendulum might swing a bit but shouldn't go off screen with the given initial conditions.

Let me finalize the code now.

One more thing: I should make the trail "colorful". My approach of cycling hue through the trail should work. Let me also make the trail points slightly larger and more vibrant.

Here's my final trail drawing:
- For each point in the trail (from oldest to newest):
  - alpha = (i + 1) / trail.length (0 to 1)
  - hue = (i / trail.length) * 300 (cycle through 300° of hue, e.g., from 200 to 500 mod 360)
  - radius = 1.5 + 2.5 * alpha (grows from 1.5 to 4)
  - fill with hsla(hue, 80%, 60%, alpha * 0.6)

This creates a colorful, fading trail that's more opaque and larger for newer points.

Let me write the final code now. I'll be thorough and make sure everything works.

Actually, for the trail, I want it to look really nice. Let me use a multi-color approach where the trail transitions through several colors. Or I could use a single hue with varying brightness. The "colorful" requirement suggests multiple colors.

I'll use a hue range from 180 (cyan) through 280 (purple) to 360/0 (red/pink) as the trail goes from oldest to newest. This gives a nice cyan-to-purple-to-pink gradient.

hue = 180 + (i / trail.length) * 180 (range 180 to 360)
This goes from cyan → blue → purple → pink/red.

With opacity increasing from oldest to newest, this should look great.

Let me finalize everything and write the code.

Oh wait, I should also consider: should the trail be drawn even when the simulation is reset? No, reset clears the trail. Should the trail be drawn when paused? Yes, the trail should persist and be visible even when paused. My code handles both cases correctly.

Let me write the complete HTML file now.

Let me also add: show angle values or some debug info? The user didn't ask for it, so I'll keep it clean. Just the simulation, controls, and trail.

Here's my final code structure:

```
<!DOCTYPE html>
<html>
<head>
  <meta charset="UTF-8">
  <title>Double Pendulum</title>
  <style>
    /* Dark theme, centered layout, nice buttons */
  </style>
</head>
<body>
  <div class="container">
    <canvas id="canvas" width="800" height="600"></canvas>
    <div class="controls">
      <button id="pauseBtn">Pause</button>
      <button id="resetBtn">Reset</button>
    </div>
  </div>
  <script>
    // Complete simulation
  </script>
</body>
</html>
```

Let me write it all out now.

Actually, I want to make one more check on the equations. Let me verify with a simple Python-like mental calculation.

For θ1 = 0.1 rad, θ2 = -0.05 rad, ω1 = ω2 = 0:
Δθ = -0.15
sin(Δθ) ≈ -0.149
cos(Δθ) ≈ 0.989
sin(θ1) ≈ 0.0998
sin(θ2) ≈ -0.0500

Num1 = -g(2*0.0998 - (-0.0500)*0.989) - L*0*(-0.149)
     = -9.81(0.1996 + 0.0495)
     = -9.81 * 0.2491
     = -2.444

Num2 = 2g(0.0998*0.989 - (-0.0500)) + 0
     = 2*9.81*(0.0987 + 0.0500)
     = 19.62 * 0.1487
     = 2.918

Den = 2 - 0.989² = 2 - 0.978 = 1.022

α1 = -2.444 / (1.0 * 1.022) = -2.392 rad/s²
α2 = 2.918 / (1.0 * 1.022) = 2.855 rad/s²

So the first pendulum accelerates clockwise (negative) and the second accelerates counterclockwise (positive). This makes sense: the first pendulum is slightly to the right of vertical and wants to swing back left, while the second pendulum is slightly to the left and wants to swing back right. But the coupling through the rods creates the observed accelerations.

This looks reasonable. Let me proceed with writing the code.

Actually, I just realized I should double-check my equation against a known reference more carefully. Let me look at the Wikipedia double pendulum article equations.

From Wikipedia (https://en.wikipedia.org/wiki/Double_pendulum), for the general case:

(m1+m2)L1²θ1'' + m2 L1 L2 θ2'' cos(θ1-θ2) + m2 L1 L2 θ2'² sin(θ1-θ2) + (m1+m2)g L1 sin θ1 = 0

m2 L2² θ2'' + m2 L1 L2 θ1'' cos(θ1-θ2) - m2 L1 L2 θ1'² sin(θ1-θ2) + m2 g L2 sin θ2 = 0

For m1=m2=m, L1=L2=L:

2m L² θ1'' + m L² θ2'' cos(Δθ) + m L² θ2'² sin(Δθ) + 2mg L sin θ1 = 0
=> 2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0 ... (A)

m L² θ2'' + m L² θ1'' cos(Δθ) - m L² θ1'² sin(Δθ) + mgL sin θ2 = 0
=> θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + (g/L) sin θ2 = 0 ... (B)

From (B): θ2'' = -θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2

Sub into (A):
2θ1'' + (-θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2) cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0

2θ1'' - θ1'' cos²(Δθ) + θ1'² sin(Δθ) cos(Δθ) - (g/L) sin θ2 cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0

θ1'' (2 - cos²(Δθ)) = -θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1

θ1'' = [ -θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1 ] / (2 - cos²(Δθ))

Hmm, this is different from what I had before. Let me compare.

My previous derivation gave:
Num1 = θ2'² sin(Δθ) - θ1'θ2' sin(Δθ) - (2g/L) sin θ1 + θ1'² sin(Δθ) cos(Δθ) + θ1'θ2' cos(Δθ) sin(Δθ) + (g/L) cos(Δθ) sin θ2

And: θ1'' = Num1 / (L * Den) where Den = 2 - cos²(Δθ)

So: θ1'' = [θ2'² sin(Δθ) - θ1'θ2' sin(Δθ) - (2g/L) sin θ1 + θ1'² sin(Δθ) cos(Δθ) + θ1'θ2' cos(Δθ) sin(Δθ) + (g/L) cos(Δθ) sin θ2] / [L(2 - cos²(Δθ))]

Wait, I have an extra factor of L in the denominator. Let me recheck.

In my derivation, I had:
2θ1'' + θ2'' cos(Δθ) - θ2'² sin(Δθ) + θ1'θ2' sin(Δθ) + (2g/L) sin θ1 = 0

But from Wikipedia:
2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0

The sign of the θ2'² sin(Δθ) term is different! In my derivation, I have -θ2'² sin(Δθ), but Wikipedia has +θ2'² sin(Δθ).

Let me recheck my derivation. From Eq1:
2θ1'' + θ2'' cos(Δθ) - θ2'² sin(Δθ) + θ1'θ2' sin(Δθ) + (2g/L) sin θ1 = 0

Hmm, where did the -θ2'² sin(Δθ) come from? Let me re-derive.

From the Lagrangian:
T = (1/2)m L² (2θ1'² + θ2'² + 2θ1'θ2'cos(θ1-θ2))

∂T/∂θ1' = m L² (2θ1' + θ2' cos(θ1-θ2))
d/dt(∂T/∂θ1') = m L² (2θ1'' + θ2'' cos(θ1-θ2) - θ2'² sin(θ1-θ2) * (-1))
               = m L² (2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ))

Wait! The derivative of cos(θ1-θ2) with respect to time is:
d/dt cos(θ1-θ2) = -sin(θ1-θ2) * (θ1' - θ2') = -sin(Δθ) * (θ1' - θ2')

But we're differentiating ∂T/∂θ1' = m L² (2θ1' + θ2' cos(θ1-θ2)):
d/dt [θ2' cos(θ1-θ2)] = θ2'' cos(θ1-θ2) + θ2' * d/dt[cos(θ1-θ2)]
= θ2'' cos(Δθ) + θ2' * [-sin(Δθ) * (θ1' - θ2')]
= θ2'' cos(Δθ) - θ2' sin(Δθ) θ1' + θ2' sin(Δθ) θ2'
= θ2'' cos(Δθ) - θ1'θ2' sin(Δθ) + θ2'² sin(Δθ)

So d/dt(∂T/∂θ1') = m L² (2θ1'' + θ2'' cos(Δθ) - θ1'θ2' sin(Δθ) + θ2'² sin(Δθ))

∂T/∂θ1 = m L² θ1'θ2' * d/dθ1 [cos(θ1-θ2)] = m L² θ1'θ2' * (-sin(θ1-θ2)) = -m L² θ1'θ2' sin(Δθ)

∂V/∂θ1 = 2mgL sin θ1

Euler-Lagrange: d/dt(∂T/∂θ1') - ∂T/∂θ1 + ∂V/∂θ1 = 0
m L² (2θ1'' + θ2'' cos(Δθ) - θ1'θ2' sin(Δθ) + θ2'² sin(Δθ)) - (-m L² θ1'θ2' sin(Δθ)) + 2mgL sin θ1 = 0
m L² (2θ1'' + θ2'' cos(Δθ) - θ1'θ2' sin(Δθ) + θ2'² sin(Δθ) + θ1'θ2' sin(Δθ)) + 2mgL sin θ1 = 0
m L² (2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ)) + 2mgL sin θ1 = 0

Divide by m L²:
2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0

This matches Wikipedia! The sign of θ2'² sin(Δθ) is POSITIVE. I had it wrong in my earlier derivation.

Let me re-derive the second equation too.

∂T/∂θ2' = m L² (θ2' + θ1' cos(θ1-θ2))
d/dt(∂T/∂θ2') = m L² (θ2'' + θ1'' cos(Δθ) + θ1' * d/dt[cos(Δθ)])
d/dt[cos(Δθ)] = -sin(Δθ) * (θ1' - θ2')
θ1' * d/dt[cos(Δθ)] = -θ1' sin(Δθ) (θ1' - θ2') = -θ1'² sin(Δθ) + θ1'θ2' sin(Δθ)

d/dt(∂T/∂θ2') = m L² (θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + θ1'θ2' sin(Δθ))

∂T/∂θ2 = m L² θ1'θ2' * d/dθ2 [cos(θ1-θ2)] = m L² θ1'θ2' * (-sin(θ1-θ2) * (-1)) = m L² θ1'θ2' sin(Δθ)

∂V/∂θ2 = mgL sin θ2

Euler-Lagrange: d/dt(∂T/∂θ2') - ∂T/∂θ2 + ∂V/∂θ2 = 0
m L² (θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + θ1'θ2' sin(Δθ)) - m L² θ1'θ2' sin(Δθ) + mgL sin θ2 = 0
m L² (θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ)) + mgL sin θ2 = 0

Divide by m L²:
θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + (g/L) sin θ2 = 0

This also matches Wikipedia!

So the correct equations are:
(A) 2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
(B) θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + (g/L) sin θ2 = 0

From (B): θ2'' = -θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2

Sub into (A):
2θ1'' + (-θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2) cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
2θ1'' - θ1'' cos²(Δθ) + θ1'² sin(Δθ) cos(Δθ) - (g/L) sin θ2 cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
θ1'' (2 - cos²(Δθ)) = -θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1
θ1'' = [ -θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1 ]

I notice there's a typo in my code. In the rk4Step function, I have `k4q1` and `k4q2` instead of `k4w1` and `k4w2`. Let me fix that.

Actually, looking more carefully at the code I wrote:
```javascript
w1: w1 + (dt/6)*(k1w1 + 2*k2w1 + 2*k3w1 + k4q1),
w2: w2 + (dt/6)*(k1w2 + 2*k2w2 + 2*k3w2 + k4q2)
```

This should be `k4w1` and `k4w2`. Let me fix this.

Also, I need to double-check the derivatives function. Let me re-derive the correct equations one more time.

From the Lagrangian for m1=m2=m, L1=L2=L:

Equations of motion (from Wikipedia and my corrected derivation):
(A) 2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
(B) θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + (g/L) sin θ2 = 0

From (B): θ2'' = -θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2

Sub into (A):
2θ1'' + (-θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2) cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0

2θ1'' - θ1'' cos²(Δθ) + θ1'² sin(Δθ) cos(Δθ) - (g/L) sin θ2 cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0

θ1'' (2 - cos²(Δθ)) = -θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1

θ1'' = [ -θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1 ] / (2 - cos²(Δθ))

Now for θ2'':
From (A): θ2'' cos(Δθ) = -2θ1'' - θ2'² sin(Δθ) - (2g/L) sin θ1
θ2'' = (-2θ1'' - θ2'² sin(Δθ) - (2g/L) sin θ1) / cos(Δθ)

But this has cos(Δθ) in the denominator, which could be zero. Let me use the other substitution.

From (A) again, solve for θ2'' differently. Actually, let me substitute θ1'' into (B):

From my expression for θ1'':
θ1'' = [ -w1² s c + (g/L) sin2 c - w2² s - 2(g/L) sin1 ] / den

where den = 2 - c², s = sin(Δθ), c = cos(Δθ), w1 = θ1', w2 = θ2', sin1 = sin(θ1), sin2 = sin(θ2)

Sub into (B):
θ2'' = -θ1'' c + w1² s - (g/L) sin2

θ2'' = -c * [ -w1² s c + (g/L) sin2 c - w2² s - 2(g/L) sin1 ] / den + w1² s - (g/L) sin2

= [ w1² s c² - (g/L) sin2 c² + w2² s c + 2(g/L) sin1 c ] / den + w1² s - (g/L) sin2

= [ w1² s c² - (g/L) sin2 c² + w2² s c + 2(g/L) sin1 c + w1² s den - (g/L) sin2 den ] / den

Now den = 2 - c², so:
w1² s den = w1² s (2 - c²) = 2w1² s - w1² s c²
-(g/L) sin2 den = -(g/L) sin2 (2 - c²) = -2(g/L) sin2 + (g/L) sin2 c²

So numerator:
w1² s c² - (g/L) sin2 c² + w2² s c + 2(g/L) sin1 c + 2w1² s - w1² s c² - 2(g/L) sin2 + (g/L) sin2 c²

The w1² s c² terms cancel: w1² s c² - w1² s c² = 0
The (g/L) sin2 c² terms cancel: -(g/L) sin2 c² + (g/L) sin2 c² = 0

Numerator = w2² s c + 2(g/L) sin1 c + 2w1² s - 2(g/L) sin2

θ2'' = [ 2w1² s + w2² s c + 2(g/L) sin1 c - 2(g/L) sin2 ] / den

= [ 2w1² sinΔ + w2² sinΔ cosΔ + 2(g/L) sin θ1 cosΔ - 2(g/L) sin θ2 ] / (2 - cos²Δ)

So:
θ1'' = [ -w1² s c + (g/L) sin2 c - w2² s - 2(g/L) sin1 ] / den
θ2'' = [ 2w1² s + w2² s c + 2(g/L) sin1 c - 2(g/L) sin2 ] / den

Let me verify with the equilibrium case (θ1=θ2=0, w1=w2=0):
s = 0, c = 1, sin1 = 0, sin2 = 0
θ1'' = [0 + 0 - 0 - 0] / 2 = 0 ✓
θ2'' = [0 + 0 + 0 - 0] / 2 = 0 ✓

Good. Now let me check my code. In my code I had:

```javascript
const a1 = (-w1*w1*s*c + (g/L)*sin2*c - w2*w2*s - 2*(g/L)*sin1) / den;
const a2 = (2*w1*w1*s*c - 2*(g/L)*sin2*c + 2*w2*w2*s + (g/L)*sin1*(4 - 2*c*c)) / den;
```

My a1 matches the correct θ1'' expression. ✓

My a2 is wrong. The correct expression is:
a2 = (2*w1*w1*s + w2*w2*s*c + 2*(g/L)*sin1*c - 2*(g/L)*sin2) / den

My code has:
a2 = (2*w1*w1*s*c - 2*(g/L)*sin2*c + 2*w2*w2*s + (g/L)*sin1*(4 - 2*c*c)) / den

These are different. Let me fix this.

Also, I have the typo k4q1/k4q2 instead of k4w1/k4w2.

Let me rewrite the corrected code.

Correct derivatives:
a1 = (-w1² s c + (g/L) sin2 c - w2² s - 2(g/L) sin1) / den
a2 = (2 w1² s + w2² s c + 2(g/L) sin1 c - 2(g/L) sin2) / den

Where:
s = sin(θ2 - θ1)
c = cos(θ2 - θ1)
den = 2 - c²

Let me fix the file.

Good, the fixes are applied. Now let me verify the file looks correct and then test it in the browser.

The code looks good. Now let me also update the comment in the derivatives function since I fixed the formula. The comment still has the old (incorrect) formula. Let me update it.

Actually, let me first test in the browser, then fix the comment if needed. Let me also update the comment to match the correct formula.

Actually, let me update the comment now and then test.

Wait, the comment says:
```
//   θ2'' = [ 2θ1'² sinΔ cosΔ - 2(g/L) sin θ2 cosΔ
//            + 2θ2'² sinΔ + (g/L) sin θ1 (4 - 2cos²Δ) ] / (2 - cos²Δ)
```

But the correct formula is:
```
θ2'' = [ 2θ1'² sinΔ + θ2'² sinΔ cosΔ + 2(g/L) sin θ1 cosΔ - 2(g/L) sin θ2 ] / (2 - cos²Δ)
```

Let me update the comment.

Now let me test in the browser.

The page loaded. Let me take a screenshot to see what it looks like.

The screenshot doesn't return an image. Let me use the motion action or try to inspect the canvas. Actually, let me try a different approach - use the browser's inspect action or try screenshot again.

Actually, the "images" field says "model does not support images". This means I can't see the screenshot. Let me try to inspect the DOM and check the canvas state programmatically.

Let me use the inspect action to see what's on the canvas.

I can't see images. Let me use the browser's console to evaluate the physics state and verify it's correct. Let me use the evaluate action or try to run some JavaScript in the page to check the state.

Actually, let me just check the DOM structure and verify the canvas has the right initial positions by evaluating JavaScript.

Let me try to use the motion action to see the pendulum in motion, or I can evaluate JS to check positions.

Actually, let me try evaluating some JS in the page context to verify the initial positions.

The physics checks out:
- Bob1 is at (547.2, 25.0) — to the right and above the pivot (as expected for θ1=120°)
- Bob2 is at (517.7, 192.4) — to the right and below the pivot (as expected)
- alpha1 = -10.02 rad/s² (negative, so θ1 decreases — swinging back toward 0) ✓
- alpha2 = -4.74 rad/s² (negative, so θ2 decreases — becoming more negative)

Wait, I expected alpha2 to be positive (θ2 should increase from -10° toward 0). But it's negative. Let me think about this...

θ2 = -10° is slightly to the left of downward. The second bob is at (0.692, 0.485) relative to pivot, which is to the right and below. But the second rod connects from Bob1 (which is at (0.866, -0.5)) to Bob2 (0.692, 0.485). 

The second rod goes from upper-right (Bob1) downward and slightly left to Bob2. The angle θ2 = -10° means the second rod is slightly to the left of straight down (from Bob1's perspective).

Gravity pulls Bob2 downward. But Bob2 is connected to Bob1, which is also moving. The torque on θ2 depends on the relative positions.

Let me think: θ2 = -10° means the second rod is slightly to the left of the downward direction from Bob1. So Bob2 is slightly to the left of where it would be if θ2 = 0 (straight down from Bob1). Gravity should pull it back toward θ2 = 0, which means θ2 should increase (become less negative). So alpha2 should be positive.

But my calculation shows alpha2 = -4.74 rad/s², which is negative. This is suspicious.

Let me recheck the derivative formula. Actually, let me reconsider.

θ2 = -10° from downward. This means the second rod points slightly to the left of straight down. To bring it back to straight down (θ2 = 0), θ2 needs to increase (become less negative). So alpha2 should be positive if gravity is the only force.

But there's also the coupling with θ1. The first pendulum is at θ1 = 120°, which is far from equilibrium. The first pendulum is swinging down (alpha1 < 0), and this affects the second pendulum through the coupling.

Let me check: at t=0, w1=w2=0. The coupling terms in a2 are:
a2 = (2*w1²*s + w2²*s*c + 2*(g/L)*sin1*c - 2*(g/L)*sin2) / den

With w1=w2=0:
a2 = (2*(g/L)*sin1*c - 2*(g/L)*sin2) / den
   = 2(g/L)(sin1*c - sin2) / den

sin1 = sin(120°) = 0.866
c = cos(θ2-θ1) = cos(-10° - 120°) = cos(-130°) = cos(130°) = -0.643
sin2 = sin(-10°) = -0.174

sin1*c - sin2 = 0.866*(-0.643) - (-0.174) = -0.557 + 0.174 = -0.383

a2 = 2*9.81*(-0.383) / den

den = 2 - c² = 2 - (-0.643)² = 2 - 0.413 = 1.587

a2 = 2*9.81*(-0.383) / 1.587 = -7.52 / 1.587 = -4.74 rad/s²

So alpha2 is negative because the coupling term (sin1*c) dominates and is negative. This makes physical sense: the first pendulum is at 120° (far from equilibrium), and its motion creates a strong coupling that initially pulls the second pendulum in the negative direction.

Actually, let me reconsider. The first pendulum is at 120°, swinging down. As it swings, it drags the second pendulum with it. The second pendulum's angle θ2 is measured from downward relative to the first bob. As the first bob swings from 120° toward smaller angles (toward the right and down), the second bob will be dragged along.

The negative alpha2 means the second pendulum initially rotates clockwise (more negative θ2), which means it swings further to the left relative to the first bob. This could be correct given the complex coupled dynamics.

Let me verify this is physically reasonable. The first bob is at upper-right (120° from downward). The second bob is at lower-right (slightly left of straight down from the first bob). As the first bob starts swinging down and to the left (toward smaller θ1), it will pull the second bob. The second bob might initially swing further left (more negative θ2) before the pendulum dynamics take over.

This seems physically plausible for a double pendulum with these initial conditions. The coupled motion is complex and non-intuitive.

So the physics is correct. Let me now check the simulation in the browser by running it for a few seconds and checking positions.

Actually, I can't see the browser visualization. Let me trust the physics and verify by running a simulation in Node.js and checking that the energy is conserved and the motion looks reasonable.

Let me run a quick simulation in Node.js to verify the physics.

I see a problem! The energy is NOT conserved. It changes from 0.149 to 4.815 to -6.16 to 9.17 to 7.33 to 21.78 to 7.92 to 2.56. This is a huge variation, indicating the integrator is unstable or there's a bug.

Wait, the "rel change" shows 0.0000% for all frames, but that's because I'm computing `(e/energy()-1)*100` where `energy()` is called again, which gives the current energy, not the initial energy. Let me fix that calculation.

But the absolute energy values clearly show it's not conserved: 0.149 → 4.815 → -6.16 → 9.17 → 7.33 → 21.78 → 7.92 → 2.56. This is a massive change.

Something is wrong. Let me check the energy calculation and the equations of motion.

Wait, the energy at t=0 is 0.149 J. But earlier I calculated the initial PE as:
V = -mgL(2cos θ1 + cos θ2) = -9.81*1*(2*cos(120°) + cos(-10°)) = -9.81*(2*(-0.5) + 0.9848) = -9.81*(-1 + 0.9848) = -9.81*(-0.0152) = 0.149 J

And T = 0 (since w1=w2=0). So E = 0.149 J. That's correct.

But after 30 frames (0.48s), the energy is 4.815 J. That's a huge increase. The pendulum shouldn't gain that much energy.

Let me check if there's a bug in my energy function or the equations of motion.

Energy function:
T = 0.5 * L² * (2w1² + w2² + 2w1*w2*cos(θ1-θ2))

Wait, I'm using m=1 implicitly. Let me check:
T = (1/2) * m * L² * (2θ1'² + θ2'² + 2θ1'θ2'cos(θ1-θ2))

With m=1, L=1:
T = 0.5 * (2w1² + w2² + 2w1*w2*cos(Δθ))

But my code has:
T = 0.5 * L*L * (2*w1*w1 + w2*w2 + 2*w1*w2*Math.cos(dTh))

That's the same thing. So the energy function is correct.

Let me check the equations of motion again. Maybe there's a sign error.

Actually, wait. Let me re-examine the equations more carefully. The energy is blowing up, which suggests the equations might be wrong.

Let me re-derive the equations very carefully.

Lagrangian for double pendulum with m1=m2=m, L1=L2=L:

Positions:
r1 = (L sin θ1, -L cos θ1)  [standard math coords, y up]
r2 = (L sin θ1 + L sin θ2, -L cos θ1 - L cos θ2)

Velocities:
v1 = (L θ1' cos θ1, L θ1' sin θ1)
v2 = (L θ1' cos θ1 + L θ2' cos θ2, L θ1' sin θ1 + L θ2' sin θ2)

v1² = L² θ1'² (cos² θ1 + sin² θ1) = L² θ1'²
v2² = L² [θ1' cos θ1 + θ2' cos θ2]² + L² [θ1' sin θ1 + θ2' sin θ2]²
    = L² [θ1'² cos² θ1 + 2θ1'θ2' cos θ1 cos θ2 + θ2'² cos² θ2
         + θ1'² sin² θ1 + 2θ1'θ2' sin θ1 sin θ2 + θ2'² sin² θ2]
    = L² [θ1'² (cos² θ1 + sin² θ1) + θ2'² (cos² θ2 + sin² θ2)
         + 2θ1'θ2' (cos θ1 cos θ2 + sin θ1 sin θ2)]
    = L² [θ1'² + θ2'² + 2θ1'θ2' cos(θ1 - θ2)]

T = (1/2) m (v1² + v2²) = (1/2) m L² [θ1'² + θ1'² + θ2'² + 2θ1'θ2' cos(θ1-θ2)]
  = (1/2) m L² [2θ1'² + θ2'² + 2θ1'θ2' cos(Δθ)]

V = m g y1 + m g y2 = m g (-L cos θ1) + m g (-L cos θ1 - L cos θ2)
  = -m g L (2 cos θ1 + cos θ2)

This matches what I had before. So the Lagrangian is correct.

Now, the Euler-Lagrange equations. Let me be very careful.

For θ1:
∂T/∂θ1' = m L² (2θ1' + θ2' cos(Δθ))
∂T/∂θ1 = m L² θ1' θ2' (-sin(Δθ)) * (1) = -m L² θ1' θ2' sin(Δθ)
  Wait, ∂T/∂θ1 = ∂/∂θ1 [0.5 m L² (2θ1'² + θ2'² + 2θ1'θ2' cos(θ1-θ2))]
  The only term with θ1 is the cos(θ1-θ2) term.
  ∂/∂θ1 [2θ1'θ2' cos(θ1-θ2)] = 2θ1'θ2' * (-sin(θ1-θ2)) * 1 = -2θ1'θ2' sin(Δθ)
  So ∂T/∂θ1 = 0.5 m L² * (-2θ1'θ2' sin(Δθ)) = -m L² θ1'θ2' sin(Δθ)

Hmm wait, I need to be careful. The Lagrangian is:
L = T - V = 0.5 m L² (2θ1'² + θ2'² + 2θ1'θ2' cos(θ1-θ2)) + m g L (2 cos θ1 + cos θ2)

∂L/∂θ1' = m L² (2θ1' + θ2' cos(θ1-θ2))
d/dt (∂L/∂θ1') = m L² (2θ1'' + θ2'' cos(θ1-θ2) + θ2' * d/dt[cos(θ1-θ2)])
d/dt[cos(θ1-θ2)] = -sin(θ1-θ2) * (θ1' - θ2')
So: d/dt (∂L/∂θ1') = m L² (2θ1'' + θ2'' cos(Δθ) - θ2' sin(Δθ)(θ1' - θ2'))
                    = m L² (2θ1'' + θ2'' cos(Δθ) - θ1'θ2' sin(Δθ) + θ2'² sin(Δθ))

∂L/∂θ1 = m L² θ1'θ2' (-sin(θ1-θ2)) + m g L (-2 sin θ1)
       = -m L² θ1'θ2' sin(Δθ) - 2 m g L sin θ1

Euler-Lagrange: d/dt(∂L/∂θ1') - ∂L/∂θ1 = 0
m L² (2θ1'' + θ2'' cos(Δθ) - θ1'θ2' sin(Δθ) + θ2'² sin(Δθ)) - (-m L² θ1'θ2' sin(Δθ) - 2 m g L sin θ1) = 0
m L² (2θ1'' + θ2'' cos(Δθ) - θ1'θ2' sin(Δθ) + θ2'² sin(Δθ) + θ1'θ2' sin(Δθ)) + 2 m g L sin θ1 = 0
m L² (2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ)) + 2 m g L sin θ1 = 0

Divide by m L²:
2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0   ... (1)

For θ2:
∂L/∂θ2' = m L² (θ2' + θ1' cos(θ1-θ2))
d/dt(∂L/∂θ2') = m L² (θ2'' + θ1'' cos(Δθ) + θ1' * d/dt[cos(Δθ)])
d/dt[cos(Δθ)] = -sin(Δθ) * (θ1' - θ2')
So: d/dt(∂L/∂θ2') = m L² (θ2'' + θ1'' cos(Δθ) - θ1' sin(Δθ)(θ1' - θ2'))
                    = m L² (θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + θ1'θ2' sin(Δθ))

∂L/∂θ2 = m L² θ1'θ2' (-sin(θ1-θ2)) * (-1) + m g L (-sin θ2)
Wait, let me redo this. ∂L/∂θ2:
The cos(θ1-θ2) term: ∂/∂θ2 [2θ1'θ2' cos(θ1-θ2)] = 2θ1'θ2' * (-sin(θ1-θ2)) * (-1) = 2θ1'θ2' sin(Δθ)
So ∂T/∂θ2 = 0.5 m L² * 2θ1'θ2' sin(Δθ) = m L² θ1'θ2' sin(Δθ)
∂V/∂θ2 = m g L (-sin θ2) = -m g L sin θ2
∂L/∂θ2 = m L² θ1'θ2' sin(Δθ) + m g L sin θ2

Wait, that's ∂L/∂θ2 = ∂T/∂θ2 - ∂V/∂θ2. But V = -mgL(2cos θ1 + cos θ2), so ∂V/∂θ2 = mgL sin θ2.
And L = T - V, so ∂L/∂θ2 = ∂T/∂θ2 - ∂V/∂θ2 = m L² θ1'θ2' sin(Δθ) - (-mgL sin θ2) = m L² θ1'θ2' sin(Δθ) + mgL sin θ2

Hmm wait. V = -mgL(2cos θ1 + cos θ2). So ∂V/∂θ2 = -mgL * (-sin θ2) = mgL sin θ2.
L = T - V, so ∂L/∂θ2 = ∂T/∂θ2 - ∂V/∂θ2.

∂T/∂θ2 = m L² θ1'θ2' sin(Δθ) (from above)
∂V/∂θ2 = mgL sin θ2

So ∂L/∂θ2 = m L² θ1'θ2' sin(Δθ) - mgL sin θ2

Euler-Lagrange: d/dt(∂L/∂θ2') - ∂L/∂θ2 = 0
m L² (θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + θ1'θ2' sin(Δθ)) - (m L² θ1'θ2' sin(Δθ) - mgL sin θ2) = 0
m L² (θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + θ1'θ2' sin(Δθ) - θ1'θ2' sin(Δθ)) + mgL sin θ2 = 0
m L² (θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ)) + mgL sin θ2 = 0

Divide by m L²:
θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + (g/L) sin θ2 = 0   ... (2)

So the equations are:
(1) 2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
(2) θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + (g/L) sin θ2 = 0

This matches what I had before. Good.

Now let me solve for θ1'' and θ2''.

From (2): θ2'' = -θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2

Sub into (1):
2θ1'' + (-θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2) cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
2θ1'' - θ1'' cos²(Δθ) + θ1'² sin(Δθ) cos(Δθ) - (g/L) sin θ2 cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
θ1'' (2 - cos²(Δθ)) = -θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1

θ1'' = [-θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1] / (2 - cos²(Δθ))

Now for θ2'':
From (1): θ2'' cos(Δθ) = -2θ1'' - θ2'² sin(Δθ) - (2g/L) sin θ1
θ2'' = [-2θ1'' - θ2'² sin(Δθ) - (2g/L) sin θ1] / cos(Δθ)

But this has cos(Δθ) in the denominator (could be zero). Let me use the other approach.

Sub θ1'' into (2):
θ2'' = -θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2

θ2'' = -cos(Δθ) * [-θ1'² s c + (g/L) sin2 c - θ2'² s - (2g/L) sin1] / den + θ1'² s - (g/L) sin2

where den = 2 - c², s = sin(Δθ), c = cos(Δθ), sin1 = sin(θ1), sin2 = sin(θ2)

= [θ1'² s c² - (g/L) sin2 c² + θ2'² s c + (2g/L) sin1 c] / den + θ1'² s - (g/L) sin2

= [θ1'² s c² - (g/L) sin2 c² + θ2'² s c + (2g/L) sin1 c + θ1'² s den - (g/L) sin2 den] / den

Now den = 2 - c²:
θ1'² s den = θ1'² s (2 - c²) = 2θ1'² s - θ1'² s c²
-(g/L) sin2 den = -(g/L) sin2 (2 - c²) = -2(g/L) sin2 + (g/L) sin2 c²

Numerator = θ1'² s c² - (g/L) sin2 c² + θ2'² s c + (2g/L) sin1 c + 2θ1'² s - θ1'² s c² - 2(g/L) sin2 + (g/L) sin2 c²

The θ1'² s c² terms: θ1'² s c² - θ1'² s c² = 0
The (g/L) sin2 c² terms: -(g/L) sin2 c² + (g/L) sin2 c² = 0

Numerator = θ2'² s c + (2g/L) sin1 c + 2θ1'² s - 2(g/L) sin2

θ2'' = [2θ1'² s + θ2'² s c + 2(g/L) sin1 c - 2(g/L) sin2] / den

So:
θ1'' = [-θ1'² s c + (g/L) sin2 c - θ2'² s - 2(g/L) sin1] / den
θ2'' = [2θ1'² s + θ2'² s c + 2(g/L) sin1 c - 2(g/L) sin2] / den

This matches my code! So the equations are correct.

But the energy is not conserved. Let me check the energy function again.

Energy = T + V
T = 0.5 m L² (2θ1'² + θ2'² + 2θ1'θ2' cos(Δθ))
V = -mgL (2 cos θ1 + cos θ2)

With m=1, L=1:
T = 0.5 (2w1² + w2² + 2w1 w2 cos(Δθ))
V = -g (2 cos θ1 + cos θ2)

Let me verify at t=0:
w1=0, w2=0 → T=0
θ1=120°, θ2=-10°
V = -9.81 * (2*cos(120°) + cos(-10°)) = -9.81 * (2*(-0.5) + 0.9848) = -9.81 * (-0.0152) = 0.149

E = 0 + 0.149 = 0.149. Correct.

After 30 frames (0.48s):
θ1=53.07°, θ2=-54.04°, w1=-4.444, w2=-3.943
T = 0.5 * (2*(-4.444)² + (-3.943)² + 2*(-4.444)*(-3.943)*cos(53.07°-(-54.04°)))
  = 0.5 * (2*19.75 + 15.55 + 2*17.52*cos(107.11°))
  = 0.5 * (39.50 + 15.55 + 35.04*(-0.294))
  = 0.5 * (55.05 - 10.30)
  = 0.5 * 44.75
  = 22.37

V = -9.81 * (2*cos(53.07°) + cos(-54.04°))
  = -9.81 * (2*0.601 + 0.588)
  = -9.81 * (1.202 + 0.588)
  = -9.81 * 1.790
  = -17.56

E = 22.37 - 17.56 = 4.81

That matches the output (4.815). But the initial energy was 0.149. So energy increased by a factor of 32!

This means the equations of motion are wrong, or the integrator is unstable.

Wait, let me double-check by computing the expected energy at t=0.48s if energy were conserved. If E = 0.149 J is conserved, then at any time T + V = 0.149.

But I'm getting E = 4.815 J. That's way off. So either:
1. The equations of motion are wrong
2. The RK4 implementation has a bug
3. The time step is too large

Let me check the RK4 implementation. Actually, let me first verify the equations by checking a simple case where I know the answer.

For a simple pendulum (θ2 = 0, w2 = 0, and the second mass doesn't move):
θ1'' = [-θ1'² sin(0-θ1) cos(0-θ1) + (g/L) sin(0) cos(0-θ1) - 0 - 2(g/L) sin θ1] / (2 - cos²(0-θ1))
     = [-θ1'² (-sin θ1) cos θ1 + 0 - 0 - 2(g/L) sin θ1] / (2 - cos² θ1)
     = [θ1'² sin θ1 cos θ1 - 2(g/L) sin θ1] / (2 - cos² θ1)

For a simple pendulum, θ1'' = -(g/L) sin θ1. But my formula gives something different because of the cos² term in the denominator. 

Hmm, for the double pendulum with θ2 = 0 and w2 = 0, the second mass is hanging straight down from the first mass. But the first mass is moving. The system is not a simple pendulum because the second mass adds inertia.

Actually, for a double pendulum with θ2 = 0 (second rod hanging straight down), the system behaves like a simple pendulum with effective mass 2m (since both masses move together). But the equation should be:
θ1'' = -(g/L) sin θ1 for a simple pendulum of length L and mass m.
For a double pendulum with θ2=0, the effective equation is:
(2m)L² θ1'' + ... = 0, which gives θ1'' = -(g/L) sin θ1 (since the second mass just adds to the inertia but doesn't change the torque).

Wait, let me check. If θ2 = 0 and w2 = 0, and we assume θ2 stays 0 (the second rod stays vertical), then:
From equation (2): θ2'' + θ1'' cos(θ1) - θ1'² sin(θ1) + (g/L) sin(0) = 0
θ2'' = -θ1'' cos(θ1) + θ1'² sin(θ1)

For θ2 to stay 0, we need θ2'' = 0, which means:
θ1'' cos(θ1) = θ1'² sin(θ1)
θ1'' = θ1'² tan(θ1)

But for a simple pendulum, θ1'' = -(g/L) sin θ1. These are different. So the assumption that θ2 stays 0 is wrong; the second pendulum will start moving.

This is correct: in a double pendulum, you can't assume one angle stays fixed. The coupling is real.

So the equations are correct. Let me check the RK4 implementation and the time step.

Actually, let me try a much smaller time step and see if energy is conserved better.

Let me also check if there's a bug in my RK4 implementation by testing with a simple harmonic oscillator.

Actually, let me first check the RK4 code for bugs.

Looking at my code:
```javascript
function rk4Step(dt) {
  const d1 = derivatives(t1, t2, w1, w2);
  const k1t1=w1, k1t2=w2, k1w1=d1.a1, k1w2=d1.a2;
  const d2 = derivatives(t1+.5*dt*k1t1, t2+.5*dt*k1t2, w1+.5*dt*k1w1, w2+.5*dt*k1w2);
  const k2t1=w1+.5*dt*k1w1, k2t2=w2+.5*dt*k1w2, k2w1=d2.a1, k2w2=d2.a2;
  const d3 = derivatives(t1+.5*dt*k2t1, t2+.5*dt*k2t2, w1+.5*dt*k2w1, w2+.5*dt*k2w2);
  const k3t1=w1+.5*dt*k2w1, k3t2=w2+.5*dt*k2w2, k3w1=d3.a1, k3w2=d3.a2;
  const d4 = derivatives(t1+dt*k3t1, t2+dt*k3t2, w1+dt*k3w1, w2+dt*k3w2);
  const k4t1=w1+dt*k3w1, k4t2=w2+dt*k3w2, k4w1=d4.a1, k4w2=d4.a2;
  t1 += (dt/6)*(k1t1+2*k2t1+2*k3t1+k4t1);
  t2 += (dt/6)*(k1t2+2*k2t2+2*k3t2+k4t2);
  w1 += (dt/6)*(k1w1+2*k2w1+2*k3w1+k4w1);
  w2 += (dt/6)*(k1w2+2*k2w2+2*k3w2+k4w2);
}
```

Wait, there's a bug! For k2, k3, k4, the "t1" and "t2" components should be the intermediate state values, not the velocity components.

For a state vector y = [θ1, θ2, ω1, ω2], the derivative is f(y) = [ω1, ω2, α1, α2].

RK4:
k1 = f(y0)
k2 = f(y0 + 0.5*dt*k1)
k3 = f(y0 + 0.5*dt*k2)
k4 = f(y0 + dt*k3)
y1 = y0 + (dt/6)*(k1 + 2*k2 + 2*k3 + k4)

So for k2:
y_temp = y0 + 0.5*dt*k1
= [θ1 + 0.5*dt*ω1, θ2 + 0.5*dt*ω2, ω1 + 0.5*dt*α1, ω2 + 0.5*dt*α2]

k2 = f(y_temp) = [ω1_temp, ω2_temp, α1_temp, α2_temp]
= [ω1 + 0.5*dt*α1, ω2 + 0.5*dt*α2, α1_temp, α2_temp]

So k2t1 = ω1 + 0.5*dt*α1 (this is the first component of k2)
k2t2 = ω2 + 0.5*dt*α2 (second component)
k2w1 = α1_temp (third component)
k2w2 = α2_temp (fourth component)

Looking at my code:
```javascript
const k2t1=w1+.5*dt*k1w1, k2t2=w2+.5*dt*k1w2, k2w1=d2.a1, k2w2=d2.a2;
```

k2t1 = w1 + 0.5*dt*k1w1 = ω1 + 0.5*dt*α1. This is correct (it's the first component of k2).
k2t2 = w2 + 0.5*dt*k1w2 = ω2 + 0.5*dt*α2. Correct.
k2w1 = d2.a1. This is α1 evaluated at the intermediate state. Correct.
k2w2 = d2.a2. Correct.

But then when I compute the final update:
```javascript
t1 += (dt/6)*(k1t1+2*k2t1+2*k3t1+k4t1);
```

k1t1 = ω1 (first component of k1)
k2t1 = ω1 + 0.5*dt*α1 (first component of k2)
k3t1 = ω1 + 0.5*dt*α1_temp2 (first component of k3, where α1_temp2 is from k2's state)
k4t1 = ω1 + dt*α1_temp3 (first component of k4)

The first component of the RK4 step is:
θ1_new = θ1 + (dt/6)*(ω1 + 2*(ω1+0.5*dt*α1) + 2*(ω1+0.5*dt*α1_temp2) + (ω1+dt*α1_temp3))

Hmm, this is... actually correct for RK4. The first component of the derivative is always the velocity, and RK4 handles this correctly.

Wait, but there's a subtlety. In the standard RK4 formulation:
y' = f(t, y)

For our system, y = [θ1, θ2, ω1, ω2] and f(y) = [ω1, ω2, α1(θ1,θ2,ω1,ω2), α2(θ1,θ2,ω1,ω2)].

This is an autonomous system (f doesn't depend on t explicitly). RK4 for autonomous systems:
k1 = f(y0)
k2 = f(y0 + (h/2)*k1)
k3 = f(y0 + (h/2)*k2)
k4 = f(y0 + h*k3)
y1 = y0 + (h/6)*(k1 + 2*k2 + 2*k3 + k4)

This is exactly what I'm doing. So the RK4 implementation looks correct.

Let me check if the issue is with the time step. DT = 0.0008, and I do 20 steps per frame, so 16ms of simulation per frame. At 60fps, that's roughly real-time.

Let me try with a much smaller DT and see if energy is conserved.

Actually, wait. Let me re-examine the energy computation in my Node.js test. I think there might be a bug in the energy function.

```javascript
function energy() {
  const dTh = t2 - t1;
  const T = 0.5 * L*L * (2*w1*w1 + w2*w2 + 2*w1*w2*Math.cos(dTh));
  const V = -g*L * (2*Math.cos(t1) + Math.cos(t2));
  return T + V;
}
```

With L=1, this is:
T = 0.5 * (2*w1² + w2² + 2*w1*w2*cos(t2-t1))
V = -g * (2*cos(t1) + cos(t2))

This looks correct.

Let me check the energy at frame 30 more carefully:
t1 = 53.07° = 0.926 rad
t2 = -54.04° = -0.943 rad
w1 = -4.444
w2 = -3.943
dTh = -54.04° - 53.07° = -107.11° = -1.869 rad
cos(dTh) = cos(-107.11°) = -0.294

T = 0.5 * (2*19.75 + 15.55 + 2*(-4.444)*(-3.943)*(-0.294))
  = 0.5 * (39.50 + 15.55 + 2*17.52*(-0.294))
  = 0.5 * (39.50 + 15.55 - 10.30)
  = 0.5 * 44.75
  = 22.375

V = -9.81 * (2*cos(53.07°) + cos(-54.04°))
  = -9.81 * (2*0.6007 + 0.5878)
  = -9.81 * (1.2014 + 0.5878)
  = -9.81 * 1.7892
  = -17.55

E = 22.375 - 17.55 = 4.825

But the output says 4.815. Close enough (rounding).

So the energy really is 4.815 J at frame 30, compared to 0.149 J at frame 0. That's a 32x increase. Something is definitely wrong.

Let me check if the equations of motion are truly correct by comparing with a known reference implementation.

Actually, let me look up the standard double pendulum equations from a reliable source. The Wikipedia article gives:

For the general case (m1, m2, L1, L2):
(m1+m2)L1²θ1'' + m2 L1 L2 θ2'' cos(θ1-θ2) + m2 L1 L2 θ2'² sin(θ1-θ2) + (m1+m2)g L1 sin θ1 = 0
m2 L2² θ2'' + m2 L1 L2 θ1'' cos(θ1-θ2) - m2 L1 L2 θ1'² sin(θ1-θ2) + m2 g L2 sin θ2 = 0

For m1=m2=m, L1=L2=L:
2mL²θ1'' + mL²θ2'' cos(Δθ) + mL²θ2'² sin(Δθ) + 2mgL sin θ1 = 0
mL²θ2'' + mL²θ1'' cos(Δθ) - mL²θ1'² sin(Δθ) + mgL sin θ2 = 0

Divide by mL²:
2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0  ... (A)
θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + (g/L) sin θ2 = 0  ... (B)

This matches my equations. So the equations are correct.

Now, the energy should be conserved. Let me check if there's a bug in the RK4 implementation by testing with a much smaller DT.

Actually, let me check the energy conservation with DT = 0.0001 (much smaller) and see if it's better.

The energy error is the same regardless of DT. This means the error is not from the time step — it's systematic. The final energy is always 2.561780 J regardless of the time step. This suggests there's a bug in the equations of motion or the energy function.

Wait, the final energy is exactly the same for all DT values. That's very suspicious. It suggests the simulation is converging to some state that has E = 2.56 J, regardless of the integration accuracy.

Let me check: after 600 frames with DT=0.0008 and 20 steps/frame, the total simulation time is 600 * 20 * 0.0008 = 9.6 seconds. With DT=0.00005 and 320 steps/frame, the total time is 600 * 320 * 0.00005 = 9.6 seconds. So the total simulation time is the same. But the energy should be conserved better with smaller DT.

But it's exactly the same! This suggests the error is not numerical but systematic — either the equations are wrong, or there's a bug in the implementation that doesn't depend on DT.

Let me check the equations one more time. Actually, let me look at the energy after just a few steps to see how quickly it diverges.

Actually, wait. Let me check if there's a bug in the energy function. The total simulation time is the same (9.6s) but the energy is the same for all DT. This is very strange.

Let me check: after 1 step with DT=0.0008 vs DT=0.00005, is the energy different?

Actually, let me check a simpler thing: let me verify the equations by computing the initial acceleration and checking if it matches a known result.

At t=0: θ1=120°, θ2=-10°, w1=0, w2=0.

θ1'' = [-0 + (g/L)*sin(-10°)*cos(-130°) - 0 - 2(g/L)*sin(120°)] / (2 - cos²(-130°))

cos(-130°) = cos(130°) = -0.6428
sin(-10°) = -0.1736
sin(120°) = 0.8660

θ1'' = [9.81*(-0.1736)*(-0.6428) - 2*9.81*0.8660] / (2 - 0.4132)
     = [9.81*0.1116 - 16.99] / 1.5868
     = [1.095 - 16.99] / 1.5868
     = -15.895 / 1.5868
     = -10.017 rad/s²

This matches my earlier calculation. So the equations give α1 = -10.017 rad/s².

For a simple pendulum of length L=1, at θ=120°, the acceleration would be:
α = -(g/L)*sin(120°) = -9.81*0.866 = -8.50 rad/s²

The double pendulum has α1 = -10.02 rad/s², which is more negative. This makes sense because the second pendulum adds coupling effects.

Now let me check the energy after 1 step with different DTs.

Actually, let me just compute the energy after a very small number of steps to see if the energy is conserved for small DT.

Interesting! After just 1 step, the energy is exactly conserved (delta = 0). But after 600 frames, it's wildly off. So the RK4 step itself is correct, but something goes wrong over many steps.

Wait, the delta after 1 step is exactly 0. That's suspicious. RK4 doesn't conserve energy exactly for a nonlinear system like the double pendulum. It should have a small error.

Let me check if there's a bug in the energy function or if the RK4 is doing something weird.

Actually, the delta being exactly 0 after 1 step might be a coincidence or there might be a bug. Let me check after more steps.

Hmm, wait. Let me re-examine. After 1 step:
- theta1 changed from 120° to 119.999816° (change of -0.000184° = -3.2e-6 rad)
- theta2 changed from -10° to -10.000087° (change of -0.000087° = -1.5e-6 rad)
- omega1 = -0.00801422
- omega2 = -0.00378866

The change in theta is w*dt = 0 * 0.0008 = 0 for the first step (since initial omega is 0). But after the step, omega is non-zero. The RK4 computes the average slope, so the theta change is (dt/6)*(k1 + 2k2 + 2k3 + k4) where k1=0 (initial omega=0).

k1t1 = w1 = 0
k2t1 = w1 + 0.5*dt*k1w1 = 0 + 0.5*0.0008*(-10.017) = -0.004007
k3t1 = w1 + 0.5*dt*k2w1 ≈ 0 + 0.5*0.0008*(-10.017) ≈ -0.004007 (roughly)
k4t1 = w1 + dt*k3w1 ≈ 0 + 0.0008*(-10.017) ≈ -0.008014

theta1 change = (0.0008/6)*(0 + 2*(-0.004007) + 2*(-0.004007) + (-0.008014))
= (0.0008/6)*(-0.008014 - 0.008014 - 0.008014)
= (0.0008/6)*(-0.024042)
= -0.00003206 rad
= -0.00184°

That matches the output (119.999816° = 120° - 0.000184°... wait, 0.00184° not 0.000184°).

Hmm, 0.00003206 rad * 180/π = 0.001837°. So theta1 should be 120 - 0.001837 = 119.998163°. But the output says 119.999816°. The difference is 0.000184°. That's off by a factor of 10.

Wait, let me recheck. 0.00003206 rad * 180/π = 0.001837°. But the output shows a change of 0.000184°. That's exactly 1/10 of what I calculated.

Hmm, let me recalculate more carefully.

k1w1 = d1.a1 = -10.017 rad/s² (from my earlier calculation)

k2t1 = w1 + 0.5*dt*k1w1 = 0 + 0.5*0.0008*(-10.017) = -0.0040068

For k2w1, I need to evaluate derivatives at the intermediate state:
t1_k2 = t1 + 0.5*dt*k1t1 = t1 + 0 = t1 (since k1t1 = w1 = 0)
t2_k2 = t2 + 0.5*dt*k1t2 = t2 + 0 = t2
w1_k2 = w1 + 0.5*dt*k1w1 = -0.0040068
w2_k2 = w2 + 0.5*dt*k1w2 = 0 + 0.5*0.0008*(-4.736) = -0.0018944

d2 = derivatives(t1, t2, w1_k2, w2_k2)
a1_k2 = (-w1_k2²*s*c + (g/L)*sin2*c - w2_k2²*s - 2(g/L)*sin1) / den

With w1_k2 = -0.004, w2_k2 = -0.002:
w1_k2² = 1.6e-5, w2_k2² = 3.6e-6
These are very small, so a1_k2 ≈ a1 = -10.017 (almost unchanged)

So k2w1 ≈ -10.017
k3t1 = w1 + 0.5*dt*k2w1 ≈ 0 + 0.5*0.0008*(-10.017) = -0.0040068
k3w1 ≈ -10.017 (similar)

k4t1 = w1 + dt*k3w1 ≈ 0 + 0.0008*(-10.017) = -0.008014
k4w1 ≈ -10.017

theta1_change = (dt/6)*(k1t1 + 2*k2t1 + 2*k3t1 + k4t1)
= (0.0008/6)*(0 + 2*(-0.004007) + 2*(-0.004007) + (-0.008014))
= (0.0008/6)*(-0.008014 - 0.008014 - 0.008014)
= (0.0008/6)*(-0.024042)
= 0.0001333 * (-0.024042)
= -0.000003206 rad

Wait, 0.0008/6 = 0.0001333. Multiply by -0.024042 = -0.000003206 rad.

In degrees: -0.000003206 * 180/π = -0.0001837°.

So theta1 changes by -0.000184°, from 120° to 119.999816°. That matches the output!

I made a calculation error earlier. The change is indeed -0.0000032 rad, not -0.000032 rad.

OK so the RK4 step is correct. The energy after 1 step is exactly conserved (to machine precision). That's actually expected for RK4 with a small enough time step for this system.

But after many steps, the energy diverges. This suggests the energy is not being conserved over long times, which is expected for RK4 (it's not symplectic). But the divergence seems too large.

Let me check after 100 steps, 1000 steps, etc.

Actually, let me check if the energy oscillates or monotonically diverges.

Wait, the initial energy shown is 0.15305453, not 0.14903594 as before. Something is different. Let me check...

Oh wait, the loop starts at step 0 and does 100 steps before printing. But the initial energy is printed at step 0 (before any steps). So the initial energy should be 0.14903594. But it shows 0.15305453.

Wait, the code does:
```
for (let step = 0; step <= 2000; step += 100) {
  for (let i = 0; i < 100; i++) rk4Step(DT);
  ...
}
```

At step=0, it does 100 RK4 steps first, then prints. So the "Step 0" output is actually after 100 steps, not before. Let me fix the test.

But more importantly, the energy is clearly not conserved. It grows from ~0.15 to ~5.0 and then becomes negative. This is a clear sign that the equations of motion are wrong or there's a bug in the implementation.

Let me check the equations one more time. I'll compare with a known correct implementation.

Actually, let me check a very simple case: the double pendulum with both angles near 0 (small oscillations). In this limit, the equations should decouple to two simple pendulums with frequencies sqrt(g/L) and sqrt(g/L * (1+1/2)) or something like that.

For small angles (θ1, θ2 ≈ 0):
sin θ ≈ θ, cos θ ≈ 1
Δθ = θ2 - θ1 ≈ small
sin(Δθ) ≈ Δθ, cos(Δθ) ≈ 1

a1 = (-w1²*Δθ*1 + (g/L)*θ2*1 - w2²*Δθ - 2(g/L)*θ1) / (2 - 1)
   = (-w1²*Δθ + (g/L)*θ2 - w2²*Δθ - 2(g/L)*θ1) / 1
   ≈ -(2g/L)*θ1 + (g/L)*θ2  (ignoring w² terms for small velocities)

a2 = (2*w1²*Δθ + w2²*Δθ*1 + 2(g/L)*θ1*1 - 2(g/L)*θ2) / (2 - 1)
   ≈ 2(g/L)*θ1 - 2(g/L)*θ2

So for small oscillations:
θ1'' = -(2g/L)*θ1 + (g/L)*θ2
θ2'' = 2(g/L)*θ1 - 2(g/L)*θ2

In matrix form:
[θ1'']   = (g/L) * [-2   1] [θ1]
[θ2'']              [ 2  -2] [θ2]

The eigenvalues of [-2 1; 2 -2] are:
det([-2-λ, 1; 2, -2-λ]) = (-2-λ)² - 2 = λ² + 4λ + 4 - 2 = λ² + 4λ + 2 = 0
λ = (-4 ± sqrt(16-8))/2 = (-4 ± sqrt(8))/2 = (-4 ± 2√2)/2 = -2 ± √2

λ1 = -2 + √2 ≈ -0.586
λ2 = -2 - √2 ≈ -3.414

So the normal mode frequencies are:
ω1 = sqrt((g/L)*|λ1|) = sqrt(9.81*0.586) = sqrt(5.75) = 2.40 rad/s
ω2 = sqrt((g/L)*|λ2|) = sqrt(9.81*3.414) = sqrt(33.49) = 5.79 rad/s

The frequencies of a double pendulum with equal masses and lengths are:
ω1 = sqrt((g/L)*(2-√2)) ≈ sqrt(9.81*0.586) = 2.40 rad/s
ω2 = sqrt((g/L)*(2+√2)) ≈ sqrt(9.81*3.414) = 5.79 rad/s

This matches! So the small-angle limit is correct.

Now, the equations look correct. So why is the energy diverging?

Let me check if there's a bug in the RK4 implementation. Actually, let me try a different integrator (like symplectic Euler or velocity Verlet) and see if energy is conserved better.

Actually, let me first check if the RK4 implementation has a bug by testing on a simple harmonic oscillator.

Wait, actually, I just realized something. Let me look at the energy function again:

```javascript
function energy() {
  const dTh = t2 - t1;
  const T = 0.5 * L*L * (2*w1*w1 + w2*w2 + 2*w1*w2*Math.cos(dTh));
  const V = -g*L * (2*Math.cos(t1) + Math.cos(t2));
  return T + V;
}
```

This is the energy for m=1, L=1. Let me verify:
T = 0.5 * m * L² * (2ω1² + ω2² + 2ω1ω2 cos(θ1-θ2))
With m=1, L=1: T = 0.5 * (2ω1² + ω2² + 2ω1ω2 cos(Δθ))

That matches.

V = -m*g*L*(2 cos θ1 + cos θ2)
With m=1, L=1: V = -g*(2 cos θ1 + cos θ2)

That matches.

So the energy function is correct.

Let me check the equations one more time by looking at the initial acceleration and comparing with a numerical derivative of the energy.

Actually, for a Hamiltonian system, dE/dt = 0. Let me check dE/dt analytically.

E = T + V
dE/dt = dT/dt + dV/dt

T = 0.5 * (2ω1² + ω2² + 2ω1ω2 cos(Δθ))
dT/dt = 2ω1*α1 + ω2*α2 + ω1*ω2*cos(Δθ)*(-sin(Δθ))*(ω1-ω2) + (ω1*ω2 + ω2*ω1)*cos(Δθ)*(ω1-ω2)... 

Hmm, this is getting complicated. Let me just trust the equations and check if there's a bug in the code.

Actually, wait. I just realized something. Let me re-examine the derivatives function:

```javascript
const a1 = (-w1*w1*s*c + (g/L)*sin2*c - w2*w2*s - 2*(g/L)*sin1) / den;
const a2 = (2*w1*w1*s + w2*w2*s*c + 2*(g/L)*sin1*c - 2*(g/L)*sin2) / den;
```

where s = sin(Δθ), c = cos(Δθ), den = 2 - c², sin1 = sin(θ1), sin2 = sin(θ2).

Let me verify this against the derived equations one more time:

θ1'' = [-θ1'² s c + (g/L) sin2 c - θ2'² s - 2(g/L) sin1] / den
θ2'' = [2θ1'² s + θ2'² s c + 2(g/L) sin1 c - 2(g/L) sin2] / den

My code:
a1 = (-w1*w1*s*c + (g/L)*sin2*c - w2*w2*s - 2*(g/L)*sin1) / den
a2 = (2*w1*w1*s + w2*w2*s*c + 2*(g/L)*sin1*c - 2*(g/L)*sin2) / den

These match. So the equations are correct.

Hmm, but the energy is clearly not conserved. Let me check if the issue is with the energy computation in the loop. Maybe the energy function is being called with the wrong state.

Actually, let me just try a completely different approach: use a symplectic integrator (like the implicit midpoint method or symplectic Euler) and see if energy is conserved better.

Or, let me try to find the bug by comparing with a known reference. Let me implement the equations from a different source.

From the paper "Double Pendulum" by various authors, or from the Wikipedia article, the equations for equal masses and lengths are:

θ1'' = ( -g*(2*sin(θ1) - sin(θ2)*cos(θ1-θ2)) - L*θ2'²*sin(θ1-θ2) ) / ( L*(2 - cos²(θ1-θ2)) )

θ2'' = ( 2*g*(sin(θ1)*cos(θ1-θ2) - sin(θ2)) + 2*L*θ1'²*sin(θ1-θ2) ) / ( L*(2 - cos²(θ1-θ2)) )

Wait, this is a different form! Let me check if this is equivalent to what I have.

From Wikipedia (https://en.wikipedia.org/wiki/Double_pendulum), the equations for the special case m1=m2, L1=L2 are sometimes given as:

θ1'' = [ -g (2 sin θ1 - sin(θ2) cos(θ1-θ2)) - L θ2'² sin(θ1-θ2) ] / [ L (2 - cos²(θ1-θ2)) ]

θ2'' = [ 2g (sin(θ1) cos(θ1-θ2) - sin(θ2)) + 2L θ1'² sin(θ1-θ2) ] / [ L (2 - cos²(θ1-θ2)) ]

Let me check if this is equivalent to my equations.

My θ1'': [-w1² s c + (g/L) sin2 c - w2² s - 2(g/L) sin1] / den

Multiply by L: L*θ1'' = [-L w1² s c + g sin2 c - L w2² s - 2g sin1] / den

Wikipedia θ1'': [-g (2 sin1 - sin2 c) - L w2² s] / (L * den)
= [-2g sin1 + g sin2 c - L w2² s] / (L * den)

Multiply by L: L*θ1''_wiki = [-2g sin1 + g sin2 c - L w2² s] / den

My L*θ1'': [-L w1² s c + g sin2 c - L w2² s - 2g sin1] / den

Comparing:
Wiki: -2g sin1 + g sin2 c - L w2² s
Mine: -L w1² s c + g sin2 c - L w2² s - 2g sin1

Difference: mine has an extra term "-L w1² s c".

So the Wikipedia formula doesn't have the -θ1'² sinΔ cosΔ term, but mine does!

Let me recheck my derivation. From the Euler-Lagrange equation:

2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0  ... (1)

Wait, this is the equation BEFORE solving for θ1'' and θ2''. Let me re-derive.

From (1): 2θ1'' + θ2'' cos(Δθ) = -θ2'² sin(Δθ) - (2g/L) sin θ1
From (2): θ1'' cos(Δθ) + θ2'' = θ1'² sin(Δθ) - (g/L) sin θ2

Let me solve this system again.

From (2): θ2'' = θ1'² sin(Δθ) - (g/L) sin θ2 - θ1'' cos(Δθ)

Sub into (1):
2θ1'' + (θ1'² sin(Δθ) - (g/L) sin θ2 - θ1'' cos(Δθ)) cos(Δθ) = -θ2'² sin(Δθ) - (2g/L) sin θ1
2θ1'' + θ1'² sin(Δθ) cos(Δθ) - (g/L) sin θ2 cos(Δθ) - θ1'' cos²(Δθ) = -θ2'² sin(Δθ) - (2g/L) sin θ1
θ1'' (2 - cos²(Δθ)) = -θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1

θ1'' = [-θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1] / (2 - cos²(Δθ))

So my equation has the -θ1'² s c term. Let me check the Wikipedia version again.

Actually, I think the Wikipedia version might be for a different convention or might be wrong. Let me check another source.

From the "Double Pendulum" Wikipedia article, the equations for general m1, m2, L1, L2 are:

(m1+m2)L1²θ1'' + m2 L1 L2 θ2'' cos(θ1-θ2) + m2 L1 L2 θ2'² sin(θ1-θ2) + (m1+m2)g L1 sin θ1 = 0

m2 L2² θ2'' + m2 L1 L2 θ1'' cos(θ1-θ2) - m2 L1 L2 θ1'² sin(θ1-θ2) + m2 g L2 sin θ2 = 0

For m1=m2=m, L1=L2=L:
2m L² θ1'' + m L² θ2'' cos(Δθ) + m L² θ2'² sin(Δθ) + 2m g L sin θ1 = 0
=> 2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0  ... (A)

m L² θ2'' + m L² θ1'' cos(Δθ) - m L² θ1'² sin(Δθ) + m g L sin θ2 = 0
=> θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + (g/L) sin θ2 = 0  ... (B)

This matches my equations. So the Wikipedia general equations are consistent with my derivation.

Now, solving (A) and (B) for θ1'' and θ2'':

From (B): θ2'' = -θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2

Sub into (A):
2θ1'' + (-θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2) cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
2θ1'' - θ1'' cos²(Δθ) + θ1'² sin(Δθ) cos(Δθ) - (g/L) sin θ2 cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
θ1'' (2 - cos²(Δθ)) = -θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1

θ1'' = [-θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1] / (2 - cos²(Δθ))

This is what I have. So my θ1'' is correct.

Now for θ2'':
From (A): θ2'' cos(Δθ) = -2θ1'' - θ2'² sin(Δθ) - (2g/L) sin θ1
θ2'' = (-2θ1'' - θ2'² sin(Δθ) - (2g/L) sin θ1) / cos(Δθ)

This has cos(Δθ) in the denominator. To avoid this, I substitute θ1'' into (B):

θ2'' = -θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2

= -cos(Δθ) * [-θ1'² s c + (g/L) sin2 c - θ2'² s - 2(g/L) sin1] / den + θ1'² s - (g/L) sin2

= [θ1'² s c² - (g/L) sin2 c² + θ2'² s c + 2(g/L) sin1 c] / den + θ1'² s - (g/L) sin2

= [θ1'² s c² - (g/L) sin2 c² + θ2'² s c + 2(g/L) sin1 c + θ1'² s (2 - c²) - (g/L) sin2 (2 - c²)] / den

= [θ1'² s c² - (g/L) sin2 c² + θ2'² s c + 2(g/L) sin1 c + 2θ1'² s - θ1'² s c² - 2(g/L) sin2 + (g/L) sin2 c²] / den

= [θ2'² s c + 2(g/L) sin1 c + 2θ1'² s - 2(g/L) sin2] / den

θ2'' = [2θ1'² s + θ2'² s c + 2(g/L) sin1 c - 2(g/L) sin2] / den

This matches my code. So both equations are correct.

But the energy is not conserved! This is very strange. Let me check if there's a bug in the RK4 code that's not immediately obvious.

Actually, let me try using a simple Euler integrator and see if the energy still diverges. If it does, the problem is with the equations. If it doesn't, the problem is with RK4.

Wait, I already tested Euler: after 1 Euler step, the energy changed by 0.00005, which is larger than RK4 (which was 0). So RK4 is better than Euler for one step.

But after many steps, RK4 energy diverges while Euler would diverge even faster. So the issue is the equations.

Actually wait, let me re-examine the energy after 1 RK4 step. It was exactly 0.14903594, same as the initial. That's suspicious. For a nonlinear system, RK4 shouldn't conserve energy exactly. Let me check with more precision.

Actually, the delta was 0.0000000000, which is 10 decimal places. That's suspicious. Let me check with more steps.

Let me try 10 steps and see if energy is still conserved.

Actually, I think the issue might be that I'm testing with the wrong energy function. Let me double-check the energy formula.

For the double pendulum with m1=m2=m, L1=L2=L:

T = 0.5*m*(v1² + v2²)

v1² = L²*θ1'²
v2² = L²*(θ1'² + θ2'² + 2*θ1'*θ2'*cos(θ1-θ2))

T = 0.5*m*L²*(θ1'² + θ1'² + θ2'² + 2*θ1'*θ2'*cos(θ1-θ2))
  = 0.5*m*L²*(2*θ1'² + θ2'² + 2*θ1'*θ2'*cos(θ1-θ2))

With m=1, L=1:
T = 0.5*(2*w1² + w2² + 2*w1*w2*cos(Δθ))

This matches my energy function. So the energy function is correct.

Hmm, let me try running the simulation with a different set of initial conditions (small angles) and see if energy is conserved.

Actually, let me try a completely different approach: let me check if the equations conserve energy by computing dE/dt analytically and checking if it's zero.

At t=0 (zero velocities), dE/dt = 0, which is correct (since w1=w2=0, the kinetic terms vanish and the potential derivative is also zero because... wait, it shouldn't be zero).

Actually, at t=0 with w1=w2=0:
dT/dt = 2*w1*a1 + w2*a2 + ... = 0 (since all terms have w1 or w2)
dV/dt = g*(2*sin(t1)*w1 + sin(t2)*w2) = 0 (since w1=w2=0)

So dE/dt = 0 at t=0 simply because the velocities are zero. This doesn't verify that the equations are correct.

At the arbitrary state (w1=2, w2=-1.5), dE/dt = 20.95, which is NOT zero. This means the equations do NOT conserve energy! There's a bug in the equations.

Let me re-derive the equations more carefully. I think there might be an error in my derivation or in the Wikipedia equations I'm using.

Let me re-derive from scratch, being very careful.

Lagrangian L = T - V

T = 0.5*m*(v1² + v2²)

Position of m1: r1 = (L*sin(θ1), -L*cos(θ1))
Velocity of m1: v1 = (L*θ1'*cos(θ1), L*θ1'*sin(θ1))
v1² = L²*θ1'²*(cos²(θ1) + sin²(θ1)) = L²*θ1'²

Position of m2: r2 = (L*sin(θ1) + L*sin(θ2), -L*cos(θ1) - L*cos(θ2))
Velocity of m2: v2 = (L*θ1'*cos(θ1) + L*θ2'*cos(θ2), L*θ1'*sin(θ1) + L*θ2'*sin(θ2))
v2² = L²*(θ1'*cos(θ1) + θ2'*cos(θ2))² + L²*(θ1'*sin(θ1) + θ2'*sin(θ2))²
    = L²*[θ1'²*cos²(θ1) + 2θ1'*θ2'*cos(θ1)*cos(θ2) + θ2'²*cos²(θ2)
         + θ1'²*sin²(θ1) + 2θ1'*θ2'*sin(θ1)*sin(θ2) + θ2'²*sin²(θ2)]
    = L²*[θ1'²*(cos²(θ1)+sin²(θ1)) + θ2'²*(cos²(θ2)+sin²(θ2))
         + 2θ1'*θ2'*(cos(θ1)*cos(θ2) + sin(θ1)*sin(θ2))]
    = L²*[θ1'² + θ2'² + 2θ1'*θ2'*cos(θ1-θ2)]

T = 0.5*m*L²*(θ1'² + θ1'² + θ2'² + 2θ1'*θ2'*cos(θ1-θ2))
  = 0.5*m*L²*(2θ1'² + θ2'² + 2θ1'*θ2'*cos(θ1-θ2))

V = m*g*y1 + m*g*y2
  = m*g*(-L*cos(θ1)) + m*g*(-L*cos(θ1) - L*cos(θ2))
  = -m*g*L*(2*cos(θ1) + cos(θ2))

L = T - V = 0.5*m*L²*(2θ1'² + θ2'² + 2θ1'*θ2'*cos(θ1-θ2)) + m*g*L*(2*cos(θ1) + cos(θ2))

Now, Euler-Lagrange equations:

d/dt(∂L/∂θ1') - ∂L/∂θ1 = 0

∂L/∂θ1' = m*L²*(2θ1' + θ2'*cos(θ1-θ2))

d/dt(∂L/∂θ1') = m*L²*(2θ1'' + θ2''*cos(θ1-θ2) + θ2'*(-sin(θ1-θ2))*(θ1'-θ2'))
               = m*L²*(2θ1'' + θ2''*cos(Δθ) - θ2'*sin(Δθ)*(θ1'-θ2'))
               = m*L²*(2θ1'' + θ2''*cos(Δθ) - θ1'*θ2'*sin(Δθ) + θ2'²*sin(Δθ))

∂L/∂θ1 = m*L²*θ1'*θ2'*(-sin(θ1-θ2))*(1) + m*g*L*(-2*sin(θ1))
        = -m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1)

Wait, I need to be careful with the sign. ∂L/∂θ1 = ∂T/∂θ1 - ∂V/∂θ1.

∂T/∂θ1 = 0.5*m*L²*2*θ1'*θ2'*(-sin(θ1-θ2)) = -m*L²*θ1'*θ2'*sin(Δθ)
∂V/∂θ1 = -m*g*L*2*(-sin(θ1)) = 2*m*g*L*sin(θ1)

So ∂L/∂θ1 = ∂T/∂θ1 - ∂V/∂θ1 = -m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1)

Wait, V = -m*g*L*(2*cos(θ1) + cos(θ2))
∂V/∂θ1 = -m*g*L*2*(-sin(θ1)) = 2*m*g*L*sin(θ1)

L = T - V, so ∂L/∂θ1 = ∂T/∂θ1 - ∂V/∂θ1 = -m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1)

Hmm, that's ∂L/∂θ1 = -m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1).

But earlier I had ∂L/∂θ1 = -m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1). Let me re-check.

Actually, V = -m*g*L*(2*cos(θ1) + cos(θ2))
∂V/∂θ1 = -m*g*L*2*(-sin(θ1)) = 2*m*g*L*sin(θ1)

So ∂L/∂θ1 = ∂T/∂θ1 - ∂V/∂θ1 = (-m*L²*θ1'*θ2'*sin(Δθ)) - (2*m*g*L*sin(θ1))
= -m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1)

Euler-Lagrange:
d/dt(∂L/∂θ1') - ∂L/∂θ1 = 0
m*L²*(2θ1'' + θ2''*cos(Δθ) - θ1'*θ2'*sin(Δθ) + θ2'²*sin(Δθ)) - (-m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1)) = 0
m*L²*(2θ1'' + θ2''*cos(Δθ) - θ1'*θ2'*sin(Δθ) + θ2'²*sin(Δθ) + θ1'*θ2'*sin(Δθ)) + 2*m*g*L*sin(θ1) = 0
m*L²*(2θ1'' + θ2''*cos(Δθ) + θ2'²*sin(Δθ)) + 2*m*g*L*sin(θ1) = 0

Divide by m*L²:
2θ1'' + θ2''*cos(Δθ) + θ2'²*sin(Δθ) + (2*g/L)*sin(θ1) = 0  ... (1)

This matches what I had before. Good.

Now for θ2:

∂L/∂θ2' = m*L²*(θ2' + θ1'*cos(θ1-θ2))

d/dt(∂L/∂θ2') = m*L²*(θ2'' + θ1''*cos(Δθ) + θ1'*(-sin(Δθ))*(θ1'-θ2'))
               = m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + θ1'*θ2'*sin(Δθ))

∂L/∂θ2 = ∂T/∂θ2 - ∂V/∂θ2

∂T/∂θ2 = 0.5*m*L²*2*θ1'*θ2'*(-sin(θ1-θ2))*(-1) = m*L²*θ1'*θ2'*sin(Δθ)
Wait, let me redo this. T = 0.5*m*L²*(2θ1'² + θ2'² + 2θ1'*θ2'*cos(θ1-θ2))
∂T/∂θ2 = 0.5*m*L²*2*θ1'*θ2'*(-sin(θ1-θ2))*(-1) = m*L²*θ1'*θ2'*sin(Δθ)
(The derivative of cos(θ1-θ2) with respect to θ2 is -sin(θ1-θ2)*(-1) = sin(θ1-θ2) = sin(Δθ))

∂V/∂θ2 = -m*g*L*(-sin(θ2)) = m*g*L*sin(θ2)

∂L/∂θ2 = m*L²*θ1'*θ2'*sin(Δθ) - m*g*L*sin(θ2)

Wait, ∂L/∂θ2 = ∂T/∂θ2 - ∂V/∂θ2. And V = -m*g*L*(2*cos(θ1) + cos(θ2)), so ∂V/∂θ2 = -m*g*L*(-sin(θ2)) = m*g*L*sin(θ2).

So ∂L/∂θ2 = m*L²*θ1'*θ2'*sin(Δθ) - m*g*L*sin(θ2)

Euler-Lagrange:
d/dt(∂L/∂θ2') - ∂L/∂θ2 = 0
m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + θ1'*θ2'*sin(Δθ)) - (m*L²*θ1'*θ2'*sin(Δθ) - m*g*L*sin(θ2)) = 0
m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + θ1'*θ2'*sin(Δθ) - θ1'*θ2'*sin(Δθ)) + m*g*L*sin(θ2) = 0
m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ)) + m*g*L*sin(θ2) = 0

Divide by m*L²:
θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + (g/L)*sin(θ2) = 0  ... (2)

This also matches what I had before. So the Euler-Lagrange equations are:

(1) 2θ1'' + θ2''*cos(Δθ) + θ2'²*sin(Δθ) + (2g/L)*sin(θ1) = 0
(2) θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + (g/L)*sin(θ2) = 0

Now, solving for θ1'' and θ2'':

From (2): θ2'' = -θ1''*cos(Δθ) + θ1'²*sin(Δθ) - (g/L)*sin(θ2)

Sub into (1):
2θ1'' + (-θ1''*cos(Δθ) + θ1'²*sin(Δθ) - (g/L)*sin(θ2))*cos(Δθ) + θ2'²*sin(Δθ) + (2g/L)*sin(θ1) = 0
2θ1'' - θ1''*cos²(Δθ) + θ1'²*sin(Δθ)*cos(Δθ) - (g/L)*sin(θ2)*cos(Δθ) + θ2'²*sin(Δθ) + (2g/L)*sin(θ1) = 0
θ1''*(2 - cos²(Δθ)) = -θ1'²*sin(Δθ)*cos(Δθ) + (g/L)*sin(θ2)*cos(Δθ) - θ2'²*sin(Δθ) - (2g/L)*sin(θ1)

θ1'' = [-θ1'²*sin(Δθ)*cos(Δθ) + (g/L)*sin(θ2)*cos(Δθ) - θ2'²*sin(Δθ) - (2g/L)*sin(θ1)] / (2 - cos²(Δθ))

Now for θ2'', substitute θ1'' into (2):
θ2'' = -θ1''*cos(Δθ) + θ1'²*sin(Δθ) - (g/L)*sin(θ2)

Let den = 2 - cos²(Δθ), s = sin(Δθ), c = cos(Δθ), sin1 = sin(θ1), sin2 = sin(θ2).

θ1'' = [-w1²*s*c + (g/L)*sin2*c - w2²*s - 2*(g/L)*sin1] / den

θ2'' = -c * [-w1²*s*c + (g/L)*sin2*c - w2²*s - 2*(g/L)*sin1] / den + w1²*s - (g/L)*sin2

= [w1²*s*c² - (g/L)*sin2*c² + w2²*s*c + 2*(g/L)*sin1*c] / den + w1²*s - (g/L)*sin2

= [w1²*s*c² - (g/L)*sin2*c² + w2²*s*c + 2*(g/L)*sin1*c + w1²*s*den - (g/L)*sin2*den] / den

den = 2 - c²
w1²*s*den = w1²*s*(2-c²) = 2*w1²*s - w1²*s*c²
-(g/L)*sin2*den = -(g/L)*sin2*(2-c²) = -2*(g/L)*sin2 + (g/L)*sin2*c²

Numerator = w1²*s*c² - (g/L)*sin2*c² + w2²*s*c + 2*(g/L)*sin1*c + 2*w1²*s - w1²*s*c² - 2*(g/L)*sin2 + (g/L)*sin2*c²

The w1²*s*c² terms cancel: w1²*s*c² - w1²*s*c² = 0
The (g/L)*sin2*c² terms cancel: -(g/L)*sin2*c² + (g/L)*sin2*c² = 0

Numerator = w2²*s*c + 2*(g/L)*sin1*c + 2*w1²*s - 2*(g/L)*sin2

θ2'' = [2*w1²*s + w2²*s*c + 2*(g/L)*sin1*c - 2*(g/L)*sin2] / den

This is exactly what I have in my code. So the equations are correct.

But dE/dt is not zero at the arbitrary state! There must be a bug in the dE/dt calculation or in the equations.

Let me recompute dE/dt more carefully.

E = T + V
T = 0.5*(2*w1² + w2² + 2*w1*w2*cD)  where cD = cos(θ2-θ1)
V = -g*(2*cos(θ1) + cos(θ2))

dE/dt = dT/dt + dV/dt

dT/dt = 0.5 * d/dt[2*w1² + w2² + 2*w1*w2*cos(θ2-θ1)]
       = 0.5 * [4*w1*w1' + 2*w2*w2' + 2*w1'*w2*cos(θ2-θ1) + 2*w1*w2'*cos(θ2-θ1) + 2*w1*w2*(-sin(θ2-θ1))*(θ2'-θ1')]
       = 2*w1*a1 + w2*a2 + w1'*w2*cos(Δθ) + w1*w2'*cos(Δθ) - w1*w2*sin(Δθ)*(w2-w1)
       = 2*w1*a1 + w2*a2 + a1*w2*cos(Δθ) + w1*a2*cos(Δθ) - w1*w2*sin(Δθ)*(w2-w1)

dV/dt = -g * [-2*sin(θ1)*θ1' - sin(θ2)*θ2']
       = g * [2*sin(θ1)*w1 + sin(θ2)*w2]

So dE/dt = 2*w1*a1 + w2*a2 + a1*w2*cos(Δθ) + w1*a2*cos(Δθ) - w1*w2*sin(Δθ)*(w2-w1) + g*(2*sin(θ1)*w1 + sin(θ2)*w2)

Let me compute this at the arbitrary state and see if it's zero.

At θ1=1, θ2=-0.5, w1=2, w2=-1.5:
Δθ = -0.5 - 1 = -1.5
s = sin(-1.5) = -0.9975
c = cos(-1.5) = 0.0707
den = 2 - c² = 2 - 0.005 = 1.995

a1 = [-4*(-0.9975)*0.0707 + 9.81*sin(-0.5)*0.0707 - 2.25*(-0.9975) - 2*9.81*sin(1)] / 1.995
   = [0.2818 + 9.81*(-0.4794)*0.0707 - (-2.244) - 19.62*0.8415] / 1.995
   = [0.2818 - 0.3325 + 2.244 - 16.509] / 1.995
   = [-14.316] / 1.995
   = -7.177

a2 = [2*4*(-0.9975) + 2.25*(-0.9975)*0.0707 + 2*9.81*sin(1)*0.0707 - 2*9.81*sin(-0.5)] / 1.995
   = [-7.98 - 0.159 + 2*9.81*0.8415*0.0707 - 19.62*(-0.4794)] / 1.995
   = [-7.98 - 0.159 + 1.171 + 9.406] / 1.995
   = [2.438] / 1.995
   = 1.222

Now dE/dt:
2*w1*a1 = 2*2*(-7.177) = -28.708
w2*a2 = (-1.5)*1.222 = -1.833
a1*w2*c = (-7.177)*(-1.5)*0.0707 = 0.759
w1*a2*c = 2*1.222*0.0707 = 0.173
-w1*w2*s*(w2-w1) = -2*(-1.5)*(-0.9975)*(-1.5-2) = 3*(-0.9975)*(-3.5) = 10.474

Wait, let me redo this carefully:
-w1*w2*sin(Δθ)*(w2-w1) = -2*(-1.5)*(-0.9975)*(-1.5-2) = 3*(-0.9975)*(-3.5) = 3*3.491 = 10.474

dT/dt = -28.708 - 1.833 + 0.759 + 0.173 + 10.474 = -19.135

dV/dt = 9.81 * (2*sin(1)*2 + sin(-0.5)*(-1.5))
       = 9.81 * (4*0.8415 + (-0.4794)*(-1.5))
       = 9.81 * (3.366 + 0.719)
       = 9.81 * 4.085
       = 40.074

dE/dt = -19.135 + 40.074 = 20.939

This matches the computed value of 20.947. So dE/dt is definitely not zero.

This means the equations of motion do NOT conserve energy. There must be an error in the equations.

Let me re-examine the Euler-Lagrange equations. Maybe I made a sign error somewhere.

Actually, let me check the Wikipedia equations again. From Wikipedia:

(m1+m2)L1²θ1'' + m2 L1 L2 θ2'' cos(θ1-θ2) + m2 L1 L2 θ2'² sin(θ1-θ2) + (m1+m2)g L1 sin θ1 = 0

For m1=m2=m, L1=L2=L:
2mL²θ1'' + mL²θ2'' cos(Δθ) + mL²θ2'² sin(Δθ) + 2mgL sin θ1 = 0
=> 2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0  ... (A)

m2 L2² θ2'' + m2 L1 L2 θ1'' cos(θ1-θ2) - m2 L1 L2 θ1'² sin(θ1-θ2) + m2 g L2 sin θ2 = 0

For m1=m2=m, L1=L2=L:
mL²θ2'' + mL²θ1'' cos(Δθ) - mL²θ1'² sin(Δθ) + mgL sin θ2 = 0
=> θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + (g/L) sin θ2 = 0  ... (B)

These match my equations. So the Wikipedia equations are the same as mine.

But the energy is not conserved! This is very strange. Let me check if the Wikipedia equations are actually correct by testing energy conservation with a known-correct implementation.

Actually, let me look for a different source. Let me check the "Double pendulum" article more carefully. Maybe the equations on Wikipedia are for a different convention.

Wait, I just realized something. The Wikipedia equations might be using a different sign convention for the angles. Let me check.

In my convention:
- θ is measured from the downward vertical
- Positive θ means the pendulum swings to the right (counterclockwise from downward in standard math, but in canvas coords with y-down, it's clockwise)
- Position: (L sin θ, L cos θ) in canvas coords (y-down)

In the Wikipedia article, they might use:
- θ measured from the downward vertical
- Position: (L sin θ, -L cos θ) in standard math coords (y-up)

The equations should be the same regardless of the coordinate system, as long as the Lagrangian is derived consistently.

Let me re-derive the Lagrangian in standard math coords (y-up) and see if I get the same equations.

Standard math coords (y-up):
r1 = (L sin θ1, -L cos θ1)
r2 = (L sin θ1 + L sin θ2, -L cos θ1 - L cos θ2)

v1 = (L θ1' cos θ1, L θ1' sin θ1)
v1² = L² θ1'²

v2 = (L θ1' cos θ1 + L θ2' cos θ2, L θ1' sin θ1 + L θ2' sin θ2)
v2² = L²(θ1' cos θ1 + θ2' cos θ2)² + L²(θ1' sin θ1 + θ2' sin θ2)²

This is the same as before. So the Lagrangian is the same.

Hmm, let me try a completely different approach. Let me look up the equations from a reliable source and compare.

Actually, I just realized I should check the equations from the "Double pendulum" Wikipedia article more carefully. The article gives:

θ1'' = ( -g (2 sin θ1 - sin(θ2) cos(θ1-θ2)) - L θ2'² sin(θ1-θ2) ) / ( L (2 - cos²(θ1-θ2)) )

θ2'' = ( 2g (sin(θ1) cos(θ1-θ2) - sin(θ2)) + 2L θ1'² sin(θ1-θ2) ) / ( L (2 - cos²(θ1-θ2)) )

Let me check if THIS form conserves energy.

θ1''_wiki = [-g(2 sin1 - sin2 c) - L w2² s] / (L * den)
           = [-2g sin1 + g sin2 c - L w2² s] / (L * den)

θ2''_wiki = [2g(sin1 c - sin2) + 2L w1² s] / (L * den)
           = [2g sin1 c - 2g sin2 + 2L w1² s] / (L * den)

My θ1'' = [-w1² s c + (g/L) sin2 c - w2² s - 2(g/L) sin1] / den
         = [-L w1² s c + g sin2 c - L w2² s - 2g sin1] / (L * den)

My θ2'' = [2 w1² s + w2² s c + 2(g/L) sin1 c - 2(g/L) sin2] / den
         = [2 L w1² s + L w2² s c + 2g sin1 c - 2g sin2] / (L * den)

Comparing:
Wiki θ1'': -2g sin1 + g sin2 c - L w2² s
My θ1'': -L w1² s c + g sin2 c - L w2² s - 2g sin1

Difference: my θ1'' has an extra term -L w1² s c.

Wiki θ2'': 2g sin1 c - 2g sin2 + 2L w1² s
My θ2'': 2 L w1² s + L w2² s c + 2g sin1 c - 2g sin2

Difference: my θ2'' has an extra term L w2² s c.

So the Wikipedia equations don't have the velocity-squared coupling terms (-L w1² s c and L w2² s c). These terms come from the d/dt of the cos(Δθ) term in the Lagrangian.

Wait, but I derived these terms from the Euler-Lagrange equations. Let me re-check.

Actually, I think the issue might be that the Wikipedia equations are for a different system or have a typo. Let me check the source of the Wikipedia equations.

Actually, let me just test the Wikipedia equations and see if they conserve energy.

Neither form conserves energy! Both have dE/dt ≈ 20. So both equations are wrong, or my dE/dt calculation is wrong.

Let me re-derive dE/dt very carefully.

E = T + V
T = 0.5 * (2*w1² + w2² + 2*w1*w2*cos(θ2-θ1))
V = -g * (2*cos(θ1) + cos(θ2))

Let me compute dE/dt step by step.

Let u = θ2 - θ1 = Δθ
du/dt = w2 - w1

T = 0.5 * (2*w1² + w2² + 2*w1*w2*cos(u))

dT/dt = 0.5 * d/dt[2*w1² + w2² + 2*w1*w2*cos(u)]
       = 0.5 * [4*w1*w1' + 2*w2*w2' + 2*(d/dt(w1*w2*cos(u)))]

d/dt(w1*w2*cos(u)) = w1'*w2*cos(u) + w1*w2'*cos(u) + w1*w2*(-sin(u))*du/dt
                    = a1*w2*cos(u) + w1*a2*cos(u) - w1*w2*sin(u)*(w2-w1)

So dT/dt = 0.5 * [4*w1*a1 + 2*w2*a2 + 2*(a1*w2*cos(u) + w1*a2*cos(u) - w1*w2*sin(u)*(w2-w1))]
         = 2*w1*a1 + w2*a2 + a1*w2*cos(u) + w1*a2*cos(u) - w1*w2*sin(u)*(w2-w1)

V = -g * (2*cos(θ1) + cos(θ2))
dV/dt = -g * (-2*sin(θ1)*w1 - sin(θ2)*w2)
       = g * (2*sin(θ1)*w1 + sin(θ2)*w2)

So dE/dt = 2*w1*a1 + w2*a2 + a1*w2*cos(u) + w1*a2*cos(u) - w1*w2*sin(u)*(w2-w1) + g*(2*sin(θ1)*w1 + sin(θ2)*w2)

Now, for the equations to conserve energy, we need dE/dt = 0. Let me substitute the equations of motion and check.

From the Lagrangian:
(1) 2θ1'' + θ2''*cos(u) + θ2'²*sin(u) + (2g/L)*sin(θ1) = 0
(2) θ2'' + θ1''*cos(u) - θ1'²*sin(u) + (g/L)*sin(θ2) = 0

Multiply (1) by w1 and (2) by w2:
2*w1*θ1'' + w1*θ2''*cos(u) + w1*θ2'²*sin(u) + (2g/L)*w1*sin(θ1) = 0
w2*θ2'' + w2*θ1''*cos(u) - w2*θ1'²*sin(u) + (g/L)*w2*sin(θ2) = 0

Add them:
2*w1*a1 + w2*a2 + (w1*a2 + w2*a1)*cos(u) + w1*w2²*sin(u) - w2*w1²*sin(u) + (2g/L)*w1*sin(θ1) + (g/L)*w2*sin(θ2) = 0

Multiply by L (since my energy has g, not g/L... wait, my energy uses g, not g/L. Let me check).

Actually, my equations use g/L for the gravitational terms. Let me rewrite.

From (1): 2*a1 + a2*c + w2²*s + (2g/L)*sin1 = 0
Multiply by w1: 2*w1*a1 + w1*a2*c + w1*w2²*s + (2g/L)*w1*sin1 = 0  ... (A)

From (2): a2 + a1*c - w1²*s + (g/L)*sin2 = 0
Multiply by w2: w2*a2 + w2*a1*c - w2*w1²*s + (g/L)*w2*sin2 = 0  ... (B)

Add (A) and (B):
2*w1*a1 + w2*a2 + w1*a2*c + w2*a1*c + w1*w2²*s - w2*w1²*s + (2g/L)*w1*sin1 + (g/L)*w2*sin2 = 0

Now, the dE/dt expression is:
dE/dt = 2*w1*a1 + w2*a2 + a1*w2*c + w1*a2*c - w1*w2*s*(w2-w1) + g*(2*sin1*w1 + sin2*w2)

Let me rearrange the terms:
dE/dt = (2*w1*a1 + w2*a2 + w1*a2*c + w2*a1*c) - w1*w2*s*(w2-w1) + g*(2*sin1*w1 + sin2*w2)

From (A)+(B):
2*w1*a1 + w2*a2 + w1*a2*c + w2*a1*c = -w1*w2²*s + w2*w1²*s - (2g/L)*w1*sin1 - (g/L)*w2*sin2

So:
dE/dt = [-w1*w2²*s + w2*w1²*s - (2g/L)*w1*sin1 - (g/L)*w2*sin2] - w1*w2*s*(w2-w1) + g*(2*sin1*w1 + sin2*w2)

Now, -w1*w2²*s + w2*w1²*s = w1*w2*s*(w1 - w2) = -w1*w2*s*(w2 - w1)

So:
dE/dt = [-w1*w2*s*(w2-w1) - (2g/L)*w1*sin1 - (g/L)*w2*sin2] - w1*w2*s*(w2-w1) + g*(2*sin1*w1 + sin2*w2)

= -2*w1*w2*s*(w2-w1) - (2g/L)*w1*sin1 - (g/L)*w2*sin2 + 2g*w1*sin1 + g*w2*sin2

Now, group the gravitational terms:
g*w1*sin1*(2 - 2/L) + g*w2*sin2*(1 - 1/L)

Wait, this doesn't look right. Let me be more careful with the units.

My equations use g/L for the gravitational acceleration term. But my energy uses g directly (since V = -g*L*(2*cos(θ1) + cos(θ2)) with L=1, so V = -g*(2*cos(θ1) + cos(θ2))).

Actually, with L=1, g/L = g. So the equations simplify to:
(1) 2*a1 + a2*c + w2²*s + 2*g*sin1 = 0
(2) a2 + a1*c - w1²*s + g*sin2 = 0

And dE/dt:
dE/dt = 2*w1*a1 + w2*a2 + a1*w2*c + w1*a2*c - w1*w2*s*(w2-w1) + g*(2*sin1*w1 + sin2*w2)

From (A)+(B) with L=1:
2*w1*a1 + w2*a2 + w1*a2*c + w2*a1*c = -w1*w2²*s + w2*w1²*s - 2*g*w1*sin1 - g*w2*sin2

dE/dt = [-w1*w2²*s + w2*w1²*s - 2*g*w1*sin1 - g*w2*sin2] - w1*w2*s*(w2-w1) + g*(2*w1*sin1 + w2*sin2)

= -w1*w2²*s + w2*w1²*s - 2*g*w1*sin1 - g*w2*sin2 - w1*w2*s*(w2-w1) + 2*g*w1*sin1 + g*w2*sin2

The gravitational terms cancel: -2*g*w1*sin1 + 2*g*w1*sin1 = 0, and -g*w2*sin2 + g*w2*sin2 = 0.

dE/dt = -w1*w2²*s + w2*w1²*s - w1*w2*s*(w2-w1)
       = w1*w2*s*(-w2 + w1 - w2 + w1)
       = w1*w2*s*(2*w1 - 2*w2)
       = 2*w1*w2*s*(w1 - w2)

This is NOT zero in general! So the equations do NOT conserve energy.

This is a major discovery. The equations I derived from the Euler-Lagrange equations do not conserve energy. But the Euler-Lagrange equations should give equations that conserve energy (since the Lagrangian has no explicit time dependence).

There must be an error in my derivation. Let me re-derive very carefully.

Actually, wait. Let me re-examine the Euler-Lagrange equations. I think the issue might be with the ∂V/∂θ terms.

V = -m*g*L*(2*cos(θ1) + cos(θ2))

∂L/∂θ1 = ∂T/∂θ1 - ∂V/∂θ1

∂T/∂θ1: 
T = 0.5*m*L²*(2*θ1'² + θ2'² + 2*θ1'*θ2'*cos(θ1-θ2))
∂T/∂θ1 = 0.5*m*L²*2*θ1'*θ2'*(-sin(θ1-θ2)) = -m*L²*θ1'*θ2'*sin(θ1-θ2) = -m*L²*θ1'*θ2'*sin(Δθ)

∂V/∂θ1 = -m*g*L*2*(-sin(θ1)) = 2*m*g*L*sin(θ1)

∂L/∂θ1 = -m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1)

Wait, ∂V/∂θ1 should be:
V = -m*g*L*(2*cos(θ1) + cos(θ2))
∂V/∂θ1 = -m*g*L*2*(-sin(θ1)) = 2*m*g*L*sin(θ1)

So ∂L/∂θ1 = ∂T/∂θ1 - ∂V/∂θ1 = -m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1)

This is correct. And the Euler-Lagrange equation gives:
d/dt(∂L/∂θ1') - ∂L/∂θ1 = 0
m*L²*(2*θ1'' + θ2''*cos(Δθ) - θ1'*θ2'*sin(Δθ) + θ2'²*sin(Δθ)) - (-m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1)) = 0
m*L²*(2*θ1'' + θ2''*cos(Δθ) - θ1'*θ2'*sin(Δθ) + θ2'²*sin(Δθ) + θ1'*θ2'*sin(Δθ)) + 2*m*g*L*sin(θ1) = 0
m*L²*(2*θ1'' + θ2''*cos(Δθ) + θ2'²*sin(Δθ)) + 2*m*g*L*sin(θ1) = 0

Divide by m*L²:
2*θ1'' + θ2''*cos(Δθ) + θ2'²*sin(Δθ) + (2*g/L)*sin(θ1) = 0

This is correct.

For θ2:
∂T/∂θ2 = 0.5*m*L²*2*θ1'*θ2'*(-sin(θ1-θ2))*(-1) = m*L²*θ1'*θ2'*sin(Δθ)
Wait, ∂T/∂θ2 = ∂/∂θ2 [0.5*m*L²*(2*θ1'² + θ2'² + 2*θ1'*θ2'*cos(θ1-θ2))]
The only term with θ2 is the cos(θ1-θ2) term.
∂/∂θ2 [2*θ1'*θ2'*cos(θ1-θ2)] = 2*θ1'*θ2'*(-sin(θ1-θ2))*(-1) = 2*θ1'*θ2'*sin(θ1-θ2) = 2*θ1'*θ2'*sin(Δθ)
So ∂T/∂θ2 = 0.5*m*L²*2*θ1'*θ2'*sin(Δθ) = m*L²*θ1'*θ2'*sin(Δθ)

∂V/∂θ2 = -m*g*L*(-sin(θ2)) = m*g*L*sin(θ2)

∂L/∂θ2 = ∂T/∂θ2 - ∂V/∂θ2 = m*L²*θ1'*θ2'*sin(Δθ) - m*g*L*sin(θ2)

d/dt(∂L/∂θ2') = m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + θ1'*θ2'*sin(Δθ))

Euler-Lagrange:
m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + θ1'*θ2'*sin(Δθ)) - (m*L²*θ1'*θ2'*sin(Δθ) - m*g*L*sin(θ2)) = 0
m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + θ1'*θ2'*sin(Δθ) - θ1'*θ2'*sin(Δθ)) + m*g*L*sin(θ2) = 0
m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ)) + m*g*L*sin(θ2) = 0

Divide by m*L²:
θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + (g/L)*sin(θ2) = 0

This is also correct.

So the Euler-Lagrange equations are:
(1) 2*θ1'' + θ2''*cos(Δθ) + θ2'²*sin(Δθ) + (2*g/L)*sin(θ1) = 0
(2) θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + (g/L)*sin(θ2) = 0

And I showed that these equations lead to dE/dt = 2*w1*w2*s*(w1-w2), which is NOT zero.

This means there's an inconsistency. The Euler-Lagrange equations should conserve energy (since ∂L/∂t = 0). But my derivation shows they don't.

Let me re-check the dE/dt derivation. Maybe I made an error.

Actually, I think I found the issue. When I compute dE/dt, I use:
dE/dt = dT/dt + dV/dt

But the Euler-Lagrange equations are:
d/dt(∂L/∂θ1') - ∂L/∂θ1 = 0
d/dt(∂L/∂θ2') - ∂L/∂θ2 = 0

These can be written as:
d/dt(∂T/∂θ1') - ∂T/∂θ1 + ∂V/∂θ1 = 0  (since L = T - V)
d/dt(∂T/∂θ2') - ∂T/∂θ2 + ∂V/∂θ2 = 0

Now, dE/dt = d/dt(T + V) = dT/dt + dV/dt

dT/dt = ∂T/∂θ1 * θ1' + ∂T/∂θ2 * θ2' + ∂T/∂θ1' * θ1'' + ∂T/∂θ2' * θ2'
       = ∂T/∂θ1 * w1 + ∂T/∂θ2 * w2 + ∂T/∂θ1' * a1 + ∂T/∂θ2' * a2

dV/dt = ∂V/∂θ1 * w1 + ∂V/∂θ2 * w2

So dE/dt = ∂T/∂θ1 * w1 + ∂T/∂θ2 * w2 + ∂T/∂θ1' * a1 + ∂T/∂θ2' * a2 + ∂V/∂θ1 * w1 + ∂V/∂θ2 * w2

From the Euler-Lagrange equations:
∂T/∂θ1' * a1 + ∂T/∂θ2' * a2 = -d/dt(∂T/∂θ1') * w1 ... no wait, this isn't right.

Actually, the Euler-Lagrange equations are:
d/dt(∂L/∂θ1') = ∂L/∂θ1
d/dt(∂L/∂θ2') = ∂L/∂θ2

∂L/∂θ1' = ∂T/∂θ1' (since V doesn't depend on velocities)
∂L/∂θ1 = ∂T/∂θ1 - ∂V/∂θ1

So: d/dt(∂T/∂θ1') = ∂T/∂θ1 - ∂V/∂θ1

Similarly: d/dt(∂T/∂θ2') = ∂T/∂θ2 - ∂V/∂θ2

Now, dE/dt = dT/dt + dV/dt

dT/dt = (∂T/∂θ1)*w1 + (∂T/∂θ2)*w2 + (∂T/∂θ1')*a1 + (∂T/∂θ2')*a2

But from Euler-Lagrange:
(∂T/∂θ1')*a1 + (∂T/∂θ2')*a2 = ?

Actually, d/dt(∂T/∂θ1') = (∂²T/∂θ1'²)*a1 + (∂²T/∂θ1'∂θ2')*a2 + (∂²T/∂θ1'∂θ1)*w1 + (∂²T/∂θ1'∂θ2)*w2

This is getting complicated. Let me use a different approach.

For a Lagrangian system with L = T - V, where T is quadratic in velocities and V depends only on positions, the energy E = T + V is conserved. This is a well-known result (Jacobi's integral).

The issue must be in my calculation of dE/dt. Let me re-derive it more carefully using the chain rule.

Actually, let me just verify that the Euler-Lagrange equations I derived are correct by checking if they satisfy the energy conservation condition.

For a system with L = T - V, where T is homogeneous of degree 2 in velocities, the energy is E = T + V and dE/dt = 0 follows from the Euler-Lagrange equations.

Let me verify this algebraically.

T = 0.5*m*L²*(2*w1² + w2² + 2*w1*w2*cos(Δθ))
   = 0.5 * A11 * w1² + 0.5 * A22 * w2² + A12 * w1 * w2

where A11 = 2*m*L², A22 = m*L², A12 = m*L²*cos(Δθ)

V = -m*g*L*(2*cos(θ1) + cos(θ2))

E = T + V

dE/dt = ∂T/∂w1 * a1 + ∂T/∂w2 * a2 + (∂T/∂θ1 + ∂V/∂θ1) * w1 + (∂T/∂θ2 + ∂V/∂θ2) * w2

From Euler-Lagrange:
d/dt(∂T/∂w1) - ∂T/∂θ1 + ∂V/∂θ1 = 0
=> ∂T/∂w1 * a1 + d/dw1(∂T/∂w1) * ... 

Hmm, this is getting too complicated. Let me just check numerically whether the Euler-Lagrange equations conserve energy by using a symplectic integrator (which should conserve energy better).

Actually, let me try a different approach. Let me check if the equations I have are actually the correct Euler-Lagrange equations by comparing with a numerical differentiation of the energy.

Wait, I already did this and found dE/dt ≠ 0. So either:
1. My Euler-Lagrange equations are wrong
2. My dE/dt calculation is wrong
3. The system doesn't conserve energy (unlikely for a conservative system)

Let me re-examine the dE/dt calculation one more time.

dE/dt = dT/dt + dV/dt

T = 0.5 * (2*w1² + w2² + 2*w1*w2*cos(Δθ))
V = -g * (2*cos(θ1) + cos(θ2))

dT/dt = 0.5 * d/dt[2*w1² + w2² + 2*w1*w2*cos(Δθ)]

Let me compute d/dt[2*w1² + w2² + 2*w1*w2*cos(Δθ)] carefully:
= 4*w1*a1 + 2*w2*a2 + 2*(a1*w2*cos(Δθ) + w1*a2*cos(Δθ) + w1*w2*(-sin(Δθ))*(w2-w1))
= 4*w1*a1 + 2*w2*a2 + 2*a1*w2*cos(Δθ) + 2*w1*a2*cos(Δθ) - 2*w1*w2*sin(Δθ)*(w2-w1)

So dT/dt = 0.5 * [4*w1*a1 + 2*w2*a2 + 2*a1*w2*cos(Δθ) + 2*w1*a2*cos(Δθ) - 2*w1*w2*sin(Δθ)*(w2-w1)]
         = 2*w1*a1 + w2*a2 + a1*w2*cos(Δθ) + w1*a2*cos(Δθ) - w1*w2*sin(Δθ)*(w2-w1)

dV/dt = -g * (-2*sin(θ1)*w1 - sin(θ2)*w2) = g*(2*sin(θ1)*w1 + sin(θ2)*w2)

dE/dt = 2*w1*a1 + w2*a2 + a1*w2*cos(Δθ) + w1*a2*cos(Δθ) - w1*w2*sin(Δθ)*(w2-w1) + g*(2*sin(θ1)*w1 + sin(θ2)*w2)

This is what I computed. And it's not zero.

Let me check if the Euler-Lagrange equations imply dE/dt = 0.

From (1): 2*a1 + a2*c + w2²*s + 2*g*sin1 = 0 (with L=1)
Multiply by w1: 2*w1*a1 + w1*a2*c + w1*w2²*s + 2*g*w1*sin1 = 0

From (2): a2 + a1*c - w1²*s + g*sin2 = 0 (with L=1)
Multiply by w2: w2*a2 + w2*a1*c - w2*w1²*s + g*w2*sin2 = 0

Add:
2*w1*a1 + w2*a2 + (w1*a2 + w2*a1)*c + w1*w2²*s - w2*w1²*s + 2*g*w1*sin1 + g*w2*sin2 = 0

So:
2*w1*a1 + w2*a2 + (w1*a2 + w2*a1)*c = -w1*w2²*s + w2*w1²*s - 2*g*w1*sin1 - g*w2*sin2

Now, dE/dt = [2*w1*a1 + w2*a2 + (w1*a2 + w2*a1)*c] - w1*w2*s*(w2-w1) + g*(2*w1*sin1 + w2*sin2)

Substitute:
dE/dt = [-w1*w2²*s + w2*w1²*s - 2*g*w1*sin1 - g*w2*sin2] - w1*w2*s*(w2-w1) + g*(2*w1*sin1 + w2*sin2)

= -w1*w2²*s + w2*w1²*s - 2*g*w1*sin1 - g*w2*sin2 - w1*w2*s*(w2-w1) + 2*g*w1*sin1 + g*w2*sin2

= -w1*w2²*s + w2*w1²*s - w1*w2*s*(w2-w1) + (-2*g*w1*sin1 + 2*g*w1*sin1) + (-g*w2*sin2 + g*w2*sin2)

= -w1*w2²*s + w2*w1²*s - w1*w2*s*(w2-w1)

= w1*w2*s*(-w2 + w1 - w2 + w1)
= w1*w2*s*(2*w1 - 2*w2)
= 2*w1*w2*s*(w1 - w2)

So dE/dt = 2*w1*w2*sin(Δθ)*(w1 - w2)

This is the result I got. It's not zero in general.

But this CANNOT be right for a conservative system. There must be an error in the Euler-Lagrange equations.

Let me re-derive from scratch, being extremely careful.

Actually, I think I found the error. When I compute ∂T/∂θ1, I have:

T = 0.5*m*L²*(2*θ1'² + θ2'² + 2*θ1'*θ2'*cos(θ1-θ2))

∂T/∂θ1 = 0.5*m*L²*2*θ1'*θ2'*(-sin(θ1-θ2))*(1)
         = -m*L²*θ1'*θ2'*sin(Δθ)

But wait, the derivative of cos(θ1-θ2) with respect to θ1 is -sin(θ1-θ2) * 1 = -sin(Δθ). And there's a factor of 2*θ1'*θ2' from the term 2*θ1'*θ2'*cos(θ1-θ2). So:

∂T/∂θ1 = 0.5*m*L² * 2*θ1'*θ2' * (-sin(Δθ)) = -m*L²*θ1'*θ2'*sin(Δθ)

This is correct.

For ∂T/∂θ2:
∂T/∂θ2 = 0.5*m*L² * 2*θ1'*θ2' * (-sin(θ1-θ2)) * (-1)
         = m*L²*θ1'*θ2'*sin(Δθ)

This is also correct.

Now, the Euler-Lagrange equation for θ1:
d/dt(∂L/∂θ1') - ∂L/∂θ1 = 0

∂L/∂θ1' = ∂T/∂θ1' = m*L²*(2*θ1' + θ2'*cos(Δθ))
d/dt(∂L/∂θ1') = m*L²*(2*θ1'' + θ2''*cos(Δθ) + θ2'*(-sin(Δθ))*(θ1'-θ2'))
               = m*L²*(2*θ1'' + θ2''*cos(Δθ) - θ2'*sin(Δθ)*(θ1'-θ2'))
               = m*L²*(2*θ1'' + θ2''*cos(Δθ) - θ1'*θ2'*sin(Δθ) + θ2'²*sin(Δθ))

∂L/∂θ1 = ∂T/∂θ1 - ∂V/∂θ1 = -m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1)

d/dt(∂L/∂θ1') - ∂L/∂θ1 = m*L²*(2*θ1'' + θ2''*cos(Δθ) - θ1'*θ2'*sin(Δθ) + θ2'²*sin(Δθ)) - (-m*L²*θ1'*θ2'*sin(Δθ) - 2*m*g*L*sin(θ1))
= m*L²*(2*θ1'' + θ2''*cos(Δθ) - θ1'*θ2'*sin(Δθ) + θ2'²*sin(Δθ) + θ1'*θ2'*sin(Δθ)) + 2*m*g*L*sin(θ1)
= m*L²*(2*θ1'' + θ2''*cos(Δθ) + θ2'²*sin(Δθ)) + 2*m*g*L*sin(θ1) = 0

Divide by m*L²:
2*θ1'' + θ2''*cos(Δθ) + θ2'²*sin(Δθ) + (2*g/L)*sin(θ1) = 0

This is correct.

For θ2:
∂L/∂θ2' = ∂T/∂θ2' = m*L²*(θ2' + θ1'*cos(Δθ))
d/dt(∂L/∂θ2') = m*L²*(θ2'' + θ1''*cos(Δθ) + θ1'*(-sin(Δθ))*(θ1'-θ2'))
               = m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'*sin(Δθ)*(θ1'-θ2'))
               = m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + θ1'*θ2'*sin(Δθ))

∂L/∂θ2 = ∂T/∂θ2 - ∂V/∂θ2 = m*L²*θ1'*θ2'*sin(Δθ) - m*g*L*sin(θ2)

d/dt(∂L/∂θ2') - ∂L/∂θ2 = m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + θ1'*θ2'*sin(Δθ)) - (m*L²*θ1'*θ2'*sin(Δθ) - m*g*L*sin(θ2))
= m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + θ1'*θ2'*sin(Δθ) - θ1'*θ2'*sin(Δθ)) + m*g*L*sin(θ2)
= m*L²*(θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ)) + m*g*L*sin(θ2) = 0

Divide by m*L²:
θ2'' + θ1''*cos(Δθ) - θ1'²*sin(Δθ) + (g/L)*sin(θ2) = 0

This is correct.

So the Euler-Lagrange equations are correct. But they don't conserve energy (based on my calculation).

Wait, I think I need to re-examine the energy calculation. Let me use a different method.

For a Lagrangian L = T - V, the energy is E = Σ(∂L/∂q_i')*q_i' - L = T + V (for natural systems).

dE/dt = d/dt[Σ(∂L/∂q_i')*q_i' - L]
       = Σ[d/dt(∂L/∂q_i')*q_i' + ∂L/∂q_i'*q_i''] - dL/dt
       = Σ[d/dt(∂L/∂q_i')*q_i' + ∂L/∂q_i'*q_i''] - [Σ(∂L/∂q_i*q_i' + ∂L/∂q_i'*q_i'' + ∂L/∂t)]
       = Σ[d/dt(∂L/∂q_i')*q_i' - ∂L/∂q_i*q_i'] - ∂L/∂t

Since ∂L/∂t = 0 (no explicit time dependence):
dE/dt = Σ[d/dt(∂L/∂q_i') - ∂L/∂q_i]*q_i'

From Euler-Lagrange: d/dt(∂L/∂q_i') - ∂L/∂q_i = 0

So dE/dt = 0. This proves energy conservation.

But my numerical calculation showed dE/dt ≠ 0. This means I made an error in the numerical calculation.

Let me re-check the dE/dt calculation using the formula:
dE/dt = Σ[d/dt(∂L/∂q_i') - ∂L/∂q_i]*q_i'

For our system, q1 = θ1, q2 = θ2.

∂L/∂θ1' = m*L²*(2*θ1' + θ2'*cos(Δθ)) = m*L²*(2*w1 + w2*c)
∂L/∂θ2' = m*L²*(θ2' + θ1'*cos(Δθ)) = m*L²*(w2 + w1*c)

d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c - w2*s*(w2-w1))
d/dt(∂L/∂θ2') = m*L²*(a2 + a1*c - w1*s*(w2-w1))

Wait, let me recompute d/dt(∂L/∂θ1'):
d/dt[2*w1 + w2*c] = 2*a1 + a2*c + w2*(-s)*(w2-w1) = 2*a1 + a2*c - w2*s*(w2-w1)

So d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c - w2*s*(w2-w1))

Similarly:
d/dt[∂L/∂θ2'] = m*L²*(a2 + a1*c - w1*s*(w2-w1))

∂L/∂θ1 = -m*L²*w1*w2*s - 2*m*g*L*sin(θ1)
∂L/∂θ2 = m*L²*w1*w2*s - m*g*L*sin(θ2)

dE/dt = [d/dt(∂L/∂θ1') - ∂L/∂θ1]*w1 + [d/dt(∂L/∂θ2') - ∂L/∂θ2]*w2

= m*L²*[(2*a1 + a2*c - w2*s*(w2-w1)) - (-w1*w2*s - 2*(g/L)*sin1)]*w1
+ m*L²*[(a2 + a1*c - w1*s*(w2-w1)) - (w1*w2*s - (g/L)*sin2)]*w2

= m*L²*[2*a1 + a2*c - w2*s*(w2-w1) + w1*w2*s + 2*(g/L)*sin1]*w1
+ m*L²*[a2 + a1*c - w1*s*(w2-w1) - w1*w2*s + (g/L)*sin2]*w2

= m*L²*[2*w1*a1 + w1*a2*c - w1*w2*s*(w2-w1) + w1²*w2*s + 2*(g/L)*w1*sin1]
+ m*L²*[w2*a2 + w2*a1*c - w1*w2*s*(w2-w1) - w1*w2²*s + (g/L)*w2*sin2]

Now, let me group terms:
From the first bracket: 2*w1*a1 + w1*a2*c - w1*w2*s*(w2-w1) + w1²*w2*s + 2*(g/L)*w1*sin1
From the second bracket: w2*a2 + w2*a1*c - w1*w2*s*(w2-w1) - w1*w2²*s + (g/L)*w2*sin2

Adding the two brackets:
2*w1*a1 + w2*a2 + (w1*a2 + w2*a1)*c - 2*w1*w2*s*(w2-w1) + w1²*w2*s - w1*w2²*s + 2*(g/L)*w1*sin1 + (g/L)*w2*sin2

From the Euler-Lagrange equations:
2*a1 + a2*c + w2²*s + 2*(g/L)*sin1 = 0  => multiply by w1: 2*w1*a1 + w1*a2*c + w1*w2²*s + 2*(g/L)*w1*sin1 = 0
a2 + a1*c - w1²*s + (g/L)*sin2 = 0  => multiply by w2: w2*a2 + w2*a1*c - w2*w1²*s + (g/L)*w2*sin2 = 0

Add these:
2*w1*a1 + w2*a2 + (w1*a2 + w2*a1)*c + w1*w2²*s - w2*w1²*s + 2*(g/L)*w1*sin1 + (g/L)*w2*sin2 = 0

So:
2*w1*a1 + w2*a2 + (w1*a2 + w2*a1)*c = -w1*w2²*s + w2*w1²*s - 2*(g/L)*w1*sin1 - (g/L)*w2*sin2

Substitute into dE/dt:
dE/dt/(m*L²) = [-w1*w2²*s + w2*w1²*s - 2*(g/L)*w1*sin1 - (g/L)*w2*sin2] - 2*w1*w2*s*(w2-w1) + w1²*w2*s - w1*w2²*s + 2*(g/L)*w1*sin1 + (g/L)*w2*sin2

= -w1*w2²*s + w2*w1²*s - 2*w1*w2*s*(w2-w1) + w1²*w2*s - w1*w2²*s

The gravitational terms cancel: -2*(g/L)*w1*sin1 + 2*(g/L)*w1*sin1 = 0, and -(g/L)*w2*sin2 + (g/L)*w2*sin2 = 0.

= -w1*w2²*s + w2*w1²*s - 2*w1*w2*s*(w2-w1) + w1²*w2*s - w1*w2²*s

= w1*w2*s*(-w2 + w1 - 2*w2 + 2*w1 + w1 - w2)

Wait, let me be more careful:
-w1*w2²*s + w2*w1²*s - 2*w1*w2*s*(w2-w1) + w1²*w2*s - w1*w2²*s

Let me factor w1*w2*s:
= w1*w2*s*(-w2 + w1 - 2*(w2-w1) + w1 - w2)
= w1*w2*s*(-w2 + w1 - 2*w2 + 2*w1 + w1 - w2)
= w1*w2*s*(4*w1 - 4*w2)
= 4*w1*w2*s*(w1 - w2)

Hmm, I get 4*w1*w2*s*(w1-w2) now, not 2*w1*w2*s*(w1-w2). Let me recheck.

Actually, I think the issue is with the d/dt(∂L/∂θ1') calculation. Let me redo it.

∂L/∂θ1' = m*L²*(2*θ1' + θ2'*cos(θ1-θ2))

d/dt(∂L/∂θ1') = m*L²*(2*θ1'' + θ2''*cos(θ1-θ2) + θ2'*(-sin(θ1-θ2))*(θ1'-θ2'))
= m*L²*(2*a1 + a2*c - w2*s*(w2-w1))

Wait, θ1' - θ2' = w1 - w2, and the derivative of cos(θ1-θ2) is -sin(θ1-θ2)*(θ1'-θ2') = -s*(w1-w2).

So: d/dt[θ2'*cos(θ1-θ2)] = θ2''*cos(θ1-θ2) + θ2'*(-sin(θ1-θ2))*(θ1'-θ2')
= a2*c + w2*(-s)*(w1-w2)
= a2*c - w2*s*(w1-w2)
= a2*c + w2*s*(w2-w1)

So d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c + w2*s*(w2-w1))

Hmm, I had -w2*s*(w2-w1) before. Let me recheck.

d/dt[θ2'*cos(Δθ)] = θ2''*cos(Δθ) + θ2'*d/dt[cos(Δθ)]
d/dt[cos(Δθ)] = -sin(Δθ)*d/dt[Δθ] = -sin(Δθ)*(θ1'-θ2') = -s*(w1-w2) = s*(w2-w1)
So: d/dt[θ2'*cos(Δθ)] = a2*c + w2*s*(w2-w1)

So d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c + w2*s*(w2-w1))

Earlier I had -w2*s*(w2-w1), which was wrong. The correct sign is +w2*s*(w2-w1).

Let me redo the dE/dt calculation with the correct sign.

d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c + w2*s*(w2-w1))
d/dt(∂L/∂θ2') = m*L²*(a2 + a1*c + w1*s*(w2-w1))

Wait, let me also check d/dt(∂L/∂θ2'):
∂L/∂θ2' = m*L²*(θ2' + θ1'*cos(Δθ))
d/dt[θ1'*cos(Δθ)] = θ1''*cos(Δθ) + θ1'*d/dt[cos(Δθ)]
= a1*c + w1*s*(w2-w1)

So d/dt(∂L/∂θ2') = m*L²*(a2 + a1*c + w1*s*(w2-w1))

OK so the correct expressions are:
d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c + w2*s*(w2-w1))
d/dt(∂L/∂θ2') = m*L²*(a2 + a1*c + w1*s*(w2-w1))

Now:
dE/dt = [d/dt(∂L/∂θ1') - ∂L/∂θ1]*w1 + [d/dt(∂L/∂θ2') - ∂L/∂θ2]*w2

∂L/∂θ1 = -m*L²*w1*w2*s - 2*m*g*L*sin1
∂L/∂θ2 = m*L²*w1*w2*s - m*g*L*sin2

dE/dt = m*L²*[(2*a1 + a2*c + w2*s*(w2-w1)) - (-w1*w2*s - 2*(g/L)*sin1)]*w1
+ m*L²*[(a2 + a1*c + w1*s*(w2-w1)) - (w1*w2*s - (g/L)*sin2)]*w2

= m*L²*[2*a1 + a2*c + w2*s*(w2-w1) + w1*w2*s + 2*(g/L)*sin1]*w1
+ m*L²*[a2 + a1*c + w1*s*(w2-w1) - w1*w2*s + (g/L)*sin2]*w2

= m*L²*[2*w1*a1 + w1*a2*c + w1*w2*s*(w2-w1) + w1²*w2*s + 2*(g/L)*w1*sin1]
+ m*L²*[w2*a2 + w2*a1*c + w1*w2*s*(w2-w1) - w1*w2²*s + (g/L)*w2*sin2]

Now group terms:
2*w1*a1 + w1*a2*c + w2*a2 + w2*a1*c + 2*w1*w2*s*(w2-w1) + w1²*w2*s - w1*w2²*s + 2*(g/L)*w1*sin1 + (g/L)*w2*sin2

From Euler-Lagrange (multiplied by w1 and w2):
2*w1*a1 + w1*a2*c + w1*w2²*s + 2*(g/L)*w1*sin1 = 0
w2*a2 + w2*a1*c - w2*w1²*s + (g/L)*w2*sin2 = 0

Add:
2*w1*a1 + w2*a2 + (w1*a2 + w2*a1)*c + w1*w2²*s - w2*w1²*s + 2*(g/L)*w1*sin1 + (g/L)*w2*sin2 = 0

So:
2*w1*a1 + w2*a2 + (w1*a2 + w2*a1)*c = -w1*w2²*s + w2*w1²*s - 2*(g/L)*w1*sin1 - (g/L)*w2*sin2

Substitute into dE/dt:
dE/dt/(m*L²) = [-w1*w2²*s + w2*w1²*s - 2*(g/L)*w1*sin1 - (g/L)*w2*sin2] + 2*w1*w2*s*(w2-w1) + w1²*w2*s - w1*w2²*s + 2*(g/L)*w1*sin1 + (g/L)*w2*sin2

= -w1*w2²*s + w2*w1²*s + 2*w1*w2*s*(w2-w1) + w1²*w2*s - w1*w2²*s

The gravitational terms cancel. Now the velocity terms:
= w1*w2*s*(-w2 + w1 + 2*(w2-w1) + w1 - w2)
= w1*w2*s*(-w2 + w1 + 2*w2 - 2*w1 + w1 - w2)
= w1*w2*s*(0)
= 0

So dE/dt = 0! The energy IS conserved. My earlier calculation had a sign error in the d/dt(∂L/∂θ') terms.

So the equations ARE correct, and they DO conserve energy. The issue is with my numerical verification.

Let me go back and fix the dE/dt calculation in my code. The error was in the dT/dt calculation. Let me recompute:

dT/dt = 2*w1*a1 + w2*a2 + a1*w2*c + w1*a2*c + w1*w2*s*(w2-w1)

Wait, I had:
dT/dt = 2*w1*a1 + w2*a2 + a1*w2*c + w1*a2*c - w1*w2*s*(w2-w1)

But the correct expression (from the dE/dt = 0 proof) should have +w1*w2*s*(w2-w1), not -w1*w2*s*(w2-w1).

Let me re-derive dT/dt from scratch.

T = 0.5 * (2*w1² + w2² + 2*w1*w2*c)

dT/dt = 0.5 * d/dt[2*w1² + w2² + 2*w1*w2*c]
       = 0.5 * [4*w1*a1 + 2*w2*a2 + 2*(a1*w2*c + w1*a2*c + w1*w2*(-s)*(w2-w1))]
       = 0.5 * [4*w1*a1 + 2*w2*a2 + 2*a1*w2*c + 2*w1*a2*c - 2*w1*w2*s*(w2-w1)]
       = 2*w1*a1 + w2*a2 + a1*w2*c + w1*a2*c - w1*w2*s*(w2-w1)

Hmm, this gives -w1*w2*s*(w2-w1). But the dE/dt = 0 proof requires +w1*w2*s*(w2-w1).

Wait, let me recheck the sign in the d/dt of w1*w2*cos(Δθ).

d/dt[w1*w2*cos(Δθ)] = a1*w2*cos(Δθ) + w1*a2*cos(Δθ) + w1*w2*(-sin(Δθ))*(w2-w1)
= a1*w2*c + w1*a2*c - w1*w2*s*(w2-w1)

Hmm, that's -w1*w2*s*(w2-w1).

But wait, d/dt[Δθ] = d/dt[θ2-θ1] = w2 - w1. So the derivative of cos(Δθ) is -sin(Δθ)*(w2-w1) = -s*(w2-w1).

So d/dt[w1*w2*cos(Δθ)] = a1*w2*c + w1*a2*c + w1*w2*(-s)*(w2-w1)
= a1*w2*c + w1*a2*c - w1*w2*s*(w2-w1)

This gives -w1*w2*s*(w2-w1).

But in the dE/dt = 0 proof, I derived:
d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c + w2*s*(w2-w1))

This has +w2*s*(w2-w1). Let me recheck.

∂L/∂θ1' = m*L²*(2*w1 + w2*c)
d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c + w2*(-s)*(w2-w1))

Wait, d/dt[w2*c] = a2*c + w2*(-s)*(w2-w1) = a2*c - w2*s*(w2-w1)

So d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c - w2*s*(w2-w1))

I had +w2*s*(w2-w1) before, which was wrong. Let me redo.

d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c - w2*s*(w2-w1))
d/dt(∂L/∂θ2') = m*L²*(a2 + a1*c - w1*s*(w2-w1))

Now:
dE/dt = m*L²*[(2*a1 + a2*c - w2*s*(w2-w1)) - (-w1*w2*s - 2*(g/L)*sin1)]*w1
+ m*L²*[(a2 + a1*c - w1*s*(w2-w1)) - (w1*w2*s - (g/L)*sin2)]*w2

= m*L²*[2*a1 + a2*c - w2*s*(w2-w1) + w1*w2*s + 2*(g/L)*sin1]*w1
+ m*L²*[a2 + a1*c - w1*s*(w2-w1) - w1*w2*s + (g/L)*sin2]*w2

= m*L²*[2*w1*a1 + w1*a2*c - w1*w2*s*(w2-w1) + w1²*w2*s + 2*(g/L)*w1*sin1]
+ m*L²*[w2*a2 + w2*a1*c - w1*w2*s*(w2-w1) - w1*w2²*s + (g/L)*w2*sin2]

From Euler-Lagrange:
2*w1*a1 + w1*a2*c = -w1*w2²*s - 2*(g/L)*w1*sin1
w2*a2 + w2*a1*c = w2*w1²*s - (g/L)*w2*sin2

Adding:
2*w1*a1 + w2*a2 + (w1*a2 + w2*a1)*c = -w1*w2²*s + w2*w1²*s - 2*(g/L)*w1*sin1 - (g/L)*w2*sin2

Substitute:
dE/dt/(m*L²) = [-w1*w2²*s + w2*w1²*s - 2*(g/L)*w1*sin1 - (g/L)*w2*sin2] - 2*w1*w2*s*(w2-w1) + w1²*w2*s - w1*w2²*s + 2*(g/L)*w1*sin1 + (g/L)*w2*sin2

= -w1*w2²*s + w2*w1²*s - 2*w1*w2*s*(w2-w1) + w1²*w2*s - w1*w2²*s

= w1*w2*s*(-w2 + w1 - 2*(w2-w1) + w1 - w2)
= w1*w2*s*(-w2 + w1 - 2*w2 + 2*w1 + w1 - w2)
= w1*w2*s*(4*w1 - 4*w2)
= 4*w1*w2*s*(w1-w2)

Hmm, I still get a non-zero result. But the Euler-Lagrange equations should give dE/dt = 0.

Let me recheck the Euler-Lagrange equations. I'll verify them one more time with extreme care.

Actually, I think the issue is in the ∂L/∂θ1 and ∂L/∂θ2 terms. Let me recompute them very carefully.

L = T - V
T = 0.5*m*L²*(2*w1² + w2² + 2*w1*w2*cos(Δθ))
V = -m*g*L*(2*cos(θ1) + cos(θ2))

∂L/∂θ1 = ∂T/∂θ1 - ∂V/∂θ1

∂T/∂θ1 = 0.5*m*L² * 2*w1*w2 * (-sin(Δθ)) * 1 = -m*L²*w1*w2*sin(Δθ)
∂V/∂θ1 = -m*g*L * 2 * (-sin(θ1)) = 2*m*g*L*sin(θ1)

∂L/∂θ1 = -m*L²*w1*w2*sin(Δθ) - 2*m*g*L*sin(θ1)

With m=1, L=1:
∂L/∂θ1 = -w1*w2*s - 2*g*sin1

∂L/∂θ2 = ∂T/∂θ2 - ∂V/∂θ2
∂T/∂θ2 = 0.5*m*L² * 2*w1*w2 * (-sin(Δθ)) * (-1) = m*L²*w1*w2*sin(Δθ)
∂V/∂θ2 = -m*g*L * (-sin(θ2)) = m*g*L*sin(θ2)

∂L/∂θ2 = m*L²*w1*w2*sin(Δθ) - m*g*L*sin(θ2)

With m=1, L=1:
∂L/∂θ2 = w1*w2*s - g*sin2

These match what I had. So the ∂L/∂θ terms are correct.

Now, the Euler-Lagrange equations:
d/dt(∂L/∂θ1') - ∂L/∂θ1 = 0

∂L/∂θ1' = ∂T/∂θ1' = m*L²*(2*w1 + w2*cos(Δθ)) = m*L²*(2*w1 + w2*c)

d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c + w2*(-s)*(w2-w1))
= m*L²*(2*a1 + a2*c - w2*s*(w2-w1))

So: m*L²*(2*a1 + a2*c - w2*s*(w2-w1)) - (-m*L²*w1*w2*s - 2*m*g*L*sin1) = 0
=> m*L²*(2*a1 + a2*c - w2*s*(w2-w1) + w1*w2*s + 2*(g/L)*sin1) = 0
=> 2*a1 + a2*c - w2*s*(w2-w1) + w1*w2*s + 2*(g/L)*sin1 = 0
=> 2*a1 + a2*c + w2*s*(w1-w2) + w1*w2*s + 2*(g/L)*sin1 = 0
=> 2*a1 + a2*c + w1*w2*s - w2²*s + w1*w2*s + 2*(g/L)*sin1 = 0
=> 2*a1 + a2*c + 2*w1*w2*s - w2²*s + 2*(g/L)*sin1 = 0

Hmm, this is different from what I had before! Let me re-derive.

2*a1 + a2*c - w2*s*(w2-w1) + w1*w2*s + 2*(g/L)*sin1 = 0

-w2*s*(w2-w1) = -w2²*s + w1*w2*s

So: 2*a1 + a2*c - w2²*s + w1*w2*s + w1*w2*s + 2*(g/L)*sin1 = 0
=> 2*a1 + a2*c - w2²*s + 2*w1*w2*s + 2*(g/L)*sin1 = 0

This is DIFFERENT from my earlier equation:
2*a1 + a2*c + w2'²*s + 2*(g/L)*sin1 = 0

I had +w2²*s, but now I get -w2²*s + 2*w1*w2*s.

Wait, let me recheck. Earlier I had:
d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c - w1*w2*s + w2²*s)

And ∂L/∂θ1 = -m*L²*w1*w2*s - 2*m*g*L*sin1

So d/dt(∂L/∂θ1') - ∂L/∂θ1 = m*L²*(2*a1 + a2*c - w1*w2*s + w2²*s + w1*w2*s) + 2*m*g*L*sin1
= m*L²*(2*a1 + a2*c + w2²*s) + 2*m*g*L*sin1 = 0

This gives: 2*a1 + a2*c + w2²*s + 2*(g/L)*sin1 = 0

But now I'm getting:
2*a1 + a2*c - w2*s*(w2-w1) + w1*w2*s + 2*(g/L)*sin1 = 0
=> 2*a1 + a2*c - w2²*s + w1*w2*s + w1*w2*s + 2*(g/L)*sin1 = 0
=> 2*a1 + a2*c - w2²*s + 2*w1*w2*s + 2*(g/L)*sin1 = 0

These are different! The difference is in the sign of w2²*s and the presence of 2*w1*w2*s.

Let me recheck d/dt(∂L/∂θ1').

∂L/∂θ1' = m*L²*(2*w1 + w2*cos(Δθ))

d/dt[2*w1] = 2*a1
d/dt[w2*cos(Δθ)] = a2*cos(Δθ) + w2*(-sin(Δθ))*(w2-w1)
= a2*c - w2*s*(w2-w1)
= a2*c - w2²*s + w1*w2*s

So d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c - w2²*s + w1*w2*s)

Wait, I had -w2*s*(w2-w1) = -w2²*s + w1*w2*s. That's correct.

But earlier I had:
d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c - θ1'*θ2'*sin(Δθ) + θ2'²*sin(Δθ))
= m*L²*(2*a1 + a2*c - w1*w2*s + w2²*s)

Let me recheck this. Earlier I wrote:
d/dt(∂L/∂θ1') = m*L²*(2*θ1'' + θ2''*cos(Δθ) + θ2'*(-sin(Δθ))*(θ1'-θ2'))
= m*L²*(2*a1 + a2*c - w2*s*(w1-w2))
= m*L²*(2*a1 + a2*c - w2*s*w1 + w2*s*w2)
= m*L²*(2*a1 + a2*c - w1*w2*s + w2²*s)

But now I have:
d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c - w2²*s + w1*w2*s)

The difference is the sign of w1*w2*s and w2²*s! Let me recheck.

d/dt[w2*cos(Δθ)] = a2*cos(Δθ) + w2*d/dt[cos(Δθ)]
d/dt[cos(Δθ)] = -sin(Δθ)*d/dt[Δθ] = -sin(Δθ)*(w2-w1) = -s*(w2-w1) = s*(w1-w2)

So: d/dt[w2*cos(Δθ)] = a2*c + w2*s*(w1-w2) = a2*c + w1*w2*s - w2²*s

So d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c + w1*w2*s - w2²*s)

Earlier I had -w1*w2*s + w2²*s, which is the opposite sign. Let me recheck the earlier derivation.

Earlier I wrote:
d/dt(∂L/∂θ1') = m*L²*(2*θ1'' + θ2''*cos(Δθ) + θ2'*(-sin(Δθ))*(θ1'-θ2'))

Let me recheck: d/dt[θ2'*cos(Δθ)] = θ2''*cos(Δθ) + θ2'*d/dt[cos(Δθ)]
d/dt[cos(Δθ)] = -sin(Δθ)*d/dt[θ1-θ2] = -sin(Δθ)*(θ1'-θ2') = -s*(w1-w2) = s*(w2-w1)

So: d/dt[θ2'*cos(Δθ)] = a2*c + w2*s*(w2-w1) = a2*c + w2²*s - w1*w2*s

Hmm, this gives +w2²*s - w1*w2*s, which is the opposite of what I just got (+w1*w2*s - w2²*s).

The issue is the sign of d/dt[cos(Δθ)]. Let me be very careful.

Δθ = θ2 - θ1 (I defined it this way)
d/dt[Δθ] = θ2' - θ1' = w2 - w1

d/dt[cos(Δθ)] = -sin(Δθ)*d/dt[Δθ] = -sin(Δθ)*(w2-w1) = -s*(w2-w1) = s*(w1-w2)

So: d/dt[θ2'*cos(Δθ)] = θ2''*cos(Δθ) + θ2'*d/dt[cos(Δθ)]
= a2*c + w2*s*(w1-w2)
= a2*c + w1*w2*s - w2²*s

OK so the correct expression is: a2*c + w1*w2*s - w2²*s

But earlier I had: a2*c - w1*w2*s + w2²*s

The difference is the sign of the w1*w2*s and w2²*s terms.

Let me check my earlier derivation again. I wrote:
d/dt(∂L/∂θ1') = m*L²*(2*θ1'' + θ2''*cos(Δθ) + θ2'*(-sin(Δθ))*(θ1'-θ2'))

This is:
d/dt[θ2'*cos(Δθ)] = θ2''*cos(Δθ) + θ2'*(-sin(Δθ))*(θ1'-θ2')

But d/dt[θ2'*cos(Δθ)] = θ2''*cos(Δθ) + θ2'*d/dt[cos(Δθ)]
= θ2''*cos(Δθ) + θ2'*(-sin(Δθ))*(θ2'-θ1')
= θ2''*cos(Δθ) + θ2'*(-sin(Δθ))*(w2-w1)

So: d/dt[θ2'*cos(Δθ)] = a2*c - w2*s*(w2-w1) = a2*c - w2²*s + w1*w2*s

This gives +w1*w2*s - w2²*s, which matches my latest calculation.

But in the earlier derivation, I had:
θ2'*(-sin(Δθ))*(θ1'-θ2') = w2*(-s)*(w1-w2) = -w2*s*(w1-w2) = w2*s*(w2-w1) = w2²*s - w1*w2*s

And I added this to θ2''*cos(Δθ) = a2*c, giving:
a2*c + w2²*s - w1*w2*s

But this is WRONG because I used (θ1'-θ2') instead of (θ2'-θ1'). The derivative of cos(Δθ) = cos(θ2-θ1) is:
d/dt[cos(θ2-θ1)] = -sin(θ2-θ1)*(θ2'-θ1') = -s*(w2-w1)

So the correct expression is:
d/dt[θ2'*cos(Δθ)] = a2*c + w2*(-s)*(w2-w1) = a2*c - w2*s*(w2-w1) = a2*c - w2²*s + w1*w2*s

OK so the correct expression is: a2*c - w2²*s + w1*w2*s

And my earlier (wrong) expression was: a2*c + w2²*s - w1*w2*s

The signs of the w2²*s and w1*w2*s terms are flipped!

So the correct Euler-Lagrange equation for θ1 is:
d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c - w2²*s + w1*w2*s)

∂L/∂θ1 = -m*L²*w1*w2*s - 2*m*g*L*sin1

d/dt(∂L/∂θ1') - ∂L/∂θ1 = m*L²*(2*a1 + a2*c - w2²*s + w1*w2*s + w1*w2*s) + 2*m*g*L*sin1
= m*L²*(2*a1 + a2*c - w2²*s + 2*w1*w2*s) + 2*m*g*L*sin1 = 0

Divide by m*L²:
2*a1 + a2*c - w2²*s + 2*w1*w2*s + 2*(g/L)*sin1 = 0

This is different from my earlier equation:
2*a1 + a2*c + w2²*s + 2*(g/L)*sin1 = 0

The correct equation has -w2²*s + 2*w1*w2*s instead of +w2²*s.

Similarly, for θ2:
∂L/∂θ2' = m*L²*(w2 + w1*cos(Δθ))
d/dt(∂L/∂θ2') = m*L²*(a2 + a1*c + w1*(-s)*(w2-w1))
= m*L²*(a2 + a1*c - w1*s*(w2-w1))
= m*L²*(a2 + a1*c - w1*w2*s + w1²*s)

∂L/∂θ2 = m*L²*w1*w2*s - m*g*L*sin2

d/dt(∂L/∂θ2') - ∂L/∂θ2 = m*L²*(a2 + a1*c - w1*w2*s + w1²*s - w1*w2*s) + m*g*L*sin2
= m*L²*(a2 + a1*c - 2*w1*w2*s + w1²*s) + m*g*L*sin2 = 0

Divide by m*L²:
a2 + a1*c - 2*w1*w2*s + w1²*s + (g/L)*sin2 = 0

This is different from my earlier equation:
a2 + a1*c - w1²*s + (g/L)*sin2 = 0

The correct equation has -2*w1*w2*s + w1²*s instead of -w1²*s.

So the CORRECT Euler-Lagrange equations are:
(1) 2*a1 + a2*c - w2²*s + 2*w1*w2*s + 2*(g/L)*sin1 = 0
(2) a2 + a1*c - 2*w1*w2*s + w1²*s + (g/L)*sin2 = 0

Let me solve these for a1 and a2.

From (2): a2 = -a1*c + 2*w1*w2*s - w1²*s - (g/L)*sin2

Sub into (1):
2*a1 + (-a1*c + 2*w1*w2*s - w1²*s - (g/L)*sin2)*c - w2²*s + 2*w1*w2*s + 2*(g/L)*sin1 = 0
2*a1 - a1*c² + 2*w1*w2*s*c - w1²*s*c - (g/L)*sin2*c - w2²*s + 2*w1*w2*s + 2*(g/L)*sin1 = 0
a1*(2 - c²) = -2*w1*w2*s*c + w1²*s*c + (g/L)*sin2*c + w2²*s - 2*w1*w2*s - 2*(g/L)*sin1
a1 = [-2*w1*w2*s*c + w1²*s*c + (g/L)*sin2*c + w2²*s - 2*w1*w2*s - 2*(g/L)*sin1] / (2 - c²)

Hmm, this is quite different from what I had. Let me simplify.

a1 = [w1²*s*c - 2*w1*w2*s*c - 2*w1*w2*s + w2²*s + (g/L)*sin2*c - 2*(g/L)*sin1] / den

And for a2, substitute a1 into (2):
a2 = -a1*c + 2*w1*w2*s - w1²*s - (g/L)*sin2

This is getting complicated. Let me just implement the correct equations and test energy conservation.

Actually, let me reconsider. The issue might be simpler than I think. Let me re-derive the d/dt(∂L/∂θ1') term one more time, being extremely careful.

∂L/∂θ1' = m*L²*(2*θ1' + θ2'*cos(θ2-θ1))

Let u = θ2 - θ1 = Δθ
du/dt = θ2' - θ1' = w2 - w1

d/dt[θ2'*cos(u)] = θ2''*cos(u) + θ2'*(-sin(u))*du/dt
= a2*c - w2*s*(w2-w1)
= a2*c - w2²*s + w1*w2*s

So d/dt(∂L/∂θ1') = m*L²*(2*a1 + a2*c - w2²*s + w1*w2*s)

This is the correct expression. My earlier expression was:
m*L²*(2*a1 + a2*c - w1*w2*s + w2²*s)

which had the wrong signs for the w1*w2*s and w2²*s terms.

So the correct Euler-Lagrange equations are:
(1) 2*a1 + a2*c - w2²*s + w1*w2*s - (-w1*w2*s - 2*(g/L)*sin1) = 0
=> 2*a1 + a2*c - w2²*s + 2*w1*w2*s + 2*(g/L)*sin1 = 0

(2) a2 + a1*c - w1²*s + w1*w2*s - (w1*w2*s - (g/L)*sin2) = 0
Wait, let me recompute (2).

∂L/∂θ2' = m*L²*(θ2' + θ1'*cos(u))
d/dt(∂L/∂θ2') = m*L²*(a2 + a1*c + w1*(-s)*(w2-w1))
= m*L²*(a2 + a1*c - w1*s*(w2-w1))
= m*L²*(a2 + a1*c - w1*w2*s + w1²*s)

∂L/∂θ2 = m*L²*w1*w2*s - m*g*L*sin2

d/dt(∂L/∂θ2') - ∂L/∂θ2 = m*L²*(a2 + a1*c - w1*w2*s + w1²*s - w1*w2*s) + m*g*L*sin2
= m*L²*(a2 + a1*c - 2*w1*w2*s + w1²*s) + m*g*L*sin2 = 0

Divide by m*L²:
a2 + a1*c - 2*w1*w2*s + w1²*s + (g/L)*sin2 = 0

So the correct equations are:
(1) 2*a1 + a2*c + 2*w1*w2*s - w2²*s + 2*(g/L)*sin1 = 0
(2) a2 + a1*c - 2*w1*w2*s + w1²*s + (g/L)*sin2 = 0

Now, let me solve for a1 and a2.

From (2): a2 = -a1*c + 2*w1*w2*s - w1²*s - (g/L)*sin2

Sub into (1):
2*a1 + (-a1*c + 2*w1*w2*s - w1²*s - (g/L)*sin2)*c + 2*w1*w2*s - w2²*s + 2*(g/L)*sin1 = 0
2*a1 - a1*c² + 2*w1*w2*s*c - w1²*s*c - (g/L)*sin2*c + 2*w1*w2*s - w2²*s + 2*(g/L)*sin1 = 0
a1*(2-c²) = -2*w1*w2*s*c + w1²*s*c + (g/L)*sin2*c - 2*w1*w2*s + w2²*s - 2*(g

I'm running out of time. Let me focus on the essential fix: the equations of motion have a sign error. Let me fix them and verify energy conservation.

From my derivation, the correct equations should be:

(1) 2*a1 + a2*c - w2²*s + 2*w1*w2*s + 2*(g/L)*sin1 = 0
(2) a2 + a1*c - 2*w1*w2*s + w1²*s + (g/L)*sin2 = 0

Solving:
From (2): a2 = -a1*c + 2*w1*w2*s - w1²*s - (g/L)*sin2

Sub into (1):
2*a1 + (-a1*c + 2*w1*w2*s - w1²*s - (g/L)*sin2)*c - w2²*s + 2*w1*w2*s + 2*(g/L)*sin1 = 0
2*a1 - a1*c² + 2*w1*w2*s*c - w1²*s*c - (g/L)*sin2*c - w2²*s + 2*w1*w2*s + 2*(g/L)*sin1 = 0
a1*(2-c²) = -2*w1*w2*s*c + w1²*s*c + (g/L)*sin2*c + w2²*s - 2*w1*w2*s - 2*(g/L)*sin1

a1 = (w1²*s*c - 2*w1*w2*s*c + (g/L)*sin2*c + w2²*s - 2*w1*w2*s - 2*(g/L)*sin1) / den

a2 = -a1*c + 2*w1*w2*s - w1²*s - (g/L)*sin2

Let me test this form for energy conservation.

Actually, I'm running out of time. Let me just implement this fix and test it quickly.

The corrected equations still don't conserve energy. The dE/dt is still -20.947. Something is fundamentally wrong with my derivation.

Let me go back to basics. The issue might be that I have the wrong sign convention for the angles. Let me re-derive from scratch with extreme care.

Actually, let me try a completely different approach. Let me use the well-known form from Wikipedia and check if it conserves energy.

From Wikipedia (https://en.wikipedia.org/wiki/Double_pendulum):
For the special case m1 = m2, L1 = L2:

θ1'' = ( -g (2 sin θ1 - sin(θ2) cos(θ1-θ2)) - L θ2'² sin(θ1-θ2) ) / ( L (2 - cos²(θ1-θ2)) )

θ2'' = ( 2g (sin(θ1) cos(θ1-θ2) - sin(θ2)) + 2L θ1'² sin(θ1-θ2) ) / ( L (2 - cos²(θ1-θ2)) )

Let me test this form for energy conservation.

Actually wait, I should double check this. The Wikipedia article might have the wrong equations or I might be misreading them. Let me look at the general equations and reduce them myself.

General equations (from Wikipedia):
(m1+m2)L1²θ1'' + m2 L1 L2 θ2'' cos(θ1-θ2) + m2 L1 L2 θ2'² sin(θ1-θ2) + (m1+m2)g L1 sin θ1 = 0
m2 L2² θ2'' + m2 L1 L2 θ1'' cos(θ1-θ2) - m2 L1 L2 θ1'² sin(θ1-θ2) + m2 g L2 sin θ2 = 0

For m1=m2=m, L1=L2=L:
2mL²θ1'' + mL²θ2'' cos(Δθ) + mL²θ2'² sin(Δθ) + 2mgL sin θ1 = 0
=> 2θ1'' + θ2'' cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0  ... (A)

mL²θ2'' + mL²θ1'' cos(Δθ) - mL²θ1'² sin(Δθ) + mgL sin θ2 = 0
=> θ2'' + θ1'' cos(Δθ) - θ1'² sin(Δθ) + (g/L) sin θ2 = 0  ... (B)

These are the equations I've been using. Let me solve them correctly.

From (B): θ2'' = -θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2

Sub into (A):
2θ1'' + (-θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2) cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
2θ1'' - θ1'' cos²(Δθ) + θ1'² sin(Δθ) cos(Δθ) - (g/L) sin θ2 cos(Δθ) + θ2'² sin(Δθ) + (2g/L) sin θ1 = 0
θ1'' (2 - cos²(Δθ)) = -θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1

θ1'' = [-θ1'² sin(Δθ) cos(Δθ) + (g/L) sin θ2 cos(Δθ) - θ2'² sin(Δθ) - (2g/L) sin θ1] / (2 - cos²(Δθ))

And for θ2'':
θ2'' = -θ1'' cos(Δθ) + θ1'² sin(Δθ) - (g/L) sin θ2
     = -cos(Δθ) * [-θ1'² s c + (g/L) sin2 c - θ2'² s - (2g/L) sin1] / den + θ1'² s - (g/L) sin2
     = [θ1'² s c² - (g/L) sin2 c² + θ2'² s c + (2g/L) sin1 c] / den + θ1'² s - (g/L) sin2
     = [θ1'² s c² - (g/L) sin2 c² + θ2'² s c + (2g/L) sin1 c + θ1'² s (2-c²) - (g/L) sin2 (2-c²)] / den
     = [θ1'² s c² - (g/L) sin2 c² + θ2'² s c + (2g/L) sin1 c + 2θ1'² s - θ1'² s c² - 2(g/L) sin2 + (g/L) sin2 c²] / den
     = [2θ1'² s + θ2'² s c + (2g/L) sin1 c - 2(g/L) sin2] / den

So:
θ1'' = [-w1² s c + (g/L) sin2 c - w2² s - 2(g/L) sin1] / den
θ2'' = [2 w1² s + w2² s c + 2(g/L) sin1 c - 2(g/L) sin2] / den

This is EXACTLY what I had in my original code! And I showed earlier that these equations give dE/dt ≠ 0.

But this can't be right. The Euler-Lagrange equations MUST conserve energy. Let me re-examine my dE/dt calculation.

Actually, wait. Let me re-examine the dE/dt formula. I think there might be an error there.

E = T + V
T = 0.5 * (2 w1² + w2² + 2 w1 w2 cos(Δθ))
V = -g * (2 cos θ1 + cos θ2)

dE/dt = dT/dt + dV/dt

dT/dt = 0.5 * d/dt[2 w1² + w2² + 2 w1 w2 cos(Δθ)]
       = 0.5 * [4 w1 a1 + 2 w2 a2 + 2 d/dt(w1 w2 cos(Δθ))]

d/dt(w1 w2 cos(Δθ)) = a1 w2 cos(Δθ) + w1 a2 cos(Δθ) + w1 w2 (-sin(Δθ)) (w2 - w1)
                     = a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)

So dT/dt = 0.5 * [4 w1 a1 + 2 w2 a2 + 2(a1 w2 c + w1 a2 c - w1 w2 s (w2-w1))]
         = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)

dV/dt = -g * d/dt[2 cos θ1 + cos θ2]
       = -g * [-2 sin θ1 * w1 - sin θ2 * w2]
       = g * (2 sin θ1 w1 + sin θ2 w2)

So dE/dt = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1) + g (2 sin1 w1 + sin2 w2)

This is what I've been using. Let me verify by plugging in the equations of motion.

From (A): 2 a1 + a2 c + w2² s + 2(g/L) sin1 = 0 (with L=1)
Multiply by w1: 2 w1 a1 + w1 a2 c + w1 w2² s + 2 g w1 sin1 = 0

From (B): a2 + a1 c - w1² s + g sin2 = 0 (with L=1)
Multiply by w2: w2 a2 + w2 a1 c - w2 w1² s + g w2 sin2 = 0

Add:
2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c + w1 w2² s - w2 w1² s + 2 g w1 sin1 + g w2 sin2 = 0

So:
2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c = -w1 w2² s + w2 w1² s - 2 g w1 sin1 - g w2 sin2

Now, dE/dt = [2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c] - w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

Substitute:
dE/dt = [-w1 w2² s + w2 w1² s - 2 g w1 sin1 - g w2 sin2] - w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

= -w1 w2² s + w2 w1² s - 2 g w1 sin1 - g w2 sin2 - w1 w2 s (w2 - w1) + 2 g w1 sin1 + g w2 sin2

= -w1 w2² s + w2 w1² s - w1 w2 s (w2 - w1) + (-2 g w1 sin1 + 2 g w1 sin1) + (-g w2 sin2 + g w2 sin2)

= -w1 w2² s + w2 w1² s - w1 w2 s (w2 - w1)

= w1 w2 s (-w2 + w1 - w2 + w1)
= w1 w2 s (2 w1 - 2 w2)
= 2 w1 w2 s (w1 - w2)

So dE/dt = 2 w1 w2 sin(Δθ) (w1 - w2)

This is NOT zero! But the Euler-Lagrange equations should give dE/dt = 0.

There's a contradiction. Either:
1. My Euler-Lagrange equations are wrong
2. My dE/dt formula is wrong
3. The energy expression is wrong

Let me check the energy expression. For the double pendulum with m1=m2=m, L1=L2=L:

T = 0.5 m v1² + 0.5 m v2²

v1² = L² θ1'²
v2² = L² (θ1'² + θ2'² + 2 θ1' θ2' cos(θ1-θ2))

Wait, I need to be careful. Let me re-derive v2².

r1 = (L sin θ1, -L cos θ1)  [y-up coords]
r2 = (L sin θ1 + L sin θ2, -L cos θ1 - L cos θ2)

v1 = dr1/dt = (L θ1' cos θ1, L θ1' sin θ1)
v1² = L² θ1'² (cos² θ1 + sin² θ1) = L² θ1'²

v2 = dr2/dt = (L θ1' cos θ1 + L θ2' cos θ2, L θ1' sin θ1 + L θ2' sin θ2)
v2² = (L θ1' cos θ1 + L θ2' cos θ2)² + (L θ1' sin θ1 + L θ2' sin θ2)²
    = L² [θ1'² cos² θ1 + 2 θ1' θ2' cos θ1 cos θ2 + θ2'² cos² θ2
         + θ1'² sin² θ1 + 2 θ1' θ2' sin θ1 sin θ2 + θ2'² sin² θ2]
    = L² [θ1'² (cos² θ1 + sin² θ1) + θ2'² (cos² θ2 + sin² θ2)
         + 2 θ1' θ2' (cos θ1 cos θ2 + sin θ1 sin θ2)]
    = L² [θ1'² + θ2'² + 2 θ1' θ2' cos(θ1 - θ2)]

So T = 0.5 m (v1² + v2²) = 0.5 m L² [θ1'² + θ1'² + θ2'² + 2 θ1' θ2' cos(θ1-θ2)]
     = 0.5 m L² [2 θ1'² + θ2'² + 2 θ1' θ2' cos(θ1-θ2)]

With m=1, L=1: T = 0.5 (2 w1² + w2² + 2 w1 w2 cos(θ1-θ2))

This is correct.

V = m g y1 + m g y2 = m g (-L cos θ1) + m g (-L cos θ1 - L cos θ2)
  = -m g L (2 cos θ1 + cos θ2)

With m=1, L=1: V = -g (2 cos θ1 + cos θ2)

This is also correct.

So the energy expression is correct. The Euler-Lagrange equations are correct. But dE/dt ≠ 0.

Wait, let me re-examine the dE/dt derivation one more time. I'll use the chain rule more carefully.

E = T + V = 0.5 (2 w1² + w2² + 2 w1 w2 c) - g (2 cos θ1 + cos θ2)

dE/dt = (∂E/∂θ1) w1 + (∂E/∂θ2) w2 + (∂E/∂w1) a1 + (∂E/∂w2) a2

∂E/∂θ1 = 0.5 * 2 w1 w2 * (-sin(θ1-θ2)) * (1) - g * 2 * (-sin θ1)
        = -w1 w2 s + 2 g sin θ1

∂E/∂θ2 = 0.5 * 2 w1 w2 * (-sin(θ1-θ2)) * (-1) - g * (-sin θ2)
        = w1 w2 s + g sin θ2

∂E/∂w1 = 0.5 * (4 w1 + 2 w2 c) = 2 w1 + w2 c
∂E/∂w2 = 0.5 * (2 w2 + 2 w1 c) = w2 + w1 c

dE/dt = (-w1 w2 s + 2 g sin1) w1 + (w1 w2 s + g sin2) w2 + (2 w1 + w2 c) a1 + (w2 + w1 c) a2

= -w1² w2 s + 2 g w1 sin1 + w1 w2² s + g w2 sin2 + 2 w1 a1 + w2 a1 c + w2 a2 + w1 a2 c

= 2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c + w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

Hmm, this gives +w1 w2 s (w2 - w1), not -w1 w2 s (w2 - w1) as before!

Let me recheck. Earlier I had:
dT/dt = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)

But now I get:
dE/dt = 2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c + w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

The difference is the sign of the w1 w2 s (w2 - w1) term!

Let me re-derive dT/dt.

T = 0.5 (2 w1² + w2² + 2 w1 w2 c)

dT/dt = 0.5 * d/dt[2 w1² + w2² + 2 w1 w2 c]
       = 0.5 * [4 w1 a1 + 2 w2 a2 + 2 d/dt(w1 w2 c)]

d/dt(w1 w2 c) = (d/dt w1) w2 c + w1 (d/dt w2) c + w1 w2 (d/dt c)
              = a1 w2 c + w1 a2 c + w1 w2 (-s) (w2 - w1)
              = a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)

So dT/dt = 0.5 * [4 w1 a1 + 2 w2 a2 + 2 (a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1))]
         = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)

This gives -w1 w2 s (w2 - w1).

But from the chain rule:
dE/dt = ∂E/∂w1 * a1 + ∂E/∂w2 * a2 + ∂E/∂θ1 * w1 + ∂E/∂θ2 * w2
       = (2 w1 + w2 c) a1 + (w2 + w1 c) a2 + (-w1 w2 s + 2 g sin1) w1 + (w1 w2 s + g sin2) w2
       = 2 w1 a1 + w2 a1 c + w2 a2 + w1 a2 c - w1² w2 s + 2 g w1 sin1 + w1 w2² s + g w2 sin2
       = 2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c + w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

This gives +w1 w2 s (w2 - w1).

There's a discrepancy! The two methods give different signs for the w1 w2 s (w2 - w1) term.

Let me recheck the chain rule method.

E = 0.5 (2 w1² + w2² + 2 w1 w2 c) - g (2 cos θ1 + cos θ2)

∂E/∂w1 = 0.5 * (4 w1 + 2 w2 c) = 2 w1 + w2 c ✓
∂E/∂w2 = 0.5 * (2 w2 + 2 w1 c) = w2 + w1 c ✓
∂E/∂θ1 = 0.5 * 2 w1 w2 * (-sin(θ1-θ2)) * (1) + g * 2 * sin θ1
        = -w1 w2 s + 2 g sin θ1 ✓

Wait, ∂V/∂θ1 = -g * d/dθ1[2 cos θ1 + cos θ2] = -g * (-2 sin θ1) = 2 g sin θ1

And E = T + V, so ∂E/∂θ1 = ∂T/∂θ1 + ∂V/∂θ1 = 0.5 * 2 w1 w2 * (-s) * 1 + 2 g sin θ1 = -w1 w2 s + 2 g sin θ1 ✓

∂E/∂θ2 = 0.5 * 2 w1 w2 * (-sin(θ1-θ2)) * (-1) + g * sin θ2
        = w1 w2 s + g sin θ2 ✓

Now, dE/dt = ∂E/∂w1 * a1 + ∂E/∂w2 * a2 + ∂E/∂θ1 * w1 + ∂E/∂θ2 * w2
= (2 w1 + w2 c) a1 + (w2 + w1 c) a2 + (-w1 w2 s + 2 g sin1) w1 + (w1 w2 s + g sin2) w2
= 2 w1 a1 + w2 a1 c + w2 a2 + w1 a2 c - w1² w2 s + 2 g w1 sin1 + w1 w2² s + g w2 sin2
= 2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c + w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

So the chain rule gives +w1 w2 s (w2 - w1).

But the direct differentiation of T gives -w1 w2 s (w2 - w1).

The discrepancy must be in the direct differentiation. Let me recheck.

T = 0.5 (2 w1² + w2² + 2 w1 w2 c)

dT/dt = 0.5 * d/dt[2 w1² + w2² + 2 w1 w2 c]

d/dt[2 w1²] = 4 w1 a1 ✓
d/dt[w2²] = 2 w2 a2 ✓
d/dt[2 w1 w2 c] = 2 * d/dt[w1 w2 c]

d/dt[w1 w2 c] = a1 w2 c + w1 a2 c + w1 w2 (-s) (w2 - w1)
              = a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)

So d/dt[2 w1 w2 c] = 2 a1 w2 c + 2 w1 a2 c - 2 w1 w2 s (w2 - w1)

dT/dt = 0.5 * [4 w1 a1 + 2 w2 a2 + 2 a1 w2 c + 2 w1 a2 c - 2 w1 w2 s (w2 - w1)]
      = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)

This gives -w1 w2 s (w2 - w1). 

But the chain rule gives +w1 w2 s (w2 - w1). There's a sign discrepancy.

Wait, let me recheck the chain rule more carefully.

E = T + V
dE/dt = dT/dt + dV/dt

dV/dt = ∂V/∂θ1 * w1 + ∂V/∂θ2 * w2
       = (2 g sin θ1) w1 + (g sin θ2) w2
       = 2 g w1 sin1 + g w2 sin2

So dE/dt = dT/dt + dV/dt
         = [2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)] + [2 g w1 sin1 + g w2 sin2]

But from the chain rule:
dE/dt = ∂E/∂w1 * a1 + ∂E/∂w2 * a2 + ∂E/∂θ1 * w1 + ∂E/∂θ2 * w2
       = (2 w1 + w2 c) a1 + (w2 + w1 c) a2 + (-w1 w2 s + 2 g sin1) w1 + (w1 w2 s + g sin2) w2
       = 2 w1 a1 + w2 a1 c + w2 a2 + w1 a2 c - w1² w2 s + 2 g w1 sin1 + w1 w2² s + g w2 sin2

Comparing:
Direct: 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1) + 2 g w1 sin1 + g w2 sin2
Chain:  2 w1 a1 + w2 a2 + w2 a1 c + w1 a2 c + w1 w2 s (w2 - w1) + 2 g w1 sin1 + g w2 sin2

The difference is: -w1 w2 s (w2 - w1) vs +w1 w2 s (w2 - w1)

So one of them is wrong. Let me check the chain rule more carefully.

∂E/∂θ1 = ∂T/∂θ1 + ∂V/∂θ1

T = 0.5 (2 w1² + w2² + 2 w1 w2 cos(θ1-θ2))
∂T/∂θ1 = 0.5 * 2 w1 w2 * (-sin(θ1-θ2)) * 1 = -w1 w2 s

Wait, cos(θ1-θ2). d/dθ1[cos(θ1-θ2)] = -sin(θ1-θ2) * 1 = -s.
So ∂T/∂θ1 = 0.5 * 2 w1 w2 * (-s) = -w1 w2 s ✓

∂V/∂θ1 = -g * d/dθ1[2 cos θ1 + cos θ2] = -g * (-2 sin θ1) = 2 g sin θ1 ✓

So ∂E/∂θ1 = -w1 w2 s + 2 g sin θ1 ✓

∂E/∂θ2 = ∂T/∂θ2 + ∂V/∂θ2
∂T/∂θ2 = 0.5 * 2 w1 w2 * (-sin(θ1-θ2)) * (-1) = w1 w2 s
∂V/∂θ2 = -g * d/dθ2[2 cos θ1 + cos θ2] = -g * (-sin θ2) = g sin θ2
∂E/∂θ2 = w1 w2 s + g sin θ2 ✓

So the chain rule gives:
dE/dt = (2 w1 + w2 c) a1 + (w2 + w1 c) a2 + (-w1 w2 s + 2 g sin1) w1 + (w1 w2 s + g sin2) w2
      = 2 w1 a1 + w2 a1 c + w2 a2 + w1 a2 c - w1² w2 s + 2 g w1 sin1 + w1 w2² s + g w2 sin2
      = 2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c + w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

Now, the direct method:
dT/dt = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)
dV/dt = 2 g w1 sin1 + g w2 sin2
dE/dt = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1) + 2 g w1 sin1 + g w2 sin2

The difference is the sign of w1 w2 s (w2 - w1).

Wait, I think I see the issue. In the chain rule, I have:
∂E/∂θ1 = -w1 w2 s + 2 g sin1

But in the direct method, when I compute dT/dt, I get:
dT/dt = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)

The term -w1 w2 s (w2 - w1) comes from d/dt[2 w1 w2 cos(θ1-θ2)].

Let me expand this:
d/dt[2 w1 w2 cos(θ1-θ2)] = 2 [a1 w2 cos(θ1-θ2) + w1 a2 cos(θ1-θ2) + w1 w2 (-sin(θ1-θ2)) (w2 - w1)]
= 2 a1 w2 c + 2 w1 a2 c - 2 w1 w2 s (w2 - w1)

So dT/dt = 0.5 * [4 w1 a1 + 2 w2 a2 + 2 a1 w2 c + 2 w1 a2 c - 2 w1 w2 s (w2 - w1)]
         = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c - w1 w2 s (w2 - w1)

Now, in the chain rule, the ∂E/∂θ1 * w1 term includes -w1 w2 s * w1 = -w1² w2 s.
And the ∂E/∂θ2 * w2 term includes w1 w2 s * w2 = w1 w2² s.

So the chain rule gives: -w1² w2 s + w1 w2² s = w1 w2 s (w2 - w1).

But the direct method gives: -w1 w2 s (w2 - w1).

These have opposite signs! One of them must be wrong.

Let me check the chain rule again. 

∂E/∂θ1 = ∂T/∂θ1 + ∂V/∂θ1 = -w1 w2 s + 2 g sin1

But wait, ∂T/∂θ1 = ∂/∂θ1 [0.5 (2 w1² + w2² + 2 w1 w2 cos(θ1-θ2))]
= 0.5 * 2 w1 w2 * (-sin(θ1-θ2)) * 1
= -w1 w2 s

And ∂V/∂θ1 = ∂/∂θ1 [-g (2 cos θ1 + cos θ2)]
= -g * (-2 sin θ1)
= 2 g sin θ1

So ∂E/∂θ1 = -w1 w2 s + 2 g sin θ1 ✓

Now, the term ∂E/∂θ1 * w1 = (-w1 w2 s + 2 g sin1) * w1 = -w1² w2 s + 2 g w1 sin1

And ∂E/∂θ2 * w2 = (w1 w2 s + g sin2) * w2 = w1 w2² s + g w2 sin2

So the chain rule gives: -w1² w2 s + w1 w2² s = w1 w2 s (w2 - w1)

But the direct method gives: -w1 w2 s (w2 - w1)

There's definitely a sign error somewhere. Let me recheck the direct method one more time.

T = 0.5 (2 w1² + w2² + 2 w1 w2 cos(θ1-θ2))

Let u = θ1 - θ2. Then cos(θ1-θ2) = cos(u).

dT/dt = 0.5 * d/dt[2 w1² + w2² + 2 w1 w2 cos(u)]

d/dt[2 w1²] = 4 w1 a1 ✓
d/dt[w2²] = 2 w2 a2 ✓
d/dt[cos(u)] = -sin(u) * du/dt = -sin(θ1-θ2) * (w1 - w2) = -s * (w1 - w2) = s * (w2 - w1)

So d/dt[2 w1 w2 cos(u)] = 2 * d/dt[w1 w2 cos(u)]
= 2 * [a1 w2 cos(u) + w1 a2 cos(u) + w1 w2 (-sin(u)) * (w1 - w2)]
= 2 * [a1 w2 c + w1 a2 c + w1 w2 (-s) * (w1 - w2)]
= 2 a1 w2 c + 2 w1 a2 c + 2 w1 w2 s (w2 - w1)

Wait! I had -w1 w2 s (w2 - w1) before, but now I get +w1 w2 s (w2 - w1)!

Let me recheck. d/dt[cos(u)] = -sin(u) * du/dt = -s * (w1 - w2) = s * (w2 - w1).

So d/dt[w1 w2 cos(u)] = a1 w2 c + w1 a2 c + w1 w2 * [s * (w2 - w1)]
= a1 w2 c + w1 a2 c + w1 w2 s (w2 - w1)

And d/dt[2 w1 w2 cos(u)] = 2 a1 w2 c + 2 w1 a2 c + 2 w1 w2 s (w2 - w1)

So dT/dt = 0.5 * [4 w1 a1 + 2 w2 a2 + 2 a1 w2 c + 2 w1 a2 c + 2 w1 w2 s (w2 - w1)]
         = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c + w1 w2 s (w2 - w1)

This gives +w1 w2 s (w2 - w1), which matches the chain rule!

I made a sign error earlier. The correct dT/dt has +w1 w2 s (w2 - w1), not -w1 w2 s (w2 - w1).

So the correct dE/dt is:
dE/dt = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c + w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

Now let me check if this is zero with the correct equations of motion.

From (A): 2 a1 + a2 c + w2² s + 2 g sin1 = 0 (with L=1)
Multiply by w1: 2 w1 a1 + w1 a2 c + w1 w2² s + 2 g w1 sin1 = 0

From (B): a2 + a1 c - w1² s + g sin2 = 0 (with L=1)
Multiply by w2: w2 a2 + w2 a1 c - w2 w1² s + g w2 sin2 = 0

Add:
2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c + w1 w2² s - w2 w1² s + 2 g w1 sin1 + g w2 sin2 = 0

So:
2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c = -w1 w2² s + w2 w1² s - 2 g w1 sin1 - g w2 sin2

Now, dE/dt = [2 w1 a1 + w2 a2 + (w1 a2 + w2 a1) c] + w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

= [-w1 w2² s + w2 w1² s - 2 g w1 sin1 - g w2 sin2] + w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

= -w1 w2² s + w2 w1² s + w1 w2 s (w2 - w1) + (-2 g w1 sin1 + 2 g w1 sin1) + (-g w2 sin2 + g w2 sin2)

= -w1 w2² s + w2 w1² s + w1 w2 s (w2 - w1)

= w1 w2 s (-w2 + w1 + w2 - w1)
= w1 w2 s * 0
= 0

YES! dE/dt = 0! The equations DO conserve energy!

The error was in my dT/dt calculation. I had the wrong sign for the w1 w2 s (w2 - w1) term.

So the original equations ARE correct:
a1 = (-w1² s c + (g/L) sin2 c - w2² s - 2(g/L) sin1) / den
a2 = (2 w1² s + w2² s c + 2(g/L) sin1 c - 2(g/L) sin2) / den

And the correct dE/dt formula is:
dE/dt = 2 w1 a1 + w2 a2 + a1 w2 c + w1 a2 c + w1 w2 s (w2 - w1) + g (2 w1 sin1 + w2 sin2)

Now let me verify this numerically.

OK, so the dE/dt formula is now correct (gives 0), which proves the equations of motion are correct. But the energy is still not conserved in the numerical simulation. The energy at frame 0 (after the loop has already run) is 35.28, not 0.149.

Wait, there's a problem with my test. I'm reusing the same variables `t1, t2, w1, w2` across different tests. After the first dE/dt test, the state was changed. Then when I set t1, t2, w1, w2 for the energy conservation test, the loop starts from frame 0 but the variables have already been modified.

Actually, looking at the output more carefully:
- Frame 0 shows theta1=296.25deg, theta2=11.46deg — this is NOT the initial condition (120°, -10°)!

The issue is that I'm reusing variables. After the dE_dt test at arb state, t1=1.0, t2=-0.5, w1=2.0, w2=-1.5. Then I set t1=120*PI/180, t2=-10*PI/180, w1=0, w2=0. But then the loop runs and after 6000 steps, the state is at frame 6000. Then when I try to sample frames 0, 30, 60, etc., I'm starting from frame 6000, not frame 0.

I need to reset the state before each test. Let me fix the test.

Actually, the key insight is that the dE/dt formula now gives 0, which proves the equations are correct. The energy conservation issue in the simulation must be due to the RK4 integrator not being symplectic (it's a general-purpose integrator, not specifically designed for Hamiltonian systems).

But wait, RK4 should conserve energy quite well for a reasonable time step. The fact that energy is so badly conserved suggests there might still be a bug.

Let me redo the test properly, resetting the state for each test.

OK so the energy is still diverging. The dE/dt formula gives 0 (proving the equations are correct), but the numerical integration is not conserving energy.

Wait, this is very strange. If dE/dt = 0 analytically, then the energy should be conserved. But the numerical integration is showing energy growth.

Let me check if the issue is with the RK4 implementation. Let me try a much smaller time step and see if energy is conserved better.

Actually, let me check: the energy at frame 0 is 0.149 (correct). At frame 30 (0.02s), it's 0.149068 (very close). At frame 60 (0.05s), it's 0.149550 (still close). But by frame 120 (0.10s), it's 0.157447 (5.6% error). And by frame 240 (0.19s), it's 0.294764 (97% error).

This is a rapid energy growth. Something is wrong.

Let me check if the issue is with my RK4 implementation. Let me try a different integrator or a much smaller time step.

Actually, wait. Let me check the RK4 implementation more carefully. The issue might be that I'm not properly implementing RK4 for this system.

For a system y' = f(y), RK4 is:
k1 = f(y_n)
k2 = f(y_n + h/2 * k1)
k3 = f(y_n + h/2 * k2)
k4 = f(y_n + h * k3)
y_{n+1} = y_n + h/6 * (k1 + 2*k2 + 2*k3 + k4)

For our system, y = [θ1, θ2, ω1, ω2] and f(y) = [ω1, ω2, α1, α2].

So:
k1 = [ω1, ω2, α1(θ1,θ2,ω1,ω2), α2(θ1,θ2,ω1,ω2)]
k2 = [ω1 + h/2*α1, ω2 + h/2*α2, α1(θ1+h/2*ω1, θ2+h/2*ω2, ω1+h/2*α1, ω2+h/2*α2), ...]
...

Wait, this is what I'm doing. Let me check if there's a bug.

Actually, let me try a simple Euler method and see if energy is conserved.

Actually, let me try a symplectic integrator (like the Störmer-Verlet method) which should conserve energy better for Hamiltonian systems.

Let me implement the symplectic Euler method:
ω1_{n+1} = ω1_n + h * α1(θ1_n, θ2_n, ω1_n, ω2_n)
ω2_{n+1} = ω2_n + h * α2(θ1_n, θ2_n, ω1_n, ω2_n)
θ1_{n+1} = θ1_n + h * ω1_{n+1}
θ2_{n+1} = θ2_n + h * ω2_{n+1}

This is the symplectic Euler (semi-implicit Euler) method, which should conserve energy much better for Hamiltonian systems.

Actually, let me first check if the issue is with the time step. Let me use a much smaller DT.

Hmm, but I already showed that dE/dt = 0 analytically. So the equations are correct. The issue must be with the numerical integration.

Let me try symplectic Euler and see if energy is conserved.