数值积分

对于某一函数$f(x)$，我们需要由 $t$ 时刻的状态求得未来某时间 $t+h$ 的状态。即由 $f(t)$求得 $f(t+h)$ 。为此，我们将$f(t)$ 泰勒展开 $f(x+\Delta x)=f(x)+f’(x)\Delta x+\frac{1}{2}f’’(x) \Delta x^2+\frac{1}{6}f’’’(x) \Delta x^3+\cdots$ 实际中不可能取到无穷项，故以下方法均为其的近似截断。

我们可以通过所取泰勒展开中最高阶数对积分方式进行分类。如欧拉积分只用到了一阶导数，所以是一阶方法。且其误差为 $O(n^2)$ ，步长 $h$ 减半后，误差变为原误差的四分之一。

欧拉积分

经典以折线逼近曲线的趋近式表达

显式欧拉(Explicit Euler)

最简单的积分方式，然而会有数值膨胀

$\vec v(t_1)=\vec v(t_0)+M^{-1}\vec F(t_0)\Delta t $

$x(t_1)=x(t_0)+\vec v(t_0)\Delta t$

隐式欧拉(Implicit Euler)

以未来的状态更新现在，需要解线性系统，且有数值耗散

$ \vec v(t_1)=\vec v(t_0)+M^{-1}\vec F(t_1)\Delta t $

$ x(t_1)=x(t_0)+\vec v(t_1)\Delta t $

半隐式欧拉(Semi_Implicit Euler)

显隐式欧拉的结合，数值稳定

$\vec v(t_1)=\vec v(t_0)+M^{-1}\vec F(t_0)\Delta t$

$x(t_1)=x(t_0)+\vec v(t_1)\Delta t$

跳蛙法(Leap Frog)

半隐式欧拉的变体，将速度和位移错开半个时步，用中点法更新，误差更小

$ \vec v(t_{0.5})=\vec v(t_{-0.5})+M^{-1}\vec F(t_0)\Delta t $

$ x(t_1)=x(t_0)+\vec v(t_{0.5})\Delta t $

速度韦尔莱（Velocity verlet）

引入速度简化韦尔莱积分，与跳蛙法相比改中点法为梯形法，使速度和位置能在时间上同步

$ x(t_1)=x(t_0)+\vec v(t_0)\Delta t+\frac{1}{2}M^{-1}\vec F(t_0)\Delta t^2 $

$ \vec v(t_1)=\vec v(t_0)+\frac{1}{2}M^{-1}(\vec F(t_0)+\vec F(t_1))\Delta t $

龙格库塔(Runge Kutta)

欧拉法其实就是一阶龙格库塔，高阶龙格库塔则是使用中点处导数的近似值，用之前的导数估算后面的导数

二阶龙格库塔(RK2)

$ \vec v(t_1)=\vec v(t_0)+M^{-1}\vec F(t_0)\Delta t $

$ x(t_{0.5})=x(t_0)+\frac{1}{2}\vec v(t_0)\Delta t $

$ \vec v(t_{0.5})=\vec v(t_0)+\frac{1}{2}M^{-1}\vec F(t_{0.5})\Delta t $

$ x(t_1) = x(t_0)+\vec v(t_{0.5})\Delta t $

线性系统

经典解 $Ax=b$问题，直接求逆过于费，所以需要各凭本事进行优化

矩阵分解

将一个复杂的矩阵拆解为数个简单的矩阵的乘积，从而方便求逆。这里以eigen库为例，最常用的有LU分解（分解为下三角阵 $L$和上三角阵 $U$ 的积），LLT分解（要求 $A$ 正定对称，分解为下三角阵 $L$ 和其转置 $L^T$ 的积，速度两倍于LU分解），具体方法为高斯消元

线性迭代

迭代一定次数或到达一定精度后即可停止

雅可比迭代(Jacobi)

对$Ax=b$，裂解$A$为$D-L-R$
$$
A= \left[ \begin{matrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{matrix} \right] \ b= \left[ \begin{matrix} b_1 \\ b_2 \\ b_3 \end{matrix} \right] \ D= \left[ \begin{matrix} a_{11} & 0 & 0 \\ 0 & a_{22} & 0 \\ 0 & 0 & a_{33} \end{matrix} \right]
$$

$$
L= \left[ \begin{matrix} 0 & 0 & 0 \\ a_{21} & 0& 0\\ a_{31} & a_{32} & 0 \end{matrix} \right] \ R= \left[ \begin{matrix} 0 & a_{12} & a_{13} \\ 0 & 0 & a_{23} \\ 0 & 0 & 0 \end{matrix} \right] \ f= \left[ \begin{matrix} \frac{b_1}{a_{11}} \\ \frac{b_2}{a_{22}} \\ \frac{b_3}{a_{33}} \end{matrix} \right]
$$
$$B=I-D^{-1}A，x_{k+1} =Bx_k+f$$

高斯赛德尔迭代(Gauss-Seidel)

雅可比迭代中每计算完一轮，才会将结果作用于下一轮。若是在当轮计算后就及时更新，直观上效率优于雅可比迭代，但雅可比迭代天然适于并行
$$
雅可比 \begin{cases} \begin{matrix} x_1^{k+1}=\frac{1}{a_{11}}(-a_{12}x_2^k-a_{13}x_3^k+b_1)\\ x_2^{k+1}=\frac{1}{a_{22}}(-a_{21}x_1^k-a_{23}x_3^k+b_2)\\ x_3^{k+1}=\frac{1}{a_{33}}(-a_{31}x_1^k-a_{32}x_2^k+b_3) \end{matrix} \end{cases}
$$

$$
高斯赛德尔\begin{cases} \begin{matrix} x_1^{k+1}=\frac{1}{a_{11}}(-a_{12}x_2^k-a_{13}x_3^k+b_1)\\ x_2^{k+1}=\frac{1}{a_{22}}(-a_{21}x_1^{k+1}-a_{23}x_3^k+b_2)\\ x_3^{k+1}=\frac{1}{a_{33}}(-a_{31}x_1^{k+1}-a_{32}x_2^{k+1}+b_3) \end{matrix} \end{cases}
$$

超松弛迭代(SOR)

高斯赛德尔的加速版，引入松弛因子 $\omega$ 调节步长。 $\omega\in(0,2)$ 时收敛， $1<\omega<2$ *为超松弛(过冲加速)，* $0<\omega<1$ *为亚松弛(给本身就快振荡的迭代降火)。最优* $\omega$ *取决于系数矩阵的谱半径，工程上没法精确算，一般在1.0~1.9里手调，PBD里那个stiffness参数本质上就是这个*
$$
x_i^{k+1}=(1-\omega)x_i^k+\frac{\omega}{a_{ii}}\Big(b_i-\sum_{j<i}a_{ij}x_j^{k+1}-\sum_{j>i}a_{ij}x_j^k\Big)
$$

共轭梯度(Conjugate Gradient)

以上方法对一般矩阵都适用，但对大规模对称正定稀疏矩阵(图形学里布料、柔体、隐式积分系统)，雅可比和高斯赛德尔的收敛实在过于慢。CG的思路是把解 $Ax=b$ 转化为最小化 $\phi(x)=\frac{1}{2}x^TAx-b^Tx$ ，沿一组 $A$ -共轭方向迭代，理论上 $n$ 步必收敛(实际浮点误差远早于此就达精度了)
$$
r_0=b-Ax_0,\ p_0=r_0
$$

$$
\alpha_k=\frac{r_k^Tr_k}{p_k^TAp_k},\ x_{k+1}=x_k+\alpha_kp_k,\ r_{k+1}=r_k-\alpha_kAp_k
$$

$$
\beta_k=\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k},\ p_{k+1}=r_{k+1}+\beta_kp_k
$$

实际实现只需存 $x,r,p$ 三个向量加一次 $Ap$ 乘法，对稀疏 $A$ 来说每步 $O(nnz)$ 而非 $O(n^2)$ 。布料/柔体的隐式求解器基本都靠这个

预处理共轭梯度(PCG)

CG的收敛速度由 $A$ 的条件数决定，条件数越大收敛越慢。预处理就是找个易解的近似 $M^{-1}\approx A^{-1}$ 缩小条件数。最简单的Jacobi预处理 $M=diag(A)$ 几乎零成本但立竿见影，更猛的有不完全Cholesky(IC)，工程上权衡构造成本和加速比

$$
z_k=M^{-1}r_k
$$

$$
\alpha_k=\frac{r_k^Tz_k}{p_k^TAp_k},\ \beta_k=\frac{r_{k+1}^Tz_{k+1}}{r_k^Tz_k}
$$

其余流程与CG相同，只是用 $z$ 替代 $r$ 参与方向更新。Eigen的ConjugateGradient默认就带Jacobi预处理，写起来也就一行

小结

这里简单拉表罗列了下，还是主要看矩阵性质和规模选择方法

方法 适用矩阵 单次复杂度 备注
LU/LLT分解 中小规模稠密 $O(n^3)$ 一次 一次分解多次求解，eigen直接调
Jacobi 对角占优 $O(n^2)$ 每步 天然并行但慢
Gauss-Seidel 对角占优 $O(n^2)$ 每步 比Jacobi快一倍，不好并行
SOR 对角占优 $O(n^2)$ 每步 调好$\omega$再快数倍，怕的就是调不好
CG 对称正定 $O(nnz)$ 每步 大规模稀疏首选
PCG 对称正定 $O(nnz)$ 每步 CG工业版，几乎无脑用

总体而言数值积分管时域演化，线性系统管空间约束，二者在隐式方法里耦合——隐式欧拉每步都要解一次形如 $(M-\Delta t^2K)v=…$ 的方程，布料/柔体每帧都要解Poisson或其变体。GAMES103里隐式布料、PBD/XPBD的相关部分基本都是这套组合拳，看似花样多，回头看就这几个数学工具排列组合罢了。