Double pendulum · minimax-m3.1-flash

Let me write the double pendulum simulation with RK4 integration.

Equations of motion for double pendulum with equal masses m and equal lengths L, from downward vertical angles θ1, θ2.

Standard formulation (from the classic Wikipedia double pendulum equations), with coordinates:
x1 = L sinθ1, y1 = -L cosθ1
x2 = x1 + L sinθ2, y2 = y1 - L cosθ2

I'll use the well-known equations (e.g. from myphysicslab / Wikipedia "Double pendulum" article's coupled equations):

Let Δ = θ1 - θ2.

Using the formulation with moments of inertia and the simplified m1=m2=m, L1=L2=L:

I'll use the standard equations (e.g., from the Wikipedia "Double pendulum" article's coupled equations):

Let Δ = θ1 - θ2.

Using the formulation with moments of inertia and the simplified m1=m2=m, L1=L2=L:

Mass matrix:
A11 = (m1+m2) L1², A12 = m2 L1 L2 cos(Δ)
A22 = m2 L2²

Right side:
b1 = -m1 g L1 sinθ1 - m2 g L2 sin(θ2) cos(Δ) - m2 L1 L2 ω2² sin(Δ)
b2 = -m2 g L2 sinθ2 + m2 L1 L2 ω1² sin(Δ)

Then solve for α1, α2. That's the standard linear system. For m1=m2=m, L1=L2=L=1: A11 = 2L², A12 = L² cosΔ, A22 = L².

Let me just use a general implementation with m1,m2,L1,L2 parameters (equal ones though).

Solve 2x2:
det = A11*A22 - A12²
α1 = (b1*A22 - A12*b2)/det
α2 = (A11*b2 - A12*b1)/det

Check sign conventions: with θ measured from downward vertical, positive counterclockwise. The equations above come from the standard derivation where angles measured from downward vertical with x = L sinθ, y = -L cosθ. Yes, this is the standard.

Let me verify with a simple case: m2=0 (point mass on rod 1 only). Then A11 = m1 L1², A12=0, A22 = 0 → det=0. Bad. So keep the general form but with nonzero m2.

Better use the alternative formulation that avoids the degeneracy? With equal masses it's fine (m2 = m1 = 1).

Actually let me use the well-known simplified equal-mass equations to double check. Alternative: use the "energy equations" form:

For equal masses and lengths, there's a known system:
θ1'' = [ -g(2m)sinθ1 - m g sin(θ1-2θ2) - 2 m sin(θ1-θ2)(L2 ω2² + g cos(θ1-θ2)) ] / (L (2m - m cos²(θ1-θ2)))
Hmm, I recall from Wikipedia:

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

That's for equal masses and equal lengths L. This is cleaner and standard. I'll use this one but implement via the general mass-matrix form? The general form is more robust. But equal-mass form is simpler and well known. Let me just implement the general mass-matrix version (handles unequal too) and verify numerically against energy conservation in console? Simpler: implement equal-length equal-mass using the Wikipedia simplified equations, plus RK4 with dt small (1/2000 s), and verify energy conservation.

Actually let's implement the general mass matrix with m1=m2=1, L1=L2=1, g=9.81. Then I can numerically check energy conservation. Let me pick the general form carefully — I'll implement both and cross-check? Time budget 15 min. I'll implement general mass matrix and also compute total energy to display in HUD, then check energy drift is small in browser console.

Let me write the general form with the derivation-safe equations:

Using generalized coordinates θ1, θ2 from downward vertical, positions:
x1 = L1 sinθ1, y1 = -L1 cosθ1
x2 = x1 + L2 sinθ2, y2 = y1 - L2 cosθ2

Kinetic energy T = ½ m1 L1² ω1² + ½ m2 [L1²ω1² + L2²ω2² + 2 L1 L2 ω1 ω2 cos(θ1-θ2)]
Potential V = -m1 g L1 cosθ1 - m2 g (L1 cosθ1 + L2 cosθ2)

Equations (standard from e.g. Nelson & McNulty / myphysicslab):
M11 = (m1+m2) L1²
M12 = m2 L1 L2 cos(θ1-θ2)
M22 = m2 L2²

f1 = -m1 g L1 sinθ1 - m2 g L2 sinθ2 cos(θ1-θ2) - m2 L1 L2 ω2² sin(θ1-θ2)
f2 = -m2 g L2 sinθ2 + m2 L1 L2 ω1² sin(θ1-θ2)

Hmm, need to check consistency with equal-mass Wikipedia form. Let's just test numerically with the general form against the Wikipedia equal form via console in browser, or offline with node.

Let me do a quick node check.

Also note that with this "from downward vertical" convention, x = L sinθ, y = -L cosθ, angle increases counterclockwise. Wikipedia's equal-mass equations use angles from downward vertical too I think. Let's numerically test: I'll write the general form in node and compare with equal-mass Wikipedia form.

Plan: write a node script computing with both, compare over time.

Trail: array of points for lower bob, bounded (e.g. max 900 points), fading alpha based on index; draw as a polyline with per-segment alpha and color hue based on... "colorful fading trail" — use a gradient of hues along the trail.

Let me also add controls: Pause/Resume, Reset, maybe speed slider (optional). Keep required: pause/resume, reset. Add HUD showing angles, energy, time.

Physics time-stepping: use fixed dt = 1/1000 s of sim time, accumulate real time. RK4 integrator.

Let's write node test first.

Also handle: pause/resume button toggling; reset sets angles 120° and -10°, zero velocities, clears trail.

Also, the "bounded" trail: limit number of points.

Let me write the code.

For the "angle measured from downward vertical" — 120° means bob1 is up-left? θ=120° → x = sin120 = 0.866 (right), y = -cos120 = +0.5 (up, if y up). So bob 1 up-right. Fine.

Now implement.

Draw:
- background dark gradient
- pivot: small circle at origin
- rods: lines with gradient or stroke
- bobs: circles with color
- trail: polyline segments with alpha increasing toward the head, hue shifting.

Let me write the node verification script.

The general mass-matrix form is wrong (energy not conserved), and also the "wiki" equal form gives wrong energy (E0 computed with my energy function should match). E wiki = 10.6 vs E0 = 0.149 → wiki form as I wrote it is wrong (probably sign/convention of angles measured from downward vs upward). Let's derive properly myself.

Let me derive from Lagrangian with my convention: x = L sinθ, y = -L cosθ (y up), θ from downward vertical, positive counterclockwise? Check: θ=0 → (0,-L) below pivot ✓. θ=90° → (L, 0) to the right. Increasing θ goes counterclockwise in standard math orientation (x right, y up): from pointing down to pointing right is counterclockwise? From (0,-1) to (1,0): rotating counterclockwise by 90° takes (0,-1) to (1,0)? Rotation by +90°: (x,y)→(-y,x): (0,-1)→(1,0). Yes counterclockwise. Good.

Derive Lagrange with q = θ1, θ2.

Positions:
r1 = (L1 s1, -L1 c1)
r2 = (L1 s1 + L2 s2, -L1 c1 - L2 c2)

Velocities: dr1/dt = L1 c1 ω1, L1 s1 ω1
dr2/dt = (L1 c1 ω1 + L2 c2 ω2, L1 s1 ω1 + L2 s2 ω2)

T = ½ m1 |dr1|² + ½ m2 |dr2|² = ½ m1 L1² ω1² + ½ m2 (L1²ω1² + L2²ω2² + 2 L1 L2 ω1ω2 (c1c2 + s1s2)) = ... c1c2+s1s2 = cos(θ1-θ2) ✓.

V = m1 g y1 + m2 g y2 = -m1 g L1 c1 - m2 g (L1 c1 + L2 c2).

Now Lagrange equations. Define
T = ½ A ω1² + ½ B ω2² + C ω1 ω2, with A = (m1+m2)L1², B = m2 L2², C = m2 L1 L2 cos Δ, Δ = θ1-θ2.

∂T/∂ω1 = A ω1 + C ω2; d/dt = A ω1' + C' ω2 + C ω2', C' = -m2 L1L2 sinΔ ω1.
∂T/∂θ1 = -C sinΔ ω1ω2 (since d/dθ1 cos Δ = -sinΔ). Also ∂T/∂2 = +C sinΔ ω1ω2 (d/dθ2 cosΔ = +sinΔ).
∂V/∂θ1 = m1 g L1 s1 + m2 g L1 s1 = g L1 (m1+m2) s1.
∂V/∂θ2 = m2 g L2 s2.

Equation 1: d/dt ∂T/∂ω1 - ∂T/∂θ1 + ∂V/∂θ1 = 0:
A ω1' + C ω2' - m2 L1L2 sinΔ ω1 ω2 + C sinΔ ω1ω2 + gL1(m1+m2) s1 = 0
→ A ω1' + C ω2' + gL1(m1+m2) s1 - m2 L1L2 sinΔ ω1 ω2 = 0.

Hmm, note A ω1' term: A = (m1+m2)L1², so gL1(m1+m2)s1 = (m1+m2) g L1 s1.

Equation 2: d/dt ∂T/∂ω2 - ∂T/∂θ2 + ∂V/∂θ2 = 0:
B ω2' + C ω1' + m2 L1L2 sinΔ ω1 ω2 - C sinΔ ω1ω2 + m2 g L2 s2 = 0
→ C ω1' + B ω2' + m2 g L2 s2 = 0.

So the equations in matrix form (M ω' = f):
A ω1' + C ω2' = m2 L1L2 sinΔ ω1ω2 - (m1+m2) g L1 s1
C ω1' + B ω2' = -m2 g L2 s2

Interesting: velocity coupling term only in eq1, and RHS of eq2 has no velocity term. Hmm, but symmetry seems broken. Actually the "m2 L1 L2 sinΔ ω1ω2" appears in eq1 only because ∂T/∂θ1 = -C sinΔ ω1ω2 contributes and ∂T/∂θ2 = +C sinΔ ω1ω2 also contributes but cancels in eq2's d/dt term? Let's redo eq2:

d/dt (B ω2 + C ω1) = B ω2' + C' ω1 + C ω1' = B ω2' - m2L1L2 sinΔ ω1 ω2 + C ω1'.
∂T/∂θ2 = +C sinΔ ω1 ω2.
So eq2: B ω2' - m2L1L2 sinΔ ω1ω2 + C ω1' - C sinΔ ω1ω2 + m2 g L2 s2 = 0.

Hmm wait, ∂V/∂θ2 = d/dθ2 [-m2 g L2 c2] = m2 g L2 s2 ✓. And in eq1: ∂V/∂1 = m1 g L1 s1 + m2 g L1 s1 = (m1+m2) g L1 s1 ✓.

So eq1: A ω1' + C ω2' + [(-m2L1L2 sinΔ) + C sinΔ] ω1ω2 + (m1+m2) g L1 s1 = 0
where the ω1' term in d/dt: d/dt(Aω1 + Cω2) = A ω1' + C' ω2 + C ω2' = Aω1' - m2L1L2 sinΔ ω1ω2 + Cω2'.
∂T/∂θ1 = -C sinΔ ω1ω2 → subtract → +C sinΔ ω1ω2.
So coefficient of ω1ω2 in eq1: -m2L1L2 sinΔ + C sinΔ = sinΔ(C - m2L1L2) = m2L1L2 sinΔ (cosΔ - 1).

eq2 coefficient: -m2L1L2 sinΔ - C sinΔ = -m2L1L2 sinΔ(1+cosΔ).

So:
A ω1' + C ω2' = -m2L1L2 sinΔ(cosΔ-1) ω1ω2 - (m1+m2) g L1 s1
C ω1' + B ω2' = m2L1L2 sinΔ(1+cosΔ) ω1ω2 - m2 g L2 s2

Hmm, that's asymmetric in a suspicious way. Let me sanity check symmetry: swap the roles? The system isn't symmetric under swapping arms in general, so maybe fine. But let me verify with an independent approach: numerically derive equations from Lagrange via autodiff-ish check? Better: verify energy conservation numerically with this derived system — if the equations are right, energy is conserved.

Let me code it and test energy conservation. Also compare against a brute-force Lagrangian computation with finite differences? Energy conservation test is a good check.

Actually there's a cleaner check: use the general form derived via symbolic differentiation of T and V with a numeric Lagrangian: d/dt(∂L/∂ω) - ∂L/∂q, where ∂/∂ω computed analytically and ∂/∂q computed numerically by finite difference on T and V. That validates my derivation. Let's do that.

Wait — but also the "cross" term sign. Let me recheck T cross term: |dr2|² = (L1c1ω1+L2c2ω2)² + (L1s1ω1+L2s2ω2)² = L1²ω1²(c1²+s1²) + L2²ω2² + 2L1L2ω1ω2(c1c2+s1s2) = ... + 2L1L2ω1ω2 cos(θ1-θ2) ✓.

C = m2 L1 L2 cosΔ. ∂T/∂ω1 = A ω1 + C ω2 where A=(m1+m2)L1² ✓.

Fine, let's test numerically.

Also note: maybe I should just use a cleaner known-good approach: compute accelerations from the Lagrangian using symbolic derivatives via a simple CAS-ish approach: define T and V as functions of θ,ω and compute ∂/∂θ numerically and ∂/∂ω analytically, then integrate the Lagrangian equations with numerical time derivative. That's messy. Energy conservation test with the derived equations is enough, plus check against the Wikipedia known-good equations for equal mass/length (which are in terms of angles measured from... the wiki double pendulum article uses angles from vertical downward? Let me recall the wiki equations:

For double pendulum with equal masses and lengths, from Wikipedia (Double pendulum), the equations of motion in terms of angles from vertical:
θ1'' = ( -g(2 sinθ1 + sin(θ1-2θ2)) - 2 sin(θ1-θ2)(L2 θ2'² + g cos(θ1-θ2)) ) / (L(2 - cos²(θ1-θ2)))
θ2'' = ( 2 sin(θ1-θ2)(L1 θ1'² + g cos(θ1-θ2)) ) / (L(2 - cos²(θ1-θ2)))

That's for angles measured from vertical (down). My test showed energy not conserved with my energy function, so either the wiki form or my energy function/angle convention is off. Note wiki's θ2'' has no gravity term?? Actually the wiki form appears to have them in a "reduced" form where... hmm, that second equation looks like it lacks gravity — likely that wiki snippet is for the version after some manipulation (the "θ1'' and θ2''" formulas in the Wikipedia article for equal masses with angles measured from the vertical). Actually I recall the wiki article's equations are for the case of angles measured from vertical and the second equation does contain gravity: θ2'' = [2 sin(θ1-θ2)(L1 θ1'² + g cos(θ1-θ2))] / (L(2-cos²(θ1-θ2))) — hmm no gravity. That corresponds to angles measured from the vertical but where θ2 is measured relative to θ1? Whatever.

Let me just test my derived system numerically for energy conservation, and also build an independent "brute force" integrator: integrate Newton's laws via constraint? Simplest independent check: numerically compute Lagrangian equations with automatic differentiation of T, V using complex-step or central differences, and check the resulting accelerations match my analytic ones. Then energy conservation with RK4 to confirm correctness overall.

Let's do it.

My analytic derivation has an error (energy blew up). Let me carefully re-derive using symbolic computation with sympy? Not available maybe. Let me use a proper numeric Lagrangian formulation instead: the generalized equations

M(q) q'' + C(q,q') + G(q) = 0.

Better approach: use the standard known-correct mass-matrix form. Let me derive it carefully.

T = ½ A ω1² + ½ B ω2² + C ω1 ω2, where A=(m1+m2)L1², B=m2L2², C=m2L1L2 cosΔ, Δ=θ1-θ2.

Note: standard textbook (Nelson & McNulty) uses:
T = ½ m1 L1² ω1² + ½ m2 (L1²ω1² + L2²ω2² + 2L1L2ω1ω2 cos(θ1-θ2))
V = -(m1+m2) g L1 cosθ1 - m2 g L2 cosθ2

Equations (Nelson/McNulty form with angles measured from vertical down):
M11 = (m1+m2)L1²
M12 = m2 L1 L2 cos(θ1-θ2)
M22 = m2 L2²

b1 = -m1 g L1 sinθ1 - m2 g L2 sinθ2 cos(θ1-θ2) - m2 L1 L2 ω2² sin(θ1-θ2)
b2 = -m2 g L2 sinθ2 + m2 L1 L2 ω1² sin(θ1-θ2)

That's the standard "inverse mass matrix" form from Wikipedia's double pendulum article (the version with matrix M and vector b). Actually Wikipedia gives:

M11 = (m1+m2) l1²
M12 = m2 l1 l2 cos(θ1-θ2)
M22 = m2 l2²
b1 = -(m1+m2) g l1 sinθ1 - m2 l2 ω2² l1 sin(θ1-θ2)
b2 = -m2 g l2 sinθ2 + m2 l1 ω1² l2 sin(θ1-θ2)

Yes! Wikipedia's "Double pendulum" article (for the case with angles measured from vertical) gives exactly:
b1 = -(m1+m2) g l1 sin θ1 - m2 l2 ω2² l1 sin(θ1-θ2)
b2 = -m2 g l2 sin θ2 + m2 l1 ω1² l2 sin(θ1-θ2)

This matches my earlier "general" test but I had an extra term -m2 g L2 sinθ2 cos(Δ) in f1 (that's from a different variant). Let me test Wikipedia's version exactly.

Wait, my earlier test had f1 = -m1 g L1 sinθ1 - m2 g L2 sinθ2 cos Δ - m2 L1L2 ω2² sinΔ, f2 = -m2 g L2 sinθ2 + m2 L1L2 ω1² sinΔ. Energy drifted a lot (max 17). Hmm, that suggests errors, or the issue was the coupling term missing. Let's check Wikipedia form: b1 = -(m1+m2) g l1 sinθ1 - m2 l1 l2 ω2² sin(θ1-θ2); b2 = -m2 g l2 sinθ2 + m2 l1 l2 ω1² sin(θ1-θ2).

Hmm, but is that right? Let's derive properly myself rather than trust memory.

Derivation (standard):
Let Δ = θ1-θ2.
Positions (x right, y up), angles from downward vertical:
x1 = l1 sinθ1, y1 = -l1 cosθ1
x2 = x1 + l2 sinθ2, y2 = y1 - l2 cosθ2

Velocities:
ẋ1 = l1 cosθ1 ω1, ẏ1 = l1 sinθ1 ω1
ẋ2 = l1 cosθ1 ω1 + l2 cosθ2 ω2, ẏ2 = l1 sinθ1 ω1 + l2 sinθ2 ω2

T = ½m1(l1²ω1²) + ½m2[ (l1 cosθ1 ω1 + l2 cosθ2 ω2)² + (l1 sinθ1 ω1 + l2 sinθ2 ω2)² ]
= ½m1 l1²ω1² + ½m2[ l1²ω1² + l2²ω2² + 2 l1 l2 ω1 ω2 cos(θ1-θ2) ]

V = m1 g y1 + m2 g y2 = -m1 g l1 cosθ1 - m2 g (l1 cosθ1 + l2 cosθ2)

Lagrangian: d/dt(∂L/∂ωi) - ∂L/∂θi = 0.

∂L/∂ω1 = (m1+m2) l1² ω1 + m2 l1 l2 cosΔ ω2
∂L/∂θ1 = -m2 l1 l2 sinΔ ω1ω2 + (m1+m2) g l1 sinθ1
  since d/dθ1 cosΔ = -sinΔ, so ∂T/∂θ1 = ½ m2 * 2 l1l2 ω1ω2 * (-sinΔ) = -m2 l1 l2 sinΔ ω1ω2. And ∂V/∂θ1 = (m1+m2) g l1 sinθ1. So ∂L/∂θ1 = ∂T/∂θ1 - ∂V/∂θ1 = -m2 l1 l2 sinΔ ω1ω2 - (m1+m2) g l1 sinθ1.

d/dt(∂L/∂ω1) = (m1+m2) l1² ω1' + m2 l1 l2 [ -sinΔ (ω1-ω2) ω2 + cosΔ ω2' ]

Equation 1: (m1+m2) l1² ω1' + m2 l1 l2 cosΔ ω2' - m2 l1 l2 sinΔ ω1ω2 + m2 l1 l2 sinΔ ω2² + m2 l1 l2 sinΔ ω1ω2 + (m1+m2) g l1 sinθ1 = 0

→ (m1+m2) l1² ω1' + m2 l1 l2 cosΔ ω2' = -(m1+m2) g l1 sinθ1 - m2 l1 l2 ω2² sinΔ

Good, matches Wikipedia b1 = -(m1+m2) g l1 sinθ1 - m2 l1 l2 ω2² sinΔ. ✓ (my earlier test used m1 instead of m1+m2 for the first term — bug!)

∂L/∂ω2 = m2 l2² ω2 + m2 l1 l2 cosΔ ω1
d/dt = m2 l2² ω2' + m2 l1 l2 [ sinΔ (ω1-ω2) ω1 + cosΔ ω1' ]
∂L/∂θ2 = +m2 l1 l2 sinΔ ω1ω2 - m2 g l2 sinθ2
 (d/dθ2 cosΔ = +sinΔ since Δ=θ1-θ2 → d cosΔ/dθ2 = -sinΔ * (-1) = +sinΔ; so ∂T/∂θ2 = +m2 l1 l2 sinΔ ω1ω2; ∂V/∂θ2 = m2 g l2 sinθ2; ∂L/∂θ2 = m2 l1 l2 sinΔ ω1ω2 - m2 g l2 sinθ2.)

Equation 2: m2 l2² ω2' + m2 l1 l2 sinΔ ω1² - m2 l1 l2 sinΔ ω1ω2 + m2 l1 l2 cosΔ ω1' - m2 l1 l2 sinΔ ω1ω2 + m2 g l2 sinθ2 = 0

→ m2 l1 l2 cosΔ ω1' + m2 l2² ω2' = m2 l1 l2 ω1² sinΔ - m2 g l2 sinθ2

Matches Wikipedia b2. ✓

So the correct matrix M and b:
M = [[(m1+m2)l1², m2 l1 l2 cosΔ],[m2 l1 l2 cosΔ, m2 l2²]]
b = [-(m1+m2) g l1 sinθ1 - m2 l1 l2 ω2² sinΔ, -m2 g l2 sinθ2 + m2 l1 l2 ω1² sinΔ]

Note the earlier "wiki" deriv2 form had gravity term in θ1'' as -g(2 sinθ1 + sin(θ1-2θ2)) — that's from a different (equivalent?) formulation. Whatever; let's use this verified matrix form and test energy conservation.

Also note signs of angles: my derived convention x=l sinθ, y=-l cosθ matches angles measured from downward vertical with positive counterclockwise. Fine.

Let metest.

Wait, this new formulation: I should double check the standard equations of motion for double pendulum (Wikipedia). The standard form for equal masses m and lengths L with angles measured from downward vertical:

M11 = 2mL², M12 = mL² cos(θ1-θ2), M22 = mL².
b1 = -2mgL sinθ1 - mL² θ2'² sin(θ1-θ2)
b2 = -mgL sinθ2 + mL² θ1'² sin(θ1-θ2)

Hmm, but Wikipedia's version uses angles measured from... Actually the classic Wikipedia "Double pendulum" equations with θ measured from vertical upward? Let me just check energy conservation numerically — that's the true test. The formulas I wrote here: I derived earlier general form and got energy drift huge, meaning sign/derivation issue. Let me test this new one.

Actually let's derive carefully.

Positions (x right, y up), angles from downward vertical, θ increases counterclockwise:
x1 = L1 sinθ1, y1 = -L1 cosθ1
x2 = x1 + L2 sinθ2, y2 = y1 - L2 cosθ2

Hmm, with y up and downward vertical at θ=0, then increasing θ moves bob to the right (x positive). So x = L sinθ, y = -L cosθ. Fine.

Velocities: dx1/dt = L1 cosθ1 θ1', dy1/dt = L1 sinθ1 θ1'. Speed1² = L1²θ1'².
dx2/dt = L1 cosθ1 θ1' + L2 cosθ2 θ2', dy2/dt = L1 sinθ1 θ1' + L2 sinθ2 θ2'.
Speed2² = L1²θ1'² + L2²θ2'² + 2L1L2θ1'θ2'(cosθ1cosθ2 + sinθ1 sinθ2) = L1²θ1'²+L2²θ2'²+2L1L2θ1'θ2'cos(θ1-θ2).

T = ½ m1 L1²θ1'² + ½ m2 [L1²θ1'² + L2²θ2'² + 2L1L2θ1'θ2'cos(θ1-θ2)]
V = m1 g y1 + m2 g y2 = -m1 g L1 cosθ1 - m2 g (L1 cosθ1 + L2 cosθ2)

Lagrange equations. Let m1=m2=m, L1=L2=L.

∂T/∂θ1' = m1 L1²θ1' + m2(L1²θ1' + L1L2θ2'cos(θ1-θ2)) = (m1+m2)L1²θ1' + m2L1L2θ2'cosΔ where Δ=θ1-θ2.
d/dt: (m1+m2)L1²θ1'' + m2L1L2[θ2''cosΔ - θ2' sinΔ (θ1'-θ2')]

∂T/∂θ1 = m2 L1L2 θ1'θ2' * (-sinΔ)
∂V/∂θ1 = (m1+m2) g L1 sinθ1
Lagrange: d/dt(∂T/∂θ1') - ∂T/∂θ1 + ∂V/∂θ1 = 0:
(m1+m2)L1²θ1'' + m2L1L2 cosΔ θ2'' - m2L1L2 sinΔ θ1'θ2' + m2L1L2 sinΔ θ1'θ2' + (m1+m2)gL1 sinθ1 = 0

Wait: -m2L1L2 sinΔ (θ1'-θ2')θ2' = -m2L1L2 sinΔ θ1'θ2' + m2L1L2 sinΔ θ2'². And -∂T/∂θ1 = +m2L1L2 sinΔ θ1'θ2'. These cancel the -m2L1L2 sinΔ θ1'θ2'. So we get:

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

So M11 θ1'' + M12 θ2'' = -m2L1L2 sinΔ θ2'² - (m1+m2) g L1 sinθ1.

Hmm, my check3 had b1 = -(m1+m2)gL1 sinθ1 - m2L1L2 w2² sinΔ — that matches! Good.

Second: ∂T/∂θ2' = m2 L2²θ2' + m2L1L2θ1'cosΔ
d/dt = m2L2²θ2'' + m2L1L2[θ1''cosΔ - θ1' sinΔ(θ1'-θ2')]
∂T/∂θ2 = m2L1L2θ1'θ2'*(+sinΔ)  [since d/dθ2 cos(θ1-θ2)= +sinΔ]
∂V/∂θ2 = m2 g L2 sinθ2
Lagrange: d/dt(∂T/∂θ2') - ∂T/∂θ2 + ∂V/∂θ2 = 0:
m2L2²θ2'' + m2L1L2cosΔθ1'' - m2L1L2 sinΔ θ1'² + m2L1L2 sinΔ θ1'θ2' - m2L1L2 sinΔ θ1'θ2' + m2gL2 sinθ2 = 0

→ m2L1L2 cosΔ θ1'' + m2L2²θ2'' = m2L1L2 sinΔ θ1'² - m2 g L2 sinθ2.

Matches check3's b2. Good, so check3 should conserve energy. Let's run.

Earlier check.js used a different (wrong) formulation with cosΔ terms in b1. Fine.

Also note the "mass matrix" determinant: M11*M22 - M12² = 2mL²·mL² - (mL²cosΔ)² > 0. Good.

Run check3.

The mass-matrix formulation with b1 = -(m1+m2) g L1 sinθ1 - m2 L1 L2 ω2² sin(Δ), b2 = -m2 g L2 sinθ2 + m2 L1 L2 ω1² sin(Δ) — wait I need to verify. Actually the standard Wikipedia double pendulum equations for equal masses/lengths... Let me recall the known form (from Wikipedia, using angles from downward vertical, x = L sinθ, y = -L cosθ):

Standard Lagrangian derivation gives:
(m1+m2) L1 θ1'' + m2 L2 θ2'' cos(θ1-θ2) = -m2 L2 θ2'² sin(θ1-θ2) - (m1+m2) g sinθ1
m2 L2 θ2'' + m2 L1 θ1'' cos(θ1-θ2) = m2 L1 θ1'² sin(θ1-θ2) - m2 g sinθ2

So b1 should be -(m1+m2) g L1 sinθ1 - m2 L2² ... hmm units. Let's check: with my check3, M11 = (m1+m2)L1², M12 = m2 L1 L2 cosΔ, M22 = m2 L2². b1 = -(m1+m2) g L1 sinθ1 - m2 L1 L2 ω2² sinΔ. Yes matches the standard form (multiply standard by L1). b2 = -m2 g L2 sinθ2 + m2 L1 L2 ω1² sinΔ. Matches. Energy drift tiny → correct.

Wait, but earlier check.js gave large drift because b1 there was wrong (had extra terms). Good, check3 is verified: energy drift 2e-9 at dt=0.0005, rod lengths constant.

Also note in check3 I used angles measured... initial 120 and -10 degrees from downward vertical. Positions: x = L sinθ, y = -L cosθ so θ=0 → bob below pivot. Good.

Now build the HTML app. Use RK4 with dt = 1/1000 s of sim time, accumulator for real-time playback. Trail: bounded array of lower bob positions, colored gradient fading.

Canvas: use devicePixelRatio scaling. Trail drawn as fading segments with hue based on index.

Let me write index.html.

Design: dark background, HUD with angle/energy readouts, Pause/Resume and Reset buttons.

Trail: store points; draw with alpha increasing toward newest, hue shifting from blue→magenta etc. Bounded to e.g. 400 points, remove oldest.

Also fixed pivot at top.

Let's write it.

Rendering scale: canvas 1200x800 area? Let's use full-window canvas with device pixel ratio, and scale = min(w,h)*0.35 / L1 (L1=1). Pivot at (w/2, h*0.2).

Let me write code carefully.

Trajectory sampling: add point each frame (or each physics substep accumulated with distance threshold). Use distance threshold to avoid duplicates when slow.

Pause/Resume toggles state; Reset restores angles and clears trail.

Also should support keyboard space.

Let me write.

For "colorful fading trail": draw segments with strokeStyle hsl(hue, 90%, 60%) alpha from 0.05 to 1.

I'll implement with a single path per segment (200 segments fine) or build path grouped by hue buckets. Simpler: loop segments, set globalAlpha and strokeStyle each, stroke separately. 300 strokes per frame is fine.

Also add a glow.

Write file.

M2 not defined — probably a hoisting issue? `const M1 = M2 = 1.0;` — M2 isn't declared. Fix: `const M1 = 1.0, M2 = 1.0;`.

Same issue with L1 = L2 = 1.0. Fix.

Wait — the "const M1 = M2 = 1.0" actually defines M1 only in module scope... yes, M2 was undeclared → ReferenceError in strict mode. Fix L too.

Also note that "L2 is not defined" means const L1 = L2 = 1.0 → L1 defined, L2 undeclared. Fix both.

Also energy formula: I used V = -(M1+M2)*G*L1*cos(th1) - M2*G*L2*cos(th2). Check: potential energy = -(m1 g L1 cos th1) - m2 g (L1 cos th1 + L2 cos th2) = -(m1+m2) g L1 cos th1 - m2 g L2 cos th2. Correct.

But my check3 script used the b1 form with -(m1+m2)G L1 sin th1 — wait, in check3 I had b1 = -(m1+m2)*g*L1*Math.sin(th1) - ... and M11=(m1+m2)L1². The standard formulation: mass matrix A11=(m1+m2)L1², A12=m2 L1 L2 cosΔ, A22 = m2 L2². And b1 = -(m1+m2) g L1 sin th1 - m2 L1 L2 ω2² sinΔ, b2 = -m2 g L2 sin th2 + m2 L1 L2 ω1² sinΔ. Yes that matches check3 which conserved energy. Good.

Hmm, but check3 used this and my check.js "general" version differed (used m1 g L1 sin th1 and extra term) which was the wrong one. check3 is correct.

Wait, is b1 correct? Standard: for the double pendulum with equal lengths, Lagrangian derivation gives:
(m1+m2)L1²θ1'' + m2 L1L2 cosΔ θ2'' = -(m1+m2) g L1 sinθ1 - m2 L1L2 sinΔ ω2²
m2 L1L2 cosΔ θ1'' + m2 L2² θ2'' = -m2 g L2 sinθ2 + m2 L1L2 sinΔ ω1²
Yes. Energy conserved per check3. Good.

But note b1 in check3 has -(m1+m2)g L1 sinθ1 — with m1+m2 both gravity terms. Fine.

Also the app's energy(): V uses -(M1+M2)G L1 cos th1 - M2 G L2 cos th2. Matches.

Now fix L declaration.

It works. Note the pivot is at x=600 (screen center) but the pendulum appears to the right — because initial θ1=120° from downward vertical: x = L sin(120°)=0.866 (right), y = -cos(120°)= +0.5 (up). So bob1 is up-right. That's correct per spec (angles from downward vertical). Fine.

Also note the sim time is very slow (0.107s after several frames?) Actually screenshots are quick. Let me wait and observe motion and trail. Let's take a motion/frames to check trail.

Wait — the screenshot shows sim time 0.107s but that render was likely captured early. Let's take another screenshot after some time.

Also potential issue: the trail's initial point is the lower bob but the trail is drawn from the lower bob; initial trail has only 4 points visible. Fine.

Let me check the motion over time.

One concern: trail min distance 0.012 world units with scale ~ 250 px/unit → 3px. Good.

Let's screenshot again after delay.

Looks great. Trail is bounded at 420 points, energy conserved. Now test pause/resume and reset.

Note: the trail's faded tail is a bit thin/dark — fine, it's fading. Colors: violet→cyan→gold. Actually hue formula 250 - 190*f gives 250 (violet) at old → 60 (gold) at new. Good.

Test pause.

Paused (button says Resume). Now verify it's frozen: take another screenshot after a moment and check sim time stayssame.

Paused holds state. Now resume, then reset.

Note: the trail head glow appears at the lower bob, but the newest trail color (gold) region overlaps the bob; looks fine.

Also note the upper bob appears at a position where the lower rod goes to the left of the lower bob? Actually the drawing shows pivot → red bob (down-right), then red bob → gold bob down-left. Correct.

Let's click Resume and verify it advances, then Reset.

Resumed (time advancing). Now reset.

Note the trail head glow circle appears at lower bob; good. Now click Reset.

One concern: the trail doesn't seem to show older colors (violet) — the oldest segments are very faint. Fine.

Click Reset at (146, 93).

Reset works: trail cleared, angles back to ~120/-10, energy back to initial 0.149036.

One issue: the trail's older portion (violet) isn't visible in the screenshot — hue goes from 250 to 60; the visible path appears teal/green. That's fine, colorful fading trail present. But the fading alpha: oldest segments nearly invisible. OK.

Also the trail head glow is at the lower bob. Good.

One concern: the trail's brightest segment is near the head — but head colors: hue = 250 - 190*f, f=1 → hue 60 (gold/yellow). The head near the bob appears gold. Good.

Let me verify rod lengths constant in-browser by measuring pixel distances? Already validated in node. Also verify coupled motion sample: the upper bob clearly moves — yes.

Minor: "Trail fades from dim violet → bright cyan → warm gold" but visible shows green/teal. Fine.

Let me check motion frames to confirm no freeze and that trail follows lower bob. Also check trail bounded at 420 — yes.

Potential improvement: trail oldest segments are nearly invisible (alpha ~0.02). Maybe raise min alpha a bit. Let's slightly increase min alpha to 0.08 and adjust hue so more colorful range visible. Also the trail currently spans a wide region.

I'll tweak alpha min to 0.06.

Also check: initial pose drawing after reset shows pivot at center and upper bob up-right at 120°, correct.

Let me also run a quick check that the trail head is exactly at the lower bob (glow circle visible at bob). Yes.

Small tweak then final verify.

Also the reset trail should be cleared — currently after reset trail has 1 point (the initial lower bob position), which is fine (cleared). The screenshot shows trail points 1-2. Good.

Also, note the lower bob trail with the hue mapping: f=1 -> hue 60 (gold), f=0 -> 268 (violet). Good.

Let me reload and check the trail colors, then also verify a keyboard shortcut. Also verify that reset restores pose: after reset, th1=120, th2=-10 — screenshot showed 119.94/-10.03 which is right after a moment of simulation. Good.

Reload and screenshot mid-motion.

Trail at start is short; fine. Let's let it run and take screenshot to check the trail colors now. Also note: the trail hue mapping — oldest violet (268), newest gold (60). Good.

Let me wait for a few seconds and screenshot.

Also the trail "fading" - older segments have low alpha. Let's check with motion frames.

Motion looks good, trail colorful and fading, energy conserved. But the trail's oldest part is barely visible (very faint) — acceptable as "fading". Actually the trail appears to only show recent ~ some portion; older segments fade to near invisible blue. That's fine.

One issue: the trail hue mapping goes 268 (violet) -> 60 (gold). But visually the older part shows as teal/blue and recent as yellow-green. Fine.

Also I notice the trail doesn't span the full path in early frames (bounded at 420). Fine.

Let me check that the trail head matches the lower bob position — yes, the lower bob is at the trail head.

Also verify reset clears trail and restores pose. Earlier reset showed trail points 2 and angles 118/-10 at t=0.065 — reset worked (already running after reset so angles advanced slightly). Good.

One more: verify initial pose at exactly 120/-10 by pausing right after reset. Let me use keyboard R while paused, and inspect immediately.

Also, the fading old segments are quite dim; maybe increase visibility slightly. Let's bump minimum alpha a bit.

Let me also double-check the trail min distance sampling: trail points 420 max bounded. Good.

Test: press R thencheck readouts.

Everything works. The trail hue goes violet→blue→cyan→green→gold but visually the earliest part (violet) appears faint blue. Fine.

One issue: the trail's oldest portion is very dim/near invisible (faint blue at right). That's fine as "fading". But maybe increase minimum visibility. Current alpha floor 0.08 — visible in screenshot as faint line. OK.

Also note the pendulum hangs mostly above pivot? In screenshot the upper bob is below pivot — fine.

Reset test: click Reset and confirm trail cleared and pose = 120/-10.

Also verify initial pose after reset visually: th1=120° means upper rod points up-left... Actually th1=120° from downward vertical → sin(120°)>0 so x1 positive (right), y1 = -cos(120°) = +0.5 up. So bob up-right. In screenshot at reset earlier, the red bob was up-right of pivot. Good.

Let's do a final reset check with a fresh reload and pause immediately to inspect initial pose.

Actually the reset earlier showed th1 118.79 after some frames — because it resumed running. Fine.

Let's confirm reset clears trail: click reset and immediately check trail points count.

Also check: after reset, does running stay as before? I set setRunning(true) on reset. Task says reset restores initial conditions and clears trail; forcing resume is acceptable but maybe better to keep paused state? Requirement: "Reset restores the stated initial conditions and clears the trail." Keeping it running is fine.

Let me click reset and screenshot.

Also clean up check js files (they're extra files; the task only requires index.html but extra files probably fine). I'll remove them to keep the dirclean.

Everything works. Reset restored initial pose (120°, -10°) and cleared trail (1 point). Energy conserved at 0.149036.

One issue: the "Reset" sets running=true — fine.

Also note: after reset, trail shows 1 point but at frame it was 2 — fine.

Minor: the trail color palette — oldest is faint violet/blue, head gold. In screenshot the visible trail was green→yellow, which is a bit of the range. Acceptable — "colorful fading trail". Let me make the hue range cover more rainbow for clarity: hue from 280 down to 40 currently = violet→blue→cyan→green→gold. That's fine.

Clean up helper js files.

Also one detail: the trail bound is 420 points, but the initial trail after reset only has 1 point — good.

Let me delete check files and finish. Actually keeping them may be fine but the task says write index.html; extra files could be noise. Removethem.