Beyond awesome | 越而胜己

本文是对 GAMES201线上课程:高级物理引擎实战指南2020 第二课的部分笔记,主要对讲义中复杂的公式做一些注释,以供后续参考。

拉格朗日视角 vs 欧拉视角(Lagrangian vs Eulerian)

简单的说,拉格朗日视角将物质视为粒子的集合,并计算每个例子的运动状态和受力情况来模拟物质的运动,老师称为“随波逐流”;欧拉视角则是将空间分割成固定位置的网格,计算每个方格中物质的流入流出,老师称为“岿kuī然不动”。

弹簧质点系统(Mass-Spring Systems)

一个弹簧质点系统有如下参数:$k$ 弹簧的劲度系数、$l_{ij}$ 连接粒子 $i$、$j$ 的弹簧的静止长度、$m_i$ 粒子质量。则根据胡克定律(Hooke's Law),粒子 $i$ 受到来自 $j$ 的弹力为

$$\mathbf{f}_{ij} = -k \left( \left| \mathbf{x}_i - \mathbf{x}_j \right|_2 - l_{ij} \right) \left( \widehat{\mathbf{x}_i - \mathbf{x}_j} \right)$$

而粒子 $i$ 受到的合力为

$$\mathbf{f}_{i} = \sum_{j}^{j \neq i} \mathbf{f}_{ij}$$

根据牛顿第二定律,

$$\frac{\partial \mathbf{v}_i}{\partial t} = \frac{1}{m_i} \mathbf{f}_i,\quad \frac{\partial \mathbf{x}_i}{\partial t} = \mathbf{v}_i$$

时间积分(Time integration)

根据分析粒子受到的合力,通过时间积分,可以计算并更新粒子的速度、位置。

这里提供了三种积分的方式:

  1. 前向欧拉方法(Forward Euler,显式)

$$\mathbf{v}_{t+1} = \mathbf{v}_t + \Delta t \frac{\mathbf{f}_t}{m}, \quad \mathbf{x}_{t+1} = \mathbf{x}_t + \Delta t \mathbf{v}_t$$

  1. 半隐式欧拉方法(Semi-implicit Euler,Symplectic Euler,仍是显式)

$$\mathbf{v}_{t+1} = \mathbf{v}_t + \Delta t \frac{\mathbf{f}_t}{m}, \quad \mathbf{x}_{t+1} = \mathbf{x}_t + \Delta t \mathbf{v}_{t+1}$$

  1. 后向欧拉方法(Backward Euler,隐式)

$$\mathbf{v}_{t+1} = \mathbf{v}_t + \Delta t \frac{\partial \mathbf{v}_{t+1}}{\partial t}, \quad \mathbf{x}_{t+1} = \mathbf{x}_t + \Delta t \mathbf{v}_{t+1}$$

显式方法简单、容易实现,但是容易发散(爆炸),一般来说要求 $\Delta t \leq c \sqrt{\frac{m}{k}} (c \sim 1)$;隐式方法用未来的速度和位置信息来更新未来的信息,数值更加稳定,但是实现复杂,需要求解线性方程。

改写一下后项欧拉的方程:

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

其中 $t+1$ 时刻的加速度 $\frac{\partial \mathbf{v}_{t+1}}{\partial t}$ 使用牛顿第二定律替代。假设粒子受到的合力是粒子位移的函数,则 $\frac{\partial \mathbf{v}_{t+1}}{\partial t} = \frac{1}{m} \mathbf{f}(\mathbf{x}_{t+1})$。

这里使用质量矩阵 $\mathbf{M}$ 来表示许多粒子的质量。$\mathbf{M}$ 是对角矩阵,对角线上每个元素都是一个粒子的质量,因此 $\mathbf{M}^{-1}$ 即为所有粒子质量的倒数的对角矩阵。

消去 $\mathbf{x}_{t+1}$:

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

一阶泰勒展开 $\mathbf{f}(\mathbf{x}_t + \Delta t \mathbf{v}_{t+1}) = \mathbf{f}(x_t) + \frac{\partial \mathbf{f}}{\partial t}(\mathbf{x}_t) \Delta t \mathbf{v}_{t+1}$:

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

化简:

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

这样就得到了一个形如 $\mathbf{A} \mathbf{x} = \mathbf{b}$ 的线性方程。
由于直接求解 $\mathbf{A}^{-1}$ 过于困难,可以使用雅可比迭代(Jacobi iterations)或者共轭梯度(conjugate gradients)。老师给出了雅可比迭代在 Taichi 中的简单实现:

A = ti.var(dt=ti.f32, shape=(n, n))
x = ti.var(dt=ti.f32, shape=n)
new_x = ti.var(dt=ti.f32, shape=n)
b = ti.var(dt=ti.f32, shape=n)

@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]

可以将隐式和显式的积分统一到一个公式中:

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

参数 $\beta$ 为 $0$ 时即为显式(前向或半隐式欧拉),为 $1$ 时即为后向欧拉,为 $\frac{1}{2}$ 即为中点(也是隐式)。

如果粒子和弹簧数量巨大,可以采用稀疏矩阵、共轭梯度、position-based dynamics 等技巧。或者采用另一种不同的快速方法: fast mass-spring system solver。

参考资料

GAMES201 Lec2 Slides: https://github.com/taichi-dev/games201/releases/download/lec2/lec2-math.pdf