Position-Based Dynamics and XPBD

Most solvers compute forces, integrate them into velocities, then integrate those into positions. Position-based dynamics throws out the middle steps: guess where the particles want to go, then push them until the constraints hold.

Skip the forces

Start from a symplectic-Euler prediction of the unconstrained positions,

$$\tilde{\mathbf{x}} = \mathbf{x}^{n} + h\,\mathbf{v}^{n} + h^{2} M^{-1}\mathbf{f}_{\text{ext}},$$

then move the predicted state so every constraint $C(\mathbf{x}) = 0$ is satisfied, and recover velocity from the position change,

$$\mathbf{v}^{n+1} = \frac{\mathbf{x}^{n+1} - \mathbf{x}^{n}}{h}.$$

PBD does the projection with a fixed number of Gauss-Seidel passes.¹ It is fast and unconditionally stable, which is why cloth and ropes in games run on it. The catch is the stiffness. A stiff constraint is one you correct hard; PBD "corrects hard" by iterating more, so the effective stiffness rides on the iteration count and the timestep. The same cloth feels limp at 5 iterations and board-stiff at 50, with no physical number you can point to.

Compliance makes stiffness a number

XPBD fixes this by giving each constraint a compliance $\alpha$, the inverse of its stiffness, and a Lagrange multiplier $\lambda$ that accumulates the constraint force.² Projecting the predicted positions against the constraint and its compliance is one Newton step on the coupled system

$$\begin{bmatrix} M & -\nabla C^{\top} \\ \nabla C & \tilde{\alpha} \end{bmatrix} \begin{bmatrix} \Delta\mathbf{x} \\ \Delta\lambda \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ -C - \tilde{\alpha}\,\lambda \end{bmatrix}, \qquad \tilde{\alpha} = \frac{\alpha}{h^{2}}.$$

Eliminating $\Delta\mathbf{x} = M^{-1}\nabla C^{\top}\,\Delta\lambda$ leaves a scalar update per constraint,

$$\Delta\lambda = \frac{-C - \tilde{\alpha}\,\lambda}{\nabla C^{\top} M^{-1}\nabla C + \tilde{\alpha}}, \qquad \Delta\mathbf{x}_{i} = w_{i}\,\nabla_{i}C\,\Delta\lambda,$$

where $w_{i}$ is the inverse mass of particle $i$. The scaling $\tilde{\alpha} = \alpha/h^{2}$ is the whole trick: it cancels the timestep out of the update, so $\alpha$ becomes a physical, step- and iteration-independent quantity. As $\alpha \to 0$ the denominator collapses to $\nabla C^{\top} M^{-1}\nabla C$ and you recover rigid PBD; a larger $\alpha$ enlarges the denominator, so the correction is deliberately partial and the constraint is allowed to stretch.

One step by hand

Take a single distance constraint between two particles. Particle $i$ is pinned, so $w_{i} = 0$, at $\mathbf{x}_{i} = (0, 0)$. Particle $j$ is free, $w_{j} = 1$, at $\mathbf{x}_{j} = (0, -1.2)$. The rest length is $L_{0} = 1.0$ and $h = 1/60$. Take $\alpha = 0$ first.

The constraint value is the current length minus the rest length,

$$C = \lVert \mathbf{x}_{i} - \mathbf{x}_{j} \rVert - L_{0} = 1.2 - 1.0 = 0.2 .$$

The gradient points along the segment,

$$\mathbf{n} = \frac{\mathbf{x}_{i} - \mathbf{x}_{j}}{\lVert \mathbf{x}_{i} - \mathbf{x}_{j} \rVert} = (0, 1).$$

With $\alpha = 0$ the scaled compliance is $\tilde{\alpha} = 0$, and since $\nabla C^{\top} M^{-1}\nabla C = w_{i} + w_{j} = 0 + 1 = 1$, the denominator is $w_{i} + w_{j} + \tilde{\alpha} = 1$. Starting from $\lambda = 0$,

$$\Delta\lambda = \frac{-C - \tilde{\alpha}\,\lambda}{w_{i} + w_{j} + \tilde{\alpha}} = \frac{-0.2 - 0}{1} = -0.2 .$$

The pinned particle does not move ($w_{i} = 0$). The free particle moves by

$$\Delta\mathbf{x}_{j} = -w_{j}\,\mathbf{n}\,\Delta\lambda = -(1)\,(0, 1)\,(-0.2) = (0,\ 0.2),$$

so $\mathbf{x}_{j}$ rises from $(0, -1.2)$ to $(0, -1.0)$. Recompute: the length is now exactly $1.0$, so $C = 0$. One step, constraint solved.

Now let $\alpha > 0$, say $\alpha = 10^{-3}$, keeping $h = 1/60$. Then $\tilde{\alpha} = \alpha/h^{2} = 10^{-3} \cdot 3600 = 3.6$, and

$$\Delta\lambda = \frac{-0.2 - 0}{0 + 1 + 3.6} = \frac{-0.2}{4.6} \approx -0.0435, \qquad \Delta\mathbf{x}_{j} = (0,\ 0.0435).$$

The particle only climbs to $\mathbf{x}_{j} \approx (0, -1.1565)$, leaving $C \approx 0.157$ of stretch after the step. That residual is not solver error; it is the compliant spring holding its length under load, which is exactly what a finite stiffness is supposed to do.

The rope, live

The demo below is that same update run on a chain of twelve particles, top one pinned, gravity pulling the rest down. Each frame runs $S$ substeps, and each substep does one projection sweep with the formula above.

XPBD rope · compliance vs substepsinteractive

compliance α3.5e-4

substeps S8

rope length1.980

stretch+0.0%

compliance α = 3.5e-4 (stiff → soft)substeps S = 8

PauseResetSwing: off

Slide S from 1 to 20: the stretch barely moves, because α, not the iteration count, sets the stiffness. Slide α right and the rope sags below the rest-length line.

The stiff end holds its length: the tip sits on the dashed rest-length line. The soft end visibly stretches under gravity, settling where the compliant tension balances the weight below. The point worth watching is the substep slider: in plain PBD, more substeps means a stiffer rope, but here $\alpha$ owns the stiffness and $S$ only buys accuracy, so the rope's droop hardly changes as you sweep it.

That link between compliance and the timestep is also the bridge to the rest of the solver story. XPBD with one iteration per substep is a quasi-Newton step on the same implicit-Euler system the constraint solvers target directly, which is why it slots so naturally into the massively parallel, substepped engines coming up next.

Footnotes

1. Müller, Heidelberger, Hennix & Ratcliff, Position Based Dynamics (2007). ↩

2. Macklin, Müller & Chentanez, XPBD: Position-Based Simulation of Compliant Constrained Dynamics (MiG 2016). ↩