The kick-drift, live¶
Velocity Verlet's first half: a half-kick from the acceleration at the position, v' = v + a dt / 2, then the drift with that half-step velocity, x' = x + v' dt. The method's second half, the closing kick with the acceleration at the new position, is lusine.integrator.kick; whole steps of the two, the method's order and an orbit are on the kick's page. Below, the library's kick-drift (lusine.integrator.kick_drift) takes a body from rest under gravity at a step length swept from none to half a second, 3.125 ms a step, twenty steps a second, on your GPU as you read: from rest under a constant force it must land on the parabola −g dt² / 2. Beside it the program steps a state that turns a little each step, its mirror, and the same state in space, each component of it stepped alone.
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.
- From rest: the chart draws how far the body falls in one step against the step's length, on a log scale, a point a step over one whole sweep, and the parabola g dt² / 2 beside it: the two lie one over the other. Below it, the step, the fall and the half-step velocity, −g dt / 2.
- The checks: a turning state stepped, and its mirror stepped, the same number with its sign turned; a parabola, the fall within 0.00005% of g dt² / 2 and the half-step velocity −g dt / 2 to a rounding; no force, a drift, a moving body with no acceleration going to x + v dt to a rounding with its velocity unchanged; the mirror, the mirrored state under the mirrored force stepping to the mirrored state, bit for bit, so nothing in the step prefers a direction; in three dimensions, a step of a body in space giving in each component exactly the scalar step of that component alone. 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 turning state and the sweep are the program's settings, which the host may write; this page offers no control to change them yet.
The program¶
The sweep, the steps and the checks, the panels the page draws, and the one procedure, main, run each step: the sweep moves on, the kick-drift is taken from each state, and the checks are taken.
// The kick-drift's live page: velocity Verlet's first half taken from rest under gravity at a step swept from none to
// half a second, beside the parabola it must land on, and from states that turn each step, held to the drift with no
// force, the mirror and the three components of a step in space, with the page's checks computed here on the GPU and
// shown by the panels the program declares.
namespace kick_drift;
lusine.version version = { { number = 1 } };
// The reader's settings: gravity; a body's position, velocity and acceleration, which the program turns each step;
// and how far each step moves the step length, from none.
host struct {
float<m/s2> gravity = 9.81 m/s2; float<m> position = 0.3 m; float<m/s> velocity = 1.1 m/s;
float<m/s2> acceleration = -0.5 m/s2; float<s> stride = 0.003125 s;
} settings;
// The sweep: the step it is at, from 0 to 160 and round again, the step length there, and an angle that turns the
// body's state.
struct { integer step = -1; float<s> dt; float turn; } sweep;
// From rest under gravity: where the step leaves the body and its half-step velocity, the fall's depth and the
// parabola's, g dt^2 / 2; with no force, where it leaves a moving body and its half-step velocity; a state and its
// mirror, each stepped; and a body in space, its components and each stepped alone.
struct {
float<m> dropped; float<m/s> half; float<m> fall; float<m> parabola; float<m> drifted; float<m/s> coasting;
float<m> x; float<m/s> v; float<m> mirror_x; float<m/s> mirror_v;
vector<float<m>, 3> space_x; vector<float<m/s>, 3> space_v; float<m> x0; float<m/s> v0; float<m> x1; float<m/s> v1;
float<m> x2; float<m/s> v2;
} step;
// The checks: from rest under gravity the step landing within 0.00005% of -g dt^2 / 2, its half-step velocity
// -g dt / 2 to a rounding; with no force a drift, x + v dt to a rounding and the velocity unchanged; the mirrored state
// under the mirrored force stepping to the mirrored state exactly; and in space each component exactly the scalar
// step of that component, each at every step.
struct { bool parabola = true; bool drift = true; bool mirror = true; bool components = 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 body's state turned by the angle, so each step steps another one; its components in space each turned again.
float turned(float turn) { return 0.5 + 0.5 * cos(turn); }
vector<float<m>, 3> in_space(float<m> position, float turn)
{
return vector<float<m>, 3>(position * turned(turn), -position * turned(turn = 2.0 * turn), position);
}
vector<float<m/s>, 3> moving(float<m/s> velocity, float turn)
{
return vector<float<m/s>, 3>(velocity, velocity * turned(turn), -velocity * turned(turn = 3.0 * turn));
}
vector<float<m/s2>, 3> pulled(float<m/s2> acceleration, float turn)
{
return vector<float<m/s2>, 3>(acceleration * turned(turn = 5.0 * turn), acceleration, -acceleration);
}
bool close(float<m> got, float<m> want, float share) { return abs(got - want) <= share * abs(want); }
output table<lusine.widget> panels = {
lusine.series { { panel = "from rest", chart = "fall in one step", quantity = step.fall, steps = 161,
against = sweep.dt, log = true, label = "kick_drift" } },
lusine.series { { panel = "from rest", chart = "fall in one step", quantity = step.parabola, steps = 161,
against = sweep.dt, log = true, label = "g dt^2 / 2" } },
lusine.value { { panel = "from rest", quantity = sweep.dt, label = "step" } },
lusine.value { { panel = "from rest", quantity = step.fall, label = "fall" } },
lusine.value { { panel = "from rest", quantity = step.half, label = "half-step velocity" } },
lusine.value { { panel = "the checks", quantity = step.x, label = "a state stepped" } },
lusine.value { { panel = "the checks", quantity = step.mirror_x, label = "its mirror stepped" } },
lusine.value { { panel = "the checks", quantity = checks.parabola, label = "a parabola" } },
lusine.value { { panel = "the checks", quantity = checks.drift, label = "no force, a drift" } },
lusine.value { { panel = "the checks", quantity = checks.mirror, label = "the mirror" } },
lusine.value { { panel = "the checks", quantity = checks.components, label = "in three dimensions" } },
};
void main()
{
sweep.step = sweep.step >= 160 ? 0 : sweep.step + 1;
sweep.dt = float(sweep.step) * settings.stride;
sweep.turn = sweep.turn + 0.37;
(step.dropped, step.half) = lusine.integrator.kick_drift(position = 0.0 m, velocity = 0.0 m/s,
acceleration = -settings.gravity, sweep.dt);
step.fall = -step.dropped;
step.parabola = 0.5 * settings.gravity * sweep.dt * sweep.dt;
(step.drifted, step.coasting) = lusine.integrator.kick_drift(settings.position, settings.velocity,
acceleration = 0.0 m/s2, sweep.dt);
(step.x, step.v) = lusine.integrator.kick_drift(position = settings.position * turned(sweep.turn),
velocity = settings.velocity * turned(turn = 2.0 * sweep.turn),
settings.acceleration, sweep.dt);
(step.mirror_x, step.mirror_v) = lusine.integrator.kick_drift(position = -settings.position * turned(sweep.turn),
velocity = -settings.velocity
* turned(turn = 2.0 * sweep.turn),
acceleration = -settings.acceleration, sweep.dt);
(step.space_x, step.space_v) = lusine.integrator.kick_drift(position = in_space(settings.position, sweep.turn),
velocity = moving(settings.velocity, sweep.turn),
acceleration = pulled(settings.acceleration,
sweep.turn),
sweep.dt);
(step.x0, step.v0) = lusine.integrator.kick_drift(position = in_space(settings.position, sweep.turn).x,
velocity = moving(settings.velocity, sweep.turn).x,
acceleration = pulled(settings.acceleration, sweep.turn).x,
sweep.dt);
(step.x1, step.v1) = lusine.integrator.kick_drift(position = in_space(settings.position, sweep.turn).y,
velocity = moving(settings.velocity, sweep.turn).y,
acceleration = pulled(settings.acceleration, sweep.turn).y,
sweep.dt);
(step.x2, step.v2) = lusine.integrator.kick_drift(position = in_space(settings.position, sweep.turn).z,
velocity = moving(settings.velocity, sweep.turn).z,
acceleration = pulled(settings.acceleration, sweep.turn).z,
sweep.dt);
checks.parabola = checks.parabola && close(got = step.fall, want = step.parabola, share = 0.0000005)
&& abs(step.half + 0.5 * settings.gravity * sweep.dt) <= 0.0000002 * settings.gravity * sweep.dt;
checks.drift = checks.drift
&& abs(step.drifted - (settings.position + settings.velocity * sweep.dt))
<= 0.0000005 * (settings.position + settings.velocity * sweep.dt)
&& step.coasting == settings.velocity;
checks.mirror = checks.mirror && step.mirror_x == -step.x && step.mirror_v == -step.v;
checks.components = checks.components && step.space_x.x == step.x0 && step.space_v.x == step.v0
&& step.space_x.y == step.x1 && step.space_v.y == step.v1 && step.space_x.z == step.x2
&& step.space_v.z == step.v2;
lamps.colour = lamp_colour(holds = checks.parabola && checks.drift && checks.mirror && checks.components);
}
The law's other proof (a drift and a parabola at four steps, and states at random in space) is on its proof page.