XPBD学习笔记

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+12xn+xn1)=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~=2xnxn1,表示惯性项,于是 (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+1x~)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ΔxC(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=xig=M2C(xi)λi

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

  1. K=M\boldsymbol{K}=\boldsymbol{M}K=M,这直接抛弃了约束刚度以及约束的海森矩阵;
  2. 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} [MC(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)M1C(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=M1C(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=CjM1CjT+α~jCj(xi)α~jλij(11)
整个 XPBD 算法的一个时间步长的算法流程如下

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值