Skip to content

The pendulum constraint, live

A point held at a fixed length from a pivot, position-based: each step it is projected back onto the sphere of that radius about the pivot, p′ = c + (p − c) L / |p − c|, the correction of a distance constraint of infinite stiffness. Below, three bobs of 0.25, 0.5 and 1 m are let go 0.05 rad out from below their pivots and swing under gravity: each step the program moves them by position Verlet, its own step, and the library's constraint (lusine.integrator.pendulum_constraint) puts each back at its length, two milliseconds a step and five hundred steps a second, on your GPU as you read. Nothing else holds them: the swing is the constraint's. The program times each bob's crossings below its pivot and measures its period, which for a small swing is 2π √(L / g).

This page runs the program on your GPU through WebGPU, so it needs a browser with WebGPU (a recent Chrome, Edge, Safari or Firefox). Without it the program does not run and nothing is drawn below.

What it shows

  • The strip above the controls is green while every check holds, and turns red for good once one fails.
  • The swing: the chart draws each bob's sideways position against the time over the last three seconds: the shortest swings twice for each swing of the longest. Below it, the time, and the shortest and the longest bob's measured periods beside 2π √(L / g), 1.003 s and 2.006 s.
  • The checks: the period, every bob's measured period within 0.5% of 2π √(L / g) once it has crossed twice (its swing of 0.05 rad lengthens it by only 0.016%); held at its length, every bob within 4 roundings of its length at every step; along its own direction, a point and a pivot that change each step, the point projected to its length from the pivot and along its own direction from it, to a few roundings; twice is once, the constraint applied to a point it has already placed leaving it where it is, within 2 roundings. Each is held at every step.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. Gravity, the bobs and the step are the program's settings and table, which the host may write; this page offers no control to change them yet.

The program

The bobs, the swing and the checks, the panels the page draws, and the two procedures: start, run once before the first step, lets the bobs go at rest; main, each step, moves them by position Verlet, puts them back at their lengths, times their crossings, and takes the checks.

// The pendulum constraint's live page: three bobs of three lengths swinging under gravity, stepped by position Verlet
// (the program's own) and held to their lengths each step by the library's constraint, their periods measured from
// their crossings against 2 pi sqrt(L / g), with the page's checks computed here on the GPU and shown by the panels
// the program declares.

namespace pendulum_constraint;

lusine.version version = { { number = 1 } };

// The reader's settings: gravity and the step.
host struct { float<m/s2> gravity = 9.80665 m/s2; float<s> dt = 0.002 s; } settings;

// The bobs: each one's number, its length, where it is, where it was a step before, where gravity alone would take it,
// how many times it has crossed below its pivot going one way, when it first and last did, and its period measured
// from them, (last - first) / (crossings - 1), with the small-swing period 2 pi sqrt(L / g); and whether its period
// and its length hold this step.
struct bob {
  integer number; float<m> length; vector<float<m>, 3> position; vector<float<m>, 3> before; vector<float<m>, 3> free;
  integer crossings; float<s> first; float<s> last; float<s> period; float<s> small_swing; bool timed; bool held;
};
table<bob> bobs(capacity = 3) = {
  { number = 0, length = 0.25 m }, { number = 1, length = 0.5 m }, { number = 2, length = 1.0 m },
};

// The swing: its time; each bob's sideways position, its measured period and the small-swing one; and a point and a
// pivot that change each step, the point projected, and projected again.
struct {
  float<s> time; float<m> short_x; float<m> middle_x; float<m> long_x; float<s> long_period; float<s> long_small;
  float<s> short_period; float<s> short_small; float turn; vector<float<m>, 3> point; vector<float<m>, 3> pivot;
  float<m> length; vector<float<m>, 3> held; vector<float<m>, 3> again;
} swing;

// The checks: every bob's measured period within 0.5% of 2 pi sqrt(L / g) once it has crossed twice (its swing of
// 0.05 rad lengthens it by 0.016%); every bob at its length within 4 roundings at every step; the changing point held
// at its length from its pivot along its own direction from it; and the constraint applied twice the same as once,
// within 2 roundings.
struct { bool period = true; bool held = true; bool direction = true; bool idempotent = true; } checks;

// A strip of lamps the page draws, lit green while every check holds and red once one fails.
struct lamp { vector<float, 4> colour; };
lattice<lamp, 2> lamps(spacing = 1 m, boundary = wall);

vector<float, 4> lamp_colour(bool holds)
{
  return holds ? vector<float, 4>(0.15, 0.6, 0.3, 1.0) : vector<float, 4>(0.75, 0.15, 0.15, 1.0);
}

vector<float<m>, 3> pivot() { return vector<float<m>, 3>(0.0 m, 0.0 m, 0.0 m); }
vector<float<m/s2>, 3> down(float<m/s2> gravity) { return vector<float<m/s2>, 3>(0.0 m/s2, -gravity, 0.0 m/s2); }

// A bob at rest, 0.05 rad out from below its pivot: (L sin 0.05, -L cos 0.05, 0).
vector<float<m>, 3> out(float<m> length)
{
  return vector<float<m>, 3>(0.04997917 * length, -0.99875026 * length, 0.0 m);
}

// Position Verlet under gravity, the program's own step: x' = x + (x - x_before) + g dt^2.
vector<float<m>, 3> verlet(vector<float<m>, 3> position, vector<float<m>, 3> before, float<m/s2> gravity, float<s> dt)
{
  return position + (position - before) + down(gravity) * (dt * dt);
}

// A crossing below the pivot from +x to -x between the last step and this, and when, between the two steps.
bool crossed(vector<float<m>, 3> before, vector<float<m>, 3> position)
{
  return before.x > 0.0 m && position.x <= 0.0 m;
}
float<s> when(vector<float<m>, 3> before, vector<float<m>, 3> position, float<s> time, float<s> dt)
{
  return time - dt + dt * (before.x / (before.x - position.x));
}

float<s> small_swing(float<m> length, float<m/s2> gravity) { return 6.2831853 * sqrt(length / gravity); }

bool period_holds(integer crossings, float<s> period, float<s> small_swing)
{
  return crossings < 2 || abs(period / small_swing - 1.0) < 0.005;
}
bool at_length(vector<float<m>, 3> position, float<m> length)
{
  return abs(length(position) - length) <= 0.00000048 * length;
}
bool all_held(table<bool> held) { return all(held); }

float<m> x_of(table<vector<float<m>, 3>> position, table<integer> number, integer which)
{
  return sum(position.x) where number == which;
}
float<s> period_of(table<float<s>> period, table<integer> number, integer which)
{
  return sum(period) where number == which;
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the swing", chart = "sideways", quantity = swing.short_x, steps = 1500,
                    against = swing.time, label = "a quarter metre" } },
  lusine.series { { panel = "the swing", chart = "sideways", quantity = swing.middle_x, steps = 1500,
                    against = swing.time, label = "half a metre" } },
  lusine.series { { panel = "the swing", chart = "sideways", quantity = swing.long_x, steps = 1500,
                    against = swing.time, label = "a metre" } },
  lusine.value { { panel = "the swing", quantity = swing.time, label = "time" } },
  lusine.value { { panel = "the swing", quantity = swing.short_period, label = "0.25 m's period" } },
  lusine.value { { panel = "the swing", quantity = swing.short_small, label = "2 pi sqrt(L / g)" } },
  lusine.value { { panel = "the swing", quantity = swing.long_period, label = "1 m's period" } },
  lusine.value { { panel = "the swing", quantity = swing.long_small, label = "2 pi sqrt(L / g)" } },
  lusine.value { { panel = "the checks", quantity = checks.period, label = "the period" } },
  lusine.value { { panel = "the checks", quantity = checks.held, label = "held at its length" } },
  lusine.value { { panel = "the checks", quantity = checks.direction, label = "along its own direction" } },
  lusine.value { { panel = "the checks", quantity = checks.idempotent, label = "twice is once" } },
};

void start()
{
  bobs.position = out(bobs.length);
  bobs.before = bobs.position;
}

void main()
{
  swing.time = swing.time + settings.dt;
  bobs.free = verlet(bobs.position, bobs.before, settings.gravity, settings.dt);
  bobs.before = bobs.position;
  bobs.position = lusine.integrator.pendulum_constraint(position = bobs.free, pivot = pivot(), bobs.length);
  bobs.first = crossed(bobs.before, bobs.position) && bobs.crossings == 0
               ? when(bobs.before, bobs.position, swing.time, settings.dt) : bobs.first;
  bobs.last = crossed(bobs.before, bobs.position) ? when(bobs.before, bobs.position, swing.time, settings.dt)
                                                  : bobs.last;
  bobs.crossings = crossed(bobs.before, bobs.position) ? bobs.crossings + 1 : bobs.crossings;
  bobs.period = bobs.crossings >= 2 ? (bobs.last - bobs.first) / float(bobs.crossings - 1) : bobs.period;
  bobs.small_swing = small_swing(bobs.length, settings.gravity);
  swing.short_x = x_of(bobs.position, bobs.number, which = 0);
  swing.middle_x = x_of(bobs.position, bobs.number, which = 1);
  swing.long_x = x_of(bobs.position, bobs.number, which = 2);
  swing.short_period = period_of(bobs.period, bobs.number, which = 0);
  swing.long_period = period_of(bobs.period, bobs.number, which = 2);
  swing.short_small = period_of(period = bobs.small_swing, bobs.number, which = 0);
  swing.long_small = period_of(period = bobs.small_swing, bobs.number, which = 2);
  bobs.timed = period_holds(bobs.crossings, bobs.period, bobs.small_swing);
  bobs.held = at_length(bobs.position, bobs.length);
  checks.period = checks.period && all_held(held = bobs.timed);
  checks.held = checks.held && all_held(bobs.held);
  // A point and a pivot that change each step, the point projected onto its sphere, and projected again.
  swing.turn = swing.turn + 0.41;
  swing.point = vector<float<m>, 3>(3.0 m * cos(swing.turn), 2.0 m * sin(1.7 * swing.turn),
                                    1.5 m * cos(0.6 * swing.turn));
  swing.pivot = vector<float<m>, 3>(0.5 m * sin(0.3 * swing.turn), 0.2 m, -0.4 m * cos(swing.turn));
  swing.length = 1.6 m + 1.4 m * sin(0.8 * swing.turn);
  swing.held = lusine.integrator.pendulum_constraint(position = swing.point, pivot = swing.pivot, swing.length);
  swing.again = lusine.integrator.pendulum_constraint(position = swing.held, pivot = swing.pivot, swing.length);
  checks.direction = checks.direction
                     && abs(length(swing.held - swing.pivot) - swing.length)
                        <= 0.00000096 * (swing.length + length(swing.pivot))
                     && dot(swing.point - swing.pivot, swing.held - swing.pivot)
                        >= (1.0 - 0.0000019) * length(swing.point - swing.pivot) * length(swing.held - swing.pivot);
  checks.idempotent = checks.idempotent
                      && length(swing.again - swing.held) <= 0.00000024 * (swing.length + length(swing.pivot));
  lamps.colour = lamp_colour(holds = checks.period && checks.held && checks.direction && checks.idempotent);
}

The law's other proof (five lengths' periods, and two hundred points at random) is on its proof page.