Semi-implicit Euler, live¶
The semi-implicit (symplectic) Euler step: the velocity first, v′ = v + a(x) dt, and then the position with the new velocity, x′ = x + v′ dt. It is first order and takes one force a step, and it keeps a conservative system's energy within a band of its start rather than letting it drift; the explicit Euler, which moves the position with the old velocity, gains energy every step. Below, a spring of period 1 s is stepped by the library's step (lusine.integrator.symplectic_euler) from x = 1 m at rest, a hundred steps a period and a hundred steps a second, on your GPU as you read, and beside it the same spring stepped by the explicit Euler, which the program writes for contrast. With them the program steps five springs at five steps for the order, a mirrored pair, a fall, and a spring in a plane, whose position is a vector of two.
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 spring: the chart draws each spring's energy's drift from its start, |E / E₀ − 1|, against the time on a log scale over a ten-second run: the symplectic step's swings within a band of about 3% and comes back to it every half period, while the explicit step's climbs without end, past its start's own size within a few seconds. Every ten seconds the springs start again. Below it, the time, the symplectic step's worst drift in the run's first five periods and in its last five, the same, and the explicit step's drift now.
- The checks: the ratio of the misses of cos(ω t) at 0.3 s at steps of 1/400 and 1/800 s, and the fall's miss of −g t² / 2 at 0.5 s; energy in a band, the symplectic step's drift below ω dt at every step and its worst in the last five periods never more than 1% past its worst in the first five, and within 1% of it at the run's end; first order, halving the step halving the miss at 0.3 s, each ratio within 0.1 of 2 from 1/50 s to 1/800 s; the mirror, the mirrored state stepping to the mirrored state, bit for bit, at every step; a fall, under a constant force the step off −g t² / 2 by the first order's g t dt / 2, at 1 s within 10⁻⁵ m of −g (1 + dt) / 2 and twice its miss at 0.5 s within 1%; in a plane, the spring in a plane within ω dt of its circle, and each axis exactly the scalar spring started on that axis's component.
Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The spring, the step, the run and gravity are the program's settings, which the host may write; this page offers no control to change them yet.
The program¶
The run, the springs and the checks, the panels the page draws, and the one procedure, main, run each step: one step of every spring and of the fall, the drifts and misses after it, and the checks.
// The semi-implicit Euler step's live page: a spring of period 1 s stepped by the library's symplectic Euler, its
// energy kept in a band, beside the explicit Euler the program writes for contrast, whose energy grows; five springs
// at five steps for the order; a mirrored pair; a fall; and a spring in a plane, each axis the scalar step's, with the
// page's checks computed here on the GPU and shown by the panels the program declares.
namespace symplectic_euler;
lusine.version version = { { number = 1 } };
// The reader's settings: the spring's angular frequency (a period of 1 s) and where it is at 0.3 s, cos(0.6 pi) of its
// start (given, since a cosine taken on the device carries its own error), the step, the run's length in steps and
// gravity for the fall.
host struct {
float<s^-1> omega = 6.2831853 s^-1; float<m> at_03 = -0.30901699 m; float<s> dt = 0.01 s; integer steps = 1000;
float<m/s2> gravity = 9.81 m/s2;
} settings;
// The run: its step, from 0 to 999 and round again, and its time.
struct { integer step = -1; float<s> time; } run;
// The spring, from x = 1 m at rest: stepped by the kernel and by the explicit Euler; each one's energy against its
// start, |E / E0 - 1|; the worst of the kernel's in the run's first five periods and in its last five.
struct {
float<m> x; float<m/s> v; float<m> ex; float<m/s> ev; float drift; float explicit_drift; float worst_first;
float worst_last;
} spring;
// The order's five springs at steps of 1/50 to 1/800 s, each one's miss of cos(omega t) at 0.3 s once it reaches it,
// and the misses' ratios; the mirrored pair; the fall from rest; and the spring in a plane with its two scalar springs.
struct {
float<m> x50; float<m/s> v50; float<m> x100; float<m/s> v100; float<m> x200; float<m/s> v200; float<m> x400;
float<m/s> v400; float<m> x800; float<m/s> v800; float<m> e50; float<m> e100; float<m> e200; float<m> e400;
float<m> e800; float ratio;
float<m> mx; float<m/s> mv; float<m> nx; float<m/s> nv; float<m> fx; float<m/s> fv; float<m> half_miss;
vector<float<m>, 2> plane_x; vector<float<m/s>, 2> plane_v; float<m> ax; float<m/s> av; float<m> bx;
float<m/s> bv;
} others;
// The checks: the kernel's energy within omega dt of its start at every step, and its worst in the last five periods
// never more than 1% past its worst in the first five, and within 1% of it at the run's end (a band, no drift); halving
// the step halving the miss at 0.3 s, each ratio within 0.1 of 2; the mirrored state stepping to the mirrored state
// exactly at every step; a fall off -g t^2 / 2 by the first order's g t dt / 2, at 1 s within 10^-5 m and twice its
// miss at 0.5 s within 1%; and the spring in a plane within omega dt of its circle, each axis exactly the scalar
// spring's.
struct { bool band = true; bool order = true; bool mirror = true; bool fall = true; bool plane = 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 spring's pull, the program's own force kernel, bound before each step as a game's are.
float<m/s2> pull(float<m> position, float<s^-1> omega) { return -omega * omega * position; }
vector<float<m/s2>, 2> pull_plane(vector<float<m>, 2> position, float<s^-1> omega)
{
return -omega * omega * position;
}
// The explicit Euler, written here for contrast: both the position and the velocity from the old state.
(float<m>, float<m/s>) explicit_euler(float<m> x, float<m/s> v, float<s^-1> omega, float<s> dt)
{
return (x + v * dt, v + pull(position = x, omega) * dt);
}
// The energy against a start of x = 1 m at rest: (v^2 / 2 + omega^2 x^2 / 2) / (omega^2 / 2), less one.
float drift(float<m> x, float<m/s> v, float<s^-1> omega)
{
return abs((v * v + omega * omega * x * x) / (omega * omega * 1.0 m2) - 1.0);
}
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 spring", chart = "energy's drift", quantity = spring.drift, steps = 1000,
against = run.time, log = true, label = "symplectic Euler" } },
lusine.series { { panel = "the spring", chart = "energy's drift", quantity = spring.explicit_drift, steps = 1000,
against = run.time, log = true, label = "explicit Euler" } },
lusine.value { { panel = "the spring", quantity = run.time, label = "time" } },
lusine.value { { panel = "the spring", quantity = spring.worst_first, label = "worst, first five periods" } },
lusine.value { { panel = "the spring", quantity = spring.worst_last, label = "worst, last five periods" } },
lusine.value { { panel = "the spring", quantity = spring.explicit_drift, label = "explicit Euler's" } },
lusine.value { { panel = "the checks", quantity = others.ratio, label = "1/400 s over 1/800 s" } },
lusine.value { { panel = "the checks", quantity = others.half_miss, label = "the fall's miss at 0.5 s" } },
lusine.value { { panel = "the checks", quantity = checks.band, label = "energy in a band" } },
lusine.value { { panel = "the checks", quantity = checks.order, label = "first order" } },
lusine.value { { panel = "the checks", quantity = checks.mirror, label = "the mirror" } },
lusine.value { { panel = "the checks", quantity = checks.fall, label = "a fall" } },
lusine.value { { panel = "the checks", quantity = checks.plane, label = "in a plane" } },
};
void main()
{
// Every run's steps every spring starts again.
run.step = run.step + 1 >= settings.steps ? 0 : run.step + 1;
spring.x = run.step == 0 ? 1.0 m : spring.x;
spring.v = run.step == 0 ? 0.0 m/s : spring.v;
spring.ex = run.step == 0 ? 1.0 m : spring.ex;
spring.ev = run.step == 0 ? 0.0 m/s : spring.ev;
spring.worst_first = run.step == 0 ? 0.0 : spring.worst_first;
spring.worst_last = run.step == 0 ? 0.0 : spring.worst_last;
others.x50 = run.step == 0 ? 1.0 m : others.x50;
others.v50 = run.step == 0 ? 0.0 m/s : others.v50;
others.x100 = run.step == 0 ? 1.0 m : others.x100;
others.v100 = run.step == 0 ? 0.0 m/s : others.v100;
others.x200 = run.step == 0 ? 1.0 m : others.x200;
others.v200 = run.step == 0 ? 0.0 m/s : others.v200;
others.x400 = run.step == 0 ? 1.0 m : others.x400;
others.v400 = run.step == 0 ? 0.0 m/s : others.v400;
others.x800 = run.step == 0 ? 1.0 m : others.x800;
others.v800 = run.step == 0 ? 0.0 m/s : others.v800;
others.mx = run.step == 0 ? 0.3 m : others.mx;
others.mv = run.step == 0 ? 1.1 m/s : others.mv;
others.nx = run.step == 0 ? -0.3 m : others.nx;
others.nv = run.step == 0 ? -1.1 m/s : others.nv;
others.fx = run.step == 0 ? 0.0 m : others.fx;
others.fv = run.step == 0 ? 0.0 m/s : others.fv;
others.plane_x = run.step == 0 ? vector<float<m>, 2>(1.0 m, 0.0 m) : others.plane_x;
others.plane_v = run.step == 0 ? vector<float<m/s>, 2>(0.0 m/s, settings.omega * 1.0 m) : others.plane_v;
others.ax = run.step == 0 ? 1.0 m : others.ax;
others.av = run.step == 0 ? 0.0 m/s : others.av;
others.bx = run.step == 0 ? 0.0 m : others.bx;
others.bv = run.step == 0 ? settings.omega * 1.0 m : others.bv;
// The explicit Euler, the step this kernel is not: the position with the old velocity.
(spring.ex, spring.ev) = explicit_euler(x = spring.ex, v = spring.ev, settings.omega, settings.dt);
(spring.x, spring.v) = lusine.integrator.symplectic_euler(position = spring.x, velocity = spring.v,
acceleration = pull(position = spring.x, settings.omega),
settings.dt);
(others.x50, others.v50) = lusine.integrator.symplectic_euler(position = others.x50, velocity = others.v50,
acceleration = pull(position = others.x50,
settings.omega), dt = 0.02 s);
(others.x100, others.v100) = lusine.integrator.symplectic_euler(position = others.x100, velocity = others.v100,
acceleration = pull(position = others.x100,
settings.omega), dt = 0.01 s);
(others.x200, others.v200) = lusine.integrator.symplectic_euler(position = others.x200, velocity = others.v200,
acceleration = pull(position = others.x200,
settings.omega), dt = 0.005 s);
(others.x400, others.v400) = lusine.integrator.symplectic_euler(position = others.x400, velocity = others.v400,
acceleration = pull(position = others.x400,
settings.omega),
dt = 0.0025 s);
(others.x800, others.v800) = lusine.integrator.symplectic_euler(position = others.x800, velocity = others.v800,
acceleration = pull(position = others.x800,
settings.omega),
dt = 0.00125 s);
(others.mx, others.mv) = lusine.integrator.symplectic_euler(position = others.mx, velocity = others.mv,
acceleration = pull(position = others.mx,
settings.omega), dt = 0.013 s);
(others.nx, others.nv) = lusine.integrator.symplectic_euler(position = others.nx, velocity = others.nv,
acceleration = pull(position = others.nx,
settings.omega), dt = 0.013 s);
(others.fx, others.fv) = lusine.integrator.symplectic_euler(position = others.fx, velocity = others.fv,
acceleration = -settings.gravity, settings.dt);
(others.plane_x, others.plane_v) = lusine.integrator.symplectic_euler(position = others.plane_x,
velocity = others.plane_v,
acceleration = pull_plane(
position = others.plane_x, settings.omega),
settings.dt);
(others.ax, others.av) = lusine.integrator.symplectic_euler(position = others.ax, velocity = others.av,
acceleration = pull(position = others.ax,
settings.omega), settings.dt);
(others.bx, others.bv) = lusine.integrator.symplectic_euler(position = others.bx, velocity = others.bv,
acceleration = pull(position = others.bx,
settings.omega), settings.dt);
run.time = float(run.step + 1) * settings.dt;
spring.drift = drift(spring.x, spring.v, settings.omega);
spring.explicit_drift = drift(x = spring.ex, v = spring.ev, settings.omega);
spring.worst_first = run.step < 500 && spring.drift > spring.worst_first ? spring.drift : spring.worst_first;
spring.worst_last = run.step >= 500 && spring.drift > spring.worst_last ? spring.drift : spring.worst_last;
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.e800 = miss_at(miss = others.e800, x = others.x800, run.step, reaches = 240, settings.at_03);
others.ratio = others.e400 / others.e800;
others.half_miss = run.step == 49 ? others.fx + 0.5 * settings.gravity * 0.25 s2 : others.half_miss;
checks.band = checks.band && spring.drift < settings.omega * settings.dt
&& spring.worst_last <= 1.01 * spring.worst_first
&& (run.step + 1 != settings.steps || abs(spring.worst_last / spring.worst_first - 1.0) < 0.01);
checks.order = checks.order
&& (run.step < 239 || abs(others.e50 / others.e100 - 2.0) < 0.1
&& abs(others.e100 / others.e200 - 2.0) < 0.1
&& abs(others.e200 / others.e400 - 2.0) < 0.1
&& abs(others.e400 / others.e800 - 2.0) < 0.1);
checks.mirror = checks.mirror && others.nx == -others.mx && others.nv == -others.mv;
checks.fall = checks.fall
&& (run.step != 99 || abs(others.fx + 0.5 * settings.gravity * 1.01 s2) < 0.00001 m
&& abs((others.fx + 0.5 * settings.gravity * 1.0 s2) / others.half_miss - 2.0)
< 0.01);
checks.plane = checks.plane
&& abs(length(others.plane_x) - 1.0 m) < settings.omega * settings.dt * 1.0 m
&& others.plane_x.x == others.ax && others.plane_v.x == others.av && others.plane_x.y == others.bx
&& others.plane_v.y == others.bv;
lamps.colour = lamp_colour(holds = checks.band && checks.order && checks.mirror && checks.fall && checks.plane);
}
The law's other proof (a hundred periods, and fifty in a plane) is on its proof page.