XPBD: Extended Position-Based Dynamics

1. Why Plain PBD Wasn't Enough

Classic Position-Based Dynamics (Müller et al., 2007) simulates constraints — distance, bending, volume — by directly projecting particle positions so that a constraint function C(x) = 0 is satisfied, then iterating this projection a fixed number of times per frame (Gauss-Seidel style). It's fast, unconditionally stable, and trivial to implement. The catch: nothing in the projection step corresponds to a real stiffness value in N/m. The "stiffness" you get is an artifact of two unrelated numbers — the number of solver iterations and the timestep — multiplied together.

Halve the timestep or double the iteration count, and a rope that used to stretch visibly suddenly behaves like a rigid rod. This made PBD unusable for anything that needed to match real material behaviour, and it meant tuning a cloth sim on one machine could look completely different on another with a different frame rate.

2. Compliance and Lagrange Multipliers

XPBD reformulates the projection as the solution of a constrained optimization problem straight out of Lagrangian mechanics. Each constraint C_j(x) = 0 gets a compliance value α_j — the inverse of physical stiffness, measured in m/N (metres per Newton, i.e. how much the constraint "gives" per unit force):

The total potential energy of a constraint is U(x) = ½ · α⁻¹ · C(x)². Minimizing this energy subject to the equations of motion introduces a Lagrange multiplier λ per constraint — physically, λ is the accumulated constraint force impulse. XPBD's key trick is updating λ incrementally alongside the position update, instead of solving the full linear system:

Set α = 0 and this reduces exactly to classic PBD's projection formula — XPBD is a strict superset. Set α > 0 and you get a spring with a physically meaningful, iteration-count-independent stiffness.

3. The XPBD Solver Loop

The full per-frame algorithm interleaves substeps, each of which predicts positions, solves all constraints once (with λ reset to zero at the start of the substep, not the frame), and derives velocity from the position delta:

function xpbdStep(particles, constraints, dt, substeps) {
  const h = dt / substeps; // substep timestep

  for (let s = 0; s < substeps; s++) {
    // 1. Predict positions with external forces (gravity, wind...)
    for (const p of particles) {
      if (p.invMass === 0) continue; // pinned/static
      p.vel.y += -9.8 * h;
      p.prevPos.x = p.pos.x; p.prevPos.y = p.pos.y;
      p.pos.x += p.vel.x * h;
      p.pos.y += p.vel.y * h;
    }

    // 2. Reset Lagrange multipliers for this substep
    for (const c of constraints) c.lambda = 0;

    // 3. Solve every constraint once (Gauss-Seidel sweep)
    for (const c of constraints) solveConstraint(c, h);

    // 4. Derive velocity from the position change (PBD's core idea)
    for (const p of particles) {
      if (p.invMass === 0) continue;
      p.vel.x = (p.pos.x - p.prevPos.x) / h;
      p.vel.y = (p.pos.y - p.prevPos.y) / h;
    }
  }
}

Notice there is no separate "velocity solve" pass like in impulse-based engines — velocity is always a byproduct of the position change divided by h. This is what makes PBD and XPBD unconditionally stable: you can never diverge to infinite velocity, because velocity is bounded by how far a position can actually move in one substep.

4. Deriving the Distance Constraint

The workhorse constraint for cloth and ropes is the distance constraint, which keeps two particles at rest length L₀:

function solveDistanceConstraint(c, h) {
  const { p1, p2, restLength, compliance } = c;
  const dx = p1.pos.x - p2.pos.x;
  const dy = p1.pos.y - p2.pos.y;
  const dist = Math.hypot(dx, dy);
  if (dist < 1e-9) return;

  const C = dist - restLength;
  const nx = dx / dist, ny = dy / dist;

  const alphaTilde = compliance / (h * h);
  const denom = p1.invMass + p2.invMass + alphaTilde;
  if (denom < 1e-9) return;

  const dLambda = (-C - alphaTilde * c.lambda) / denom;
  c.lambda += dLambda;

  p1.pos.x +=  p1.invMass * dLambda * nx;
  p1.pos.y +=  p1.invMass * dLambda * ny;
  p2.pos.x -=  p2.invMass * dLambda * nx;
  p2.pos.y -=  p2.invMass * dLambda * ny;
}

The exact same skeleton — compute C, compute ∇C per particle, solve for Δλ, apply weighted position deltas — works for bending constraints (dihedral angle between two triangles), volume constraints (tetrahedron signed volume), and area constraints. Only C and its gradient change.

5. Why Substeps, Not Solver Iterations

The original PBD paper fixed the timestep at the frame rate (e.g. 1/60 s) and increased the iteration count to improve convergence and stiffness. XPBD's authors found this converges poorly — a Gauss-Seidel sweep over many constraints propagates corrections only one link per iteration, so a long cloth strip needs dozens of iterations to feel taut.

Instead, XPBD favors many small substeps with a single constraint solve each, over one big frame with many iterations. For the same total compute budget, substepping converges dramatically faster because the compliance term α̃ = α / h² shrinks as h shrinks — smaller substeps make every constraint stiffer and more accurate automatically, without retuning any parameter.

6. Damping Without Energy-Loss Bugs

Naively damping velocity after the solve (v *= 0.98) removes energy non-uniformly and can make simulations frame-rate-dependent again. XPBD instead recommends constraint-based damping: add a small extra position correction proportional to the relative velocity along the constraint gradient, scaled by a damping compliance β and the substep h:

This keeps damping proportional to the actual relative velocity at the constraint (like a physical dashpot/dampener), so a cloth sim run at 20 substeps and one run at 30 substeps settle to the same resting shape at the same rate.

7. Full Implementation: A Hanging Chain

Putting it together — a chain of particles linked by distance constraints, with the first particle pinned (invMass = 0):

class Particle {
  constructor(x, y, invMass = 1) {
    this.pos = { x, y };
    this.prevPos = { x, y };
    this.vel = { x: 0, y: 0 };
    this.invMass = invMass;
  }
}

class DistanceConstraint {
  constructor(p1, p2, compliance = 0) {
    this.p1 = p1; this.p2 = p2;
    this.restLength = Math.hypot(p1.pos.x - p2.pos.x, p1.pos.y - p2.pos.y);
    this.compliance = compliance; // 0 = rigid, higher = stretchier
    this.lambda = 0;
  }
}

// Build a 12-link chain, first link pinned in place
const chain = [];
for (let i = 0; i < 12; i++)
  chain.push(new Particle(i * 0.2, 2, i === 0 ? 0 : 1));

const links = [];
for (let i = 0; i < chain.length - 1; i++)
  links.push(new DistanceConstraint(chain[i], chain[i + 1], 0.0001)); // slightly compliant steel

// Main loop (call once per rendered frame, dt ≈ 1/60)
function simulate(dt) {
  xpbdStep(chain, links, dt, 24); // 24 substeps
}

8. Applications: Cloth, Rods, and Soft Bodies

• Cloth: a grid of distance constraints (structural + shear + bend) with per-constraint compliance gives realistic drape without the "PBD balloon" over-stiffening artifact seen in classic implementations.
• Rods and hair: XPBD extends naturally to orientation-based constraints (using quaternions as generalized coordinates) for bend/twist stiffness in Cosserat rod models — this is how Blender simulates hair strands.
• Soft bodies: tetrahedral meshes with volume- preservation constraints and neo-Hookean strain-energy constraints, solved with the same Δλ update, reproduce FEM-like elastic behaviour at a fraction of the cost.
• Rigid body joints: the compliance formulation generalizes to 6-DOF joints (hinge, ball socket, prismatic), which is why PhysX 5 and Unity's physics stack both migrated their joint solvers to XPBD-style updates.

9. Common Pitfalls

• Resetting λ at the wrong scope: λ must reset to zero once per substep, not once per frame and not once per Gauss-Seidel iteration inside a substep — getting this wrong silently breaks the compliance math.
• Forgetting α̃ = α / h²: using the raw compliance instead of the substep-rescaled value makes stiffness depend on substep count again, defeating the entire point of XPBD.
• Zero compliance with zero denom guard: when both particles in a constraint are pinned (invMass = 0 for both) the denominator is zero — always guard against division by zero before computing Δλ.
• Order-dependent artifacts: Gauss-Seidel sweeps are sequential, so constraint solve order affects the result slightly; for parallel (GPU) solvers, use graph-coloring to group non-conflicting constraints so they can be solved simultaneously without race conditions.