Rolling resistance: slowed to rest, never reversed¶
A body rolling along the ground over a step, slowed by its rolling resistance, the deceleration mu_r g against its motion: v' = v - sign(v) mu_r g dt, never past zero, so a body slowed to rest stops rather than reversing.
The coefficient is the ground's (the game's: the library states no rolling_friction), the gravity the world's magnitude, both passed in. A cut-off speed below which a game stops a body is the game's rule, after this kernel. Source: Gillespie, Fundamentals of Vehicle Dynamics, 4 (rolling resistance); Hibbeler, Engineering Mechanics: Statics, 8.8. Checked against a reference in numpy, which has its own tests.
- Source
- Gillespie, Fundamentals of Vehicle Dynamics, 4; Hibbeler, Engineering Mechanics: Statics, 8.8.
- Proof
- curve, order, invariant, limit.
rolling_resistance.trait
namespace lusine.integrator;
float<m/s> rolling_resistance(float<m/s> velocity, float friction, float<m/s2> gravity, float<s> dt)
{
let slowing = friction * gravity * dt;
if (abs(velocity) <= slowing) return velocity - velocity;
return velocity > 0.0 * velocity ? velocity - slowing : velocity + slowing;
}
A body rolling along the ground over a step, slowed by the deceleration μ_r g against its motion, v' = v − sign(v) μ_r g dt, stopped at zero rather than reversed. The coefficient is the ground's and the gravity the world's, both passed in.
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
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
from lusine.kernels._si import G0
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
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
from lusine.kernels._si import G0
The kernel¶
The law as the library states it.
In [2]:
Copied!
print(trait_run.source('integrator/rolling_resistance.trait'))
print(trait_run.source('integrator/rolling_resistance.trait'))
// lusine.integrator.rolling_resistance: a body rolling along the ground over a step, slowed by its rolling resistance,
// the deceleration mu_r g against its motion: v' = v - sign(v) mu_r g dt, never past zero, so a body slowed to rest
// stops rather than reversing. The coefficient is the ground's (the game's: the
// library states no rolling_friction), the gravity the world's magnitude, both passed in. A cut-off speed below which a game stops a body is
// the game's rule, after this kernel.
// Source: Gillespie, Fundamentals of Vehicle Dynamics, 4 (rolling resistance); Hibbeler, Engineering Mechanics:
// Statics, 8.8.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.integrator;
float<m/s> rolling_resistance(float<m/s> velocity, float friction, float<m/s2> gravity, float<s> dt)
{
let slowing = friction * gravity * dt;
if (abs(velocity) <= slowing) return velocity - velocity;
return velocity > 0.0 * velocity ? velocity - slowing : velocity + slowing;
}
The program¶
A lattice of balls, each with its own coefficient and step, so one run rolls every case of a section side by side; main() is one step of the kernel. The distance a ball runs is the page's sum of each step's new velocity times the step.
In [3]:
Copied!
ROLL = """namespace roll;
lusine.version version = { { number = 1 } };
struct ball { float<m/s> velocity; host float friction; host float<s> dt; };
lattice<ball, 1> balls(spacing = 1 m, boundary = wall);
host struct { float<m/s2> gravity; } settings;
void main() { balls.velocity = lusine.integrator.rolling_resistance(balls.velocity, balls.friction, settings.gravity,
balls.dt); }
"""
print(ROLL)
def roll(steps, v0, friction, dt, gravity):
"""Each ball's velocity after every step, one row a step."""
n = len(dt)
run = trait_run.Program(['integrator/rolling_resistance.trait'], ROLL, {'roll.balls': [n]})
run.store('roll.settings', 'gravity', gravity)
run.store('roll.balls', 'friction', np.broadcast_to(friction, (n,)))
run.store('roll.balls', 'dt', dt)
run.store('roll.balls', 'velocity', np.broadcast_to(v0, (n,)))
return run.steps('roll.main', steps, [('roll.balls', 'velocity')])[0]
ROLL = """namespace roll;
lusine.version version = { { number = 1 } };
struct ball { float<m/s> velocity; host float friction; host float<s> dt; };
lattice<ball, 1> balls(spacing = 1 m, boundary = wall);
host struct { float<m/s2> gravity; } settings;
void main() { balls.velocity = lusine.integrator.rolling_resistance(balls.velocity, balls.friction, settings.gravity,
balls.dt); }
"""
print(ROLL)
def roll(steps, v0, friction, dt, gravity):
"""Each ball's velocity after every step, one row a step."""
n = len(dt)
run = trait_run.Program(['integrator/rolling_resistance.trait'], ROLL, {'roll.balls': [n]})
run.store('roll.settings', 'gravity', gravity)
run.store('roll.balls', 'friction', np.broadcast_to(friction, (n,)))
run.store('roll.balls', 'dt', dt)
run.store('roll.balls', 'velocity', np.broadcast_to(v0, (n,)))
return run.steps('roll.main', steps, [('roll.balls', 'velocity')])[0]
namespace roll;
lusine.version version = { { number = 1 } };
struct ball { float<m/s> velocity; host float friction; host float<s> dt; };
lattice<ball, 1> balls(spacing = 1 m, boundary = wall);
host struct { float<m/s2> gravity; } settings;
void main() { balls.velocity = lusine.integrator.rolling_resistance(balls.velocity, balls.friction, settings.gravity,
balls.dt); }
Curve: it stops after v0 / (μ g), v0² / (2 μ g) along¶
Balls from 3 m/s on Earth's gravity, at three coefficients and dt = 10 ms: each speed falls in a line and stops within two steps of v0 / (μ g); the distance it runs is v0² / (2 μ g) less the summed steps' v0 dt / 2, within that again and a thousandth for binary32's rounding of a speed taken down in a thousand small subtractions.
In [4]:
Copied!
mu = np.array([0.02, 0.05, 0.1])
dt = 0.01
v = roll(1600, 3.0, mu, np.full(3, dt), G0) # v[k, i]: ball i after k + 1 steps
for i, m in enumerate(mu):
a = float(np.float32(m) * np.float32(G0))
stop = (int(np.argmax(v[:, i] == 0.0)) + 1) * dt
run = float(np.sum(v[:, i]) * dt)
assert v[-1, i] == 0.0 and abs(stop - 3.0 / a) <= 2 * dt, (stop, 3.0 / a)
assert abs(run - 9.0 / (2 * a)) <= 3.0 * dt / 2 + 1e-3 * run, (run, 9.0 / (2 * a))
t = np.arange(1, 1601) * dt
fig, ax = plt.subplots(figsize=(7, 3.2))
for i, m in enumerate(mu):
ax.plot(t, v[:, i], label=f'μ = {m}')
ax.set_xlabel('t (s)')
ax.set_ylabel('v (m/s)')
ax.legend()
plt.show()
mu = np.array([0.02, 0.05, 0.1])
dt = 0.01
v = roll(1600, 3.0, mu, np.full(3, dt), G0) # v[k, i]: ball i after k + 1 steps
for i, m in enumerate(mu):
a = float(np.float32(m) * np.float32(G0))
stop = (int(np.argmax(v[:, i] == 0.0)) + 1) * dt
run = float(np.sum(v[:, i]) * dt)
assert v[-1, i] == 0.0 and abs(stop - 3.0 / a) <= 2 * dt, (stop, 3.0 / a)
assert abs(run - 9.0 / (2 * a)) <= 3.0 * dt / 2 + 1e-3 * run, (run, 9.0 / (2 * a))
t = np.arange(1, 1601) * dt
fig, ax = plt.subplots(figsize=(7, 3.2))
for i, m in enumerate(mu):
ax.plot(t, v[:, i], label=f'μ = {m}')
ax.set_xlabel('t (s)')
ax.set_ylabel('v (m/s)')
ax.legend()
plt.show()
Order: the distance is first order in the step¶
On a world whose gravity is 8 m/s², with μ = 1/8, a ball from 2 m/s stops after 2 s and 2 m. Every value of this roll is exact in binary32, so the distance summed step by step falls short by exactly v0 dt / 2: halving the step halves the error.
In [5]:
Copied!
steps = np.array([16, 32, 64, 128])
v = roll(256, 2.0, 0.125, 1.0 / steps, 8.0)
errs = np.array([2.0 - np.sum(v[:2 * n, i]) / n for i, n in enumerate(steps)])
assert np.allclose(errs, 2.0 / steps / 2, rtol=0, atol=1e-12), errs
print('distance errors', errs, 'ratios', errs[:-1] / errs[1:])
steps = np.array([16, 32, 64, 128])
v = roll(256, 2.0, 0.125, 1.0 / steps, 8.0)
errs = np.array([2.0 - np.sum(v[:2 * n, i]) / n for i, n in enumerate(steps)])
assert np.allclose(errs, 2.0 / steps / 2, rtol=0, atol=1e-12), errs
print('distance errors', errs, 'ratios', errs[:-1] / errs[1:])
distance errors [0.0625 0.03125 0.015625 0.0078125] ratios [2. 2. 2.]
Invariant: it never reverses, either way¶
Two hundred balls at random speeds of both signs, coefficients and steps: after each step a ball moves the same way it did, or is at rest, never faster; a ball rolling the other way is the mirror image.
In [6]:
Copied!
rng = np.random.default_rng(300)
v0 = rng.normal(size=200) * 3
mu, dt = rng.uniform(0.0, 0.3, 200), rng.uniform(0.001, 0.1, 200)
run = trait_run.Program(['integrator/rolling_resistance.trait'], ROLL, {'roll.balls': [200]})
run.store('roll.settings', 'gravity', G0)
run.store('roll.balls', 'friction', mu)
run.store('roll.balls', 'dt', dt)
run.store('roll.balls', 'velocity', v0)
v = run.steps('roll.main', 100, [('roll.balls', 'velocity')])[0]
v = np.vstack([np.float32(v0).astype(float), v])
assert np.all(v[1:] * np.sign(v[:-1]) >= 0) and np.all(np.abs(v[1:]) <= np.abs(v[:-1]))
w = roll(100, -np.float32(v0).astype(float), mu, dt, G0)
assert np.array_equal(w, -v[1:])
print('no ball reversed; the mirror rolls the same')
rng = np.random.default_rng(300)
v0 = rng.normal(size=200) * 3
mu, dt = rng.uniform(0.0, 0.3, 200), rng.uniform(0.001, 0.1, 200)
run = trait_run.Program(['integrator/rolling_resistance.trait'], ROLL, {'roll.balls': [200]})
run.store('roll.settings', 'gravity', G0)
run.store('roll.balls', 'friction', mu)
run.store('roll.balls', 'dt', dt)
run.store('roll.balls', 'velocity', v0)
v = run.steps('roll.main', 100, [('roll.balls', 'velocity')])[0]
v = np.vstack([np.float32(v0).astype(float), v])
assert np.all(v[1:] * np.sign(v[:-1]) >= 0) and np.all(np.abs(v[1:]) <= np.abs(v[:-1]))
w = roll(100, -np.float32(v0).astype(float), mu, dt, G0)
assert np.array_equal(w, -v[1:])
print('no ball reversed; the mirror rolls the same')
no ball reversed; the mirror rolls the same
Limit: no resistance, and a slow ball¶
No resistance leaves the speed alone, bit for bit; a ball slower than one step's slowing stops at zero, either way, as does one exactly at it.
In [7]:
Copied!
law = lambda **kw: trait_run.law('integrator/rolling_resistance.trait', 'rolling_resistance', **kw)
assert np.array_equal(law(velocity=[1.5, -0.25], friction=0.0, gravity=G0, dt=0.1), [1.5, -0.25])
assert np.all(law(velocity=[0.01, -0.01, 0.5, -0.5], friction=[0.05, 0.05, 0.125, 0.125], gravity=[G0, G0, 8.0, 8.0],
dt=0.5) == 0.0)
print('unchanged at mu = 0; stopped below a step')
law = lambda **kw: trait_run.law('integrator/rolling_resistance.trait', 'rolling_resistance', **kw)
assert np.array_equal(law(velocity=[1.5, -0.25], friction=0.0, gravity=G0, dt=0.1), [1.5, -0.25])
assert np.all(law(velocity=[0.01, -0.01, 0.5, -0.5], friction=[0.05, 0.05, 0.125, 0.125], gravity=[G0, G0, 8.0, 8.0],
dt=0.5) == 0.0)
print('unchanged at mu = 0; stopped below a step')
unchanged at mu = 0; stopped below a step