The kick: velocity Verlet, second order, an orbit kept, and reversible¶

Lusine.integrator.kick and kick_drift: velocity Verlet's halves.

kick_drift takes the first half kick and the drift, kick the second half kick with the acceleration at the new position. Both are written once over a vector of any length, componentwise: a binding's first vector argument fixes the length, a single value the length 1, so one pair of halves steps a rod, a plane and a space. Checked against references in numpy, which have their own tests.

Source

Swope, Andersen, Berens & Wilson, J. Chem. Phys. 76, 637 (1982); Hairer, Lubich & Wanner, Geometric Numerical Integration, I.1.4.

Proof

order, invariant, limit.

kick.trait

namespace lusine.integrator;

vector<float<m/s>, _> kick(vector<float<m/s>, _> velocity, vector<float<m/s2>, _> acceleration, float<s> dt)
{
  return velocity + acceleration * (0.5 * dt);
}

(vector<float<m>, _>, vector<float<m/s>, _>) kick_drift(vector<float<m>, _> position, vector<float<m/s>, _> velocity,
                                                    vector<float<m/s2>, _> acceleration, float<s> dt)
{
  let half = velocity + acceleration * (0.5 * dt);
  return (position + half * dt, half);
}

Uses, without proving it here: integrator.kick_drift.

Velocity Verlet's closing half-kick, and with it the method, each step kick_drift, the force at the new position, then the kick. Both halves are written once over a vector of any length (vector<float<m/s>, _>): each binding takes the instance at its first vector argument's length, kick@1 the scalar form a lattice of single values binds (the spring and the fall below), kick@2 a plane's (the orbit), kick@3 a space's (the projectile); as standalone functions they are kick_1, kick_2 and kick_3. Values are binary32, as Trait's float is.

In [1]:

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¶

Both halves, as the library states them and this page runs them.

In [2]:

print(trait_run.source('integrator/kick.trait'))

print(trait_run.source('integrator/kick.trait'))

// lusine.integrator.kick and kick_drift: velocity Verlet's halves. kick_drift takes the first half kick and the drift,
// kick the second half kick with the acceleration at the new position. Both are written once over a vector of any
// length, componentwise: a binding's first vector argument fixes the length, a single value the
// length 1, so one pair of halves steps a rod, a plane and a space.
// Checked against references in numpy, which have their own tests.

namespace lusine.integrator;

vector<float<m/s>, _> kick(vector<float<m/s>, _> velocity, vector<float<m/s2>, _> acceleration, float<s> dt)
{
  return velocity + acceleration * (0.5 * dt);
}

(vector<float<m>, _>, vector<float<m/s>, _>) kick_drift(vector<float<m>, _> position, vector<float<m/s>, _> velocity,
                                                    vector<float<m/s2>, _> acceleration, float<s> dt)
{
  let half = velocity + acceleration * (0.5 * dt);
  return (position + half * dt, half);
}

Order: second¶

A unit spring (ω = 2π) from x = 1: halving the step quarters the error at t = 0.3 s, against cos(ω t).

In [3]:

SPRING = '''namespace spring;

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

struct body { float<m> position; float<m/s> velocity; float<m/s2> acceleration; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { float<s^-1> omega; float<s> dt; } settings;

float<m/s2> pull(float<m> position, float<s^-1> omega) { return -omega * omega * position; }

void start() { bodies.acceleration = pull(bodies.position, settings.omega); }

void main()
{
  (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity,
                                                                    bodies.acceleration, settings.dt);
  bodies.acceleration = pull(bodies.position, settings.omega);
  bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt);
}
'''

w = 2 * np.pi
steps = np.array([50, 100, 200, 400])
errs = []
for n in steps:
    run = trait_run.Program(['integrator/kick.trait'], SPRING, {'spring.bodies': [1]})
    run.store('spring.bodies', 'position', [1.0])
    run.store('spring.settings', 'omega', w)
    run.store('spring.settings', 'dt', 1.0 / n)
    run.step('spring.start')
    x = run.steps('spring.main', int(0.3 * n), [('spring.bodies', 'position')])[0].ravel()
    errs.append(abs(x[-1] - np.cos(w * 0.3)))
ratios = [errs[i] / errs[i + 1] for i in range(len(errs) - 1)]
print('error ratios', np.round(ratios, 3))
assert all(abs(r - 4) < 0.1 for r in ratios)
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.loglog(1 / steps, errs, 'o-', label='error at t = 0.3 s')
ax.loglog(1 / steps, errs[0] * (steps[0] / steps) ** 2, '--', label='slope 2')
ax.set_xlabel('dt (s)')
ax.legend()
plt.show()

SPRING = '''namespace spring; lusine.version version = { { number = 1 } }; struct body { float<m> position; float<m/s> velocity; float<m/s2> acceleration; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { float<s^-1> omega; float<s> dt; } settings; float<m/s2> pull(float<m> position, float<s^-1> omega) { return -omega * omega * position; } void start() { bodies.acceleration = pull(bodies.position, settings.omega); } void main() { (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity, bodies.acceleration, settings.dt); bodies.acceleration = pull(bodies.position, settings.omega); bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt); } ''' w = 2 * np.pi steps = np.array([50, 100, 200, 400]) errs = [] for n in steps: run = trait_run.Program(['integrator/kick.trait'], SPRING, {'spring.bodies': [1]}) run.store('spring.bodies', 'position', [1.0]) run.store('spring.settings', 'omega', w) run.store('spring.settings', 'dt', 1.0 / n) run.step('spring.start') x = run.steps('spring.main', int(0.3 * n), [('spring.bodies', 'position')])[0].ravel() errs.append(abs(x[-1] - np.cos(w * 0.3))) ratios = [errs[i] / errs[i + 1] for i in range(len(errs) - 1)] print('error ratios', np.round(ratios, 3)) assert all(abs(r - 4) < 0.1 for r in ratios) fig, ax = plt.subplots(figsize=(7, 3.5)) ax.loglog(1 / steps, errs, 'o-', label='error at t = 0.3 s') ax.loglog(1 / steps, errs[0] * (steps[0] / steps) ** 2, '--', label='slope 2') ax.set_xlabel('dt (s)') ax.legend() plt.show()

error ratios [4.001 3.995 3.986]

Invariant: an orbit's energy, and the way back¶

A circular Kepler orbit (GM = 1, 200 steps a period), its position, velocity and acceleration 2-vectors, so one binding of each half steps the plane (kick_drift@2, kick@2): its energy stays within binary32's reach of -1/2 for a hundred periods, with no drift; a run turned round comes back to its start within a rounding a step.

In [4]:

ORBIT = '''namespace orbit;

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

struct body { vector<float<m>, 2> position; vector<float<m/s>, 2> velocity; vector<float<m/s2>, 2> acceleration; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { float<s> dt; } settings;

vector<float<m/s2>, 2> kepler(vector<float<m>, 2> position)
{
  let r = sqrt(position.x * position.x + position.y * position.y);
  let k = -1.0 m3/s2 / (r * r * r);
  return k * position;
}

void start() { bodies.acceleration = kepler(bodies.position); }

void main()
{
  (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity,
                                                                    bodies.acceleration, settings.dt);
  bodies.acceleration = kepler(bodies.position);
  bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt);
}
'''

def orbit():
    run = trait_run.Program(['integrator/kick.trait'], ORBIT, {'orbit.bodies': [1]})
    run.store('orbit.bodies', 'position', [[1.0, 0.0]])
    run.store('orbit.bodies', 'velocity', [[0.0, 1.0]])
    run.store('orbit.settings', 'dt', 2 * np.pi / 200)
    run.step('orbit.start')
    return run

p, v = orbit().steps('orbit.main', 200 * 100, [('orbit.bodies', 'position'), ('orbit.bodies', 'velocity')])
(x, y), (vx, vy) = p[:, 0].T, v[:, 0].T
E = 0.5 * (vx * vx + vy * vy) - 1 / np.sqrt(x * x + y * y)
worst = np.abs(E + 0.5).reshape(10, -1).max(axis=1)   # worst in each ten periods
print('worst |E - E0| in each ten periods', worst)
assert worst.max() < 1e-5 and worst[-1] < 3 * worst[0]
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.semilogy(np.arange(10, 101, 10), worst, 'o-', label='velocity Verlet, binary32')
ax.set_xlabel('t (periods)')
ax.set_ylabel('|E - E0|, worst in ten periods')
ax.legend()
plt.show()

run = orbit()
v = run.steps('orbit.main', 1000, [('orbit.bodies', 'velocity')])[0]
run.store('orbit.bodies', 'velocity', [-v[-1][0]])
p, v = (o[-1][0] for o in run.steps('orbit.main', 1000, [('orbit.bodies', 'position'), ('orbit.bodies', 'velocity')]))
back = max(abs(p[0] - 1), abs(p[1]), abs(v[0]), abs(v[1] + 1))
assert back <= 2000 * EPS
print('back within', int(back / EPS), 'binary32 ulps after 2000 steps')

ORBIT = '''namespace orbit; lusine.version version = { { number = 1 } }; struct body { vector<float<m>, 2> position; vector<float<m/s>, 2> velocity; vector<float<m/s2>, 2> acceleration; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { float<s> dt; } settings; vector<float<m/s2>, 2> kepler(vector<float<m>, 2> position) { let r = sqrt(position.x * position.x + position.y * position.y); let k = -1.0 m3/s2 / (r * r * r); return k * position; } void start() { bodies.acceleration = kepler(bodies.position); } void main() { (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity, bodies.acceleration, settings.dt); bodies.acceleration = kepler(bodies.position); bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt); } ''' def orbit(): run = trait_run.Program(['integrator/kick.trait'], ORBIT, {'orbit.bodies': [1]}) run.store('orbit.bodies', 'position', [[1.0, 0.0]]) run.store('orbit.bodies', 'velocity', [[0.0, 1.0]]) run.store('orbit.settings', 'dt', 2 * np.pi / 200) run.step('orbit.start') return run p, v = orbit().steps('orbit.main', 200 * 100, [('orbit.bodies', 'position'), ('orbit.bodies', 'velocity')]) (x, y), (vx, vy) = p[:, 0].T, v[:, 0].T E = 0.5 * (vx * vx + vy * vy) - 1 / np.sqrt(x * x + y * y) worst = np.abs(E + 0.5).reshape(10, -1).max(axis=1) # worst in each ten periods print('worst |E - E0| in each ten periods', worst) assert worst.max() < 1e-5 and worst[-1] < 3 * worst[0] fig, ax = plt.subplots(figsize=(7, 3.5)) ax.semilogy(np.arange(10, 101, 10), worst, 'o-', label='velocity Verlet, binary32') ax.set_xlabel('t (periods)') ax.set_ylabel('|E - E0|, worst in ten periods') ax.legend() plt.show() run = orbit() v = run.steps('orbit.main', 1000, [('orbit.bodies', 'velocity')])[0] run.store('orbit.bodies', 'velocity', [-v[-1][0]]) p, v = (o[-1][0] for o in run.steps('orbit.main', 1000, [('orbit.bodies', 'position'), ('orbit.bodies', 'velocity')])) back = max(abs(p[0] - 1), abs(p[1]), abs(v[0]), abs(v[1] + 1)) assert back <= 2000 * EPS print('back within', int(back / EPS), 'binary32 ulps after 2000 steps')

worst |E - E0| in each ten periods [3.30726012e-06 2.90975933e-06 2.58169480e-06 2.75876771e-06
 1.91461810e-06 2.11137837e-06 2.71893364e-06 4.74049869e-06
 5.37441592e-06 6.79962426e-06]

back within 1120 binary32 ulps after 2000 steps

Limit: a fall is exact¶

Under a constant force the two half-kicks integrate exactly: the fall lands on -g t²/2 but for rounding.

In [5]:

FALL = '''namespace fall;

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

struct body { float<m> position; float<m/s> velocity; host float<m/s2> acceleration; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { float<s> dt; } settings;

void main()
{
  (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity,
                                                                    bodies.acceleration, settings.dt);
  bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt);
}
'''

run = trait_run.Program(['integrator/kick.trait'], FALL, {'fall.bodies': [1]})
run.store('fall.bodies', 'acceleration', [-G0])
run.store('fall.settings', 'dt', 0.01)
x = run.steps('fall.main', 100, [('fall.bodies', 'position')])[0][-1][0]
assert abs(x + 0.5 * G0) < 1e-5
print('fall error', x + 0.5 * G0)

FALL = '''namespace fall; lusine.version version = { { number = 1 } }; struct body { float<m> position; float<m/s> velocity; host float<m/s2> acceleration; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { float<s> dt; } settings; void main() { (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity, bodies.acceleration, settings.dt); bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt); } ''' run = trait_run.Program(['integrator/kick.trait'], FALL, {'fall.bodies': [1]}) run.store('fall.bodies', 'acceleration', [-G0]) run.store('fall.settings', 'dt', 0.01) x = run.steps('fall.main', 100, [('fall.bodies', 'position')])[0][-1][0] assert abs(x + 0.5 * G0) < 1e-5 print('fall error', x + 0.5 * G0)

fall error 4.68749999971152e-06

Limit: a projectile in three dimensions¶

A body in space under a constant gravity, its position, velocity and acceleration 3-vectors, so the program binds the halves' instances at 3 (kick_drift@3, kick@3): 64 steps of 1/64 s land where v t + g t² / 2 says, but for rounding, and each component is the scalar form's run on that component alone, bit for bit.

In [6]:

SPACE = FALL.replace('namespace fall;', 'namespace space;').replace(
    'struct body { float<m> position; float<m/s> velocity; host float<m/s2> acceleration; };',
    'struct body { vector<float<m>, 3> position; vector<float<m/s>, 3> velocity; host vector<float<m/s2>, 3> acceleration; };')
run = trait_run.Program(['integrator/kick.trait'], SPACE, {'space.bodies': [1]})
run.store('space.bodies', 'velocity', [[3.0, 1.0, 4.0]])
run.store('space.bodies', 'acceleration', [[0.0, 0.0, -9.81]])
run.store('space.settings', 'dt', 1 / 64)
x = run.steps('space.main', 64, [('space.bodies', 'position')])[0][:, 0]
assert np.allclose(x[-1], np.array([3.0, 1.0, 4.0 - 9.81 / 2]), rtol=0, atol=1e-5)
run = trait_run.Program(['integrator/kick.trait'], FALL, {'fall.bodies': [3]})
run.store('fall.bodies', 'velocity', [3.0, 1.0, 4.0])
run.store('fall.bodies', 'acceleration', [0.0, 0.0, -9.81])
run.store('fall.settings', 'dt', 1 / 64)
assert np.array_equal(run.steps('fall.main', 64, [('fall.bodies', 'position')])[0], x)
print('landed at', x[-1], '; each component the scalar form\'s, bit for bit, every step')

SPACE = FALL.replace('namespace fall;', 'namespace space;').replace( 'struct body { float<m> position; float<m/s> velocity; host float<m/s2> acceleration; };', 'struct body { vector<float<m>, 3> position; vector<float<m/s>, 3> velocity; host vector<float<m/s2>, 3> acceleration; };') run = trait_run.Program(['integrator/kick.trait'], SPACE, {'space.bodies': [1]}) run.store('space.bodies', 'velocity', [[3.0, 1.0, 4.0]]) run.store('space.bodies', 'acceleration', [[0.0, 0.0, -9.81]]) run.store('space.settings', 'dt', 1 / 64) x = run.steps('space.main', 64, [('space.bodies', 'position')])[0][:, 0] assert np.allclose(x[-1], np.array([3.0, 1.0, 4.0 - 9.81 / 2]), rtol=0, atol=1e-5) run = trait_run.Program(['integrator/kick.trait'], FALL, {'fall.bodies': [3]}) run.store('fall.bodies', 'velocity', [3.0, 1.0, 4.0]) run.store('fall.bodies', 'acceleration', [0.0, 0.0, -9.81]) run.store('fall.settings', 'dt', 1 / 64) assert np.array_equal(run.steps('fall.main', 64, [('fall.bodies', 'position')])[0], x) print('landed at', x[-1], '; each component the scalar form\'s, bit for bit, every step')

landed at [ 3.          1.         -0.90499908] ; each component the scalar form's, bit for bit, every step