Drag and gravity: the terminal velocity, at any stiffness¶
A velocity along gravity over a step under gravity and a linear (Stokes) drag, dw/dt = g - b w (lusine.laws.drag's acceleration), integrated implicitly: w' = (w + g dt) / (1 + b dt).
Stable however stiff (b dt far above 1, where an explicit step oscillates and past 2 blows up), never overshooting, at rest exactly at the terminal velocity g / b; first order in dt. Written once over a vector of any length: the drag is linear, so each component is the law along its own axis, a single value the length 1. Source: Hairer & Wanner, Solving Ordinary Differential Equations II, IV.3 (the implicit Euler method on a stiff decay). Checked against a reference in numpy, which has its own tests.
- Source
- Hairer & Wanner, Solving ODEs II, IV.3 (the implicit Euler method on a stiff decay).
- Proof
- curve, order, limit, invariant.
drag_gravity.trait
namespace lusine.integrator;
vector<float<m/s>, _> drag_gravity(vector<float<m/s>, _> velocity, vector<float<m/s2>, _> gravity, float<s^-1> drag_rate,
float<s> dt)
{
return (velocity + gravity * dt) / (1.0 + drag_rate * dt);
}
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): stable however stiff, tending to g / b. It is written once over a vector of any length: a binding over single speeds takes drag_gravity@1, the scalar form, and one over a space's velocities drag_gravity@3 (as standalone functions, drag_gravity_1 to drag_gravity_3).
In [1]:
Copied!
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
from lusine.kernels._si import G0
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
from lusine.kernels._si import G0
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
The kernel¶
The law as the library states it, and as this page runs it.
In [2]:
Copied!
print(trait_run.source('integrator/drag_gravity.trait'))
print(trait_run.source('integrator/drag_gravity.trait'))
// lusine.integrator.drag_gravity: a velocity along gravity over a step under gravity and a linear (Stokes) drag,
// dw/dt = g - b w (lusine.laws.drag's acceleration), integrated implicitly: w' = (w + g dt) / (1 + b dt). Stable
// however stiff (b dt far above 1, where an explicit step oscillates and past 2 blows up), never overshooting, at rest
// exactly at the terminal velocity g / b; first order in dt. Written once over a vector of any length:
// the drag is linear, so each component is the law along its own axis, a single value the length 1.
// Source: Hairer & Wanner, Solving Ordinary Differential Equations II, IV.3 (the implicit Euler method on a stiff
// decay).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.integrator;
vector<float<m/s>, _> drag_gravity(vector<float<m/s>, _> velocity, vector<float<m/s2>, _> gravity, float<s^-1> drag_rate,
float<s> dt)
{
return (velocity + gravity * dt) / (1.0 + drag_rate * dt);
}
The program¶
A lattice of bodies, each with its own drag rate and step, so one run steps every case of a section side by side; main() is one step of the kernel.
In [3]:
Copied!
FALL = '''namespace fall;
lusine.version version = { { number = 1 } };
struct body { float<m/s> velocity; host float<s^-1> drag_rate; host float<s> dt; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { float<m/s2> gravity; } settings;
void main() { bodies.velocity = lusine.integrator.drag_gravity(bodies.velocity, settings.gravity, bodies.drag_rate,
bodies.dt); }
'''
print(FALL)
def fall(steps, drag_rate, dt, velocity=0.0):
"""Each body's velocity after every step, one row a step: bodies of the given drag rates and steps."""
n = len(dt)
run = trait_run.Program(['integrator/drag_gravity.trait'], FALL, {'fall.bodies': [n]})
run.store('fall.settings', 'gravity', G0)
run.store('fall.bodies', 'drag_rate', np.broadcast_to(drag_rate, (n,)))
run.store('fall.bodies', 'dt', dt)
run.store('fall.bodies', 'velocity', np.broadcast_to(velocity, (n,)))
return run.steps('fall.main', steps, [('fall.bodies', 'velocity')])[0]
FALL = '''namespace fall;
lusine.version version = { { number = 1 } };
struct body { float<m/s> velocity; host float<s^-1> drag_rate; host float<s> dt; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { float<m/s2> gravity; } settings;
void main() { bodies.velocity = lusine.integrator.drag_gravity(bodies.velocity, settings.gravity, bodies.drag_rate,
bodies.dt); }
'''
print(FALL)
def fall(steps, drag_rate, dt, velocity=0.0):
"""Each body's velocity after every step, one row a step: bodies of the given drag rates and steps."""
n = len(dt)
run = trait_run.Program(['integrator/drag_gravity.trait'], FALL, {'fall.bodies': [n]})
run.store('fall.settings', 'gravity', G0)
run.store('fall.bodies', 'drag_rate', np.broadcast_to(drag_rate, (n,)))
run.store('fall.bodies', 'dt', dt)
run.store('fall.bodies', 'velocity', np.broadcast_to(velocity, (n,)))
return run.steps('fall.main', steps, [('fall.bodies', 'velocity')])[0]
namespace fall;
lusine.version version = { { number = 1 } };
struct body { float<m/s> velocity; host float<s^-1> drag_rate; host float<s> dt; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { float<m/s2> gravity; } settings;
void main() { bodies.velocity = lusine.integrator.drag_gravity(bodies.velocity, settings.gravity, bodies.drag_rate,
bodies.dt); }
Curve: w = g/b (1 - e^(-bt))¶
At a fine step and at steps of b dt = 1 and 5, side by side, against the analytic approach: the fine step lands on it at t = 12 s within 0.01 m/s, and the coarse ones rise to g / b without passing it.
In [4]:
Copied!
b = 0.5
dts = np.array([0.01, 2.0, 10.0])
w = fall(1200, b, dts) # w[k, i]: body i after k + 1 steps
exact = lambda t: G0 / b * (1 - np.exp(-b * np.asarray(t)))
assert abs(w[-1, 0] - exact(12.0)) < 0.01
for i, dt in enumerate(dts[1:], 1):
ws = w[:int(12 / dt), i]
assert np.all(np.diff(ws) >= 0) and np.all(ws <= np.float32(G0 / b) * (1 + 2 * EPS))
t = np.linspace(0, 12, 200)
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(t, exact(t), label='g/b (1 - exp(-b t))')
for i, (dt, style) in enumerate(zip(dts, ('-', 's', '^'))):
k = int(12 / dt)
ax.plot(np.arange(k + 1) * dt, np.concatenate([[0.0], w[:k, i]]), style, ms=3, label=f'b dt = {b * dt:g}')
ax.set_xlabel('t (s)')
ax.set_ylabel('w (m/s)')
ax.legend()
plt.show()
b = 0.5
dts = np.array([0.01, 2.0, 10.0])
w = fall(1200, b, dts) # w[k, i]: body i after k + 1 steps
exact = lambda t: G0 / b * (1 - np.exp(-b * np.asarray(t)))
assert abs(w[-1, 0] - exact(12.0)) < 0.01
for i, dt in enumerate(dts[1:], 1):
ws = w[:int(12 / dt), i]
assert np.all(np.diff(ws) >= 0) and np.all(ws <= np.float32(G0 / b) * (1 + 2 * EPS))
t = np.linspace(0, 12, 200)
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(t, exact(t), label='g/b (1 - exp(-b t))')
for i, (dt, style) in enumerate(zip(dts, ('-', 's', '^'))):
k = int(12 / dt)
ax.plot(np.arange(k + 1) * dt, np.concatenate([[0.0], w[:k, i]]), style, ms=3, label=f'b dt = {b * dt:g}')
ax.set_xlabel('t (s)')
ax.set_ylabel('w (m/s)')
ax.legend()
plt.show()
Order: first¶
At b = 2, halving the step halves the error at t = 1 s: the implicit Euler step is first order. The four step sizes are four bodies of one run, each read at its own step count.
In [5]:
Copied!
steps = np.array([25, 50, 100, 200])
w = fall(200, 2.0, 1.0 / steps)
errs = np.array([abs(w[n - 1, i] - G0 / 2.0 * (1 - np.exp(-2.0))) for i, n in enumerate(steps)])
ratios = errs[:-1] / errs[1:]
print('error ratios', np.round(ratios, 3))
assert all(abs(r - 2) < 0.1 for r in ratios)
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.loglog(1 / steps, errs, 'o-', label='error at t = 1 s')
ax.loglog(1 / steps, errs[0] * (steps[0] / steps), '--', label='slope 1')
ax.set_xlabel('dt (s)')
ax.legend()
plt.show()
steps = np.array([25, 50, 100, 200])
w = fall(200, 2.0, 1.0 / steps)
errs = np.array([abs(w[n - 1, i] - G0 / 2.0 * (1 - np.exp(-2.0))) for i, n in enumerate(steps)])
ratios = errs[:-1] / errs[1:]
print('error ratios', np.round(ratios, 3))
assert all(abs(r - 2) < 0.1 for r in ratios)
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.loglog(1 / steps, errs, 'o-', label='error at t = 1 s')
ax.loglog(1 / steps, errs[0] * (steps[0] / steps), '--', label='slope 1')
ax.set_xlabel('dt (s)')
ax.legend()
plt.show()
error ratios [1.987 1.993 1.995]
Limit: no drag is free fall; at g/b it stays¶
With b = 0 the step is w + g dt, to binary32's rounding of the sum; a body at the terminal velocity g / b is left there, to an ulp.
In [6]:
Copied!
free = trait_run.law('integrator/drag_gravity.trait', 'drag_gravity', velocity=3.0, gravity=G0, drag_rate=0.0,
dt=0.1)[0]
assert abs(free - (3.0 + G0 * 0.1)) <= 2 * EPS * (3.0 + G0 * 0.1)
rest = float(np.float32(G0 / 0.5))
w = trait_run.law('integrator/drag_gravity.trait', 'drag_gravity', velocity=rest, gravity=G0, drag_rate=0.5,
dt=[0.01, 0.3, 10.0])
assert np.all(np.abs(w - rest) <= EPS * rest), w - rest
print('free fall', free, 'and rest', w)
free = trait_run.law('integrator/drag_gravity.trait', 'drag_gravity', velocity=3.0, gravity=G0, drag_rate=0.0,
dt=0.1)[0]
assert abs(free - (3.0 + G0 * 0.1)) <= 2 * EPS * (3.0 + G0 * 0.1)
rest = float(np.float32(G0 / 0.5))
w = trait_run.law('integrator/drag_gravity.trait', 'drag_gravity', velocity=rest, gravity=G0, drag_rate=0.5,
dt=[0.01, 0.3, 10.0])
assert np.all(np.abs(w - rest) <= EPS * rest), w - rest
print('free fall', free, 'and rest', w)
free fall 3.9806649684906006 and rest [19.61330032 19.61330032 19.61329842]
Invariant: monotone at any stiffness¶
At b dt = 100 an explicit step would oscillate and blow up; the implicit one rises to g / b without overshoot, step by step.
In [7]:
Copied!
w = fall(50, 100.0, np.array([1.0]))[:, 0]
assert np.all(np.diff(w) >= 0) and w[0] > 0 and np.all(w <= np.float32(G0 / 100.0) * (1 + 2 * EPS))
assert abs(w[-1] - G0 / 100.0) <= 2 * EPS * G0 / 100.0
print('no overshoot at b dt = 100:', w[:4], '...', w[-1])
w = fall(50, 100.0, np.array([1.0]))[:, 0]
assert np.all(np.diff(w) >= 0) and w[0] > 0 and np.all(w <= np.float32(G0 / 100.0) * (1 + 2 * EPS))
assert abs(w[-1] - G0 / 100.0) <= 2 * EPS * G0 / 100.0
print('no overshoot at b dt = 100:', w[:4], '...', w[-1])
no overshoot at b dt = 100: [0.09709555 0.09805689 0.0980664 0.0980665 ] ... 0.09806650131940842
Limit: a velocity in space¶
The law at length 3 (drag_gravity@3): over random velocities and gravities, every component of a step is the scalar form's on that component, bit for bit; and a body thrown sideways settles to (0, 0, −g / b), its level speed decaying.
In [8]:
Copied!
rng = np.random.default_rng(90)
w, g = rng.uniform(-60, 60, (200, 3)), rng.uniform(-12, 12, (200, 3))
b, dt = rng.uniform(0, 50, 200), rng.uniform(1e-3, 1.0, 200)
three = trait_run.law('integrator/drag_gravity.trait', 'drag_gravity', vector_length=3, velocity=w, gravity=g, drag_rate=b, dt=dt)
for c in range(3):
one = trait_run.law('integrator/drag_gravity.trait', 'drag_gravity', velocity=w[:, c], gravity=g[:, c],
drag_rate=b, dt=dt)
assert np.array_equal(three[:, c], one)
SPACE = '''namespace space;
lusine.version version = { { number = 1 } };
struct body { vector<float<m/s>, 3> velocity; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { vector<float<m/s2>, 3> gravity; float<s^-1> drag_rate; float<s> dt; } settings;
void main() { bodies.velocity = lusine.integrator.drag_gravity(bodies.velocity, settings.gravity, settings.drag_rate,
settings.dt); }
'''
run = trait_run.Program(['integrator/drag_gravity.trait'], SPACE, {'space.bodies': [1]})
run.store('space.bodies', 'velocity', [[8.0, -3.0, 0.0]])
run.store('space.settings', 'gravity', [0.0, 0.0, -G0])
run.store('space.settings', 'drag_rate', 0.5)
run.store('space.settings', 'dt', 0.05)
v = run.steps('space.main', 600, [('space.bodies', 'velocity')])[0][:, 0]
assert np.all(np.diff(np.abs(v[:, 0])) <= 0) and np.all(np.diff(v[:, 2]) <= 0)
assert np.max(np.abs(v[-1] - [0.0, 0.0, -G0 / 0.5])) < 1e-3
print('200 random steps in space, each component the scalar form; after 30 s', v[-1])
rng = np.random.default_rng(90)
w, g = rng.uniform(-60, 60, (200, 3)), rng.uniform(-12, 12, (200, 3))
b, dt = rng.uniform(0, 50, 200), rng.uniform(1e-3, 1.0, 200)
three = trait_run.law('integrator/drag_gravity.trait', 'drag_gravity', vector_length=3, velocity=w, gravity=g, drag_rate=b, dt=dt)
for c in range(3):
one = trait_run.law('integrator/drag_gravity.trait', 'drag_gravity', velocity=w[:, c], gravity=g[:, c],
drag_rate=b, dt=dt)
assert np.array_equal(three[:, c], one)
SPACE = '''namespace space;
lusine.version version = { { number = 1 } };
struct body { vector<float<m/s>, 3> velocity; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { vector<float<m/s2>, 3> gravity; float<s^-1> drag_rate; float<s> dt; } settings;
void main() { bodies.velocity = lusine.integrator.drag_gravity(bodies.velocity, settings.gravity, settings.drag_rate,
settings.dt); }
'''
run = trait_run.Program(['integrator/drag_gravity.trait'], SPACE, {'space.bodies': [1]})
run.store('space.bodies', 'velocity', [[8.0, -3.0, 0.0]])
run.store('space.settings', 'gravity', [0.0, 0.0, -G0])
run.store('space.settings', 'drag_rate', 0.5)
run.store('space.settings', 'dt', 0.05)
v = run.steps('space.main', 600, [('space.bodies', 'velocity')])[0][:, 0]
assert np.all(np.diff(np.abs(v[:, 0])) <= 0) and np.all(np.diff(v[:, 2]) <= 0)
assert np.max(np.abs(v[-1] - [0.0, 0.0, -G0 / 0.5])) < 1e-3
print('200 random steps in space, each component the scalar form; after 30 s', v[-1])
200 random steps in space, each component the scalar form; after 30 s [ 2.94290976e-06 -1.10359031e-06 -1.96132469e+01]