L22 Simulation

1 单粒子模拟 Single Particle Simulation

首先研究单个粒子的运动。之后，推广到大量粒子
初始假设：粒子的运动由一个速度矢量场决定，该矢量场是位置和时间的函数：\(v(x,t)\)

计算粒子随时间变化的位置需要求解一阶常微分方程(first-order ordinary differential equation)：

\[\frac{dx}{dt} = \dot{x} = v(x,t)\]

“一阶”(First-order) 指的是取一阶导数。
“常”(Ordinary) 指的是没有偏导数，即 x 只是 t 的函数。

2 欧拉方法(Euler's Methord)

欧拉方法（又称前向欧拉法、显式欧拉法）Euler’s Method (a.k.a. Forward Euler, Explicit Euler)

简单的迭代方法(Simple iterative method)
常用(Commonly used)
精度非常低(inaccurate)
通常会变得不稳定(unstable)

\[ \begin{align} \boldsymbol{x}^{t+\Delta t} = \boldsymbol{x}^t + \Delta t \dot{\boldsymbol{x}}^t \\ \dot{\boldsymbol{x}}^{t+\Delta t} = \dot{\boldsymbol{x}}^t + \Delta t \ddot{\boldsymbol{x}}^t \end{align} \]

数值积分(numerical integration)会导致误差累积(accumulate)，欧拉积分的表现尤其糟糕
示例：

\[\boldsymbol{x}^{t+\Delta t} = \boldsymbol{x}^t + \Delta t \boldsymbol{v}(\boldsymbol{x}, t)\]

两个核心问题：

时间步长 \(Δt\) 越大，误差越大
不稳定性是一个常见且严重的问题，会导致模拟发散
在下图左边的速度场中，使用欧拉方法一定会导致粒子不断远离中心。这种现象被称为正反馈。

误差(Errors)

每一时间步的误差会累积。模拟进行得越久，精度越低
在图形学应用中，精度可能并非关键
不稳定性(Instability)
误差会不断叠加(compound)，即使原系统(underlying system)本身不会发散，也会导致模拟发散(diverge)
缺乏稳定性是模拟中的一个根本性问题，无法被忽视

3 解决不稳定性的方法

一些解决不稳定性的方法

中点法 / 改进欧拉法：取起点和终点的速度平均值
自适应步长：递归比较单步和两个半步的结果，直到误差可接受
隐式方法：使用下一时刻的速度（难度较高）
基于位置的积分 / 韦莱积分：在时间步结束后约束粒子的位置和速度

3.1 中点法(Midpoint Method)和改进欧拉方法(Modified Euler)

计算欧拉步（步骤 a）
计算欧拉步中点处的导数（步骤 b）
使用中点导数更新位置（步骤 c）

\[ \begin{align} x_{\text{中点}} &= x(t) + \Delta t/2 \cdot v(x(t), t) \\ x(t + \Delta t) &= x(t) + \Delta t \cdot v(x_{\text{中点}}, t) \end{align} \]

用两次欧拉方法，使用第一次求得的中点位置的速度作为第二次欧拉方法的速度。
与之类似的还有一种方法是改进的欧拉方法(Modified Euler)。
改进欧拉法取时间步起点和终点的速度平均值,效果更好

\[ \begin{align} \boldsymbol{x}^{t+\Delta t} &= \boldsymbol{x}^t + \frac{\Delta t}{2} \left( \dot{\boldsymbol{x}}^t + \dot{\boldsymbol{x}}^{t+\Delta t} \right) \\ \dot{\boldsymbol{x}}^{t+\Delta t} &= \dot{\boldsymbol{x}}^t + \Delta t \ddot{\boldsymbol{x}}^t \\ \boldsymbol{x}^{t+\Delta t} &= \boldsymbol{x}^t + \Delta t \dot{\boldsymbol{x}}^t + \frac{(\Delta t)^2}{2} \ddot{\boldsymbol{x}}^t \end{align} \]

与普通的欧拉方法相比，这种方法能够使用二次的信息，因此更加准确。

3.2 自适应步长(Adaptive Step Size)

自适应步长

基于误差估计选择步长的技术
非常实用的技术, 但可能需要非常小的步长！
重复执行直到误差低于阈值：
计算 \(x_T\)：一个步长为 \(T\) 的欧拉步
计算 \(x_{T/2}\)：两个步长为 \(T/2\) 的欧拉步
计算误差 \(\| x_T - x_{T/2} \|\)
如果（误差 > 阈值），则减小步长并重试

3.3 隐式/后向欧拉方法(Implicit Euler Methord)

隐式欧拉方法(Implicit Euler Methord)，非正式的称为后向欧拉方法。核心思想是为当前步使用未来时刻的导数。

\[ \begin{align} \boldsymbol{x}^{t+\Delta t} = \boldsymbol{x}^t + \Delta t \dot{\boldsymbol{x}}^{t+\Delta t} \\ \dot{\boldsymbol{x}}^{t+\Delta t} = \dot{\boldsymbol{x}}^t + \Delta t \ddot{\boldsymbol{x}}^{t+\Delta t} \end{align} \]

求解关于\(\boldsymbol{x}^{t+\Delta t}\)和\(\dot{\boldsymbol{x}}^{t+\Delta t}\)的非线性问题，可以使用求根算法，例如牛顿法
这种方法能够提供更好的稳定性

如何定义 / 量化 “稳定性”？

我们使用局部截断误差（每一步）(local truncation error (every step))/ 总累积误差（整体）(total accumulated error (overall))。绝对值无关紧要，重要的是相对于步长的阶数它的大小。
隐式欧拉法为 1 阶，这意味着
- 局部截断误差：\(O (h²)\)
- 全局截断误差：\(O (h)\) （\(h\) 为步长，即 \(Δt\)）

对 \(O (h)\) 的理解
- 如果将 \(h\) 减半，预计误差也会减半

3.4 龙格库塔方法(Runge-Kutta Families)

这是求解常微分方程的一类高级方法(A family of advanced methods for solving ODEs)

尤其擅长处理非线性问题
其四阶版本应用最广，又称 RK4

初始条件(Initial condition)：

\[\frac{dy}{dt} = f(t,y),\quad y(t_0) = y_0\]

RK4 解：

\[ \begin{align} y_{n+1} &= y_n + \frac{1}{6}h\left(k_1 + 2k_2 + 2k_3 + k_4\right) \\ t_{n+1} &= t_n + h \end{align} \]

其中

\[ \begin{align} k_1 &= f(t_n, y_n),\\ k_2 &= f\left(t_n + \frac{h}{2}, y_n + h\frac{k_1}{2}\right), \\ k_3 &= f\left(t_n + \frac{h}{2}, y_n + h\frac{k_2}{2}\right),\\ k_4 &= f\left(t_n + h, y_n + hk_3\right) \end{align} \]

更多相关知识在《数值分析(numarical analysis)》中会详细讲解

3.5 基于位置的积分/韦莱积分(Position-Based/Verlet Intergration)

基于位置的积分思路：

在改进欧拉前向步之后，约束粒子的位置以防止发散和不稳定行为
使用约束后的位置计算速度
这两种思路都会耗散能量，实现稳定
比如，要仿真弹簧系数无限大的弹簧，不用再考虑弹簧力计算加速度 \(a\)，而是用解约束的方法来更新质点位置：只要简单的移动每个质点的位置使得弹簧的长度保持原长。修正向量应该和两个质点之间的位移成比例，方向为一个质点指向另一质点。每个质点应该移动位移的一半。
只要对每个弹簧执行这样的操作，我们就可以得到稳定的仿真。
优缺点
快速且简单
不基于物理原理，会耗散能量（误差）

Verlet积分也属于一种基于位置的积分。这种方法的优点是只处理仿真中顶点的位置并且保证四阶精度。和欧拉法不同，Verlet积分按如下的方式来更新下一步位置：

\[x(t+1)=x(t)+[x(t)-x(t-1)]+a(t)\cdot \Delta t^2\]

公式中的 \(a(t)\) 也可以仅由位置信息得到，因此说Verlet积分是一种基于位置的积分(Position-Based Intergration)。

Verlet积分公式的推导(泰勒展开法)
对位置 \(x (t)\) 进行向前和向后泰勒展开（保留到 \(Δt^4\) 项）：

\[ \begin{align} x(t+\Delta t)=x(t)+v(t)\Delta t+\frac{1}{2}a(t)\Delta t^2 +\frac{1}{6}b(t)\Delta t^3 +O(\Delta t^4) \\ x(t-\Delta t)=x(t)-v(t)\Delta t+\frac{1}{2}a(t)\Delta t^2 -\frac{1}{6}b(t)\Delta t^3 +O(\Delta t^4) \\ \end{align} \]

将两式相加，消去速度项 \(v (t)\) 和三阶导数项 \(b (t)\)，得到：

\[x(t+\Delta t)=2x(t)-x(t-\Delta t)+a(t)\Delta t^2 + O(\Delta t^4)\]

4 刚体模拟(Rigid Body Simulation)

简单情况与粒子模拟类似，只需多考虑一些属性

\[\frac{d}{dt} \begin{pmatrix} X \\ \theta \\ \dot{X} \\ \omega \end{pmatrix} = \begin{pmatrix} \dot{X} \\ \omega \\ \mathbf{F}/M \\ \Gamma/I \end{pmatrix}\]

\(X\)：位置
\(\theta\)：旋转角度
\(\omega\)：角速度
\(F\)：力
\(\Gamma\)：扭矩
\(I\)：转动惯量

5 流体模拟(Fluid Simulation)

5.1 A Simple Position-Based Method

核心思路

假设水由小型刚体球体组成(small rigid-body spheres )
假设水不可压缩(cannot be compressed)（即密度恒定）
因此，只要某处密度发生变化，就应通过改变粒子位置来 “修正”
任意位置的密度都是相对于每个粒子位置的一个函数。你需要知道这个函数在某个位置处的梯度
更新方法？直接用梯度下降(gradient descent)！

5.2 Eulerian vs. Lagrangian

（“质点法”）拉格朗日法
摄影师全程跟踪同一只鸟。
（“网格法”）欧拉法
摄影师保持静止，只能拍摄所有经过某一画面的鸟（比如在时间‘t’时）。

5.3 Material Point Method (MPM)

这是一种混合方法，它结合了欧拉视角与拉格朗日视角

拉格朗日：考虑携带物质属性的粒子(particles carrying material properties)
欧拉：使用网格(grid)进行数值积分
交互方式(Interaction)：粒子将属性传递给网格，网格执行更新，然后插值(interpolate)回粒子