XPBD 公式推导

第一步：问题设定

我们有一个物理系统，包含 $N$ 个粒子。每个粒子的位置为 $x_i\in {\mathbb{R}}^3$，质量为 $m_i$。系统的总位置向量记作：

$$x=(x_1,x_2,\dots ,x_N)\in {\mathbb{R}}^{3N}$$

质量矩阵为对角矩阵：

$$M=\text{diag}(m_1{I}_3,m_2{I}_3,\dots ,m_N{I}_3)$$

系统受到一组约束，但我们不要求约束被严格满足。相反，我们给每个约束赋予一个弹性势能。对于单个约束 (C(x) \in \mathbb{R})，其势能定义为：

$$U(x)=\frac{1}{2}C(x{)}^T{\alpha }^{-1}C(x)$$

其中 $\alpha$ 是柔度（compliance），物理意义是刚度的倒数，单位是米/牛顿。对于多个约束，总势能是各约束势能之和。

---

第二步：隐式时间积分

我们采用隐式欧拉法进行时间步进。从时间 $t$ 到 $t+\Delta t$，运动方程为：

$$M\frac{x^{n+1}-2x^n+x^{n-1}}{\Delta t^2}=-\nabla U(x^{n+1})$$

为了方便，定义预测位置（也叫惯性位置）：

$$\tilde{x}=2x^n-x^{n-1}=x^n+\Delta t{v}^n$$

（这里假设速度是 $v^n=(x^n-x^{n-1})/\Delta t$，所以 $\tilde{x}=x^n+\Delta t{v}^n$。）

代入运动方程，整理得：

$$M(x^{n+1}-\tilde{x})=-\Delta t^2\nabla U(x^{n+1})$$

令 $\Delta x=x^{n+1}-\tilde{x}$，则方程变为：

$$M\Delta x=-\Delta t^2\nabla U(x^{n+1})$$

---

第三步：引入拉格朗日乘子

对于每个约束 $C(x)$，其势能梯度为：

$$\nabla U(x)=\nabla C(x{)}^T{\alpha }^{-1}C(x)$$

代入运动方程：

$$M\Delta x=-\Delta t^2\nabla C(x^{n+1}{)}^T{\alpha }^{-1}C(x^{n+1})$$

令

$$\tilde{\alpha }=\frac{\alpha }{\Delta t^2}$$

注意：$\tilde{\alpha }$ 是时间尺度调整后的柔度，它的引入是 XPBD 能够解耦时间步长的关键。

现在，我们定义一个拉格朗日乘子向量（对于多约束系统）：

$$\lambda =-{\tilde{\alpha }}^{-1}C(x^{n+1})$$

注意：这是定义，不是推导。它的作用是把非线性方程拆成两个联立的线性方程。

根据这个定义，我们有：

$$C(x^{n+1})+\tilde{\alpha }\lambda =0$$

并且运动方程可以改写为：

$$M\Delta x-\nabla C(x^{n+1}{)}^T\lambda =0$$

---

第四步：联立方程组

我们现在有了一个关于 $\Delta x$ 和 $\lambda$ 的联立方程组：

$$\{\begin{array}{l}M\Delta x-\nabla C(x^{n+1}{)}^T\lambda =0\\ C(x^{n+1})+\tilde{\alpha }\lambda =0\end{array}$$

注意这里 $x^{n+1}=\tilde{x}+\Delta x$。

---

第五步：线性化（牛顿迭代的一步）

方程组是非线性的，因为 $C(x^{n+1})$ 和 $\nabla C(x^{n+1})$ 都依赖于未知的 $x^{n+1}$。我们采用单步牛顿法来近似求解。

将 $C(x^{n+1})$ 在预测位置 $\tilde{x}$ 处一阶泰勒展开：

$$C(x^{n+1})\approx C(\tilde{x})+\nabla C(\tilde{x})\Delta x$$

将梯度也近似为在 $\tilde{x}$ 处的值：

$$\nabla C(x^{n+1})\approx \nabla C(\tilde{x})$$

代入联立方程组：

$$\{\begin{array}{l}M\Delta x-\nabla C^T\lambda =0\\ C(\tilde{x})+\nabla C\Delta x+\tilde{\alpha }\lambda =0\end{array}$$

其中 $\nabla C=\nabla C(\tilde{x})$，$C=C(\tilde{x})$。

---

第六步：求解 (\Delta \lambda) 和 (\Delta x)

从第一个方程解出 $\Delta x$：

$$\Delta x=M^{-1}\nabla C^T\lambda$$

代入第二个方程：

$$C+\nabla C(M^{-1}\nabla C^T\lambda )+\tilde{\alpha }\lambda =0$$ $$C+(\nabla C{M}^{-1}\nabla C^T+\tilde{\alpha })\lambda =0$$

解得：

$$\lambda =-\frac{C}{\nabla C{M}^{-1}\nabla C^T+\tilde{\alpha }}$$

然后代回求 $\Delta x$：

$$\Delta x=M^{-1}\nabla C^T\lambda$$

---

第七步：迭代形式（XPBD 的核心）

上面的推导给出的是单次牛顿步的解。但在实际模拟中，我们需要多次迭代来收敛到非线性方程组的真解。更重要的是，XPBD 的 $\lambda$ 是累积的，它代表从时间步开始以来累积的总约束冲量。

因此，我们采用迭代更新的形式。设当前迭代步的乘子为 $\lambda$（初始为上一帧结束时保留的值，或零），我们计算增量 $\Delta \lambda$：

由约束方程的一阶展开（注意此时当前位置是 $x$，而不是 $\tilde{x}$）：

$$C(x)+\nabla C\Delta x+\tilde{\alpha }(\lambda +\Delta \lambda )=0$$

代入 $\Delta x=M^{-1}\nabla C^T\Delta \lambda$：

$$C(x)+\nabla C{M}^{-1}\nabla C^T\Delta \lambda +\tilde{\alpha }\lambda +\tilde{\alpha }\Delta \lambda =0$$

整理得：

$$C(x)+\tilde{\alpha }\lambda +(\nabla C{M}^{-1}\nabla C^T+\tilde{\alpha })\Delta \lambda =0$$

解得：

$$\Delta \lambda =\frac{-C(x)-\tilde{\alpha }\lambda }{\nabla C{M}^{-1}\nabla C^T+\tilde{\alpha }}$$

然后更新：

$$\lambda \leftarrow \lambda +\Delta \lambda$$ $$\Delta x=M^{-1}\nabla C^T\Delta \lambda$$ $$x\leftarrow x+\Delta x$$

---

第八步：完整算法流程

对于每一帧：

1. 对于所有粒子，计算预测位置 $\tilde{x}=x+\Delta tv+\Delta t^2{M}^{-1}{f}_{\text{ext}}$
2. 初始化所有约束的 $\lambda =0$（或者保留上一帧的值）
3. 重复若干次（求解迭代）：
   • 对于每个约束：
     • 计算 $C(x)$ 和 $\nabla C(x)$
     • 计算 $\Delta \lambda =\frac{-C(x)-\tilde{\alpha }\lambda }{\nabla C{M}^{-1}\nabla C^T+\tilde{\alpha }}$
     • 更新 $\lambda \leftarrow \lambda +\Delta \lambda$
     • 计算 $\Delta x=M^{-1}\nabla C^T\Delta \lambda$
     • 更新 $x\leftarrow x+\Delta x$
4. 更新速度：$v=(x-x_{\text{old}})/\Delta t$
5. 更新位置：$x_{\text{old}}\leftarrow x$

---

推导中的关键点总结

步骤	关键操作	物理/数学意义
1	$U=\frac{1}{2}{C}^T{\alpha }^{-1}C$	定义弹性势能，赋予约束物理刚度
2	$\tilde{\alpha }=\alpha /\Delta t^2$	时间步长解耦，保证材料属性一致
3	$\lambda =-{\tilde{\alpha }}^{-1}C(x)$	引入对偶变量，将非线性方程拆为联立方程
4	泰勒展开 $C(x+\Delta x)\approx C(x)+\nabla C\Delta x$	线性化，使问题可解
5	迭代更新 $\lambda \leftarrow \lambda +\Delta \lambda$	累积乘子，解耦迭代次数与刚度
6	分母加 $\tilde{\alpha }$	正则化项，防止柔度为零时除零

这就是 XPBD 从物理原理到数值算法的完整推导链条。