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 runs at a fixed rate (typically 60 Hz for gameplay, 120 to 240 Hz for vehicle dynamics or VR hand tracking) 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 deceptively small. A minimal rigid body engine for convex shapes fits in roughly 3,000 lines of C++. Most of those lines are in collision detection and constraint solving; integration is about 10 lines. The hard part is not writing it, but making it stable, deterministic, and fast on the 50,000th body.

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:

The critical insight: 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: it preserves the phase-space volume of the system, which means it does not systematically pump energy 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:

The key: position is updated with the new velocity, not the old one. For rotation the pattern is the same. Angular velocity omega is updated by torque divided by the inertia tensor, and the orientation quaternion is updated by omega.

// 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 ---
    // torque / inertia = angular acceleration
    Vec3 angularAccel = body.inverseInertiaTensorWorld * body.torque;
    body.angularVelocity += angularAccel * dt;

    // integrate orientation quaternion
    // dq/dt = 0.5 * omega * q  (quaternion multiplication)
    Quaternion spin(0.0f,
        body.angularVelocity.x,
        body.angularVelocity.y,
        body.angularVelocity.z);
    body.orientation += (spin * body.orientation) * (0.5f * dt);
    body.orientation.normalize();

    // 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 ---
    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 ---
    let angular_accel = body.inverse_inertia_world * body.torque;
    body.angular_velocity += angular_accel * dt;

    // integrate orientation quaternion
    // dq/dt = 0.5 * omega * q
    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();

    // 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, O(n) incremental if bodies move slowly (the list stays nearly sorted). Used by PhysX, Bullet, and Jolt for the primary axis.
• 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). Box2D uses a dynamic AABB tree. Good for scenes with mixed sizes and frequent insertion/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: nearly sorted each frame, so O(n) amortized
    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: nearly sorted each frame, so O(n) amortized
    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.

// GJK intersection test for two convex shapes in 2D.
// Returns true if the shapes overlap.
// 'support' returns the Minkowski-difference support point for a direction.
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;  // should not reach here for convex shapes
}

// Minkowski-difference support: farthest point on A in dir d, minus
// farthest point on B in dir -d.
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.
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 generates up to 4 points (the intersection polygon of the two faces). 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 to the 4 points that best represent the polygon (typically the farthest points in the tangent 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.

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 penetration depth.

10Sequential impulses

The sequential impulse solver (Catto 2005^[5]) is a 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.
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. Without warm starting, a stack of 10 boxes needs roughly 30 iterations to look stable. With warm starting, 8 iterations is usually enough^[5].

// 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 contact point along the normal
            Vec2 relVel = bodyB.linearVelocity - bodyA.linearVelocity;
            float normalVel = dot(relVel, contact.normal);

            // Baumgarte bias: push apart proportional to penetration
            float bias = (0.2f / dt) * std::max(0.0f, contact.penetration - 0.01f);

            // impulse magnitude
            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.
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
                    + contact.tangent * contact.accumulated_tangent;
        bodies[contact.body_a].linear_velocity -= impulse * bodies[contact.body_a].inv_mass;
        bodies[contact.body_b].linear_velocity += impulse * bodies[contact.body_b].inv_mass;
    }

    for _iter in 0..iterations {
        for contact in contacts.iter_mut() {
            let vel_a = bodies[contact.body_a].linear_velocity;
            let vel_b = bodies[contact.body_b].linear_velocity;

            let rel_vel = vel_b - vel_a;
            let normal_vel = rel_vel.dot(contact.normal);

            let bias = (0.2 / dt) * (contact.penetration - 0.01).max(0.0);
            let lambda = -(normal_vel + bias) * contact.eff_mass_normal;

            let old = contact.accumulated_normal;
            contact.accumulated_normal = (old + lambda).max(0.0);
            let lambda = contact.accumulated_normal - old;

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

            // friction
            let rel_vel = bodies[contact.body_b].linear_velocity
                        - bodies[contact.body_a].linear_velocity;
            let tan_vel = rel_vel.dot(contact.tangent);
            let fric_lambda = -tan_vel * contact.eff_mass_tangent;
            let max_fric = contact.friction * contact.accumulated_normal;
            let old_f = contact.accumulated_tangent;
            contact.accumulated_tangent = (old_f + fric_lambda).clamp(-max_fric, max_fric);
            let fric_lambda = contact.accumulated_tangent - old_f;

            let fric_impulse = contact.tangent * fric_lambda;
            bodies[contact.body_a].linear_velocity -= fric_impulse * bodies[contact.body_a].inv_mass;
            bodies[contact.body_b].linear_velocity += fric_impulse * bodies[contact.body_b].inv_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). J = [n, r_A x n, -n, -(r_B x n)] where n is the contact normal and r is the vector from center of mass to contact point. The impulse must be non-negative (push only). This is the constraint the solver code above implements.
• 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. J = [d, r_A x d, -d, -(r_B x d)] where d is the unit vector between anchors. C = |p_B - p_A| - L where L is the rest length.
• Revolute joint (hinge). Forces two anchor points to coincide. Two (2D) or three (3D) constraint rows, one per axis of the position difference.
• Prismatic joint (slider). Allows motion along one axis, constrains all others. The Jacobian zeroes out the allowed-axis component.

12Substeps and fixed timestep

Physics must run at a fixed timestep for determinism and stability. The standard pattern: 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 visually stutter at the physics rate.

Variable dt breaks determinism. Floating-point arithmetic is not associative: (a + b) + c is not necessarily equal to a + (b + c). If the timestep varies frame-to-frame, the same initial conditions produce different results on different machines (or even different runs on the same machine). Replays desync, networked physics diverges, and QA cannot reproduce bugs. Fixed dt is non-negotiable for any simulation that needs reproducibility.

Substepping (running multiple physics ticks per frame with a smaller dt) improves stability for stiff constraints and high-speed objects. Jolt Physics defaults to 1 substep but recommends 2 for vehicle simulations. PhysX's substepping reduces the dt seen by the solver, which tightens constraint convergence.

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 handle most tunneling cases at a fraction of the cost of full swept-volume CCD. Jolt Physics uses speculative contacts as the default, with swept CCD reserved for user-flagged high-speed bodies^[9].

14Case studies

1 · Box2D (Erin Catto)

Box2D^[6] is a 2D rigid-body engine that serves as the reference implementation for sequential-impulse dynamics. The architecture is clean 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. Box2D's constraint formulation directly follows Catto's GDC presentations. Shipped in Angry Birds, Limbo, Crayon Physics Deluxe, and hundreds of other 2D titles.

2 · Jolt Physics

Jolt Physics^[9] is a modern C++ rigid-body engine released in 2022 by Jorrit Rouwe. Architecture: a lock-free broad phase (layered with configurable object layers), GJK/EPA for narrow phase, a SIMD-optimized sequential-impulse solver, and a job system for multi-threaded island solving. Jolt's performance is competitive with PhysX on standard benchmarks while being a single-developer open-source project. Adopted by Guerrilla Games for Horizon and by the Godot community.

3 · PhysX GPU solver

NVIDIA PhysX (now open-source) added a GPU-accelerated TGS solver in PhysX 4.0 (2018), available on both CPU and GPU. The GPU pipeline batches constraints into independent groups (coloring), solves each group in parallel, and uses temporal coherence (warm starting across frames on the GPU) to converge in fewer iterations. PhysX 5.0 extended this to soft bodies and cloth. The GPU path is fastest for scenes with thousands of rigid bodies; below roughly 500 bodies the CPU path wins due to launch overhead.

4 · Havok

Havok (now owned by Microsoft) pioneered multi-threaded constraint solving for consoles. The solver partitions the constraint graph into islands, distributes islands across worker threads, and uses a constraint-graph coloring to avoid write conflicts within an island. Havok's simulation islands (groups of interacting bodies) are the standard approach for limiting solver scope: sleeping islands are skipped entirely. Used in Halo, Assassin's Creed, The Elder Scrolls, and most AAA titles from 2005 to 2020.

15Pitfalls

• Jitter from insufficient iterations. The most visible bug. A stack of boxes vibrates because the solver did not converge. Diagnosis: increase iterations or enable warm starting. If warm starting is already on, check contact persistence (contacts may be flickering in and out of existence).
• 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 drops below a threshold. If the threshold is too aggressive, a body sleeps in mid-air or mid-slide. If too conservative, CPU is wasted on stationary bodies. Common symptom: objects that hover in air after a gentle toss, or objects that vibrate indefinitely because they never cross the sleep threshold.
• 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 physics solver uses floating-point arithmetic, which is not bit-identical across CPU architectures (x86 vs ARM), compiler optimizations (-ffast-math), or instruction reordering. Networked physics that relies on deterministic lockstep will diverge. Fix: use IEEE 754 strict mode, disable fast-math, and validate cross-platform bitwise identity on test inputs.
• Tunneling. Covered in section 13. Fast objects pass through thin geometry. Fix: CCD for flagged bodies.
• Bad inertia tensors. An inertia tensor that is too large (body resists rotation too much) or too small (body spins freely) makes the solver unstable. Common source: importing mesh data and computing inertia from the mesh without accounting for density or scale. Jolt and PhysX both clamp inertia ratios to prevent degenerate tensors.

16What's next

• Position-Based Dynamics (PBD/XPBD). Muller et al. 2007^[7]. Solves constraints by directly moving positions instead of computing impulses. Simpler, more stable for cloth and soft bodies. Extended PBD (XPBD, Macklin et al. 2016) adds compliance for correct stiffness behavior. The dominant approach for real-time cloth in Unreal (Chaos Cloth) and Unity.
• GPU physics. PhysX 5.0's GPU rigid-body and soft-body pipelines. Flex (NVIDIA) for unified particle-based simulation. GPU physics scales to millions of particles but has launch-overhead costs that make it slower than CPU for small scenes.
• Deformable bodies. Finite Element Method (FEM) for volumetric deformation: fracture, bending, tearing. Too expensive for real-time on general meshes, but the coarsened variants (tetrahedral FEM with 50 to 200 tets per object) ship in destruction systems.
• 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, not yet shipping in general-purpose game engines.
• Fluid simulation. SPH (Smoothed Particle Hydrodynamics) for Lagrangian fluids. FLIP/APIC for hybrid particle-grid methods. PIC/FLIP is the standard for offline VFX (Houdini); real-time variants are used for water effects in some titles (Sea of Thieves uses a custom wave simulation, not SPH).

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 the Stormer 1907 integration method for molecular dynamics. The velocity-free position update that preserves energy in 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. Frames contact resolution as a linear complementarity problem. The foundational paper for constraint-based rigid body dynamics.
4. David Baraff. "Physically Based Modeling: Rigid Body Simulation." SIGGRAPH 1997 Course Notes. cs.cmu.edu/~baraff/sigcourse. Comprehensive course notes covering equations of motion, constraint formulation, and contact resolution. The standard academic reference for rigid body dynamics.
5. Erin Catto. "Iterative Dynamics with Temporal Coherence." Game Developers Conference 2005. box2d.org/publications. Introduces the sequential-impulse solver with accumulated clamping and warm starting. The paper that defines the modern game-physics constraint solver.
6. Erin Catto. "Box2D: A 2D Physics Engine for Games." box2d.org. Open-source 2D rigid-body engine. Reference implementation of sequential-impulse dynamics. Used in Angry Birds, Limbo, and hundreds of other titles.
7. Matthias Muller, Bruno Heidelberger, Marcus Hennix, John Ratcliff. "Position Based Dynamics." Proc. Virtual Reality Interactions and Physical Simulations (VRIPhys), 2007. matthias-research.github.io/pages/publications/posBasedDyn.pdf. The PBD paper: constraint solving by direct position modification. Widely adopted for real-time cloth and soft bodies.
8. Erin Catto. "Understanding Constraints." Game Developers Conference 2014. box2d.org/publications. Comprehensive walkthrough of the Jacobian formulation, constraint derivation, and sequential impulse solver. The clearest single source on constraint-based dynamics.
9. Jorrit Rouwe. "Jolt Physics." GitHub repository, 2022. github.com/jrouwe/JoltPhysics. Modern C++ rigid-body engine with lock-free broad phase, SIMD solver, and multi-threaded island processing. Adopted by Guerrilla Games and the Godot community.