A harmonic bond, live

The force of a spring between two bonded particles, Hooke's law F = −k (r − r₀), positive apart: pulled together when stretched, pushed apart when compressed, nothing at the rest length. Below, a mass on such a bond is stretched and let go, stepped by the library's velocity Verlet (lusine.integrator.kick_drift, the bond's force lusine.laws.bond_force over the mass, lusine.integrator.kick), a hundred steps a second of a hundredth of a second each, on your GPU as you read. Beside it the same program computes the closed form of the oscillation, r(t) = r₀ + A cos(√(k / m) t), and the oscillation's energy, and checks both at every step.

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 bond: the mass starts 10 cm past the rest length of 1 m and swings for ten seconds, then is let go again. The chart draws the stepped length and the closed form against the mass's time, a point a step over the last thousand steps, the two one over the other. Below it, the time, the stepped length, the closed form and the force.
• The checks: the chart draws each check's miss at every step as a share of what that check allows, on a log scale, so a check holds while its line stays below 1. agrees holds while the stepped length has stayed within 0.2% of the stretch of the closed form at every step; keeps its energy while the energy, the mass's kinetic energy ½ m v² and the bond's ½ k (r − r₀)², has stayed within 0.05% of what it started with; restoring while the bond gives no force at all at its rest length and, at every step, a force pointing back toward it.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The stiffness (20 N/m), the rest length (1 m), the mass (2 kg), the stretch (10 cm) and the step (0.01 s) are the program's settings, which the host may write; this page offers no control to change them yet.

The program

The body, the closed form and the checks, the panels the page draws, and the two procedures: start, run once before the first step, stretches the bond and takes the acceleration there; main, each step, is one velocity Verlet step and the checks after it.

// The harmonic bond's live page: a mass on a bond, stretched and let go, stepped by the library's velocity Verlet under
// the bond's force, beside the closed form of the same oscillation, with its energy kept and the page's checks
// computed here on the GPU and shown by the panels the program declares.

namespace bond_force;

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

// The reader's settings: the bond's stiffness and rest length, the mass on it, how far it starts stretched, the step.
host struct {
  float<N/m> stiffness = 20.0 N/m; float<m> rest = 1.0 m; float<kg> mass = 2.0 kg; float<m> stretch = 0.1 m;
  float<s> dt = 0.01 s;
} settings;

// The mass, its length along the bond stepped from rest at the stretch; its time the sum of the steps it has taken.
struct { float<m> r; float<m/s> velocity; float<m/s2> acceleration; float<N> force; float<s> time; } body;

// The closed form at the body's time, r(t) = r0 + A cos(sqrt(k / m) t), the oscillation's energy, kinetic and
// stored, against what it started with, and the checks: the stepped length within 0.2% of the stretch of the closed
// form at every step, the energy within 0.05% of its start at every step (each miss also as a share of what its check
// allows, which the page charts), and the force restoring: none at the rest length, and at every step pointing back
// toward it.
struct {
  float<m> closed_form; float<m> error; float error_share; float<m> worst; float<J> energy; float<J> drift;
  float drift_share; float<J> worst_drift;
  bool agrees = true; bool keeps_energy = true; bool restoring = 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);
}

// Newton's second law: the force over the mass.
float<m/s2> newton(float<N> force, float<kg> mass) { return force / mass; }

float<m> swing(float<m> rest, float<m> stretch, float<N/m> stiffness, float<kg> mass, float<s> time)
{
  return rest + stretch * cos(sqrt(stiffness / mass) * time);
}

// The energy of the oscillation: the mass's kinetic energy and the bond's, half the stiffness times the stretch
// squared.
float<J> energy(float<kg> mass, float<m/s> velocity, float<N/m> stiffness, float<m> r, float<m> rest)
{
  let x = r - rest;
  return 0.5 * mass * velocity * velocity + 0.5 * stiffness * x * x;
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the bond", chart = "length", quantity = body.r, steps = 1000, against = body.time,
                    label = "stepped" } },
  lusine.series { { panel = "the bond", chart = "length", quantity = checks.closed_form, steps = 1000,
                    against = body.time, label = "closed form" } },
  lusine.value { { panel = "the bond", quantity = body.time, label = "time" } },
  lusine.value { { panel = "the bond", quantity = body.r, label = "stepped length" } },
  lusine.value { { panel = "the bond", quantity = checks.closed_form, label = "closed form" } },
  lusine.value { { panel = "the bond", quantity = body.force, label = "force" } },
  lusine.series { { panel = "the checks", chart = "share of the tolerance", quantity = checks.error_share,
                    steps = 1000, log = true, label = "length error" } },
  lusine.series { { panel = "the checks", chart = "share of the tolerance", quantity = checks.drift_share,
                    steps = 1000, log = true, label = "energy drift" } },
  lusine.value { { panel = "the checks", quantity = checks.worst, label = "worst length error" } },
  lusine.value { { panel = "the checks", quantity = checks.worst_drift, label = "worst energy drift" } },
  lusine.value { { panel = "the checks", quantity = checks.agrees, label = "agrees" } },
  lusine.value { { panel = "the checks", quantity = checks.keeps_energy, label = "keeps its energy" } },
  lusine.value { { panel = "the checks", quantity = checks.restoring, label = "restoring" } },
};

void start()
{
  body.r = settings.rest + settings.stretch;
  body.force = lusine.laws.bond_force(body.r, settings.stiffness, settings.rest);
  body.acceleration = newton(body.force, settings.mass);
}

void main()
{
  // Every ten seconds the mass is let go again from the stretch, so the chart's last thousand steps hold the run since.
  body.r = body.time >= 10.0 s ? settings.rest + settings.stretch : body.r;
  body.velocity = body.time >= 10.0 s ? 0.0 m/s : body.velocity;
  body.time = body.time >= 10.0 s ? 0.0 s : body.time;
  body.force = lusine.laws.bond_force(body.r, settings.stiffness, settings.rest);
  body.acceleration = newton(body.force, settings.mass);
  (body.r, body.velocity) = lusine.integrator.kick_drift(position = body.r, velocity = body.velocity,
                                                         acceleration = body.acceleration, dt = settings.dt);
  body.force = lusine.laws.bond_force(body.r, settings.stiffness, settings.rest);
  body.acceleration = newton(body.force, settings.mass);
  body.velocity = lusine.integrator.kick(body.velocity, body.acceleration, settings.dt);
  body.time = body.time + settings.dt;
  checks.closed_form = swing(settings.rest, settings.stretch, settings.stiffness, settings.mass, body.time);
  checks.error = abs(body.r - checks.closed_form);
  checks.worst = checks.error > checks.worst ? checks.error : checks.worst;
  checks.error_share = checks.error / (0.002 * settings.stretch);
  checks.agrees = checks.agrees && checks.error_share <= 1.0;
  checks.energy = energy(settings.mass, body.velocity, settings.stiffness, body.r, settings.rest);
  checks.drift = abs(checks.energy - 0.5 * settings.stiffness * settings.stretch * settings.stretch);
  checks.worst_drift = checks.drift > checks.worst_drift ? checks.drift : checks.worst_drift;
  checks.drift_share = checks.drift / (0.0005 * 0.5 * settings.stiffness * settings.stretch * settings.stretch);
  checks.keeps_energy = checks.keeps_energy && checks.drift_share <= 1.0;
  checks.restoring = checks.restoring
                     && lusine.laws.bond_force(r = settings.rest, settings.stiffness, settings.rest) == 0.0 N
                     && body.force * (body.r - settings.rest) <= 0.0 N*m;
  lamps.colour = lamp_colour(holds = checks.agrees && checks.keeps_energy && checks.restoring);
}

The law's other checks (the force against −k (r − r₀) along a sweep of lengths) are on its proof page.