Gravity and linear drag, live¶

The acceleration a body takes under gravity and a drag opposing its motion, a = g − b v, which is zero at the terminal velocity g / b. Below, a body falls from rest under this law, stepped by the library's velocity Verlet (lusine.integrator.kick_drift, the law, lusine.integrator.kick), a hundred steps a second of a hundredth of a second each, on your GPU as you read. Beside the stepped fall the same program computes the closed form of it, v(t) = (g / b)(1 − e^(−b t)), as a quantity of its own, and checks the one against the other 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 both checks hold, and turns red for good once either fails.
The fall: the body falls from rest for ten seconds and then starts again. The chart draws the stepped speed and the closed form against the body's time, a point a step over the last thousand steps (ten seconds), so it always holds one whole fall, the two one over the other: at a hundredth of a second a step they differ by less than a point's width. Below it, the body's time, its speed and the closed form at that time.
The checks: the chart draws the difference between the stepped speed and the closed form at each of the last thousand steps on a log scale, where its rise after each restart and its slow decay toward the terminal speed read plainly though they span orders of magnitude. Below it, the worst difference so far; agrees, true while that difference has stayed within 0.2% of the terminal speed g / b at every step; and free fall, the limit with no drag: with b = 0 the law gives gravity's acceleration whatever the speed, held exactly at the body's speed 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 (9.81 m/s²), the drag rate (0.5 per second) 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 (each a row naming the quantity it shows), and the two procedures: start, run once before the first step, takes the acceleration from the body's speed at rest; main, each step, is one velocity Verlet step and the checks after it.

// The drag law's live page: a body falling from rest under gravity and a linear drag, stepped by the library's
// velocity Verlet, beside the closed form of the same fall, with the page's checks computed here on the GPU and shown
// by the panels the program declares.

namespace drag;

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

// The reader's settings: gravity, the drag rate and the step.
host struct { float<m/s2> gravity = 9.81 m/s2; float<s^-1> drag_rate = 0.5 s^-1; float<s> dt = 0.01 s; } settings;

// The body, stepped from rest; its time the sum of the steps it has taken.
struct { float<m> position; float<m/s> velocity; float<m/s2> acceleration; float<s> time; } body;

// The closed form at the body's time, v(t) = (g / b)(1 - e^(-b t)), and the checks: the stepped speed within 0.2% of
// the terminal speed g / b of the closed form at every step so far (the worst miss kept), and no drag giving free fall
// at the body's speed.
struct { float<m/s> closed_form; float<m/s> error; float<m/s> worst; bool agrees = true; bool free_fall = true; } checks;

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

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

float<m/s> fall_speed(float<m/s2> gravity, float<s^-1> drag_rate, float<s> time)
{
  return gravity / drag_rate * (1.0 - exp(-drag_rate * time));
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the fall", chart = "speed", quantity = body.velocity, steps = 1000, against = body.time,
                    label = "stepped" } },
  lusine.series { { panel = "the fall", chart = "speed", quantity = checks.closed_form, steps = 1000,
                    against = body.time, label = "closed form" } },
  lusine.value { { panel = "the fall", quantity = body.time, label = "time" } },
  lusine.value { { panel = "the fall", quantity = body.velocity, label = "stepped speed" } },
  lusine.value { { panel = "the fall", quantity = checks.closed_form, label = "closed form" } },
  lusine.series { { panel = "the checks", chart = "error", quantity = checks.error, steps = 1000, log = true,
                    label = "step error" } },
  lusine.value { { panel = "the checks", quantity = checks.worst, label = "worst error" } },
  lusine.value { { panel = "the checks", quantity = checks.agrees, label = "agrees" } },
  lusine.value { { panel = "the checks", quantity = checks.free_fall, label = "free fall" } },
};
void start() { body.acceleration = lusine.laws.drag(body.velocity, settings.gravity, settings.drag_rate); }

void main()
{
  // Every ten seconds the fall starts again from rest, so the chart's last thousand steps always hold one whole fall.
  body.position = body.time >= 10.0 s ? 0.0 m : body.position;
  body.velocity = body.time >= 10.0 s ? 0.0 m/s : body.velocity;
  body.acceleration = body.time >= 10.0 s ? lusine.laws.drag(velocity = 0.0 m/s, gravity = settings.gravity,
                                                           drag_rate = settings.drag_rate)
                                          : body.acceleration;
  body.time = body.time >= 10.0 s ? 0.0 s : body.time;
  (body.position, body.velocity) = lusine.integrator.kick_drift(body.position, body.velocity, body.acceleration,
                                                                settings.dt);
  body.acceleration = lusine.laws.drag(body.velocity, settings.gravity, settings.drag_rate);
  body.velocity = lusine.integrator.kick(body.velocity, body.acceleration, settings.dt);
  body.time = body.time + settings.dt;
  checks.closed_form = fall_speed(settings.gravity, settings.drag_rate, body.time);
  checks.error = abs(body.velocity - checks.closed_form);
  checks.worst = checks.error > checks.worst ? checks.error : checks.worst;
  checks.agrees = checks.agrees && checks.error <= 0.002 * settings.gravity / settings.drag_rate;
  checks.free_fall = lusine.laws.drag(body.velocity, settings.gravity, drag_rate = 0.0 s^-1) == settings.gravity;
  lamps.colour = lamp_colour(checks.agrees, checks.free_fall);
}

The law itself and its other checks (the law at given speeds, the terminal velocity as a rest, a body thrown in a plane) are on its proof page.