Skip to content

Drag and gravity, stepped implicitly, live

A velocity along gravity over a step under gravity and a linear drag, dw/dt = g − b w, integrated implicitly: w′ = (w + g dt) / (1 + b dt). It is stable however stiff the drag (an explicit step oscillates once b dt passes 1 and blows up past 2), never overshoots, rests exactly at the terminal velocity g / b, and is first order in the step. Below, the library's step (lusine.integrator.drag_gravity) takes a body from rest at b = 0.5 per second and a step of a hundredth of a second, a hundred steps a second, on your GPU as you read, beside the exact fall w(t) = (g / b)(1 − e^(−b t)). With it the program steps bodies at b dt of 1, 5 and 100, four bodies at four steps for the order, and a body thrown sideways in space, whose velocity is 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 fall: the first chart draws the stepped speed and the exact fall against the time over a twelve-second run, the two one over the other; the second draws the stepped speed's miss of the exact fall on a log scale, rising while the speed changes fastest and falling away as it settles. Every twelve seconds the bodies start again. Below them, the time, the speed and the miss at 12 s.
  • Stiff steps: the bodies at b dt = 1, 5 and 100, each at g / b already or rising to it.
  • The order: at b = 2, the miss of the exact fall at 1 s at steps of 1/25, 1/50, 1/100 and 1/200 s, each over the next: 2 each time, the first order.
  • The checks: the exact fall, the fine body within g dt / 2 of it at every step and within 0.01 m/s at 12 s; stable however stiff, every stiff body rising at every step and never past g / b; first order, each ratio of the misses within 0.1 of 2; free fall, and at rest, no drag giving w + g dt and a body at g / b staying there at steps of 0.01, 0.3 and 10 s, each to a rounding; in space, each component of the thrown body's step exactly the scalar step of that component alone, and the body settled to (0, 0, −g / b) within 1 mm/s after a minute.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. Gravity, the drag rates 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 run, the bodies and the checks, the panels the page draws, and the one procedure, main, run each step: one step of every body, the exact fall and the misses after it, and the checks.

// The implicit drag-and-gravity step's live page: a body falling from rest at a fine step beside the exact fall, bodies
// at steps of every stiffness rising to the terminal velocity without passing it, four steps showing the first order,
// and a body thrown sideways in space settling, 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 drag_gravity;

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

// The reader's settings: gravity, the drag rate, the fine step and the run's length in fine steps; the stiff bodies'
// steps (b dt of 1 and 5); the order's drag rate; and the stiffest body's drag rate and step (b dt of 100).
host struct {
  float<m/s2> gravity = 9.80665 m/s2; float<s^-1> drag_rate = 0.5 s^-1; float<s> dt = 0.01 s; integer steps = 1200;
  float<s> coarse = 2.0 s; float<s> coarser = 10.0 s; float<s^-1> order_rate = 2.0 s^-1;
  float<s^-1> stiff_rate = 100.0 s^-1; float<s> stiff_dt = 1.0 s;
} settings;

// The run: its step, from 0 to 1199 and round again, and the fine body's time.
struct { integer step = -1; float<s> time; } run;

// The bodies: the fine one, the exact fall at its time, its miss, and its miss at 12 s; the two stiff ones and the
// stiffest; the order's four, at steps of 1/25, 1/50, 1/100 and 1/200 s, each one's miss of the exact fall at 1 s once
// it reaches it, and the misses' ratios; and the body in space, thrown at (8, -3, 0) m/s under gravity along -z.
struct {
  float<m/s> fine; float<m/s> closed_form; float<m/s> miss; float<m/s> miss_at_12; float<m/s> coarse;
  float<m/s> coarser; float<m/s> stiff; float<m/s> coarse_before; float<m/s> coarser_before; float<m/s> stiff_before;
  float<m/s> w25; float<m/s> w50; float<m/s> w100; float<m/s> w200; float<m/s> e25; float<m/s> e50; float<m/s> e100;
  float<m/s> e200; float ratio_a; float ratio_b; float ratio_c; vector<float<m/s>, 3> space;
} body;

// The checks: the fine body within g dt / 2 of the exact fall at every step and within 0.01 m/s at 12 s; every stiff
// body rising at every step and never past g / b; halving the step halving the miss at 1 s, each ratio within 0.1 of
// 2; no drag a free fall and a body at g / b staying there, each to a rounding; and in space each component exactly the
// scalar step's, the thrown body settling to (0, 0, -g / b) within 1 mm/s.
struct { bool exact_fall = true; bool stiff = true; bool order = true; bool limits = 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);
}

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

float<m/s> step(float<m/s> velocity, float<m/s2> gravity, float<s^-1> drag_rate, float<s> dt)
{
  return lusine.integrator.drag_gravity(velocity, gravity, drag_rate, dt);
}

// A miss of the exact fall at 1 s, taken at the step that reaches it.
float<m/s> miss_at(float<m/s> miss, float<m/s> velocity, integer step, integer reaches, float<m/s2> gravity,
                   float<s^-1> drag_rate)
{
  return step == reaches - 1 ? abs(velocity - fall(gravity, drag_rate, time = 1.0 s)) : miss;
}

// Rising, and never past the terminal velocity g / b by more than two roundings.
bool rises(float<m/s> velocity, float<m/s> before, float<m/s2> gravity, float<s^-1> drag_rate)
{
  return velocity >= before && velocity <= gravity / drag_rate * 1.00000024;
}

vector<float<m/s>, 3> thrown() { return vector<float<m/s>, 3>(8.0 m/s, -3.0 m/s, 0.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); }

// In space each component the scalar step of that component alone.
bool componentwise(vector<float<m/s>, 3> velocity, float<m/s2> gravity, float<s^-1> drag_rate, float<s> dt)
{
  let three = lusine.integrator.drag_gravity(velocity, gravity = down(gravity), drag_rate, dt);
  return three.x == step(velocity = velocity.x, gravity = 0.0 m/s2, drag_rate, dt)
      && three.y == step(velocity = velocity.y, gravity = 0.0 m/s2, drag_rate, dt)
      && three.z == step(velocity = velocity.z, gravity = -gravity, drag_rate, dt);
}

bool limits(float<m/s2> gravity, float<s^-1> drag_rate)
{
  let free = step(velocity = 3.0 m/s, gravity, drag_rate = 0.0 s^-1, dt = 0.1 s);
  let rest = gravity / drag_rate;
  return abs(free - (3.0 m/s + gravity * 0.1 s)) <= 0.00000024 * free
      && abs(step(velocity = rest, gravity, drag_rate, dt = 0.01 s) - rest) <= 0.00000012 * rest
      && abs(step(velocity = rest, gravity, drag_rate, dt = 0.3 s) - rest) <= 0.00000012 * rest
      && abs(step(velocity = rest, gravity, drag_rate, dt = 10.0 s) - rest) <= 0.00000012 * rest;
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the fall", chart = "speed", quantity = body.fine, steps = 1200, against = run.time,
                    label = "stepped" } },
  lusine.series { { panel = "the fall", chart = "speed", quantity = body.closed_form, steps = 1200, against = run.time,
                    label = "g / b (1 - exp(-b t))" } },
  lusine.series { { panel = "the fall", chart = "miss", quantity = body.miss, steps = 1200, against = run.time,
                    log = true, label = "stepped - exact" } },
  lusine.value { { panel = "the fall", quantity = run.time, label = "time" } },
  lusine.value { { panel = "the fall", quantity = body.fine, label = "speed" } },
  lusine.value { { panel = "the fall", quantity = body.miss_at_12, label = "miss at 12 s" } },
  lusine.value { { panel = "stiff steps", quantity = body.coarse, label = "b dt = 1" } },
  lusine.value { { panel = "stiff steps", quantity = body.coarser, label = "b dt = 5" } },
  lusine.value { { panel = "stiff steps", quantity = body.stiff, label = "b dt = 100" } },
  lusine.value { { panel = "the order", quantity = body.ratio_a, label = "1/25 s over 1/50 s" } },
  lusine.value { { panel = "the order", quantity = body.ratio_b, label = "1/50 s over 1/100 s" } },
  lusine.value { { panel = "the order", quantity = body.ratio_c, label = "1/100 s over 1/200 s" } },
  lusine.value { { panel = "the checks", quantity = checks.exact_fall, label = "the exact fall" } },
  lusine.value { { panel = "the checks", quantity = checks.stiff, label = "stable however stiff" } },
  lusine.value { { panel = "the checks", quantity = checks.order, label = "first order" } },
  lusine.value { { panel = "the checks", quantity = checks.limits, label = "free fall, and at rest" } },
  lusine.value { { panel = "the checks", quantity = checks.space, label = "in space" } },
};

void main()
{
  // Every run's steps every body starts again from rest, the thrown one from its throw.
  run.step = run.step + 1 >= settings.steps ? 0 : run.step + 1;
  body.fine = run.step == 0 ? 0.0 m/s : body.fine;
  body.coarse = run.step == 0 ? 0.0 m/s : body.coarse;
  body.coarser = run.step == 0 ? 0.0 m/s : body.coarser;
  body.stiff = run.step == 0 ? 0.0 m/s : body.stiff;
  body.w25 = run.step == 0 ? 0.0 m/s : body.w25;
  body.w50 = run.step == 0 ? 0.0 m/s : body.w50;
  body.w100 = run.step == 0 ? 0.0 m/s : body.w100;
  body.w200 = run.step == 0 ? 0.0 m/s : body.w200;
  body.space = run.step == 0 ? thrown() : body.space;
  body.coarse_before = body.coarse;
  body.coarser_before = body.coarser;
  body.stiff_before = body.stiff;
  checks.space = checks.space && componentwise(velocity = body.space, settings.gravity, settings.drag_rate,
                                               dt = 0.05 s);
  body.fine = step(velocity = body.fine, settings.gravity, settings.drag_rate, settings.dt);
  body.coarse = step(velocity = body.coarse, settings.gravity, settings.drag_rate, dt = settings.coarse);
  body.coarser = step(velocity = body.coarser, settings.gravity, settings.drag_rate, dt = settings.coarser);
  body.stiff = step(velocity = body.stiff, settings.gravity, drag_rate = settings.stiff_rate, dt = settings.stiff_dt);
  body.w25 = step(velocity = body.w25, settings.gravity, drag_rate = settings.order_rate, dt = 0.04 s);
  body.w50 = step(velocity = body.w50, settings.gravity, drag_rate = settings.order_rate, dt = 0.02 s);
  body.w100 = step(velocity = body.w100, settings.gravity, drag_rate = settings.order_rate, dt = 0.01 s);
  body.w200 = step(velocity = body.w200, settings.gravity, drag_rate = settings.order_rate, dt = 0.005 s);
  body.space = lusine.integrator.drag_gravity(velocity = body.space, gravity = down(settings.gravity),
                                              settings.drag_rate, dt = 0.05 s);
  run.time = float(run.step + 1) * settings.dt;
  body.closed_form = fall(settings.gravity, settings.drag_rate, run.time);
  body.miss = abs(body.fine - body.closed_form);
  body.miss_at_12 = run.step + 1 == settings.steps ? body.miss : body.miss_at_12;
  body.e25 = miss_at(miss = body.e25, velocity = body.w25, run.step, reaches = 25, settings.gravity,
                     drag_rate = settings.order_rate);
  body.e50 = miss_at(miss = body.e50, velocity = body.w50, run.step, reaches = 50, settings.gravity,
                     drag_rate = settings.order_rate);
  body.e100 = miss_at(miss = body.e100, velocity = body.w100, run.step, reaches = 100, settings.gravity,
                      drag_rate = settings.order_rate);
  body.e200 = miss_at(miss = body.e200, velocity = body.w200, run.step, reaches = 200, settings.gravity,
                      drag_rate = settings.order_rate);
  body.ratio_a = body.e25 / body.e50;
  body.ratio_b = body.e50 / body.e100;
  body.ratio_c = body.e100 / body.e200;
  checks.exact_fall = checks.exact_fall && body.miss <= 0.5 * settings.gravity * settings.dt
                      && (run.step + 1 != settings.steps || body.miss <= 0.01 m/s);
  checks.stiff = checks.stiff
                 && rises(velocity = body.coarse, before = body.coarse_before, settings.gravity, settings.drag_rate)
                 && rises(velocity = body.coarser, before = body.coarser_before, settings.gravity, settings.drag_rate)
                 && rises(velocity = body.stiff, before = body.stiff_before, settings.gravity,
                          drag_rate = settings.stiff_rate);
  checks.order = checks.order
                 && (run.step < 199 || abs(body.ratio_a - 2.0) < 0.1 && abs(body.ratio_b - 2.0) < 0.1
                                       && abs(body.ratio_c - 2.0) < 0.1);
  checks.limits = checks.limits && limits(settings.gravity, settings.drag_rate);
  checks.space = checks.space
                 && (run.step + 1 != settings.steps
                     || length(body.space - vector<float<m/s>, 3>(0.0 m/s, 0.0 m/s, -settings.gravity
                                                                                    / settings.drag_rate))
                        < 0.001 m/s);
  lamps.colour = lamp_colour(holds = checks.exact_fall && checks.stiff && checks.order && checks.limits
                                     && checks.space);
}

The law's other proof (the curves at three stiffnesses, and two hundred steps in space at random) is on its proof page.