Double Pendulum

DEFINITION

two rods, two bobs, one pivot
θ̈₁, θ̈₂ from the Euler–Lagrange equations
no closed form — the motion must be integrated

NOTES

One pendulum is the clock; two, hinged together, are chaos you can build with a hacksaw. The equations fall straight out of the Lagrangian and cannot be solved — only integrated — and two starts a hair apart diverge until they share nothing but energy. This is the demonstration piece of sensitive dependence: the trail the lower bob leaves is a chaotic attractor being drawn by actual mechanics, live.

PROVENANCE

Origin

The compound pendulum goes back to the Bernoullis and Euler; the Lagrangian formulation is J.-L. Lagrange, "Mécanique analytique", 1788

Standing

Public domain — classical mechanics

Constants

PHASE varies the launch angle, so no two plates swing alike

HOUDINI · VEX

The same published mathematics as a Detail Wrangle body. Paste it into a Wrangle with Run Over set to Detail; every constant is the published value plus a tweak channel, so Create Spare Parameters gives a slider that starts where the paper does.

// FORMA — PL. 60 · DOUBLE PENDULUM — Classical Lagrangian mechanics · after Euler and Lagrange
//   two rods, two bobs, one pivot
//   θ̈₁, θ̈₂ from the Euler–Lagrange equations
//   no closed form — the motion must be integrated
// Paste into a Detail Wrangle (Run Over: Detail), no inputs needed.
// Written from the published mathematics, not adapted from any code.
// Constants arrive at their published values. Press the node's Create
// Spare Parameters button and every tweak becomes a slider — starting
// at 0, the published figure, and moving in the constant's own units.
// https://forma-gen.com/#plate=pendulum

float p_mass   = 1 + chf('mass_tweak');         // m₂ / m₁ · live 0.3 .. 2
float p_len    = 1 + chf('len_tweak');          // l₂ / l₁ · live 0.5 .. 1.5
float p_launch = 0.55 + chf('launch_tweak');    // launch angle (× π) · live 0.2 .. 0.95

// The plate's own colour: FORMA's CURVES accent as a cosine ramp,
// brightest near t = 0 and t = 1, near-black around t = 0.5.
vector forma_ramp(float t){
  return set(
    0.11 + 0.1196 * cos(6.28318530718 * (t + 0)),
    0.46 + 0.5 * cos(6.28318530718 * (t + 0.05)),
    0.2453 + 0.2667 * cos(6.28318530718 * (t + 0.1)));
}

// The Euler–Lagrange equations for two rods under RK4, m₁ = l₁ = 1, and the
// second bob's trail emitted as the path. In a chaotic system the launch
// angle is the parameter: low energy gives a lazy lissajous trail, high
// enough gives full flips.
int   forma_steps = 6000;    // dt = 0.01, recording every second step
float forma_dt    = 0.01;
float forma_g     = 9.81;

vector4 forma_deriv(vector4 st; float m2, l2, g){
    float a1 = st.x, a2 = st.y, w1 = st.z, w2 = st.w;
    float d = a1 - a2;
    float den = 2.0 + m2 - m2 * cos(2.0 * d);
    float dw1 = (-g * (2.0 + m2) * sin(a1) - m2 * g * sin(a1 - 2.0 * a2)
                 - 2.0 * sin(d) * m2 * (w2 * w2 * l2 + w1 * w1 * cos(d))) / den;
    float dw2 = (2.0 * sin(d) * (w1 * w1 * (1.0 + m2) + g * (1.0 + m2) * cos(a1)
                 + w2 * w2 * l2 * m2 * cos(d))) / (l2 * den);
    return set(w1, w2, dw1, dw2);
}

vector4 st = set(3.14159265359 * p_launch, 3.14159265359 * 0.95, 0.0, 0.0);
int prim = addprim(0, "polyline");
for (int i = 0; i < forma_steps; i++){
    vector4 k1 = forma_deriv(st, p_mass, p_len, forma_g);
    vector4 k2 = forma_deriv(st + k1 * forma_dt / 2.0, p_mass, p_len, forma_g);
    vector4 k3 = forma_deriv(st + k2 * forma_dt / 2.0, p_mass, p_len, forma_g);
    vector4 k4 = forma_deriv(st + k3 * forma_dt, p_mass, p_len, forma_g);
    st += (k1 + 2.0 * k2 + 2.0 * k3 + k4) * forma_dt / 6.0;
    if (!isfinite(st.x) || !isfinite(st.y)) break;
    if (i % 2 == 0) continue;
    // the tip in world terms: rods hang downward, Houdini y-up
    float tx = sin(st.x) + p_len * sin(st.y);
    float ty = -(cos(st.x) + p_len * cos(st.y));
    int pt = addpoint(0, set(tx, ty, 0.0));
    // colour sweeps the bright lobe of the ramp along the trail
    float u = float(i) / float(forma_steps);
    setpointattrib(0, "Cd", pt, forma_ramp(0.8 + 0.3 * u));
    addvertex(0, prim, pt);
}

AFTER EFFECTS · EXPRESSION

The same published mathematics as a Shape Layer path expression. Paste it onto a Path property; every constant is the published value plus a Slider Control named _tweak, so a bare paste already draws the figure and each slider moves one constant in its own units. Trim Paths is the comet.

// FORMA — PL. 60 · DOUBLE PENDULUM — Classical Lagrangian mechanics · after Euler and Lagrange
//   two rods, two bobs, one pivot
//   θ̈₁, θ̈₂ from the Euler–Lagrange equations
//   no closed form — the motion must be integrated
// After Effects port — paste onto a Shape Layer's Path property
// (Contents › Shape › Path). Written from the published mathematics, not
// adapted from any code. Constants arrive at their published values; add a
// Slider Control (Effect › Expression Controls) named <k>_tweak and that
// constant moves in its own units, starting at 0 — the published figure.
// The plate's comet and its reveal are Trim Paths; the stroke colour is
// FORMA's CURVES accent, #3DFF88. Animation runs on time.
// This plate draws 2 separate paths at its published constants:
// duplicate the group (Contents › Group) that many times and each copy draws
// its own part, read from its position in the layer. A Slider Control named
// "part" on the layer pins one instead.
// https://forma-gen.com/#plate=pendulum

// A missing slider reads 0, so a bare paste already draws the figure.
function forma_tweak(n){ try { return effect(n)("Slider"); } catch (e){ return 0; } }
var p_mass   = 1 + forma_tweak("mass_tweak");       // m₂ / m₁ · live 0.3 .. 2
var p_len    = 1 + forma_tweak("len_tweak");        // l₂ / l₁ · live 0.5 .. 1.5
var p_launch = 0.55 + forma_tweak("launch_tweak");  // launch angle (× π) · live 0.2 .. 0.95

// The frame: the plate's W × H canvas is this comp, with the origin at the
// layer's anchor; canvas y already runs down, as After Effects' does.
var forma_W = thisComp.width, forma_H = thisComp.height, forma_t = time;
var forma_phase = 0.8438958800397813;   // this plate's own fixed phase, as the page has it
function forma_pt(x, y){ return [x - forma_W / 2, y - forma_H / 2]; }
function forma_partIndex(){
  try { return Math.round(effect("part")("Slider")); } catch (e){}
  try { return thisProperty.propertyGroup(3).propertyIndex - 1; } catch (e){ return 0; }
}
var forma_part = forma_partIndex();

// The double pendulum, integrated by RK4 at dt = 0.01 from the page's own
// launch: the upper arm at π·launch (staggered by this plate's phase), the
// lower at 0.95π, both at rest. Expressions keep no state between frames, so
// every frame integrates from the launch — 260 warm-up steps, then six steps
// for every sixtieth of a second of comp time, which is the page's own pace
// at sixty paints a second. Part 0 is the trail of the lower bob, the last
// 1,500 of every second step; part 1 is the pendulum itself, pivot to bob
// to bob. The cost grows with time — at sixty seconds a frame is ~22,000
// RK4 steps, still well inside what the expression engine does in a frame.
// parts: 2
var g = 9.81, m2 = p_mass, l2 = p_len;
function forma_deriv(st){
  var a1 = st[0], a2 = st[1], w1 = st[2], w2 = st[3];
  var d = a1 - a2, den = 2 + m2 - m2 * Math.cos(2 * d);
  var dw1 = (-g * (2 + m2) * Math.sin(a1) - m2 * g * Math.sin(a1 - 2 * a2)
             - 2 * Math.sin(d) * m2 * (w2 * w2 * l2 + w1 * w1 * Math.cos(d))) / den;
  var dw2 = (2 * Math.sin(d) * (w1 * w1 * (1 + m2) + g * (1 + m2) * Math.cos(a1)
             + w2 * w2 * l2 * m2 * Math.cos(d))) / (l2 * den);
  return [w1, w2, dw1, dw2];
}
function forma_add(a, b, s){ return [a[0] + b[0] * s, a[1] + b[1] * s, a[2] + b[2] * s, a[3] + b[3] * s]; }
function forma_rk4(st, dt){
  var k1 = forma_deriv(st), k2 = forma_deriv(forma_add(st, k1, dt / 2));
  var k3 = forma_deriv(forma_add(st, k2, dt / 2)), k4 = forma_deriv(forma_add(st, k3, dt));
  var out = [];
  for (var i = 0; i < 4; i++) out.push(st[i] + (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) * dt / 6);
  return out;
}
var px = forma_W / 2, py = forma_H * 0.36;
var s = Math.min(forma_W, forma_H) * 0.42 / (1 + l2);
function forma_tip(st){
  var x1 = px + s * Math.sin(st[0]), y1 = py + s * Math.cos(st[0]);
  return [x1 + s * l2 * Math.sin(st[1]), y1 + s * l2 * Math.cos(st[1]), x1, y1];
}
var st = [Math.PI * p_launch + (forma_phase - 0.5) * 0.6, Math.PI * 0.95, 0, 0];
var steps = 260 + 6 * Math.floor(Math.max(0, forma_t) * 60);
var trail = [];
for (var i = 0; i < steps; i++){
  st = forma_rk4(st, 0.01);
  if (i % 2) trail.push(forma_tip(st));
}
if (trail.length > 1500) trail = trail.slice(trail.length - 1500);
var pts = [];
if (forma_part === 1){
  var tip = forma_tip(st);
  pts.push(forma_pt(px, py), forma_pt(tip[2], tip[3]), forma_pt(tip[0], tip[1]));
} else {
  for (var j = 0; j < trail.length; j++) pts.push(forma_pt(trail[j][0], trail[j][1]));
}
createPath(pts, [], [], false);