Physics Simulation

01Why simulate physics

A physics engine is the runtime that turns artist-placed geometry into something that responds to forces: crates that stack, ragdolls that collapse, vehicles that corner, destruction that propagates. The simulation usually runs at a fixed rate (commonly 50 or 60 Hz for gameplay, higher for vehicle dynamics and VR hand interaction) inside an outer variable-rate render loop. Every frame the engine (1) integrates forces into velocities and positions, (2) detects collisions between all pairs of overlapping shapes, (3) generates contact points where shapes touch, (4) solves a system of constraints that prevent penetration and enforce joints, and (5) writes the resulting transforms back for rendering.

The workload is dominated by step 4. A stack of 20 boxes generates roughly 80 contact constraints (4 contact points per face-to-face contact, 20 contacts), and the solver iterates over every constraint multiple times per frame. At 8 iterations and 120 Hz that is 76,800 constraint solves per second for a single stack. Scale to a destruction scene with thousands of fragments and the constraint solver becomes the frame-rate limiter.

The core loop is small. A minimal rigid body engine for convex shapes is on the order of a few thousand lines of C++. Most of those lines are in collision detection and constraint solving; integration is about 10 lines. Most of the engineering effort goes into stability, determinism, and performance at scale.

02A short history

03Integration methods

Integration is the step that turns forces into motion. Given a body with mass m, position x, velocity v, and net force F, Newton's second law gives the acceleration:

The question is how to advance v and x by a discrete timestep dt. Four methods cover the space:

Order of operations matters. Explicit Euler updates position with the old velocity, then updates velocity. Semi-implicit Euler updates velocity first, then uses the new velocity to update position. That single reorder makes the integrator symplectic: among other things it preserves phase-space area, and for a conservative system it keeps the energy bounded instead of pumping it into oscillations. The widget below makes the difference visible.

04The semi-implicit Euler step

The complete integration step for a single rigid body, including rotation, is:

Position is updated with the new velocity. Rotation follows the same pattern: angular velocity omega is updated by the torque transformed by the inverse inertia tensor, and the orientation quaternion is advanced by omega. The code below drops the gyroscopic term ω × (Iω), as game engines commonly do by default; PhysX, for example, only adds it when a flag asks for it^[12].

// Semi-implicit Euler integration for a rigid body.
// Velocity is updated FIRST, then position uses the new velocity.
void integrate(RigidBody& body, float dt) {
    if (body.inverseMass == 0.0f) return;  // static body, skip

    // --- Linear ---
    // accumulate gravity and any applied forces into acceleration
    Vec3 linearAccel = body.force * body.inverseMass;
    body.linearVelocity += linearAccel * dt;

    // position advances using the NEW velocity (symplectic)
    body.position += body.linearVelocity * dt;

    // --- Angular ---
    // world-space inverse inertia times torque = angular acceleration.
    // The gyroscopic term omega x (I * omega) is left out, as many game
    // engines do by default; spinning bodies lose a little realism.
    Vec3 angularAccel = body.inverseInertiaTensorWorld * body.torque;
    body.angularVelocity += angularAccel * dt;

    // integrate orientation quaternion
    // dq/dt = 0.5 * omega * q  (quaternion multiplication)
    // spin is the pure quaternion (w = 0, xyz = omega)
    Quaternion spin(0.0f,
        body.angularVelocity.x,
        body.angularVelocity.y,
        body.angularVelocity.z);
    body.orientation += (spin * body.orientation) * (0.5f * dt);
    body.orientation.normalize();  // the additive update drifts off unit length

    // clear per-frame accumulators
    body.force  = Vec3(0.0f);
    body.torque = Vec3(0.0f);
}

/// Semi-implicit Euler integration for a rigid body.
/// Velocity is updated FIRST, then position uses the new velocity.
fn integrate(body: &mut RigidBody, dt: f32) {
    if body.inverse_mass == 0.0 { return; }  // static body, skip

    // --- Linear ---
    // accumulate gravity and any applied forces into acceleration
    let linear_accel = body.force * body.inverse_mass;
    body.linear_velocity += linear_accel * dt;

    // position advances using the NEW velocity (symplectic)
    body.position += body.linear_velocity * dt;

    // --- Angular ---
    // Same simplification as the C++ pane: no gyroscopic term.
    let angular_accel = body.inverse_inertia_tensor_world * body.torque;
    body.angular_velocity += angular_accel * dt;

    // integrate orientation quaternion
    // dq/dt = 0.5 * omega * q  (quaternion multiplication)
    let spin = Quat::from_xyzw(
        body.angular_velocity.x,
        body.angular_velocity.y,
        body.angular_velocity.z,
        0.0,
    );
    body.orientation = body.orientation + (spin * body.orientation) * (0.5 * dt);
    body.orientation = body.orientation.normalize();  // the additive update drifts off unit length

    // clear per-frame accumulators
    body.force  = Vec3::ZERO;
    body.torque = Vec3::ZERO;
}

05Collision detection: broad phase

With n bodies, testing every pair for collision is O(n²). At 5,000 bodies that is 12.5 million pairs per frame. The broad phase prunes the pair count to the small set of bodies whose bounding volumes actually overlap.

Three data structures dominate:

• Sweep-and-prune (SAP). Project every AABB onto one axis. Sort the endpoints. Walk the sorted list: when a "begin" endpoint is encountered, the body overlaps with every currently-open body. O(n log n) for the initial sort, close to O(n) per frame afterwards if bodies move slowly, because the list stays nearly sorted and insertion sort only pays for the endpoints that moved past each other. PhysX offers it as eSAP^[12], and Bullet as btAxisSweep3.
• Spatial hashing. Hash each AABB into grid cells. Bodies in the same cell are candidate pairs. O(n) expected time. Good for uniform-size objects (particles, bullet casings). Poor when objects vary wildly in size (a character standing on terrain).
• Dynamic AABB tree. A bounding-volume hierarchy where each leaf is a body's AABB and internal nodes are the union AABB. Queries are tree traversals. Insertion and removal are O(log n) in a balanced tree. Box2D and Bullet's default broadphase use dynamic AABB trees, and Jolt uses a four-way variant (a quadtree of boxes)^[9]. Good for scenes with mixed sizes and frequent insertion and removal.

// Sweep-and-prune on one axis.
// Returns candidate pairs whose AABBs overlap on the X axis.
// A second pass on Y (and Z for 3D) prunes further.
struct Endpoint {
    float value;       // x-coordinate of the AABB edge
    int   bodyIndex;   // which body this endpoint belongs to
    bool  isMin;       // true = left edge, false = right edge
};

std::vector<std::pair<int,int>> sweepAndPrune(
    std::vector<Endpoint>& endpoints)
{
    // insertion sort: O(n + swaps), so close to linear when the list is
    // still nearly sorted from last frame
    for (int i = 1; i < (int)endpoints.size(); ++i) {
        Endpoint key = endpoints[i];
        int j = i - 1;
        while (j >= 0 && endpoints[j].value > key.value) {
            endpoints[j + 1] = endpoints[j];
            --j;
        }
        endpoints[j + 1] = key;
    }

    // sweep: track which bodies are "open" (left edge seen, right edge not yet)
    std::vector<std::pair<int,int>> pairs;
    std::unordered_set<int> active;

    for (const auto& ep : endpoints) {
        if (ep.isMin) {
            // this body's AABB starts; it overlaps with every currently-active body
            for (int other : active) {
                int a = std::min(ep.bodyIndex, other);
                int b = std::max(ep.bodyIndex, other);
                pairs.push_back({a, b});
            }
            active.insert(ep.bodyIndex);
        } else {
            // this body's AABB ends
            active.erase(ep.bodyIndex);
        }
    }
    return pairs;
}

/// Sweep-and-prune on one axis.
/// Returns candidate pairs whose AABBs overlap on the X axis.
#[derive(Clone, Copy)]
struct Endpoint {
    value: f32,
    body_index: usize,
    is_min: bool,
}

fn sweep_and_prune(endpoints: &mut Vec<Endpoint>) -> Vec<(usize, usize)> {
    // insertion sort: O(n + swaps), so close to linear when the list is
    // still nearly sorted from last frame
    for i in 1..endpoints.len() {
        let key = endpoints[i];
        let mut j = i;
        while j > 0 && endpoints[j - 1].value > key.value {
            endpoints[j] = endpoints[j - 1];
            j -= 1;
        }
        endpoints[j] = key;
    }

    let mut pairs = Vec::new();
    let mut active = std::collections::HashSet::new();

    for ep in endpoints.iter() {
        if ep.is_min {
            for &other in &active {
                let a = ep.body_index.min(other);
                let b = ep.body_index.max(other);
                pairs.push((a, b));
            }
            active.insert(ep.body_index);
        } else {
            active.remove(&ep.body_index);
        }
    }
    pairs
}

06Narrow phase: GJK

The narrow phase tests the actual geometry of a candidate pair. For convex shapes, the dominant algorithm is GJK (Gilbert, Johnson, Keerthi, 1988^[2]).

The idea: two convex shapes overlap if and only if their Minkowski difference contains the origin. GJK searches the Minkowski difference without constructing it explicitly, using support functions to probe specific directions. It builds a simplex (up to a tetrahedron in 3D) and checks whether it contains the origin.

If GJK reports overlap, EPA (Expanding Polytope Algorithm) takes the final simplex and expands it to find the penetration depth and contact normal.

// Minkowski-difference support: farthest point on A in dir d, minus
// farthest point on B in dir -d. Declared first so gjkIntersect can call it.
Vec2 minkowskiSupport(const Shape& a, const Shape& b, Vec2 d) {
    return a.support(d) - b.support(-d);
}

// GJK intersection test for two convex shapes in 2D.
// Returns true if the shapes overlap.
bool gjkIntersect(const Shape& shapeA, const Shape& shapeB) {
    // pick an initial search direction (arbitrary; toward B's center works)
    Vec2 direction = shapeB.center() - shapeA.center();
    if (direction.lengthSquared() < 1e-8f) direction = Vec2(1.0f, 0.0f);

    // first support point on the Minkowski difference
    Vec2 support = minkowskiSupport(shapeA, shapeB, direction);
    Simplex simplex;
    simplex.add(support);

    // new search direction: toward the origin from the support point
    direction = -support;

    for (int iter = 0; iter < 32; ++iter) {
        Vec2 newPoint = minkowskiSupport(shapeA, shapeB, direction);

        // if the new point did not pass the origin, shapes are separated
        if (dot(newPoint, direction) < 0.0f) return false;

        simplex.add(newPoint);

        // check if the simplex contains the origin; if so, update direction
        if (simplex.containsOrigin(direction)) return true;
        // 'containsOrigin' also prunes the simplex to the closest feature
        // and sets 'direction' toward the origin from that feature
    }
    return false;  // iteration cap: only degenerate or touching cases get here
}

/// GJK intersection test for two convex shapes in 2D.
/// Returns true if the shapes overlap.
fn gjk_intersect(shape_a: &impl ConvexShape, shape_b: &impl ConvexShape) -> bool {
    let mut direction = shape_b.center() - shape_a.center();
    if direction.length_squared() < 1e-8 {
        direction = Vec2::X;
    }

    let support = minkowski_support(shape_a, shape_b, direction);
    let mut simplex = Simplex::new();
    simplex.add(support);
    direction = -support;

    for _iter in 0..32 {
        let new_point = minkowski_support(shape_a, shape_b, direction);

        // new point did not pass the origin: shapes are separated
        if new_point.dot(direction) < 0.0 { return false; }

        simplex.add(new_point);

        // simplex contains origin: shapes overlap
        if simplex.contains_origin(&mut direction) { return true; }
    }
    false
}

fn minkowski_support(a: &impl ConvexShape, b: &impl ConvexShape, d: Vec2) -> Vec2 {
    a.support(d) - b.support(-d)
}

07SAT (Separating Axis Theorem)

The SAT is the other major narrow-phase algorithm. For convex polytopes, it is often faster than GJK for simple shapes (boxes, triangles) because it directly produces the contact normal and penetration depth without a separate EPA pass.

The theorem: two convex shapes are separated if and only if there exists a line (an "axis") on which their projections do not overlap. For polygons, the candidate axes are the face normals of both shapes. For polyhedra in 3D, add the cross products of every edge pair between the two shapes. If all candidate axes show overlap, the shapes are colliding, and the axis with minimum overlap is the contact normal.

08Contact points and manifolds

A single collision between two boxes produces a contact manifold: a set of contact points, each with a position, normal, and penetration depth. Face-to-face contact between two boxes produces the intersection polygon of the two faces, up to eight vertices for two rectangles, which engines usually reduce to four points. Edge-to-edge generates 1 point.

Stable stacking requires persistent contact manifolds. The solver needs contacts to survive across frames so it can warm-start from the previous frame's solution. Contact point matching (by feature ID or by proximity) keeps the manifold stable as bodies shift slightly. Without persistence, every frame starts from scratch, and the solver needs many more iterations to converge.

Contact point reduction matters for performance. A face-to-face contact in 3D can produce an arbitrary polygon of contact points. Reducing it to 4 points, typically the deepest point plus points chosen to span the largest area in the contact plane, keeps the solver's constraint count bounded without losing stability.

09Constraint-based dynamics

A constraint is a restriction on how bodies may move. "These two bodies must not interpenetrate" is a contact constraint. "This body is attached to that body at a hinge" is a joint constraint. The constraint solver's job: find impulses that satisfy all active constraints simultaneously.

The formal framework (Baraff 1994^[3], Catto 2005^[5], Catto 2014^[8]): each constraint is a scalar function C of the positions and orientations of the involved bodies.

Differentiating C gives the velocity constraint; the solver enforces it with an added bias term b that corrects drift. The J matrix (the Jacobian) maps body velocities to the rate of constraint violation. The solver computes an impulse lambda along each constraint's Jacobian direction such that the velocity-level constraint is satisfied.

Position drift (bodies slowly sinking into each other) is corrected by Baumgarte stabilization: add a bias term to the velocity constraint that is proportional to the current position error. For a contact, C is the signed separation along the normal, negative while the bodies overlap, so b comes out negative and J v = −b asks for a separating velocity.

10Sequential impulses

The sequential impulse solver (Catto 2005^[5], 2006^[10]) is a projected Gauss-Seidel iteration over the constraints. For each constraint:

1. Compute the relative velocity at the contact point along the constraint direction: relVel = J * v.
2. Compute the impulse magnitude: lambda = -(relVel + bias) * effectiveMass, where effectiveMass = 1 / (J M⁻¹ Jᵀ) is the mass the constraint direction sees.
3. Clamp: for a contact, the accumulated impulse must be non-negative (you can push but not pull). This is accumulated impulse clamping, not per-iteration clamping.
4. Apply the impulse to both bodies' velocities immediately.

Warm starting reuses the previous frame's accumulated impulses as the initial guess. Because the configuration changes little between frames, the warm-started guess is close to the converged answer, and the solver needs far fewer iterations. In Catto's 2005 tests at 10 iterations, a stack of three boxes without the contact cache slid apart and fell, while a stack of ten boxes with it settled within a few seconds and stayed up^[5]. Matching contacts from one frame to the next is what makes this possible, which is why persistent manifolds (§8) matter.

// Sequential impulse solver for contact constraints.
// Iterates 'iterations' times over all contacts, accumulating impulses.
void solveContacts(
    std::vector<ContactConstraint>& contacts,
    std::vector<RigidBody>& bodies,
    int iterations, float dt)
{
    // warm start: apply last frame's accumulated impulses
    for (auto& contact : contacts) {
        Vec2 impulse = contact.normal * contact.accumulatedNormalImpulse
                     + contact.tangent * contact.accumulatedTangentImpulse;
        bodies[contact.bodyA].linearVelocity -= impulse * bodies[contact.bodyA].inverseMass;
        bodies[contact.bodyB].linearVelocity += impulse * bodies[contact.bodyB].inverseMass;
    }

    for (int iter = 0; iter < iterations; ++iter) {
        for (auto& contact : contacts) {
            RigidBody& bodyA = bodies[contact.bodyA];
            RigidBody& bodyB = bodies[contact.bodyB];

            // relative velocity at the contact along the normal (A to B);
            // positive means separating
            Vec2 relVel = bodyB.linearVelocity - bodyA.linearVelocity;
            float normalVel = dot(relVel, contact.normal);

            // Baumgarte bias b = (beta / dt) * C, where C is the signed
            // separation (negative while overlapping). The 0.01 slop lets
            // resting contacts overlap slightly instead of jittering. b <= 0,
            // so the impulse below asks for a separating velocity of -b.
            float separation = -contact.penetration;
            float bias = (0.2f / dt) * std::min(0.0f, separation + 0.01f);

            // impulse that drives J v + b to zero
            float lambda = -(normalVel + bias) * contact.effectiveMassNormal;

            // accumulated clamping: total impulse must be non-negative
            float oldAccum = contact.accumulatedNormalImpulse;
            contact.accumulatedNormalImpulse = std::max(0.0f, oldAccum + lambda);
            lambda = contact.accumulatedNormalImpulse - oldAccum;

            // apply impulse to both bodies
            Vec2 impulse = contact.normal * lambda;
            bodyA.linearVelocity -= impulse * bodyA.inverseMass;
            bodyB.linearVelocity += impulse * bodyB.inverseMass;

            // --- friction (Coulomb approximation) ---
            relVel = bodyB.linearVelocity - bodyA.linearVelocity;
            float tangentVel = dot(relVel, contact.tangent);
            float frictionLambda = -tangentVel * contact.effectiveMassTangent;

            // clamp friction to Coulomb cone
            float maxFriction = contact.friction * contact.accumulatedNormalImpulse;
            float oldFric = contact.accumulatedTangentImpulse;
            contact.accumulatedTangentImpulse = std::clamp(
                oldFric + frictionLambda, -maxFriction, maxFriction);
            frictionLambda = contact.accumulatedTangentImpulse - oldFric;

            Vec2 frictionImpulse = contact.tangent * frictionLambda;
            bodyA.linearVelocity -= frictionImpulse * bodyA.inverseMass;
            bodyB.linearVelocity += frictionImpulse * bodyB.inverseMass;
        }
    }
}

/// Sequential impulse solver for contact constraints.
/// Iterates `iterations` times over all contacts, accumulating impulses.
fn solve_contacts(
    contacts: &mut [ContactConstraint],
    bodies: &mut [RigidBody],
    iterations: usize, dt: f32,
) {
    // warm start: apply last frame's accumulated impulses
    for contact in contacts.iter() {
        let impulse = contact.normal * contact.accumulated_normal_impulse
                    + contact.tangent * contact.accumulated_tangent_impulse;
        bodies[contact.body_a].linear_velocity -= impulse * bodies[contact.body_a].inverse_mass;
        bodies[contact.body_b].linear_velocity += impulse * bodies[contact.body_b].inverse_mass;
    }

    for _iter in 0..iterations {
        for contact in contacts.iter_mut() {
            // relative velocity at the contact along the normal (A to B)
            let rel_vel = bodies[contact.body_b].linear_velocity
                        - bodies[contact.body_a].linear_velocity;
            let normal_vel = rel_vel.dot(contact.normal);

            // Baumgarte bias with the same sign convention as the C++ pane
            let separation = -contact.penetration;
            let bias = (0.2 / dt) * (separation + 0.01).min(0.0);
            let lambda = -(normal_vel + bias) * contact.effective_mass_normal;

            // accumulated clamping: total impulse must be non-negative
            let old_accum = contact.accumulated_normal_impulse;
            contact.accumulated_normal_impulse = (old_accum + lambda).max(0.0);
            let lambda = contact.accumulated_normal_impulse - old_accum;

            let impulse = contact.normal * lambda;
            bodies[contact.body_a].linear_velocity -= impulse * bodies[contact.body_a].inverse_mass;
            bodies[contact.body_b].linear_velocity += impulse * bodies[contact.body_b].inverse_mass;

            // --- friction (Coulomb approximation) ---
            let rel_vel = bodies[contact.body_b].linear_velocity
                        - bodies[contact.body_a].linear_velocity;
            let tangent_vel = rel_vel.dot(contact.tangent);
            let friction_lambda = -tangent_vel * contact.effective_mass_tangent;

            // clamp friction to Coulomb cone
            let max_friction = contact.friction * contact.accumulated_normal_impulse;
            let old_friction = contact.accumulated_tangent_impulse;
            contact.accumulated_tangent_impulse =
                (old_friction + friction_lambda).clamp(-max_friction, max_friction);
            let friction_lambda = contact.accumulated_tangent_impulse - old_friction;

            let friction_impulse = contact.tangent * friction_lambda;
            bodies[contact.body_a].linear_velocity -= friction_impulse * bodies[contact.body_a].inverse_mass;
            bodies[contact.body_b].linear_velocity += friction_impulse * bodies[contact.body_b].inverse_mass;
        }
    }
}

11Common constraints

Every constraint is a Jacobian row. The solver does not care what physical meaning the constraint has; it computes impulses that zero the velocity-level violation. The common constraints in a game engine:

• Contact (non-penetration). With velocities ordered (v_A, ω_A, v_B, ω_B), J = [−n, −(r_A × n), n, r_B × n], where n is the contact normal pointing from A to B and r is the vector from each centre of mass to the contact point, so J v is the separating speed. The impulse must be non-negative (push only). The solver code above implements the linear half of this row.
• Friction. Tangent-plane constraint with impulse clamped to the Coulomb cone. The maximum friction impulse is mu times the accumulated normal impulse.
• Distance joint. Keeps two anchor points at a fixed distance. C = |p_B − p_A| − L, where L is the rest length, and J = [−d, −(r_A × d), d, r_B × d], where d is the unit vector from anchor A to anchor B.
• Revolute joint (hinge). Forces two anchor points to coincide: two rows in 2D. In 3D the shared point takes three rows, and two more angular rows remove rotation about every axis except the hinge's.
• Prismatic joint (slider). Allows motion along one axis, constrains all others. The Jacobian zeroes out the allowed-axis component.

12Substeps and fixed timestep

Physics usually runs at a fixed timestep, for reproducibility and because solver tuning and stability limits depend on dt. The standard pattern, popularized by Glenn Fiedler's "Fix Your Timestep!"^[13]: an accumulator absorbs real elapsed time, and the physics loop ticks in fixed-size bites.

// Fixed-timestep accumulator pattern.
// 'elapsed' is wall-clock time since last render frame (variable).
// 'physDt' is the fixed physics timestep (e.g. 1/120 s).
accumulator += elapsed;
while (accumulator >= physDt) {
    stepPhysics(physDt);
    accumulator -= physDt;
}
// interpolation factor for rendering between the last two physics states
float alpha = accumulator / physDt;
renderState = lerp(previousPhysicsState, currentPhysicsState, alpha);

The rendering interpolation (the alpha lerp) smooths visual motion between discrete physics ticks. Without it, objects stutter whenever the render rate and the physics rate do not line up, because some frames get two physics steps and others none.

Variable dt breaks determinism. A variable dt comes from measured frame times, which differ from run to run and machine to machine, so the same inputs step through different sequences of dt values. Those sequences do not land on the same state: the integrator's error depends on the step size, and floating-point rounding means even one 1/30 s step and two 1/60 s steps disagree. Replays desync, networked physics diverges, and QA cannot reproduce bugs. A fixed dt is a prerequisite for reproducibility, though not sufficient on its own (see the floating-point pitfall in §15).

Substepping (running multiple physics ticks per frame with a smaller dt) improves stability for stiff constraints and high-speed objects. Jolt's PhysicsSystem::Update takes a number of collision steps that splits each update into smaller collision-plus-integration steps, and its documentation says the system is generally stable at 60 Hz with one^[9]. PhysX's TGS solver substeps internally: it splits the step into as many substeps as there are position iterations and integrates after each, which its documentation calls conceptually similar to running PGS several times at a smaller timestep with one iteration each^[12].

13Continuous collision detection

Discrete collision detection tests shapes at their positions at the end of the timestep. A fast-moving object can pass entirely through a thin wall between two frames. This is the tunneling problem: the object was on one side at t, on the other side at t + dt, and the collision detector never saw an overlap.

CCD fixes this by computing the TOI (time of impact): the earliest time during the timestep at which two swept shapes first touch. The simulation is advanced to the TOI, the collision is resolved, and the remainder of the timestep continues.

Speculative contacts are a cheaper alternative. Instead of computing exact TOI, the solver adds a constraint for any pair that will be within contact distance at the end of the step (based on current velocity). The bias velocity is set to prevent penetration at the predicted position. Speculative contacts catch much of the tunneling that discrete detection misses, at a fraction of the cost of swept-volume CCD. Jolt creates speculative contacts within a small distance (2 cm by default) during its normal discrete collision pass, and runs a swept cast only for bodies whose motion quality is set to LinearCast^[9].

14Case studies

1 · Box2D (Erin Catto)

Box2D^[6] is the most widely read implementation of sequential-impulse dynamics. The v2.4 architecture is small enough to read end to end: a dynamic AABB tree for broad phase, GJK for distance queries, SAT-based contact generation for polygons, and the sequential-impulse solver with warm starting and accumulated clamping, run for a suggested 8 velocity and 3 position iterations per step^[14]. Box2D's constraint formulation follows Catto's GDC presentations. Version 3 (2024) replaced the iteration counts with substepping and soft constraints. Shipped in 2D games such as Angry Birds, Limbo, and Crayon Physics Deluxe.

2 · Jolt Physics

Jolt Physics^[9] is an open-source C++ rigid-body engine by Jorrit Rouwe (v1.0 in 2022). Its broad phase is a quadtree of AABBs, split into configurable layers, that bodies can update from several threads without locks. Narrow phase uses GJK and EPA. The solver is a SIMD-optimized sequential-impulse solver with warm starting (10 velocity and 2 position iterations by default), run per simulation island as jobs, with large islands split into groups that can be solved in parallel. It offers opt-in cross-platform determinism at roughly an 8% cost. Its README lists Horizon Forbidden West and Death Stranding 2 among the games that use it.

3 · PhysX: PGS and TGS

NVIDIA PhysX (open source) offers two iterative solvers, projected Gauss-Seidel (PGS) and TGS, both available on CPU and GPU^[12]. PGS runs all its iterations on the same constraints and then integrates; TGS substeps, integrating after every position iteration, and the documentation calls it the generally recommended solver. The two also differ in friction: PGS clamps the two tangent directions separately, while TGS clamps their combined impulse. Gyroscopic forces are off unless a per-body flag enables them.

4 · Havok

Havok, now part of Microsoft, is a long-running commercial physics middleware whose solver internals are not publicly documented. Its site lists titles such as Assassin's Creed Odyssey, The Legend of Zelda: Breath of the Wild, and The Elder Scrolls IV: Oblivion Remastered among the games built on Havok Physics^[15].

15Pitfalls

• Jitter from an unconverged solver. A stack of boxes vibrates or sinks because the solver ran out of iterations before converging. Enable warm starting before raising the iteration count. If warm starting is already on, check contact persistence: contacts that flicker in and out of existence lose their cached impulses.
• Energy gain from explicit Euler. Covered in detail in sections 3 and 4. Springs and oscillatory forces gain energy every step, eventually causing objects to fly off. Fix: use semi-implicit Euler.
• Sleeping heuristics gone wrong. Bodies "fall asleep" (excluded from the solver) when their velocity or kinetic energy stays below a threshold for some time (half a second in Box2D, a wake counter in PhysX). Set the threshold too high or the time too short and a slowly sliding or tipping body freezes in place; set it too low and CPU goes to stationary bodies that never qualify, often because solver jitter keeps them moving.
• Constraint fight (over-constrained systems). When constraints contradict each other (a body is being pulled in opposite directions by two joints, or a body is wedged between two static walls), the solver oscillates. Each iteration pushes in one direction, the next iteration pushes back. The result is jitter or explosions. Fix: reduce over-constraint by adding compliance (softness) to joints, or redesign the constraint graph.
• Floating-point non-determinism across platforms. The basic IEEE 754 operations are bit-exact on every compliant CPU, but compiled physics code often is not: compilers fuse multiply-adds, reorder arithmetic under fast-math, x87 builds keep extra precision, and sin, cos and sqrt-based helpers come from different math libraries per platform. Networked physics that relies on deterministic lockstep then diverges. Fix: disable fast-math and FMA contraction, avoid platform math-library transcendentals in the simulation, and verify bitwise identity on test scenes. Jolt ships this as an option that keeps results identical across compilers, operating systems and CPU architectures for about an 8% cost^[9].
• Tunneling. Covered in section 13. Fast objects pass through thin geometry. Fix: CCD for flagged bodies.
• Bad inertia tensors. An inertia tensor that does not match the body's mass distribution misbehaves in both directions: too small and the body spins up from tiny off-centre impulses and destabilizes its contacts; too large and it refuses to tip over. Common source: computing inertia from imported mesh data without accounting for density or scale, or from a non-closed mesh.

16What's next

• Position-Based Dynamics (PBD/XPBD). Müller et al. 2006^[7]. Solves constraints by directly moving positions instead of computing impulses, which makes cloth and soft bodies easy to keep stable. Extended PBD (XPBD, Macklin et al. 2016) adds compliance so stiffness no longer depends on the iteration count and timestep; the PhysX documentation cites it to justify the approximations in its own solvers^[12].
• GPU physics. PhysX runs both of its rigid-body solvers on the GPU as well as the CPU^[12], and NVIDIA's Flex did unified particle-based simulation. GPU pipelines scale to very large particle and body counts, but kernel-launch and transfer overheads make them a poor fit for small scenes.
• Deformable bodies. The Finite Element Method (FEM) for volumetric deformation, fracture, bending, and tearing. Real-time use depends on coarse tetrahedral meshes and a small number of simulated objects.
• MPM (Material Point Method). Hybrid particle-grid method for snow, mud, sand, and granular materials. Disney used MPM for snow in Frozen (Stomakhin et al. 2013). Real-time MPM is an active research area.
• Fluid simulation. SPH (Smoothed Particle Hydrodynamics) for Lagrangian fluids, and FLIP/APIC for hybrid particle-grid methods. FLIP is the usual choice for offline VFX liquids (Houdini's FLIP solver); games mostly use cheaper height-field and wave models for water surfaces.

Sources

1. Loup Verlet. "Computer 'Experiments' on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules." Physical Review 159(1):98-103, 1967. journals.aps.org/pr/abstract/10.1103/PhysRev.159.98. Rediscovery of Størmer's integration method for molecular dynamics: a velocity-free position update whose energy stays bounded over long simulations.
2. E. G. Gilbert, D. W. Johnson, S. S. Keerthi. "A fast procedure for computing the distance between complex objects in three-dimensional space." IEEE Journal of Robotics and Automation 4(2):193-203, 1988. ieeexplore.ieee.org/document/2083. The GJK algorithm for convex distance computation via iterative simplex construction on the Minkowski difference.
3. David Baraff. "Fast Contact Force Computation for Nonpenetrating Rigid Bodies." Proc. SIGGRAPH 1994, pp. 23-34. dl.acm.org/doi/10.1145/192161.192168. A pivoting algorithm, derived from Dantzig's method for linear complementarity problems, that computes contact forces with static and dynamic friction without general optimization software, and that the paper reports as faster, simpler and more reliable than earlier approaches.
4. David Baraff. "Physically Based Modeling: Rigid Body Simulation." SIGGRAPH 1997 Course Notes. cs.cmu.edu/~baraff/sigcourse. Course notes covering rigid-body equations of motion, collision response, and resting contact. A standard introduction to rigid-body simulation.
5. Erin Catto. "Iterative Dynamics with Temporal Coherence." Game Developers Conference 2005. box2d.org (PDF). A linear-time projected Gauss-Seidel solver for contact and friction plus a contact cache for warm starting. Section 9 has the stacking results quoted in §10: with 10 iterations, three boxes without caching fell while ten boxes with caching stayed up.
6. Erin Catto. "Box2D: A 2D Physics Engine for Games." box2d.org. Open-source 2D rigid-body engine, the most widely read implementation of sequential-impulse dynamics.
7. Matthias Müller, Bruno Heidelberger, Marcus Hennix, John Ratcliff. "Position Based Dynamics." 3rd Workshop in Virtual Reality Interactions and Physical Simulation (VRIPHYS), 2006; journal version in J. Visual Communication and Image Representation 18(2), 2007. matthias-research.github.io/pages/publications/posBasedDyn.pdf. The PBD paper: constraint solving by direct position modification. The abstract notes the authors built a real-time cloth simulator for a games physics library on it.
8. Erin Catto. "Understanding Constraints." Game Developers Conference 2014. box2d.org (PDF). A walkthrough of the Jacobian formulation, constraint derivation, and the sequential impulse solver.
9. Jorrit Rouwe. "Jolt Physics" repository and Architecture documentation. github.com/jrouwe/JoltPhysics. Source for the Jolt claims on this page: the quadtree broad phase updated without locks, sequential impulses with warm starting, collision steps ("stable when running at 60 Hz with 1 collision step"), the Discrete and LinearCast motion qualities, the 2 cm default speculative contact distance and 10/2 default iteration counts (PhysicsSettings.h), the cross-platform determinism option at about 8% cost, and the README's list of games.
10. Erin Catto. "Fast and Simple Physics using Sequential Impulses." Game Developers Conference 2006. box2d.org (PDF). The sequential-impulse presentation: clamp the accumulated impulse rather than each increment, and "apply old accumulated impulses at the beginning of the step" for "less iterations and greater stability."
11. Erin Catto. "Numerical Methods." Game Developers Conference 2009. box2d.org (PDF). Explicit, symplectic and Verlet integration for games, with the advice that "first-order accuracy is usually sufficient for games" and that "you can safely ignore RK4."
12. NVIDIA. "Rigid Body Dynamics," PhysX 5.4 documentation. nvidia-omniverse.github.io, with the broad-phase list in Rigid Body Collision. PGS versus TGS (TGS substeps once per position iteration; both run on CPU and GPU), per-patch friction handling, the gyroscopic-forces flag, sleeping, and the eSAP broad phase.
13. Glenn Fiedler. "Fix Your Timestep!" Gaffer On Games, 2004. gafferongames.com. The accumulator loop and the render interpolation between the last two physics states used in §12.
14. Erin Catto. Box2D v2.4.1 manual, "Hello Box2D." github.com/erincatto/box2d. "The suggested iteration count for Box2D is 8 for velocity and 3 for position."
15. Havok. "Havok Physics." havok.com/havok-physics. Product page listing games built on Havok Physics.
16. Leonhard Euler. Theoria motus corporum solidorum seu rigidorum. A. F. Röse, 1765 (Eneström index E289). Euler Archive. The treatise on rigid-body rotation behind Euler's equations of motion.