XPBD学习笔记

Newton’s Method

用于求解方程 $f(\mathbf{x})=0$，其中 $f$ 是一个非线性函数，很难直接解出零点，于是可以用牛顿法，假设当前的初值为 ${\mathbf{x}}_0$，对 $f$ 泰勒展开可得：
$$f(\mathbf{x})=f({\mathbf{x}}_0)+\nabla f({\mathbf{x}}_0)\cdot \Delta \mathbf{x}=0$$
于是就可以用解线性方程组的方法解出 $\Delta \mathbf{x}$，然后更新 ${\mathbf{x}}_1={\mathbf{x}}_0+\Delta \mathbf{x}$，这个 ${\mathbf{x}}_1$ 一定比 ${\mathbf{x}}_0$ 更加接近零点。如此迭代若干次之后就可以无限逼近精确解。

XPBD相对PBD做出的改进

仿真物体的刚度（Stiffness）与时间步长和迭代次数有关，迭代次数越多，或者时间步长越短，约束的刚度就会越强。所以在求解约束的时候有必要考虑时间步长的影响。

从约束（Constraint）的弹性势能（Elastic Energy Potentials）的角度入手，结合牛顿第二定律求解约束，这样可以方便计算每个位置的受力。

XPBD算法流程

首先定义约束函数向量 $\mathbf{C}=[C_1(\mathbf{x}),...,C_m(\mathbf{x}){]}^T$，约束刚度系数（逆）矩阵 $\mathbf{\alpha }$ 为对角线为约束的刚度的倒数的对角矩阵，那么系统的弹性势能表示为
$$U(\mathbf{x})=\frac{1}{2}\mathbf{C}(\mathbf{x}{)}^T{\mathbf{\alpha }}^{-1}\mathbf{C}(\mathbf{x}),\text{\hspace{1em}(1)}$$
这里的弹性势能类似弹簧的弹性势能，$\mathbf{C}(\mathbf{x})$可以看成是偏离平衡位置的距离。然后约定函数梯度为行向量，那么系统受力
$$\mathbf{f}=-{\nabla }_{\mathbf{x}}{U}^T=-\nabla {\mathbf{C}}^T{\mathbf{\alpha }}^{-1}\mathbf{C},\text{\hspace{1em}(2)}$$
根据牛顿第二定律的离散形式，我们可以得到
$$\mathbf{M}(\frac{{\mathbf{x}}^{n+1}-2{\mathbf{x}}^n+{\mathbf{x}}^{n-1}}{\Delta t^2})=-\nabla U^T({\mathbf{x}}^{n+1})\text{\hspace{1em}(3)}$$
这里 $\mathbf{M}$ 是质量矩阵。然后引入拉格朗日乘子 $\mathbf{\lambda }\in {\mathbb{R}}^{m\times 1}$ ，并且令
$$\mathbf{\lambda }=-{\tilde{\mathbf{\alpha }}}^{-1}\mathbf{C}(\mathbf{x})\text{\hspace{1em}(4)}$$
其中 $\tilde{\mathbf{\alpha }}=\frac{\mathbf{\alpha }}{\Delta t^2}$ 。再令 $\tilde{\mathbf{x}}=2{\mathbf{x}}^n-{\mathbf{x}}^{n-1}$，表示惯性项，于是 (1)(3) 式转变为了求解
$$\begin{array}{rlr}\mathbf{g}& =\mathbf{M}({\mathbf{x}}^{n+1}-\tilde{\mathbf{x}})-\nabla \mathbf{C}({\mathbf{x}}^{n+1}{)}^T{\mathbf{\lambda }}^{n+1}& =0\\ \mathbf{h}& =\mathbf{C}({\mathbf{x}}^{n+1})+\tilde{\mathbf{\alpha }}{\mathbf{\lambda }}^{n+1}& =0\end{array}\text{\hspace{1em}(5)}$$
上面两个等式不是线性的，因为 $\mathbf{C}$ 是个未知的函数，所以需要用牛顿法来逼近真实解。我们将上面两式一阶泰勒展开得到
$$\begin{array}{rl}\mathbf{g}({\mathbf{x}}_i,{\mathbf{\lambda }}_i)+\mathbf{K}\Delta \mathbf{x}-\nabla \mathbf{C}({\mathbf{x}}_i{)}^T\Delta \mathbf{\lambda }& =0\\ \mathbf{h}({\mathbf{x}}_i,{\mathbf{\lambda }}_i)+\nabla \mathbf{C}({\mathbf{x}}_i)\Delta \mathbf{x}+\tilde{\mathbf{\alpha }}\Delta \mathbf{\lambda }& =0\end{array}\text{\hspace{1em}(6)}$$
其中 $\mathbf{K}=\frac{\partial \mathbf{g}}{\partial {\mathbf{x}}_i}=\mathbf{M}-{\nabla }^2\mathbf{C}({\mathbf{x}}_i){\mathbf{\lambda }}_i$ 。

然后论文里给出了两个化简：

1. 令 $\mathbf{K}=\mathbf{M}$，这直接抛弃了约束刚度以及约束的海森矩阵；
2. 令 $\mathbf{g}({\mathbf{x}}_i,{\mathbf{\lambda }}_i)=0$，在第一轮迭代（${\mathbf{x}}_0=\tilde{\mathbf{x}}$, ${\mathbf{\lambda }}_0=0$）是准确的，后面的话好像说这个 $\boldsymbol{g} $ 的变化幅度很小，所以可以认为等于 0 。

化简后，(6)式演变成
$$\left[\begin{array}{cc}\mathbf{M}& -\nabla {\mathbf{C}}^T({\mathbf{x}}_i)\\ \nabla \mathbf{C}({\mathbf{x}}_i)& \tilde{\mathbf{\alpha }}\end{array}\right]\left[\begin{array}{c}\Delta \mathbf{x}\\ \Delta \mathbf{\lambda }\end{array}\right]=-\left[\begin{array}{c}0\\ \mathbf{h}({\mathbf{x}}_i,{\mathbf{\lambda }}_i)\end{array}\right]\text{\hspace{1em}(7)}$$
然后可以用对 $\mathbf{M}$ 求舒尔补的方法求解下面的线性方程组
$$[\nabla \mathbf{C}({\mathbf{x}}_i){\mathbf{M}}^{-1}\nabla \mathbf{C}({\mathbf{x}}_i{)}^T+\tilde{\mathbf{\alpha }}]\Delta \mathbf{\lambda }=-\mathbf{C}({\mathbf{x}}_i)-\tilde{\mathbf{\alpha }}{\mathbf{\lambda }}_i\text{\hspace{1em}(8)}$$
得到增量 $\Delta \mathbf{\lambda }$，然后计算
$$\Delta \mathbf{x}={\mathbf{M}}^{-1}\nabla \mathbf{C}({\mathbf{x}}_i{)}^T\Delta \mathbf{\lambda }\text{\hspace{1em}(9)}$$
最终完成一次迭代
$$\begin{array}{rl}{\mathbf{x}}_{i+1}& ={\mathbf{x}}_i+\Delta \mathbf{x}\\ {\mathbf{\lambda }}_{i+1}& ={\mathbf{\lambda }}_i+\Delta \mathbf{\lambda }\end{array}\text{\hspace{1em}(10)}$$
(8) 式仍然是一个线性方程组，论文里仿照PBD论文中的 Gauss-Seidel 方法，一个一个地求解约束，求解完一个约束就更新一次位置，于是对于第 $j$ 个约束需要计算
$$\Delta {\lambda }_j=\frac{-C_j({\mathbf{x}}_i)-{\tilde{\alpha }}_j{\lambda }_{ij}}{\nabla C_j{\mathbf{M}}^{-1}\nabla C_j^T+{\tilde{\alpha }}_j}\text{\hspace{1em}(11)}$$
整个 XPBD 算法的一个时间步长的算法流程如下