Velocity Verlet, live¶
Velocity Verlet, a whole step of it: the kick-drift (lusine.integrator.kick_drift), a half-kick from the acceleration at the position and the drift with that half-step velocity, then the force at the new position, then the closing half-kick (lusine.integrator.kick). It is second order, it takes one force a step, it keeps a conservative system's energy in a band with no drift, and it is time-reversible: a run turned round comes back to where it began. Below, a body is stepped round a circular Kepler orbit (GM = 1 m³/s², radius 1 m), two hundred steps a period and two hundred steps a second, on your GPU as you read, and the program adds its energy up at every step. Beside it another orbit is run out for 450 steps, turned round and run back; four springs are stepped at four steps for the order; a body falls under a constant force; and a projectile flies in space, its position a vector of three.
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 orbit: the first chart draws the body's path, y against x, a point a step over the last period: a circle. The second draws its energy's miss of the circle's −1/2 J/kg, in millionths of a J/kg, at each of the last thousand steps on a log scale, a rounding's wobble with no drift. Below them, the energy and its worst miss so far.
- The checks: how near its start the turned-back orbit came, and the ratio of the springs' misses of cos(ω t) at 0.3 s at steps of 1/200 and 1/400 s; energy kept, the orbit's energy within 10⁻⁵ J/kg of −1/2 at every step; the way back, the turned-back orbit within 2000 roundings of its start, its velocity the start's turned round; second order, halving the step quartering the miss, each ratio within 0.1 of 4 from 1/50 s to 1/400 s; a fall is exact, under a constant force the two half-kicks landing on −g t² / 2 within 10⁻⁵ m at 1 s; a projectile in space, launched at (3, 1, 4) m/s and landing after a second where v t + g t² / 2 says within 10⁻⁵ m, each component exactly the scalar step's of that component alone at every step.
Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The orbit, the springs, gravity and the steps are the program's settings, which the host may write; this page offers no control to change them yet.
The program¶
The orbits, the springs, the fall, the projectile and the checks, the panels the page draws, and the two procedures: start, run once before the first step, takes the orbit's acceleration at its start; main, each step, is one velocity Verlet step of every body and the checks after it.
// The kick's live page: whole velocity Verlet steps, the library's kick-drift and kick, round a circular Kepler orbit
// whose energy stays in its band, an orbit run out and turned back to its start, four springs at four steps for the
// second order, a fall that is exact, and a projectile in space, each component the scalar step's, with the page's
// checks computed here on the GPU and shown by the panels the program declares.
namespace kick;
lusine.version version = { { number = 1 } };
// The reader's settings: the orbit's GM and its step (two hundred a period of 2 pi s); the run's length in steps, and
// the steps out before the way back; the springs' angular frequency and where they are at 0.3 s, cos(0.6 pi) of their
// start (given, since a cosine taken on the device carries its own error); gravity; and the projectile's step.
host struct {
float<m3/s2> gm = 1.0 m3/s2; float<s> dt = 0.031415927 s; integer steps = 1000; integer out = 450;
float<s^-1> omega = 6.2831853 s^-1; float<m> at_03 = -0.30901699 m;
float<m/s2> gravity = 9.81 m/s2; float<s> projectile_dt = 0.015625 s;
} settings;
// The run: its step, from 0 to 999 and round again.
struct { integer step = -1; } run;
// The orbit, a body in a plane from (0, 1) m at (-1, 0) m/s, never restarted: its position, velocity and acceleration,
// its energy per unit mass, v^2 / 2 - GM / r, its miss of the circle's -1/2 (also in millionths of a J/kg, for the
// chart) and the worst so far.
struct {
vector<float<m>, 2> position; vector<float<m/s>, 2> velocity; vector<float<m/s2>, 2> acceleration;
float<J/kg> energy; float<J/kg> miss; float millionths; float<J/kg> worst; float<m> x; float<m> y;
} orbit;
// The way back: another orbit from the same start, run out for 450 steps, its velocity turned round, and run back;
// how far from its start it came back.
struct {
vector<float<m>, 2> position; vector<float<m/s>, 2> velocity; vector<float<m/s2>, 2> acceleration; float<m> back;
} returning;
// The order's four springs at steps of 1/50 to 1/400 s, each one's position, velocity and acceleration, its miss of
// cos(omega t) at 0.3 s once it reaches it, and the misses' ratio of the last two; the fall from rest; the projectile
// in space and its components each stepped alone.
struct {
float<m> x50; float<m/s> v50; float<m/s2> a50; float<m> x100; float<m/s> v100; float<m/s2> a100; float<m> x200;
float<m/s> v200; float<m/s2> a200; float<m> x400; float<m/s> v400; float<m/s2> a400; float<m> e50; float<m> e100;
float<m> e200; float<m> e400; float ratio; float<m> fx; float<m/s> fv;
vector<float<m>, 3> px; vector<float<m/s>, 3> pv; float<m> p0; float<m/s> q0; float<m> p1; float<m/s> q1; float<m> p2;
float<m/s> q2;
} others;
// The checks: the orbit's energy within 10^-5 J/kg of -1/2 at every step, with no drift; the way back within 2000
// roundings of the start; halving the step quartering the miss at 0.3 s, each ratio within 0.1 of 4; a fall landing on
// -g t^2 / 2 within 10^-5 m at 1 s; and the projectile landing where v t + g t^2 / 2 says within 10^-5 m, each
// component exactly the scalar step's at every step.
struct { bool energy = true; bool reversible = true; bool order = true; bool fall = true; bool space = 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);
}
// The program's own force kernels: Kepler's pull toward the origin, -GM x / r^3, and the spring's, -omega^2 x.
vector<float<m/s2>, 2> kepler(vector<float<m>, 2> position, float<m3/s2> gm)
{
let r = length(position);
return -gm / (r * r * r) * position;
}
float<m/s2> pull(float<m> position, float<s^-1> omega) { return -omega * omega * position; }
vector<float<m>, 2> top() { return vector<float<m>, 2>(0.0 m, 1.0 m); }
vector<float<m/s>, 2> leftward() { return vector<float<m/s>, 2>(-1.0 m/s, 0.0 m/s); }
vector<float<m/s>, 3> launch() { return vector<float<m/s>, 3>(3.0 m/s, 1.0 m/s, 4.0 m/s); }
vector<float<m/s2>, 3> down(float<m/s2> gravity) { return vector<float<m/s2>, 3>(0.0 m/s2, 0.0 m/s2, -gravity); }
float<J/kg> energy(vector<float<m>, 2> position, vector<float<m/s>, 2> velocity, float<m3/s2> gm)
{
return 0.5 * dot(velocity, velocity) - gm / length(position);
}
float<m> miss_at(float<m> miss, float<m> x, integer step, integer reaches, float<m> at_03)
{
return step == reaches - 1 ? abs(x - at_03) : miss;
}
output table<lusine.widget> panels = {
lusine.series { { panel = "the orbit", chart = "the orbit", quantity = orbit.y, steps = 200, against = orbit.x,
label = "y against x" } },
lusine.series { { panel = "the orbit", chart = "energy's miss", quantity = orbit.millionths, steps = 1000,
log = true, label = "|E + 1/2|, millionths of a J/kg" } },
lusine.value { { panel = "the orbit", quantity = orbit.energy, label = "energy" } },
lusine.value { { panel = "the orbit", quantity = orbit.worst, label = "worst miss" } },
lusine.value { { panel = "the checks", quantity = returning.back, label = "came back within" } },
lusine.value { { panel = "the checks", quantity = others.ratio, label = "1/200 s over 1/400 s" } },
lusine.value { { panel = "the checks", quantity = checks.energy, label = "energy kept" } },
lusine.value { { panel = "the checks", quantity = checks.reversible, label = "the way back" } },
lusine.value { { panel = "the checks", quantity = checks.order, label = "second order" } },
lusine.value { { panel = "the checks", quantity = checks.fall, label = "a fall is exact" } },
lusine.value { { panel = "the checks", quantity = checks.space, label = "a projectile in space" } },
};
void start()
{
orbit.position = top();
orbit.velocity = leftward();
orbit.acceleration = kepler(position = orbit.position, settings.gm);
}
void main()
{
// Every run's steps everything but the orbit starts again.
run.step = run.step + 1 >= settings.steps ? 0 : run.step + 1;
returning.position = run.step == 0 ? top() : returning.position;
returning.velocity = run.step == 0 ? leftward() : returning.velocity;
returning.acceleration = run.step == 0 ? kepler(position = top(), settings.gm) : returning.acceleration;
returning.velocity = run.step == settings.out ? -returning.velocity : returning.velocity;
others.x50 = run.step == 0 ? 1.0 m : others.x50;
others.v50 = run.step == 0 ? 0.0 m/s : others.v50;
others.a50 = run.step == 0 ? pull(position = 1.0 m, settings.omega) : others.a50;
others.x100 = run.step == 0 ? 1.0 m : others.x100;
others.v100 = run.step == 0 ? 0.0 m/s : others.v100;
others.a100 = run.step == 0 ? pull(position = 1.0 m, settings.omega) : others.a100;
others.x200 = run.step == 0 ? 1.0 m : others.x200;
others.v200 = run.step == 0 ? 0.0 m/s : others.v200;
others.a200 = run.step == 0 ? pull(position = 1.0 m, settings.omega) : others.a200;
others.x400 = run.step == 0 ? 1.0 m : others.x400;
others.v400 = run.step == 0 ? 0.0 m/s : others.v400;
others.a400 = run.step == 0 ? pull(position = 1.0 m, settings.omega) : others.a400;
others.fx = run.step == 0 ? 0.0 m : others.fx;
others.fv = run.step == 0 ? 0.0 m/s : others.fv;
others.px = run.step == 0 ? vector<float<m>, 3>(0.0 m, 0.0 m, 0.0 m) : others.px;
others.pv = run.step == 0 ? launch() : others.pv;
others.p0 = run.step == 0 ? 0.0 m : others.p0;
others.q0 = run.step == 0 ? 3.0 m/s : others.q0;
others.p1 = run.step == 0 ? 0.0 m : others.p1;
others.q1 = run.step == 0 ? 1.0 m/s : others.q1;
others.p2 = run.step == 0 ? 0.0 m : others.p2;
others.q2 = run.step == 0 ? 4.0 m/s : others.q2;
// One velocity Verlet step of each: the kick-drift, the force at the new position, the kick.
(orbit.position, orbit.velocity) = lusine.integrator.kick_drift(orbit.position, orbit.velocity, orbit.acceleration,
settings.dt);
orbit.acceleration = kepler(position = orbit.position, settings.gm);
orbit.velocity = lusine.integrator.kick(orbit.velocity, orbit.acceleration, settings.dt);
(returning.position, returning.velocity) = lusine.integrator.kick_drift(returning.position, returning.velocity,
returning.acceleration, settings.dt);
returning.acceleration = kepler(position = returning.position, settings.gm);
returning.velocity = lusine.integrator.kick(returning.velocity, returning.acceleration, settings.dt);
(others.x50, others.v50) = lusine.integrator.kick_drift(position = others.x50, velocity = others.v50,
acceleration = others.a50, dt = 0.02 s);
others.a50 = pull(position = others.x50, settings.omega);
others.v50 = lusine.integrator.kick(velocity = others.v50, acceleration = others.a50, dt = 0.02 s);
(others.x100, others.v100) = lusine.integrator.kick_drift(position = others.x100, velocity = others.v100,
acceleration = others.a100, dt = 0.01 s);
others.a100 = pull(position = others.x100, settings.omega);
others.v100 = lusine.integrator.kick(velocity = others.v100, acceleration = others.a100, dt = 0.01 s);
(others.x200, others.v200) = lusine.integrator.kick_drift(position = others.x200, velocity = others.v200,
acceleration = others.a200, dt = 0.005 s);
others.a200 = pull(position = others.x200, settings.omega);
others.v200 = lusine.integrator.kick(velocity = others.v200, acceleration = others.a200, dt = 0.005 s);
(others.x400, others.v400) = lusine.integrator.kick_drift(position = others.x400, velocity = others.v400,
acceleration = others.a400, dt = 0.0025 s);
others.a400 = pull(position = others.x400, settings.omega);
others.v400 = lusine.integrator.kick(velocity = others.v400, acceleration = others.a400, dt = 0.0025 s);
(others.fx, others.fv) = lusine.integrator.kick_drift(position = others.fx, velocity = others.fv,
acceleration = -settings.gravity, dt = 0.01 s);
others.fv = lusine.integrator.kick(velocity = others.fv, acceleration = -settings.gravity, dt = 0.01 s);
(others.px, others.pv) = lusine.integrator.kick_drift(position = others.px, velocity = others.pv,
acceleration = down(settings.gravity),
dt = settings.projectile_dt);
others.pv = lusine.integrator.kick(velocity = others.pv, acceleration = down(settings.gravity),
dt = settings.projectile_dt);
(others.p0, others.q0) = lusine.integrator.kick_drift(position = others.p0, velocity = others.q0,
acceleration = 0.0 m/s2, dt = settings.projectile_dt);
others.q0 = lusine.integrator.kick(velocity = others.q0, acceleration = 0.0 m/s2, dt = settings.projectile_dt);
(others.p1, others.q1) = lusine.integrator.kick_drift(position = others.p1, velocity = others.q1,
acceleration = 0.0 m/s2, dt = settings.projectile_dt);
others.q1 = lusine.integrator.kick(velocity = others.q1, acceleration = 0.0 m/s2, dt = settings.projectile_dt);
(others.p2, others.q2) = lusine.integrator.kick_drift(position = others.p2, velocity = others.q2,
acceleration = -settings.gravity, dt = settings.projectile_dt);
others.q2 = lusine.integrator.kick(velocity = others.q2, acceleration = -settings.gravity,
dt = settings.projectile_dt);
orbit.x = orbit.position.x;
orbit.y = orbit.position.y;
orbit.energy = energy(orbit.position, orbit.velocity, settings.gm);
orbit.miss = abs(orbit.energy + 0.5 J/kg);
orbit.millionths = orbit.miss / 0.000001 J/kg;
orbit.worst = orbit.miss > orbit.worst ? orbit.miss : orbit.worst;
returning.back = run.step + 1 == 2 * settings.out ? length(returning.position - top()) : returning.back;
others.e50 = miss_at(miss = others.e50, x = others.x50, run.step, reaches = 15, settings.at_03);
others.e100 = miss_at(miss = others.e100, x = others.x100, run.step, reaches = 30, settings.at_03);
others.e200 = miss_at(miss = others.e200, x = others.x200, run.step, reaches = 60, settings.at_03);
others.e400 = miss_at(miss = others.e400, x = others.x400, run.step, reaches = 120, settings.at_03);
others.ratio = others.e200 / others.e400;
checks.energy = checks.energy && orbit.miss < 0.00001 J/kg;
checks.reversible = checks.reversible
&& (run.step + 1 != 2 * settings.out
|| returning.back <= 0.00024 m
&& length(returning.velocity + leftward()) <= 0.00024 m/s);
checks.order = checks.order
&& (run.step < 119 || abs(others.e50 / others.e100 - 4.0) < 0.1
&& abs(others.e100 / others.e200 - 4.0) < 0.1
&& abs(others.e200 / others.e400 - 4.0) < 0.1);
checks.fall = checks.fall && (run.step != 99 || abs(others.fx + 0.5 * settings.gravity * 1.0 s2) < 0.00001 m);
checks.space = checks.space && others.px.x == others.p0 && others.pv.x == others.q0 && others.px.y == others.p1
&& others.pv.y == others.q1 && others.px.z == others.p2 && others.pv.z == others.q2
&& (run.step != 63 || abs(others.px.x - 3.0 m) < 0.00001 m && abs(others.px.y - 1.0 m) < 0.00001 m
&& abs(others.px.z - (4.0 m - 0.5 * settings.gravity * 1.0 s2)) < 0.00001 m);
lamps.colour = lamp_colour(holds = checks.energy && checks.reversible && checks.order && checks.fall
&& checks.space);
}
The law's other proof (a hundred periods of the orbit, and the order at four steps) is on its proof page.