物理模拟：拉格朗日视角（一）
本文是对 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）
根据分析粒子受到的合力，通过时间积分，可以计算并更新粒子的速度、位置。
这里提供了三种积分的方式：
- 前向欧拉方法（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$$
- 半隐式欧拉方法（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}$$
- 后向欧拉方法（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