1. Improvement Over PBD
As discussed in this tutorial, the biggest issue with PBD is that its apparent stiffness depends on both the timestep and the number of solver iterations. This happens because the stiffness parameter in classic PBD is not physically meaningful. XPBD addresses this by replacing purely hard geometric constraints with compliant constraints. A scalar compliance parameter is introduced to model an object's finite stiffness. The constraint is then solved by minimizing an energy written in terms of the constraint functions. This makes XPBD behave similarly to an implicit FEM method, which has two main benefits. First, the simulated stiffness no longer depends directly on the timestep size or the number of iterations per timestep, making it possible to simulate elastic and dissipative materials more consistently. Compared with PBD, the stiffness parameter has a more physical meaning. Second, the energy-based formulation provides an estimate of the elastic force induced by the constraints, allowing the user to simulate a wider range of phenomena.
2. Algorithm
Starting from Newton's second law, XPBD writes the internal force as the negative gradient of an elastic potential energy:
\[ \mathbf{M}\ddot{\mathbf{x}} = -\nabla U^T(\mathbf{x}) \]Here, $\mathbf{M}$ is the mass matrix, $\mathbf{x}$ is the position, and $U(\mathbf{x})$ denotes the constraint-based elastic potential energy.
In XPBD, this potential energy is defined as:
\[ U(\mathbf{x}) = \frac{1}{2} C(\mathbf{x})^T \mathbf{\alpha}^{-1} C(\mathbf{x}) \]Using a central-difference discretization, the acceleration at the next timestep is approximated as:
\[ \ddot{\mathbf{x}}^{n+1} \approx \frac{\mathbf{x}^{n+1} - 2\mathbf{x}^n + \mathbf{x}^{n-1}}{\Delta t^2} \]2.1 Governing Equations
The elastic force is computed by taking the derivative of the potential energy with respect to position $\mathbf{x}$:
\[ \begin{aligned} \mathbf{f} &= -\nabla_{\mathbf{x}} U(\mathbf{x}) \\ &= -\nabla \mathbf{C}^{T}(\mathbf{x})\boldsymbol{\alpha}^{-1}\mathbf{C}(\mathbf{x}) \\ &= -\nabla \mathbf{C}^{T}(\mathbf{x})\boldsymbol{\lambda}. \end{aligned} \]where
\[ \begin{aligned} \boldsymbol{\lambda} &= \tilde{\boldsymbol{\alpha}}^{-1}\mathbf{C}(\mathbf{x}) \\ &= -\frac{\boldsymbol{\alpha}}{\Delta t^2}\mathbf{C}(\mathbf{x}) \end{aligned} \]is the vector of constraint multipliers.
The force-equilibrium equation is
\[ \begin{aligned} \mathbf{M}(\mathbf{x}^{n+1}-\tilde{\mathbf{x}}) - \nabla \mathbf{C}^{T}(\mathbf{x}^{n+1})\boldsymbol{\lambda}^{n+1} &= \mathbf{0}. \end{aligned} \]where
\[ \begin{aligned} \tilde{\mathbf{x}} &= 2\mathbf{x}^{n} - \mathbf{x}^{n-1}. \end{aligned} \]Rearranging the constraint-multiplier equation gives
\[ \begin{aligned} \tilde{\boldsymbol{\alpha}}\boldsymbol{\lambda} &= -\mathbf{C}(\mathbf{x}^{n+1}) \\ \Rightarrow \mathbf{C}(\mathbf{x}^{n+1}) + \tilde{\boldsymbol{\alpha}}\boldsymbol{\lambda} &= \mathbf{0}. \end{aligned} \]Therefore, we have the following set of constraint equations, omitting the superscript $n+1$ for simplicity:
\[ \left\{ \begin{aligned} \mathbf{M}(\mathbf{x}-\tilde{\mathbf{x}}) - \nabla \mathbf{C}^{T}(\mathbf{x})\boldsymbol{\lambda} &= \mathbf{0}, \\ \mathbf{C}(\mathbf{x}) + \tilde{\boldsymbol{\alpha}}\boldsymbol{\lambda} &= \mathbf{0}. \end{aligned} \right. \]We can write the equations above as
\[ \left\{ \begin{aligned} \mathbf{g}(\mathbf{x},\boldsymbol{\lambda}) &= \mathbf{0}, \\ \mathbf{h}(\mathbf{x},\boldsymbol{\lambda}) &= \mathbf{0}. \end{aligned} \right. \]At iteration $i$, we approximate the system using a first-order Taylor expansion:
\[ \left\{ \begin{aligned} \mathbf{g}(\mathbf{x}_i+\Delta \mathbf{x},\boldsymbol{\lambda}_i+\Delta \boldsymbol{\lambda}) &= \mathbf{g}(\mathbf{x}_i,\boldsymbol{\lambda}_i) + \frac{\partial \mathbf{g}}{\partial \mathbf{x}}\Delta \mathbf{x} + \frac{\partial \mathbf{g}}{\partial \boldsymbol{\lambda}}\Delta \boldsymbol{\lambda}=\mathbf{0}, \\ \mathbf{h}(\mathbf{x}_i+\Delta \mathbf{x},\boldsymbol{\lambda}_i+\Delta \boldsymbol{\lambda}) &= \mathbf{h}(\mathbf{x}_i,\boldsymbol{\lambda}_i) + \frac{\partial \mathbf{h}}{\partial \mathbf{x}}\Delta \mathbf{x} + \frac{\partial \mathbf{h}}{\partial \boldsymbol{\lambda}}\Delta \boldsymbol{\lambda}=\mathbf{0}. \end{aligned} \right. \]Rearranging this linearized system gives
\[ \begin{aligned} \begin{bmatrix} \frac{\partial \mathbf{g}}{\partial \mathbf{x}} & \frac{\partial \mathbf{g}}{\partial \boldsymbol{\lambda}} \\ \frac{\partial \mathbf{h}}{\partial \mathbf{x}} & \frac{\partial \mathbf{h}}{\partial \boldsymbol{\lambda}} \end{bmatrix} \begin{bmatrix} \Delta \mathbf{x} \\ \Delta \boldsymbol{\lambda} \end{bmatrix} &= \begin{bmatrix} -\mathbf{g}(\mathbf{x}_i,\boldsymbol{\lambda}_i) \\ -\mathbf{h}(\mathbf{x}_i,\boldsymbol{\lambda}_i) \end{bmatrix}. \end{aligned} \]2.2 Simplifications
Solving this system with a fixed-point iteration is equivalent to an implicit FEM method. A large linear system would need to be built and solved every iteration, which is very expensive. However, we can make two simplifications to avoid an implicit solve.
First, notice that
\[ \begin{aligned} \frac{\partial \mathbf{g}}{\partial \mathbf{x}} &= \mathbf{M} - \sum_j \lambda_j \nabla^2 C_j(\mathbf{x}), \\ \frac{\partial \mathbf{g}}{\partial \boldsymbol{\lambda}} &= -\nabla \mathbf{C}^{T}(\mathbf{x}), \\ \frac{\partial \mathbf{h}}{\partial \mathbf{x}} &= \nabla \mathbf{C}(\mathbf{x}), \\ \frac{\partial \mathbf{h}}{\partial \boldsymbol{\lambda}} &= \tilde{\boldsymbol{\alpha}}. \end{aligned} \]$\frac{\partial \mathbf{g}}{\partial \mathbf{x}}$ requires computing the second-order derivative of the constraint $\mathbf{C}(\mathbf{x})$. We approximate it by dropping the higher-order term:
\[ \begin{aligned} \frac{\partial \mathbf{g}}{\partial \mathbf{x}} &\approx \mathbf{M}. \end{aligned} \]Second, we assume that
\[ \begin{aligned} \mathbf{g}(\mathbf{x}_i,\boldsymbol{\lambda}_i) &= \mathbf{0}. \end{aligned} \]Let $\mathbf{x}_0=\tilde{\mathbf{x}}$ and $\boldsymbol{\lambda}_0=\mathbf{0}$ be the initial guesses for the linear system. Then the assumption holds automatically:
\[ \begin{aligned} \mathbf{g}(\mathbf{x}_0,\boldsymbol{\lambda}_0) &= \mathbf{M}(\mathbf{x}_0-\tilde{\mathbf{x}}) - \nabla \mathbf{C}^{T}(\mathbf{x}_0)\boldsymbol{\lambda}_0 \\ &= \mathbf{0}. \end{aligned} \]Suppose at iteration $i$,
\[ \begin{aligned} \mathbf{g}(\mathbf{x}_i,\boldsymbol{\lambda}_i) &= \mathbf{0}, \end{aligned} \]which is
\[ \begin{aligned} \mathbf{M}(\mathbf{x}_i-\tilde{\mathbf{x}}) &= \nabla \mathbf{C}^{T}(\mathbf{x}_i)\boldsymbol{\lambda}_i. \end{aligned} \]By the first approximation, we know
\[ \begin{aligned} \mathbf{M}\Delta \mathbf{x} &= \nabla \mathbf{C}^{T}(\mathbf{x}_i)\Delta \boldsymbol{\lambda}. \end{aligned} \]After one iteration update,
\[ \left\{ \begin{aligned} \mathbf{x}_{i+1} &= \mathbf{x}_i+\Delta \mathbf{x}, \\ \boldsymbol{\lambda}_{i+1} &= \boldsymbol{\lambda}_i+\Delta \boldsymbol{\lambda}. \end{aligned} \right. \]The residual after this update is
\[ \begin{aligned} \mathbf{g}(\mathbf{x}_{i+1},\boldsymbol{\lambda}_{i+1}) &= \mathbf{M}(\mathbf{x}_{i+1}-\tilde{\mathbf{x}}) - \nabla \mathbf{C}^{T}(\mathbf{x}_{i+1})\boldsymbol{\lambda}_{i+1} \\ &= \mathbf{M}(\mathbf{x}_i+\Delta \mathbf{x}-\tilde{\mathbf{x}}) - \nabla \mathbf{C}^{T}(\mathbf{x}_{i+1}) (\boldsymbol{\lambda}_i+\Delta \boldsymbol{\lambda}) \\ &= \mathbf{M}(\mathbf{x}_i-\tilde{\mathbf{x}}) + \mathbf{M}\Delta \mathbf{x} - \nabla \mathbf{C}^{T}(\mathbf{x}_{i+1}) (\boldsymbol{\lambda}_i+\Delta \boldsymbol{\lambda}) \\ &= \nabla \mathbf{C}^{T}(\mathbf{x}_i)\boldsymbol{\lambda}_i + \nabla \mathbf{C}^{T}(\mathbf{x}_i)\Delta \boldsymbol{\lambda} - \nabla \mathbf{C}^{T}(\mathbf{x}_{i+1}) (\boldsymbol{\lambda}_i+\Delta \boldsymbol{\lambda}) \\ &= \left[ \nabla \mathbf{C}^{T}(\mathbf{x}_i) - \nabla \mathbf{C}^{T}(\mathbf{x}_{i+1}) \right] (\boldsymbol{\lambda}_i+\Delta \boldsymbol{\lambda}). \end{aligned} \]If the constraint gradient is constant,
\[ \begin{aligned} \nabla \mathbf{C}(\mathbf{x}_i) - \nabla \mathbf{C}(\mathbf{x}_{i+1}) &= \mathbf{0}, \end{aligned} \]or changes only slightly,
\[ \begin{aligned} \nabla \mathbf{C}(\mathbf{x}_i) - \nabla \mathbf{C}(\mathbf{x}_{i+1}) &\approx \mathbf{0}, \end{aligned} \]then
\[ \begin{aligned} \mathbf{g}(\mathbf{x}_{i+1},\boldsymbol{\lambda}_{i+1}) &= \mathbf{0}. \end{aligned} \]2.3 Update
Given these two assumptions, we have
\[ \begin{aligned} \begin{bmatrix} \mathbf{M} & -\nabla \mathbf{C}^{T}(\mathbf{x}_i) \\ \nabla \mathbf{C}(\mathbf{x}_i) & \tilde{\boldsymbol{\alpha}} \end{bmatrix} \begin{bmatrix} \Delta \mathbf{x} \\ \Delta \boldsymbol{\lambda} \end{bmatrix} &= - \begin{bmatrix} \mathbf{0} \\ \mathbf{h}(\mathbf{x}_i,\boldsymbol{\lambda}_i) \end{bmatrix}. \end{aligned} \]Solving the above equation yields
\[ \begin{aligned} \left[ \nabla \mathbf{C}(\mathbf{x}_i) \mathbf{M}^{-1} \nabla \mathbf{C}^{T}(\mathbf{x}_i) + \tilde{\boldsymbol{\alpha}} \right] \Delta \boldsymbol{\lambda} &= - \mathbf{C}(\mathbf{x}_i) - \tilde{\boldsymbol{\alpha}}\boldsymbol{\lambda}_i. \end{aligned} \]For a single constraint, the XPBD multiplier update is
\[ \begin{aligned} \Delta \lambda_j &= \frac{ -C_j(\mathbf{x}_i) - \tilde{\alpha}_j \lambda_j }{ \nabla C_j \mathbf{M}^{-1} \nabla C_j^{T} + \tilde{\alpha}_j }. \end{aligned} \]The position update is
\[ \begin{aligned} \Delta \mathbf{x} &= \mathbf{M}^{-1}\nabla C_j^{T}\Delta \lambda_j. \end{aligned} \]3. Interactive Preview
Here is an interactive demo running an XPBD simulation in the browser. You can adjust the compliance, resolution, timestep, iteration count, and total simulation time. Press Run to compute the simulation, then use the progress slider to inspect the generated frames.