Newton’s Method
用于求解方程 f(x)=0f(\boldsymbol{x})=0f(x)=0,其中 fff 是一个非线性函数,很难直接解出零点,于是可以用牛顿法,假设当前的初值为 x0\boldsymbol{x}_0x0,对 fff 泰勒展开可得:
f(x)=f(x0)+∇f(x0)⋅Δx=0
f(\boldsymbol{x})=f(\boldsymbol{x}_0)+\nabla f(\boldsymbol{x}_0)\cdot\Delta\boldsymbol{x}=0
f(x)=f(x0)+∇f(x0)⋅Δx=0
于是就可以用解线性方程组的方法解出 Δx\Delta\boldsymbol{x}Δx,然后更新 x1=x0+Δx\boldsymbol{x}_1=\boldsymbol{x}_0+\Delta\boldsymbol{x}x1=x0+Δx,这个 x1\boldsymbol{x}_1x1 一定比 x0\boldsymbol{x}_0x0 更加接近零点。如此迭代若干次之后就可以无限逼近精确解。
XPBD相对PBD做出的改进
仿真物体的刚度(Stiffness)与时间步长和迭代次数有关,迭代次数越多,或者时间步长越短,约束的刚度就会越强。所以在求解约束的时候有必要考虑时间步长的影响。
从约束(Constraint)的弹性势能(Elastic Energy Potentials)的角度入手,结合牛顿第二定律求解约束,这样可以方便计算每个位置的受力。
XPBD算法流程
首先定义约束函数向量 C=[C1(x),...,Cm(x)]T\boldsymbol{C}=[C_1(\boldsymbol{x}),...,C_m(\boldsymbol{x})]^TC=[C1(x),...,Cm(x)]T,约束刚度系数(逆)矩阵 α\boldsymbol{\alpha}α 为对角线为约束的刚度的倒数的对角矩阵,那么系统的弹性势能表示为
U(x)=12C(x)Tα−1C(x),(1)
U(\boldsymbol{x}) = \frac{1}{2}\boldsymbol{C}(\boldsymbol{x})^T \boldsymbol{\alpha}^{-1} \boldsymbol{C}(\boldsymbol{x}),
\tag{1}
U(x)=21C(x)Tα−1C(x),(1)
这里的弹性势能类似弹簧的弹性势能,C(x)\boldsymbol{C}(\boldsymbol{x})C(x)可以看成是偏离平衡位置的距离。然后约定函数梯度为行向量,那么系统受力
f=−∇xUT=−∇CTα−1C,(2)
\boldsymbol{f}=-\nabla_\boldsymbol{x}U^T
= -\nabla \boldsymbol{C}^T\boldsymbol{\alpha}^{-1}\boldsymbol{C},
\tag{2}
f=−∇xUT=−∇CTα−1C,(2)
根据牛顿第二定律的离散形式,我们可以得到
M(xn+1−2xn+xn−1Δt2)=−∇UT(xn+1)(3)
\boldsymbol{M}(\frac{\boldsymbol{x}^{n+1}-2\boldsymbol{x}^n+\boldsymbol{x}^{n-1}}{\Delta t^2}) = -\nabla U^T(\boldsymbol{x}^{n+1})
\tag{3}
M(Δt2xn+1−2xn+xn−1)=−∇UT(xn+1)(3)
这里 M\boldsymbol{M}M 是质量矩阵。然后引入拉格朗日乘子 λ∈Rm×1\boldsymbol{\lambda}\in \mathbb{R}^{m\times 1}λ∈Rm×1 ,并且令
λ=−α~−1C(x)(4)
\boldsymbol{\lambda}=-\boldsymbol{\tilde \alpha}^{-1}\boldsymbol{C}(\boldsymbol{x})
\tag{4}
λ=−α~−1C(x)(4)
其中 α~=αΔt2\boldsymbol{\tilde \alpha}=\frac{\boldsymbol{\alpha}}{\Delta t^2}α~=Δt2α 。再令 x~=2xn−xn−1\boldsymbol{\tilde x}=2\boldsymbol{x}^n-\boldsymbol{x}^{n-1}x~=2xn−xn−1,表示惯性项,于是 (1)(3) 式转变为了求解
g=M(xn+1−x~)−∇C(xn+1)Tλn+1=0h=C(xn+1)+α~λn+1=0(5)
\begin{aligned}
\boldsymbol{g} &= \boldsymbol{M}(\boldsymbol{x}^{n+1}-\boldsymbol{\tilde x})-\nabla\boldsymbol{C}(\boldsymbol{x}^{n+1})^T\boldsymbol{\lambda}^{n+1} &= \boldsymbol{0} \\
\boldsymbol{h} &= \boldsymbol{C}(\boldsymbol{x}^{n+1})+\boldsymbol{\tilde \alpha}\boldsymbol{\lambda}^{n+1} &= \boldsymbol{0}
\end{aligned}
\tag{5}
gh=M(xn+1−x~)−∇C(xn+1)Tλn+1=C(xn+1)+α~λn+1=0=0(5)
上面两个等式不是线性的,因为 C\boldsymbol{C}C 是个未知的函数,所以需要用牛顿法来逼近真实解。我们将上面两式一阶泰勒展开得到
g(xi,λi)+KΔx−∇C(xi)TΔλ=0h(xi,λi)+∇C(xi)Δx+α~Δλ=0(6)
\begin{aligned}
\boldsymbol{g}(\boldsymbol{x}_i,\boldsymbol{\lambda}_i) + \boldsymbol{K}\Delta\boldsymbol{x} - \nabla\boldsymbol{C}(\boldsymbol{x}_i)^T \Delta\boldsymbol{\lambda} &= \boldsymbol{0}
\\
\boldsymbol{h}(\boldsymbol{x}_i, \boldsymbol{\lambda}_i) + \nabla\boldsymbol{C}(\boldsymbol{x}_i)\Delta\boldsymbol{x}+\boldsymbol{\tilde \alpha}\Delta\boldsymbol{\lambda} &= \boldsymbol{0}
\end{aligned}
\tag{6}
g(xi,λi)+KΔx−∇C(xi)TΔλh(xi,λi)+∇C(xi)Δx+α~Δλ=0=0(6)
其中 K=∂g∂xi=M−∇2C(xi)λi\boldsymbol{K}=\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{x}_i}=\boldsymbol{M}-\nabla^2\boldsymbol{C}(\boldsymbol{x}_i)\boldsymbol{\lambda}_iK=∂xi∂g=M−∇2C(xi)λi 。
然后论文里给出了两个化简:
- 令 K=M\boldsymbol{K}=\boldsymbol{M}K=M,这直接抛弃了约束刚度以及约束的海森矩阵;
- 令 g(xi,λi)=0\boldsymbol{g}(\boldsymbol{x}_i,\boldsymbol{\lambda}_i)=\boldsymbol{0}g(xi,λi)=0,在第一轮迭代(x0=x~\boldsymbol{x}_0=\boldsymbol{\tilde x}x0=x~, λ0=0\boldsymbol{\lambda}_0=\boldsymbol{0}λ0=0)是准确的,后面的话好像说这个 $\boldsymbol{g} $ 的变化幅度很小,所以可以认为等于 0 。
化简后,(6)式演变成
[M−∇CT(xi)∇C(xi)α~][ΔxΔλ]=−[0h(xi,λi)](7)
\left[ \begin{matrix}
\boldsymbol{M} & -\nabla\boldsymbol{C}^T(\boldsymbol{x}_i) \\
\nabla\boldsymbol{C}(\boldsymbol{x}_i) & \boldsymbol{\tilde \alpha}
\end{matrix} \right]
\left[ \begin{matrix}
\Delta\boldsymbol{x} \\
\Delta\boldsymbol{\lambda}
\end{matrix} \right]
= -
\left[ \begin{matrix}
\boldsymbol{0} \\
\boldsymbol{h}(\boldsymbol{x}_i, \boldsymbol{\lambda}_i)
\end{matrix} \right]
\tag{7}
[M∇C(xi)−∇CT(xi)α~][ΔxΔλ]=−[0h(xi,λi)](7)
然后可以用对 M\boldsymbol{M}M 求舒尔补的方法求解下面的线性方程组
[∇C(xi)M−1∇C(xi)T+α~]Δλ=−C(xi)−α~λi(8)
[\nabla\boldsymbol{C}(\boldsymbol{x}_i)\boldsymbol{M}^{-1}\nabla\boldsymbol{C}(\boldsymbol{x}_i)^T + \boldsymbol{\tilde \alpha}]\Delta\boldsymbol{\lambda} = -\boldsymbol{C}(\boldsymbol{x}_i)-\boldsymbol{\tilde\alpha}\boldsymbol{\lambda}_i
\tag{8}
[∇C(xi)M−1∇C(xi)T+α~]Δλ=−C(xi)−α~λi(8)
得到增量 Δλ\Delta\boldsymbol{\lambda}Δλ,然后计算
Δx=M−1∇C(xi)TΔλ(9)
\Delta\boldsymbol{x} = \boldsymbol{M}^{-1}\nabla\boldsymbol{C}(\boldsymbol{x}_i)^T\Delta\boldsymbol{\lambda}
\tag{9}
Δx=M−1∇C(xi)TΔλ(9)
最终完成一次迭代
xi+1=xi+Δxλi+1=λi+Δλ(10)
\begin{aligned}
\boldsymbol{x}_{i+1} &= \boldsymbol{x}_i+\Delta\boldsymbol{x} \\
\boldsymbol{\lambda}_{i+1} &= \boldsymbol{\lambda}_i+\Delta\boldsymbol{\lambda}
\end{aligned}
\tag{10}
xi+1λi+1=xi+Δx=λi+Δλ(10)
(8) 式仍然是一个线性方程组,论文里仿照PBD论文中的 Gauss-Seidel 方法,一个一个地求解约束,求解完一个约束就更新一次位置,于是对于第 jjj 个约束需要计算
Δλj=−Cj(xi)−α~jλij∇CjM−1∇CjT+α~j(11)
\Delta\lambda_j = \frac{-C_j(\boldsymbol{x}_i)-\tilde\alpha_j\lambda_{ij}}{\nabla C_j\boldsymbol{M}^{-1}\nabla C_j^T+\tilde\alpha_j}
\tag{11}
Δλj=∇CjM−1∇CjT+α~j−Cj(xi)−α~jλij(11)
整个 XPBD 算法的一个时间步长的算法流程如下
