物理引擎中的时间积分方法及求解

介绍常用的时间积分方法，及最终的求解过程。

---

0 物理系统描述

在物理引擎中，借助牛顿第二运动定律对系统进行描述，即

$$\begin{array}{rl}\mathit{f}& =\mathit{f}(\mathit{x})\\ \frac{\partial {\mathit{v}}_i}{\partial t}& =\frac{{\mathit{f}}_i}{m_i}\\ \frac{\partial {\mathit{x}}_i}{\partial t}& ={\mathit{v}}_i\end{array}$$

有时候也会用 $\mathit{q}$ 来表示模型中节点的位置，那么系统描述即为：

$$\begin{array}{rl}\mathit{f}& =\mathit{f}(\mathit{x})\\ \ddot{\mathit{q}}& =\frac{\mathit{f}}{\mathit{M}}\end{array}$$

上述方程就是物理引擎中的 ODE 部分。该方程组的解法主要有显式（Forward Euler、Semi-implicit Euler）、隐式（Backward Euler）等。

---

1、时间积分方法

在仿真计算过程中，已知模型中节点在 $t$ 时刻的位置、速度等信息，进一步求解其在 $t+1$ 时刻的速度、位置。

1.1 显式时间积分（Explicit / Forward Euler）

计算方式如下：

$$\begin{array}{rl}{\mathit{v}}_{t+1}& ={\mathit{v}}_t+\mathrm{\Delta }t\frac{{\mathit{f}}_t}{m}\\ {\mathit{x}}_{t+1}& ={\mathit{x}}_t+\mathrm{\Delta }t{\mathit{v}}_t\end{array}$$

在该方法中，由模型中节点的位置 ${\mathit{x}}_t$ 直接计算得到受力 ${\mathit{f}}_t$ ，进而可直接计算得到 $t+1$ 时刻模型中节点的速度、位置。

1.2 半隐式积分（Explicit / Semi-implicit Euler aka. Symplectic Euler）

计算方式如下：

$$\begin{array}{rl}{\mathit{v}}_{t+1}& ={\mathit{v}}_t+\mathrm{\Delta }t\frac{{\mathit{f}}_t}{m}\\ {\mathit{x}}_{t+1}& ={\mathit{x}}_t+\mathrm{\Delta }t{\mathit{v}}_{t+1}\end{array}$$

在该方法中，同样由模型中节点的位置 ${\mathit{x}}_t$ 直接计算得到受力 ${\mathit{f}}_t$ ，进而可直接计算得到 $t+1$ 时刻模型中节点的速度、位置。

Tips：Forward Euler 和 Semi-implicit Euler 略微不同，在计算 ${\mathit{x}}_{t+1}$ 的时候，一个是用了 ${\mathit{v}}_t$，另一个是用了 ${\mathit{v}}_{t+1}$。现在，Forward Euler 用的较少，Semi-implicit Euler 用的较多一些。虽然，Semi-implicit Euler 只差别了一点点，但是准确性上会有本质性的提升。

1.3 仿真流程（显式积分）

在使用显式（Forward Euler or Semi-implicit Euler）进行仿真的时候，仿真流程有如下几个步骤：

• 计算节点受力 ${\mathit{f}}_t=\mathit{f}({\mathit{x}}_t)$
• 计算新的速度 ${\mathit{v}}_{t+1}={\mathit{v}}_t+\mathrm{\Delta }t\frac{{\mathit{f}}_t}{m}$
• 碰撞检测（此时会更正速度）
• 计算新的位置 ${\mathit{x}}_{t+1}={\mathit{x}}_t+\mathrm{\Delta }t{\mathit{v}}_{t+1}$ （Semi-implicit Euler）

显式时间积分器的性能缺陷： Easy to explore

$$\mathrm{\Delta }t\le c\sqrt{\frac{m}{k}}(c\sim 1)$$

关于稳定性 Stability 和爆炸 Explode 问题：

（略）

1.4 隐式积分（implicit Euler）

计算方式如下：

$$\begin{array}{rl}{\mathit{v}}_{t+1}& ={\mathit{v}}_t+\mathrm{\Delta }t\frac{{\mathit{f}}_{t+1}}{m}\\ {\mathit{x}}_{t+1}& ={\mathit{x}}_t+\mathrm{\Delta }t{\mathit{v}}_{t+1}\end{array}$$

亦或者记作如下的形式

$$\begin{array}{rl}{\mathit{x}}_{t+1}& ={\mathit{x}}_t+\mathrm{\Delta }t{\mathit{v}}_{t+1}\\ {\mathit{v}}_{t+1}& ={\mathit{v}}_t+\mathrm{\Delta }t{\mathbf{M}}^{-1}\mathit{f}({\mathit{x}}_{t+1})\end{array}$$

上述系统方程，可以转化成一个非线性的偏微分方程（PDE），通常有两种思路：（1）化简消去 ${\mathit{x}}_{t+1}$；（2）化简消去 ${\mathit{v}}_{t+1}$

（1）隐式积分求解 - 消去 ${\mathit{x}}_{t+1}$

化简消去 ${\mathit{x}}_{t+1}$ 后，得到系统方程，即

$${\mathit{v}}_{t+1}={\mathit{v}}_t+\mathrm{\Delta }t{\mathbf{M}}^{-1}\mathit{f}({\mathit{x}}_t+\mathrm{\Delta }t{\mathit{v}}_{t+1})$$

这是一个关于 ${\mathit{v}}_{t+1}$ 的偏微分方程 PDE 。

（2）隐式积分求解 - 消去 ${\mathit{v}}_{t+1}$

化简消去 ${\mathit{v}}_{t+1}$ 后，得到系统方程，即

$${\mathit{x}}_{t+1}={\mathit{x}}_t+\mathrm{\Delta }t{\mathit{v}}_t+\mathrm{\Delta }{t}^2{\mathbf{M}}^{-1}\mathit{f}({\mathit{x}}_{t+1})$$

这是一个关于 ${\mathit{x}}_{t+1}$ 的偏微分方程 PDE 。

Tips：在这里，消去 ${\mathit{x}}_{t+1}$ 而保留 ${\mathit{v}}_{t+1}$ 的用意，应该是便于进行碰撞处理。在一些物理引擎中，会先采用 ${\mathit{x}}_t$ 或 ${\mathit{x}}_t+\mathrm{\Delta }t{\mathit{v}}_t$ 作为模型中节点的位置，进行碰撞检测。得到碰撞结果后，进一步更新/限制节点的速度 ${\mathit{v}}_{t+1}$ ，这样的好处好象是稳定性会好一些，不会出现节点的位置穿过碰撞界限等。

---

2 物理引擎中 PDE 的求解

如第 1 节中看到，在物理仿真中，通过空间上的离散化，计算得到了模型中节点上的受力 $\mathit{f}$ ；通过运动方程，描述模型的运动规律，得到了一组 ODE；对模型的运动状态在时间上进行离散，并通过显式/隐式积分器，可以得到一组 PDE。那么，最终，就需要进行 PDE 的求解。（更为复杂的，比如带约束等情况，后面会再整理。这里只描述最基本的模型运动仿真）

（在物理引擎中，会得到各种 ODE、PDE、线性系统、非线性系统，名称上比较混乱。）其中涉及的求解方法也非常多，比如，线性化、牛顿法、Jacobin、CG（Conjugate Gradient）、优化隐式方法等等，非常多样

2.1 线性化（one step of Newton's method）

简单地求解，可以实现一步牛顿法，将 $\mathit{f}$ 在 ${\mathit{x}}_t$ 处一阶泰勒展开，得到如下：

（1）速度上的求解

$${\mathit{v}}_{t+1}={\mathit{v}}_t+\mathrm{\Delta }t{\mathbf{M}}^{-1}[\mathit{f}({\mathit{x}}_t)+\frac{\partial \mathit{f}}{\partial \mathit{x}}({\mathit{x}}_t)\mathrm{\Delta }t{\mathit{v}}_{t+1}]$$

操作之后，系统变成了线性系统。整理之后，为

$$[\mathbf{I}-\mathrm{\Delta }{t}^2{\mathbf{M}}^{-1}\frac{\partial \mathit{f}}{\partial \mathit{x}}({\mathit{x}}_t)]{\mathit{v}}_{t+1}={\mathit{v}}_t+\mathrm{\Delta }t{\mathbf{M}}^{-1}\mathit{f}({\mathit{x}}_t)$$

（2）位置上的求解

$${\mathit{x}}_{t+1}={\mathit{x}}_t+\mathrm{\Delta }t{\mathit{v}}_t+\mathrm{\Delta }{t}^2{\mathbf{M}}^{-1}[\mathit{f}({\mathit{x}}_t)+\frac{\partial \mathit{f}}{\partial \mathit{x}}({\mathit{x}}_t)\mathrm{\Delta }t{\mathit{v}}_t]$$

这里不是很确定，还需要再对照资料确认一下。

如上所述，线性化之后会得到如下的线性系统：

$$\begin{array}{rl}\mathbf{A}& =[\mathbf{I}-\mathrm{\Delta }{t}^2{\mathbf{M}}^{-1}\frac{\partial \mathit{f}}{\partial \mathit{x}}({\mathit{x}}_t)]\\ \mathbf{b}& ={\mathit{v}}_t+\mathrm{\Delta }t{\mathbf{M}}^{-1}\mathit{f}({\mathit{x}}_t)\\ \mathbf{A}{\mathit{v}}_{t+1}& =\mathbf{b}\end{array}$$

对该线性系统的求解，可以采用 Jacobin / Gauss-Seidel iterations，或者 Conjugate gradients 等方法进行求解。

Jacobin iteration 求解线性系统的方法

单步的 Jacobin 迭代过程如下所示：

A = []
x = []
new_x = []
b = []

@ti.kernel
def iterate():
  for i in range(n):
    r = b[i]
    for j in range(n):
      if i != j:
        r -= A[i,j] * x[j]

    new_x[i] = r / A[i,i]

  for i in range(n):
    x[i] = new_x[i]

Jacobin 迭代的思想就是，每次只让一行（一个点）满足该方程，依次计算完成，就能让所有的点都（曾经）满足方程。多迭代几次，就能接近线性方程组的解了。（大概是这么个意思吧）

2.2 线性化（one step of Newton's method） - 进阶

在上节所述的线性系统中，加入一个系数 $\beta$ ，那么，就会得到如下线性系统：

$$[\mathbf{I}-\beta \mathrm{\Delta }{t}^2{\mathbf{M}}^{-1}\frac{\partial \mathit{f}}{\partial \mathit{x}}({\mathit{x}}_t)]{\mathit{v}}_{t+1}={\mathit{v}}_t+\mathrm{\Delta }t{\mathbf{M}}^{-1}\mathit{f}({\mathit{x}}_t)$$

其中，当 $\beta =0$ 时，为 forward / semi-implicit Euler (explicit integrator)
当 $\beta =1/2$ 时，为 middle-point (implicit integrator)
当 $\beta =1$ 时，为 backward Euler (implicit integrator)