Semi-implicit Euler

Deep dive

How it works

In: positions x, velocities v, force function F(x, v), masses, fixed step h · Out: x and v after one step · State: x and v themselves (no history) · Cost: one force evaluation plus 4 multiply-adds per degree of freedom

1. Accumulate forces F(x, v) for the current state (gravity, springs, user forces).
2. Velocity first: v ← v + h·F/m. Apply damping: v ← v / (1 + h·c).
3. (Constraint solvers run here, on velocities.)
4. Then position: x ← x + h·v, with the new v.
5. Fixed timestep: real frame time feeds an accumulator; run as many fixed steps as fit, render with α = accumulator/h interpolation between the previous and current state.

For the oscillator x'' = −ω²x, one step is the matrix [[1 − ω²h², h], [−ω²h, 1]] acting on (x, v): determinant 1 (area preserving, the symplectic property), trace 2 − ω²h². Eigenvalues stay on the unit circle iff |trace| < 2 ⇔ ωh < 2. Explicit Euler's matrix has determinant 1 + ω²h² > 1: every step inflates phase-space area, i.e. energy.

real time ──|──frame──|────frame────|─frame─|
sim ticks   |h|h|h|h|h|h|h|h|h|h|h|h|h|h|
render      interpolate between the last two ticks with α = acc / h

Code

#[derive(Clone, Copy, Default)]
pub struct Body { pub x: [f32; 2], pub v: [f32; 2], pub inv_mass: f32, pub damping: f32 }

pub struct Spring { pub a: usize, pub b: usize, pub rest: f32, pub k: f32, pub c: f32 }

pub struct World {
    pub bodies: Vec<Body>, pub springs: Vec<Spring>, pub gravity: [f32; 2], pub h: f32,
    prev_x: Vec<[f32; 2]>, force: Vec<[f32; 2]>, acc: f32,
}

impl World {
    pub fn new(bodies: Vec<Body>, springs: Vec<Spring>, gravity: [f32; 2], h: f32) -> Self {
        let n = bodies.len();
        World { prev_x: bodies.iter().map(|b| b.x).collect(), force: vec![[0.0; 2]; n], bodies, springs, gravity, h, acc: 0.0 }
    }

    fn compute_forces(&mut self) {
        for (f, b) in self.force.iter_mut().zip(&self.bodies) {
            let m = if b.inv_mass > 0.0 { 1.0 / b.inv_mass } else { 0.0 };
            *f = [self.gravity[0] * m, self.gravity[1] * m];
        }
        for s in &self.springs {
            let (pa, pb) = (self.bodies[s.a], self.bodies[s.b]);
            let d = [pb.x[0] - pa.x[0], pb.x[1] - pa.x[1]];
            let len = (d[0] * d[0] + d[1] * d[1]).sqrt().max(1e-6);
            let n = [d[0] / len, d[1] / len];
            let rel_v = (pb.v[0] - pa.v[0]) * n[0] + (pb.v[1] - pa.v[1]) * n[1];
            let mag = s.k * (len - s.rest) + s.c * rel_v; // Hooke + damping along the spring
            for k in 0..2 {
                self.force[s.a][k] += mag * n[k];
                self.force[s.b][k] -= mag * n[k];
            }
        }
    }

    pub fn step(&mut self) {
        let h = self.h;
        self.compute_forces();
        for (b, f) in self.bodies.iter_mut().zip(&self.force) {
            if b.inv_mass == 0.0 { continue; } // static body
            for k in 0..2 { b.v[k] += h * f[k] * b.inv_mass; }      // 1. velocity
            let damp = 1.0 / (1.0 + h * b.damping);                  // Pade damping, stable for any h
            for k in 0..2 { b.v[k] *= damp; b.x[k] += h * b.v[k]; } // 2. position with new velocity
        }
    }

    /// Consume frame time in fixed steps; returns alpha for render interpolation.
    pub fn advance(&mut self, frame_dt: f32) -> f32 {
        self.acc += frame_dt.min(0.25); // clamp: avoid the "spiral of death" after a stall
        while self.acc >= self.h {
            for (p, b) in self.prev_x.iter_mut().zip(&self.bodies) { *p = b.x; }
            self.step();
            self.acc -= self.h;
        }
        self.acc / self.h
    }

    pub fn render_pos(&self, i: usize, alpha: f32) -> [f32; 2] {
        let (p, x) = (self.prev_x[i], self.bodies[i].x);
        [p[0] + (x[0] - p[0]) * alpha, p[1] + (x[1] - p[1]) * alpha]
    }

    /// Largest stable step for the stiffest undamped spring: h < 2 / omega.
    pub fn max_stable_h(&self) -> f32 {
        self.springs.iter().map(|s| {
            let w2 = s.k * (self.bodies[s.a].inv_mass + self.bodies[s.b].inv_mass);
            if w2 > 0.0 { 2.0 / w2.sqrt() } else { f32::INFINITY }
        }).fold(f32::INFINITY, f32::min)
    }
}

Complexity and constants

• Same cost as explicit Euler: one force pass plus 4 multiply-adds per degree of freedom.
• Measured with the code above, unit spring (ω = 1), 10⁶ steps: at ωh = 0.1 the energy stays within ±5 % of its start (the O(h) error band, about ωh/2); at ωh = 1.9 it stays bounded but swings widely; at ωh = 2.1 it overflows f32 after about 60 steps.
• Typical game ticks: 60 Hz (h ≈ 16.7 ms) or 120-240 Hz substeps. At 60 Hz, a spring is stable up to ω < 120 rad/s, about 19 Hz natural frequency; stiffer ropes or suspensions need substeps or implicit springs.
• f32 is standard in engines; positions far from the origin (> order of 10⁴-10⁵ units) lose sub-millimetre precision (f32 has 24 significand bits), a common source of jitter.

Efficiency

• Cost: same as explicit Euler: one force pass plus 4 multiply-adds per degree of freedom, no stored history. Measured with the code above (unit spring, 10⁶ steps): at ωh = 0.1 the energy stays within ±5 % of its start; at ωh = 2.1 it overflows f32 after about 60 steps. Typical ticks are 60 Hz (h ≈ 16.7 ms) to 240 Hz substeps, stable up to ω < 120 rad/s at 60 Hz. Bound: the force evaluation and constraint solving, not the integrator.
• Levers:
  • fixed tick with an accumulator and render interpolation — deterministic, frame-rate independent (Fiedler) — all real-time loops
  • substeps k — stability limit rises k-fold, cost grows k-fold (page above) — stiff content
  • sleeping bodies — skips resting objects, a big win in large scenes (page above) — scenes with many idle bodies
  • implicit damping v ← v/(1 + h·c) — stable for any h, same cost — damped systems
  • f32 with positions near the origin — keeps sub-millimetre precision (24 significand bits, page above) — large worlds
• Naive → standard → state of the art:
  • explicit Euler (position first) · energy grows every step
  • semi-implicit Euler, fixed tick · bounded energy error O(h), same cost
  • Verlet or XPBD with substeps, implicit springs for stiff content · second order, or stable at large ωh

Implementation notes

• Order inside a step (Box2D style): integrate velocities (forces, gravity, damping) → solve velocity constraints (contacts, joints) → integrate positions → position correction.
• Damping: v *= 1/(1 + h·c) is the Padé form of e^(−hc); v *= 1 − h·c reverses velocity when h·c > 1 and explodes when h·c > 2.
• Fixed tick + interpolation decouples simulation from frame rate, gives determinism, and makes client-side prediction and rollback replay identical steps.
• Substepping: for stiff content, run k substeps of h/k; stability grows linearly with k, cost too.
• Sleeping: skip bodies whose velocity stays below a threshold for a while; big win in large scenes.
• Determinism: same binary, same inputs, same order of force accumulation; avoid f32::mul_add on some targets and plain multiply-add on others, and avoid transcendental functions whose implementations differ by platform.

Pitfalls

• Updating position before velocity: that is explicit Euler, springs gain energy and explode.
• Variable dt = frame time: frame-rate-dependent behaviour, non-determinism, instability on a slow frame.
• No clamp on accumulated time: after a stall the loop tries to catch up with hundreds of steps and stalls again.
• Rendering the raw simulation state without interpolation: visible stutter when tick and display rates differ.
• Forgetting that two light bodies on one spring have a higher ω than one body against a wall (reduced mass).

Variants and alternatives

In production

Games: Box2D (b2Island::Solve: integrate velocities with Padé damping, solve constraints, integrate positions) and Jolt Physics, which uses a symplectic Euler integrator. Robotics research: MuJoCo, whose default Euler integrator is the semi-implicit Euler method. Bullet and most game physics engines follow the same velocity-then-position order; many particle systems in game engines.

Implementations

Sources

• Ernst Hairer, Christian Lubich, Gerhard Wanner, Geometric Numerical Integration, 2nd ed., Springer, 2006. book
• Erin Catto, Iterative Dynamics with Temporal Coherence, Game Developers Conference, 2005. pdf
• Glenn Fiedler, Fix Your Timestep!, gafferongames.com, 2004. post
• Box2D source (v2.4), b2_island.cpp. code
• David Baraff, Andrew Witkin, Large Steps in Cloth Simulation, SIGGRAPH 1998. paper
• MuJoCo documentation, Computation: numerical integration (Euler, RK4, implicit integrators). docs