Semi-implicit Euler: first order, and the energy kept¶

The semi-implicit (symplectic) Euler step, the velocity first and then the position with the new velocity: v' = v + a(x) dt, x' = x + v' dt.

First order, one force evaluation a step, and a conservative system's energy kept within a bound of its start rather than drifting (the explicit Euler, with the old velocity, gains energy every step). The acceleration at the position is the game's force kernels', bound before it. 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 law steps a rod, a plane and a space. Source: Hairer, Lubich & Wanner, Geometric Numerical Integration (2nd ed.), I.1.2 and VI.3. Checked against a reference in numpy, which has its own tests.

Source

Hairer, Lubich & Wanner, Geometric Numerical Integration (2nd ed.), I.1.2 and VI.3.

Proof

order, invariant, limit.

symplectic_euler.trait

namespace lusine.integrator;

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

The velocity first, v' = v + a(x) dt, then the position with the new velocity, x' = x + v' dt, here on a harmonic oscillator of period 1 s. It is written once over a vector of any length: a binding over single values takes symplectic_euler@1, the scalar form, and one over a plane's symplectic_euler@2 (as standalone functions, symplectic_euler_1 to symplectic_euler_3).

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¶

The law as the library states it, and as this page runs it.

In [2]:

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

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

// lusine.integrator.symplectic_euler: the semi-implicit (symplectic) Euler step, the velocity first and then the
// position with the new velocity: v' = v + a(x) dt, x' = x + v' dt. First order, one force evaluation a step, and a
// conservative system's energy kept within a bound of its start rather than drifting (the explicit Euler, with the old
// velocity, gains energy every step). The acceleration at the position is the game's force kernels', bound before it.
// 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 law steps a rod, a plane and a space.
// Source: Hairer, Lubich & Wanner, Geometric Numerical Integration (2nd ed.), I.1.2 and VI.3.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.integrator;

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

The program¶

A lattice of oscillators, each with its own step, so one run steps every case of a section side by side. The spring's pull, -omega^2 x, is the program's own kernel, bound before the step as a game's force kernels are; main() is one step.

In [3]:

SPRING = '''namespace spring;

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

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

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

void main()
{
  bodies.acceleration = pull(bodies.position, settings.omega);
  (bodies.position, bodies.velocity) = lusine.integrator.symplectic_euler(bodies.position, bodies.velocity,
                                                                          bodies.acceleration, bodies.dt);
}
'''
print(SPRING)
w = 2 * np.pi

def spring(steps, dt, position, velocity):
    """Each oscillator's position and velocity after every step, one row a step."""
    n = len(dt)
    run = trait_run.Program(['integrator/symplectic_euler.trait'], SPRING, {'spring.bodies': [n]})
    run.store('spring.settings', 'omega', w)
    run.store('spring.bodies', 'dt', dt)
    run.store('spring.bodies', 'position', position)
    run.store('spring.bodies', 'velocity', velocity)
    return run.steps('spring.main', steps, [('spring.bodies', 'position'), ('spring.bodies', 'velocity')])

SPRING = '''namespace spring; lusine.version version = { { number = 1 } }; struct body { float<m> position; float<m/s> velocity; float<m/s2> acceleration; host float<s> dt; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { float<s^-1> omega; } settings; float<m/s2> pull(float<m> position, float<s^-1> omega) { return -omega * omega * position; } void main() { bodies.acceleration = pull(bodies.position, settings.omega); (bodies.position, bodies.velocity) = lusine.integrator.symplectic_euler(bodies.position, bodies.velocity, bodies.acceleration, bodies.dt); } ''' print(SPRING) w = 2 * np.pi def spring(steps, dt, position, velocity): """Each oscillator's position and velocity after every step, one row a step.""" n = len(dt) run = trait_run.Program(['integrator/symplectic_euler.trait'], SPRING, {'spring.bodies': [n]}) run.store('spring.settings', 'omega', w) run.store('spring.bodies', 'dt', dt) run.store('spring.bodies', 'position', position) run.store('spring.bodies', 'velocity', velocity) return run.steps('spring.main', steps, [('spring.bodies', 'position'), ('spring.bodies', 'velocity')])

namespace spring;

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

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

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

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

Order: first¶

Halving the step halves the error at t = 0.3 s against cos(ω t): the first order. The five step sizes are five oscillators of one run, each read at its own step count.

In [4]:

steps = np.array([50, 100, 200, 400, 800])
x, _ = spring(240, 1.0 / steps, np.ones(5), np.zeros(5))
errs = np.array([abs(x[int(0.3 * n) - 1, i] - np.cos(w * 0.3)) 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 = 0.3 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([50, 100, 200, 400, 800]) x, _ = spring(240, 1.0 / steps, np.ones(5), np.zeros(5)) errs = np.array([abs(x[int(0.3 * n) - 1, i] - np.cos(w * 0.3)) 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 = 0.3 s') ax.loglog(1 / steps, errs[0] * (steps[0] / steps), '--', label='slope 1') ax.set_xlabel('dt (s)') ax.legend() plt.show()

error ratios [2.022 2.01  2.005 2.003]

Invariant: the energy stays in its band¶

Over a hundred periods at a hundred steps a period, the energy's worst deviation in the last ten periods is the first ten's: bounded, no drift. The explicit Euler (the position with the old velocity, the step this kernel is not), written here in numpy for contrast, gains energy every step.

In [5]:

dt = 0.01
x, v = (a[:, 0] for a in spring(100 * 100, np.array([dt]), [1.0], [0.0]))
energy = lambda x, v: (0.5 * v * v + 0.5 * w * w * x * x) / (0.5 * w * w)
worst = np.abs(energy(x, v) - 1).reshape(10, -1).max(axis=1)      # worst in each ten periods
print('worst |E/E0 - 1| in each ten periods', np.round(worst, 5))
assert worst[0] < w * dt and abs(worst[-1] / worst[0] - 1) < 0.01
ex, ev, explicit = 1.0, 0.0, []
for _ in range(3000):
    ex, ev = ex + ev * dt, ev - w * w * ex * dt
    explicit.append(abs(energy(ex, ev) - 1))
assert explicit[-1] > 1
t = dt * np.arange(1, len(x) + 1)
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.semilogy(t, np.abs(energy(x, v) - 1) + 1e-9, lw=0.5, label='semi-implicit Euler')
ax.semilogy(t[:3000], explicit, label='explicit Euler (numpy, for contrast)')
ax.set_xlabel('t (periods)')
ax.set_ylabel('|E/E0 - 1|')
ax.legend()
plt.show()

dt = 0.01 x, v = (a[:, 0] for a in spring(100 * 100, np.array([dt]), [1.0], [0.0])) energy = lambda x, v: (0.5 * v * v + 0.5 * w * w * x * x) / (0.5 * w * w) worst = np.abs(energy(x, v) - 1).reshape(10, -1).max(axis=1) # worst in each ten periods print('worst |E/E0 - 1| in each ten periods', np.round(worst, 5)) assert worst[0] < w * dt and abs(worst[-1] / worst[0] - 1) < 0.01 ex, ev, explicit = 1.0, 0.0, [] for _ in range(3000): ex, ev = ex + ev * dt, ev - w * w * ex * dt explicit.append(abs(energy(ex, ev) - 1)) assert explicit[-1] > 1 t = dt * np.arange(1, len(x) + 1) fig, ax = plt.subplots(figsize=(7, 3.5)) ax.semilogy(t, np.abs(energy(x, v) - 1) + 1e-9, lw=0.5, label='semi-implicit Euler') ax.semilogy(t[:3000], explicit, label='explicit Euler (numpy, for contrast)') ax.set_xlabel('t (periods)') ax.set_ylabel('|E/E0 - 1|') ax.legend() plt.show()

worst |E/E0 - 1| in each ten periods [0.03243 0.03243 0.03243 0.03241 0.03239 0.03241 0.03242 0.03243 0.03243
 0.03241]

Invariant: the mirror¶

The spring's pull is odd, so a mirrored state steps to the mirrored state, bit for bit, step after step.

In [6]:

x, v = spring(50, np.full(4, 0.013), [0.3, -0.3, -0.7, 0.7], [1.1, -1.1, 0.2, -0.2])
assert np.array_equal(x[:, 0], -x[:, 1]) and np.array_equal(v[:, 0], -v[:, 1])
assert np.array_equal(x[:, 2], -x[:, 3]) and np.array_equal(v[:, 2], -v[:, 3])
print('mirrored exactly for 50 steps')

x, v = spring(50, np.full(4, 0.013), [0.3, -0.3, -0.7, 0.7], [1.1, -1.1, 0.2, -0.2]) assert np.array_equal(x[:, 0], -x[:, 1]) and np.array_equal(v[:, 0], -v[:, 1]) assert np.array_equal(x[:, 2], -x[:, 3]) and np.array_equal(v[:, 2], -v[:, 3]) print('mirrored exactly for 50 steps')

mirrored exactly for 50 steps

Limit: a fall¶

Under a constant force the step is off -g t²/2 by g t dt / 2, the first-order term: after a second of 0.01 s steps, x = -g (1 + 0.01) / 2, to binary32's rounding of a hundred steps.

In [7]:

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.symplectic_euler(bodies.position, bodies.velocity,
                                                                          bodies.acceleration, settings.dt);
}
'''
run = trait_run.Program(['integrator/symplectic_euler.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][:, 0]
assert abs(x[-1] + 0.5 * G0 * 1.01) < 1e-5
assert abs((x[-1] + 0.5 * G0) / (x[49] + 0.5 * G0 * 0.25) - 2) < 0.01     # the term g t dt / 2: twice at t = 1 s
print('fall error', x[-1] + 0.5 * G0, 'against -g t dt / 2 =', -0.5 * G0 * 0.01)

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.symplectic_euler(bodies.position, bodies.velocity, bodies.acceleration, settings.dt); } ''' run = trait_run.Program(['integrator/symplectic_euler.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][:, 0] assert abs(x[-1] + 0.5 * G0 * 1.01) < 1e-5 assert abs((x[-1] + 0.5 * G0) / (x[49] + 0.5 * G0 * 0.25) - 2) < 0.01 # the term g t dt / 2: twice at t = 1 s print('fall error', x[-1] + 0.5 * G0, 'against -g t dt / 2 =', -0.5 * G0 * 0.01)

fall error -0.04903086166381865 against -g t dt / 2 = -0.04903325

Invariant: an oscillator in a plane¶

An isotropic oscillator (the pull −ω² x on a 2-vector position, so the program binds symplectic_euler@2) started on a circle of radius 1: over fifty periods its radius stays within the scalar oscillator's energy band, and each axis is the scalar oscillator started on that axis's component, bit for bit, step after step.

In [8]:

PLANE = '''namespace plane;

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

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

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

void main()
{
  bodies.acceleration = pull(bodies.position, settings.omega);
  (bodies.position, bodies.velocity) = lusine.integrator.symplectic_euler(bodies.position, bodies.velocity,
                                                                          bodies.acceleration, bodies.dt);
}
'''
run = trait_run.Program(['integrator/symplectic_euler.trait'], PLANE, {'plane.bodies': [1]})
run.store('plane.settings', 'omega', w)
run.store('plane.bodies', 'dt', [0.01])
run.store('plane.bodies', 'position', [[1.0, 0.0]])
run.store('plane.bodies', 'velocity', [[0.0, w]])
p, v = (a[:, 0] for a in run.steps('plane.main', 5000, [('plane.bodies', 'position'), ('plane.bodies', 'velocity')]))
radius = np.sqrt(p[:, 0] ** 2 + p[:, 1] ** 2)
assert np.max(np.abs(radius - 1)) < w * 0.01
x, vx = spring(5000, np.array([0.01, 0.01]), [1.0, 0.0], [0.0, w])
assert np.array_equal(p[:, 0], x[:, 0]) and np.array_equal(v[:, 0], vx[:, 0])
assert np.array_equal(p[:, 1], x[:, 1]) and np.array_equal(v[:, 1], vx[:, 1])
print(f'radius within {np.max(np.abs(radius - 1)):.4f} of 1 over 50 periods; each axis the scalar oscillator, bit for bit')
fig, ax = plt.subplots(figsize=(4, 4))
ax.plot(p[:, 0], p[:, 1], lw=0.3)
ax.set_aspect('equal')
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
plt.show()

PLANE = '''namespace plane; lusine.version version = { { number = 1 } }; struct body { vector<float<m>, 2> position; vector<float<m/s>, 2> velocity; vector<float<m/s2>, 2> acceleration; host float<s> dt; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { float<s^-1> omega; } settings; vector<float<m/s2>, 2> pull(vector<float<m>, 2> position, float<s^-1> omega) { return -omega * omega * position; } void main() { bodies.acceleration = pull(bodies.position, settings.omega); (bodies.position, bodies.velocity) = lusine.integrator.symplectic_euler(bodies.position, bodies.velocity, bodies.acceleration, bodies.dt); } ''' run = trait_run.Program(['integrator/symplectic_euler.trait'], PLANE, {'plane.bodies': [1]}) run.store('plane.settings', 'omega', w) run.store('plane.bodies', 'dt', [0.01]) run.store('plane.bodies', 'position', [[1.0, 0.0]]) run.store('plane.bodies', 'velocity', [[0.0, w]]) p, v = (a[:, 0] for a in run.steps('plane.main', 5000, [('plane.bodies', 'position'), ('plane.bodies', 'velocity')])) radius = np.sqrt(p[:, 0] ** 2 + p[:, 1] ** 2) assert np.max(np.abs(radius - 1)) < w * 0.01 x, vx = spring(5000, np.array([0.01, 0.01]), [1.0, 0.0], [0.0, w]) assert np.array_equal(p[:, 0], x[:, 0]) and np.array_equal(v[:, 0], vx[:, 0]) assert np.array_equal(p[:, 1], x[:, 1]) and np.array_equal(v[:, 1], vx[:, 1]) print(f'radius within {np.max(np.abs(radius - 1)):.4f} of 1 over 50 periods; each axis the scalar oscillator, bit for bit') fig, ax = plt.subplots(figsize=(4, 4)) ax.plot(p[:, 0], p[:, 1], lw=0.3) ax.set_aspect('equal') ax.set_xlabel('x (m)') ax.set_ylabel('y (m)') plt.show()

radius within 0.0161 of 1 over 50 periods; each axis the scalar oscillator, bit for bit