物理模拟:拉格朗日视角(一)
本文是对 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