Sequential impulse solver

Deep dive

How it works

Parts:

• contact constraint — anchors rA, rB, normal n, tangent t, effective masses, bias
• accumulated impulses Pn, Pt — warm-start state, kept across steps by feature id
• solver bodies — compact v, ω, inverse mass and inertia of one island
• iteration loop — k sweeps over all contacts, friction before normal
• position correction — Baumgarte bias, slop, restitution threshold

Operations: prepare, warm start, iterate, integrate, store (steps below).

Per contact point between bodies A and B with normal n (A → B), offsets rA, rB from the centers of mass:

1. Relative velocity: dv = vB + ωB × rB − vA − ωA × rA; normal speed vn = dv·n.
2. Effective mass along a direction d: m_eff = 1 / (1/mA + 1/mB + (rA×d)ᵀ IA⁻¹ (rA×d) + (rB×d)ᵀ IB⁻¹ (rB×d)). In 2D the inertia terms are (rA×d)²/IA. Precompute once per step.
3. Bias (target normal velocity): speculative gap s > 0 → −s/dt (may approach just enough to touch); penetration → β/dt · max(−s − slop, 0); restitution → max(bias, −e·vn) when approaching faster than the threshold.
4. Warm start: apply last step's accumulated impulses Pn·n + Pt·t.
5. Iterate k times over all contacts:
   • friction: λ = −m_t · (dv·t); Pt_new = clamp(Pt + λ, −μPn, μPn); apply Pt_new − Pt.
   • normal: λ = m_n · (bias − dv·n); Pn_new = max(Pn + λ, 0); apply Pn_new − Pn.
6. Integrate velocities into positions (semi-implicit Euler); store Pn, Pt keyed by contact feature id for the next step.

Applying impulse P at the contact: vA −= P/mA, ωA −= IA⁻¹(rA × P), vB += P/mB, ωB += IB⁻¹(rB × P).

step: gravity → broadphase/narrowphase → prepare (m_eff, bias)
      → warm start → k × [friction, normal per contact] → integrate → store impulses

Code

#[derive(Clone, Copy, Default, Debug)]
pub struct V2 { pub x: f32, pub y: f32 }
fn v(x: f32, y: f32) -> V2 { V2 { x, y } }
fn dot(a: V2, b: V2) -> f32 { a.x * b.x + a.y * b.y }
fn cross(a: V2, b: V2) -> f32 { a.x * b.y - a.y * b.x }
fn cross_sv(w: f32, r: V2) -> V2 { v(-w * r.y, w * r.x) } // ω × r in 2D
fn add(a: V2, b: V2) -> V2 { v(a.x + b.x, a.y + b.y) }
fn sub(a: V2, b: V2) -> V2 { v(a.x - b.x, a.y - b.y) }
fn scale(a: V2, s: f32) -> V2 { v(a.x * s, a.y * s) }

pub struct Body { pub vel: V2, pub ang_vel: f32, pub inv_mass: f32, pub inv_inertia: f32 }

/// One contact point. `normal` points from body a to body b; `separation` < 0 = penetration.
/// Accumulated impulses persist across frames (matched by feature id) for warm starting.
pub struct ContactPoint {
    pub a: usize, pub b: usize,
    pub ra: V2, pub rb: V2,            // contact point relative to each body's center of mass
    pub normal: V2, pub separation: f32, pub friction: f32, pub restitution: f32,
    pub pn: f32, pub pt: f32,          // accumulated normal / tangent impulse
    mass_n: f32, mass_t: f32, bias: f32,
}

impl ContactPoint {
    pub fn new(a: usize, b: usize, ra: V2, rb: V2, normal: V2, separation: f32) -> Self {
        ContactPoint { a, b, ra, rb, normal, separation, friction: 0.6, restitution: 0.0,
            pn: 0.0, pt: 0.0, mass_n: 0.0, mass_t: 0.0, bias: 0.0 }
    }
}

const BAUMGARTE: f32 = 0.2;           // fraction of penetration removed per step
const LINEAR_SLOP: f32 = 0.005;       // allowed penetration (m): keeps contacts persistent
const RESTITUTION_THRESHOLD: f32 = 1.0; // m/s: no bounce below this approach speed

fn rel_vel(bodies: &[Body], c: &ContactPoint) -> V2 {
    let (a, b) = (&bodies[c.a], &bodies[c.b]);
    sub(add(b.vel, cross_sv(b.ang_vel, c.rb)), add(a.vel, cross_sv(a.ang_vel, c.ra)))
}

fn apply(bodies: &mut [Body], c: &ContactPoint, p: V2) {
    let (ra, rb) = (c.ra, c.rb);
    let a = &mut bodies[c.a];
    a.vel = sub(a.vel, scale(p, a.inv_mass)); a.ang_vel -= a.inv_inertia * cross(ra, p);
    let b = &mut bodies[c.b];
    b.vel = add(b.vel, scale(p, b.inv_mass)); b.ang_vel += b.inv_inertia * cross(rb, p);
}

/// Velocity solve for one step. Positions are integrated afterwards with the new velocities.
pub fn solve(bodies: &mut [Body], contacts: &mut [ContactPoint], dt: f32, iterations: usize) {
    for c in contacts.iter_mut() {
        let (a, b) = (&bodies[c.a], &bodies[c.b]);
        let t = v(c.normal.y, -c.normal.x);
        let (ra, rb) = (c.ra, c.rb);
        let k = |d: V2| a.inv_mass + b.inv_mass + a.inv_inertia * cross(ra, d).powi(2)
                        + b.inv_inertia * cross(rb, d).powi(2);
        c.mass_n = 1.0 / k(c.normal);   // effective mass along the normal
        c.mass_t = 1.0 / k(t);
        // bias = target normal velocity. Gap: may close it within this step (speculative).
        // Penetration: Baumgarte pushes out a fraction of the depth beyond the slop.
        c.bias = if c.separation > 0.0 { -c.separation / dt }
                 else { BAUMGARTE / dt * (-c.separation - LINEAR_SLOP).max(0.0) };
        let vn = dot(rel_vel(bodies, c), c.normal);
        if c.separation <= 0.0 && vn < -RESTITUTION_THRESHOLD { // bounce only when touching
            c.bias = c.bias.max(-c.restitution * vn);
        }
    }
    for c in contacts.iter() { // warm start: re-apply last step's impulses
        let p = add(scale(c.normal, c.pn), scale(v(c.normal.y, -c.normal.x), c.pt));
        apply(bodies, c, p);
    }
    for _ in 0..iterations {
        for c in contacts.iter_mut() {
            let t = v(c.normal.y, -c.normal.x);
            // Friction first: the normal constraint, solved last, matters more.
            let lambda = -c.mass_t * dot(rel_vel(bodies, c), t);
            let max_f = c.friction * c.pn;
            let new_pt = (c.pt + lambda).clamp(-max_f, max_f); // Coulomb cone on the total
            let dp = new_pt - c.pt; c.pt = new_pt;
            apply(bodies, c, scale(t, dp));

            let lambda = c.mass_n * (c.bias - dot(rel_vel(bodies, c), c.normal));
            let new_pn = (c.pn + lambda).max(0.0); // clamp the accumulated impulse, not lambda
            let dp = new_pn - c.pn; c.pn = new_pn;
            apply(bodies, c, scale(c.normal, dp));
        }
    }
}

A 1 kg box resting on static ground at 60 Hz settles with the two normal impulses summing to m·g·dt = 0.1667 N·s per step and zero velocity.

Complexity and constants

• O(k · c) per step for c contact points and k iterations; each constraint update is order of 20-40 flops in 2D plus two body loads/stores.
• Memory per contact point: order of 64-100 bytes (anchors, normal, masses, accumulated impulses, bias).
• Error decays geometrically per iteration, with a rate that worsens with mass ratio and chain length; that is why warm starting (start near the answer) matters more than raising k.
• Typical budgets: Box2D v2 8 velocity / 3 position iterations; Jolt defaults to 10 velocity and 2 position steps; Box2D v3 recommends 4 substeps.

Efficiency

• Cost: O(k·c) per step; one constraint update is order of 20-40 flops in 2D plus two body loads and stores, and a contact point stores order of 64-100 bytes (this page). The bound is memory access to the bodies and the serial Gauss-Seidel dependency, not arithmetic.
• Levers:
  • warm start — same stability with fewer iterations; starting near the answer beats raising k — persistent contacts
  • graph colouring plus SIMD — constraints of one colour solved in parallel (Box2D v3) — many contacts, many cores
  • compact solver bodies — inner loop touches only hot data — large islands
  • soft constraints with substeps — replaces extra iterations and Baumgarte bias (Box2D v3, Solver2D) — stiff stacks
  • sleeping islands — resting bodies cost nothing (Box2D sleeps after 0.5 s) — settled scenes
• Naive → standard → state of the art:
  • solve all contacts at once as one linear complementarity problem · dense direct solve, order of c³
  • sequential impulses with warm start, 8 velocity + 3 position iterations (Box2D v2) · O(k·c)
  • soft step with 4 substeps and graph-coloured SIMD sweeps (Box2D v3) · O(k·c), parallel across colours

Implementation notes

• Accumulated clamping: clamp the running total (Pn ≥ 0, |Pt| ≤ μPn), never the per-iteration λ. Clamping λ alone makes the solver unable to take back impulse and stacks never settle.
• Contact matching: warm starting needs the same physical contact across frames. Key points by feature ids from the narrowphase (reference edge, incident vertex) and copy Pn, Pt when ids match.
• Solver bodies in SoA: copy v, ω, inverse mass and inverse inertia of the bodies in the island into a compact array; the inner loop then touches only hot data. Static bodies get inverse mass 0 and are never written.
• Position correction: Baumgarte bias injects energy (bouncy pops after deep penetration). Alternatives: a separate split-impulse / pseudo-velocity pass (Bullet), nonlinear Gauss–Seidel on positions (Box2D v2's position iterations), or soft constraints with substepping (Box2D v3).
• Block solver: solve the two normal constraints of a 2D manifold together as a 2×2 LCP (Box2D v2 does this); fixes rocking boxes.
• Parallelism: Gauss–Seidel is sequential per body. Colour the constraint graph so no two constraints in a colour share a dynamic body, then solve each colour in parallel with SIMD (Box2D v3 graph colouring), or use mass splitting (Tonge et al.).
• Islands and sleeping: union-find over contacts and joints; an island whose bodies stay below a velocity threshold for about half a second goes to sleep (Box2D's b2_timeToSleep is 0.5 s).

Pitfalls

• No warm starting, or warm starting with impulses from mismatched contacts: jitter and sinking stacks.
• Friction solved with the current iteration's normal impulse instead of the accumulated one: friction cone wrong, objects slide on slopes.
• Restitution applied to slow contacts: resting objects vibrate forever. Use a threshold (1 m/s in Box2D).
• Restitution with speculative contacts: the gap is closed without a bounce; handle restitution in a separate pass using the pre-solve approach speed (Box2D v3's relax/restitution phase).
• Baumgarte factor near 1 or no slop: contacts flicker between touching and separated, and stacks explode.
• Variable time step: impulses and bias scale with dt; run physics at a fixed step and interpolate rendering.

Variants and alternatives

In production

Games: Box2D (v2 sequential impulses with block solver; v3 soft step with graph colouring), Jolt Physics, Chipmunk2D. Robotics, VR and visual effects: Bullet (btSequentialImpulseConstraintSolver), used for simulation in robotics and machine learning per its README. PhysX offers PGS by default and TGS optionally.

Implementations

Sources

• Erin Catto, Iterative Dynamics with Temporal Coherence, GDC 2005. pdf
• Erin Catto, Fast and Simple Physics using Sequential Impulses, GDC 2006. pdf
• Erin Catto, Modeling and Solving Constraints, GDC 2009. pdf
• Erin Catto, Solver2D, box2d.org blog, 2024. post
• Richard Tonge, Feodor Benevolenski, Andrey Voroshilov, Mass Splitting for Jitter-Free Parallel Rigid Body Simulation, SIGGRAPH 2012. paper
• Kenny Erleben, Stable, Robust, and Versatile Multibody Dynamics Animation, PhD thesis, University of Copenhagen, 2004.