计算机图形学-15：数值积分与物理模拟

初值问题

已知随位置和时间变化的速度场 $\mathbf v(\mathbf x,t)$，粒子轨迹满足一阶 常微分方程 Ordinary Differential Equation：

$$\frac{\mathrm d\mathbf x}{\mathrm dt} = \mathbf v(\mathbf x,t), \qquad \mathbf x(t_0)=\mathbf x_0.$$

精确解满足积分形式

$$\mathbf x(t+h) = \mathbf x(t) + \int_t^{t+h} \mathbf v(\mathbf x(\tau),\tau)\,\mathrm d\tau.$$

给定初始状态后，这构成 初值问题 Initial Value Problem。数值积分的任务，是只在有限时刻查询速度或加速度，近似这段连续积分。多个粒子、质点—弹簧系统和刚体都可写成统一状态方程

$$\dot{\mathbf y} = \mathbf f(\mathbf y,t),$$

其中 $\mathbf y$ 可以同时包含位置、速度、姿态和其他状态。

图中每个位置的箭头给出当地速度，粒子沿着经过自身位置的速度方向连续移动，轨迹不是把相邻箭头直接连起来。数值积分正是在有限时刻查询这个场，再近似真实轨迹上的积分。

显式欧拉

显式 Euler Explicit Euler 用区间起点的导数代表整个时间步：

$$\mathbf y_{n+1} = \mathbf y_n +h\mathbf f(\mathbf y_n,t_n).$$

它也可以从带积分余项的 Taylor 展开得到：

$$\mathbf y(t_n+h) = \mathbf y(t_n) +h\dot{\mathbf y}(t_n) +\int_0^h (h-s)\ddot{\mathbf y}(t_n+s)\,\mathrm ds = \mathbf y(t_n) +h\dot{\mathbf y}(t_n) +O(h^2).$$

舍去积分余项，就得到 Euler 更新。积分形式对向量状态仍然严格成立，不需要假设所有分量共享同一个中值点。若解足够光滑，每一步的 局部截断误差 Local Truncation Error 为 $O(h^2)$；固定总模拟时间时大约需要 $1/h$ 步，误差累积后的 全局误差 Global Error 为 $O(h)$。

对粒子状态 $(\mathbf x,\mathbf v)$，常见写法为

$$\mathbf x_{n+1} = \mathbf x_n+h\mathbf v_n,$$

$$\mathbf v_{n+1} = \mathbf v_n +h\mathbf a(\mathbf x_n,\mathbf v_n,t_n).$$

公式简单、每步便宜，但刚性弹簧、高频振动和大时间步很容易让它失稳。

半隐式 Euler Semi-Implicit Euler 仍用当前状态计算加速度，却先更新速度，再用新速度更新位置：

$$\mathbf v_{n+1} = \mathbf v_n+h\mathbf a(\mathbf x_n,\mathbf v_n,t_n),$$

$$\mathbf x_{n+1} = \mathbf x_n+h\mathbf v_{n+1}.$$

无阻尼简谐振子满足 $\dot x=v$、$\dot v=-\omega^2x$。半隐式更新可写成

$$\begin{pmatrix} x_{n+1}\\v_{n+1} \end{pmatrix} = \underbrace{ \begin{pmatrix} 1-h^2\omega^2&h\\ -h\omega^2&1 \end{pmatrix} }_{A_h} \begin{pmatrix} x_n\\v_n \end{pmatrix}.$$

$A_h$ 的行列式为 1，迹为 $2-h^2\omega^2$，特征值满足

$$\mu^2-(2-h^2\omega^2)\mu+1=0.$$

两根乘积为 1。要让它们成为单位圆上的共轭复数，必须有

$$|2-h^2\omega^2|<2,$$

也就是 $0<h\omega<2$。$h\omega=2$ 时出现重根 $-1$，一般初值仍可能产生线性增长；超过 2 后会出现模长大于 1 的实根。半隐式 Euler 在有效区间内通常比显式 Euler 更好地控制能量，但仍是条件稳定，不能用大步长解决任意刚性问题。

稳定性

精度描述数值解离真实解有多远，稳定性 Stability 描述误差在反复更新中是否被放大。即使每一步误差不大，只要更新规则持续放大扰动，轨迹仍会发散。

用标量测试方程

$$\dot y=\lambda y$$

分析积分器。精确一步的放大因子为 $e^{h\lambda}$。显式 Euler 给出

$$y_{n+1} =(1+h\lambda)y_n,$$

数值放大因子为

$$g_E=1+h\lambda.$$

当真实系统应衰减，即 $\operatorname{Re}\lambda<0$ 时，只有 $|1+h\lambda|\le1$ 的步长才稳定。$\lambda$ 的负实部绝对值越大，允许的 $h$ 越小，这就是 刚性系统 Stiff System 对显式积分器的限制。

简谐运动的一阶状态矩阵具有特征值 $\lambda=\pm i\omega$。此时

$$|g_E| = |1\pm ih\omega| = \sqrt{1+h^2\omega^2}>1.$$

每一步都会放大振幅，所以本应封闭的圆形相轨迹会向外螺旋。减小步长只能减慢能量增长，无法让显式 Euler 对纯虚特征值真正稳定。

截断误差决定“算得准不准”，稳定性决定“误差会不会失控”；高阶方法也不自动拥有更大的稳定区域。

阻尼应作为动力学中的力项加入，例如 $\mathbf f_d=-c\mathbf v$；它会耗散能量、改变系统特征值，但并不等于把弹簧刚度 $k$ 调小。刚度决定最快的恢复模态，积分器仍需按该时间尺度选择步长。

图中数值轨迹从本应封闭的轨道逐步向外螺旋，正是放大因子超过单位模的表现。减小步长只能让发散变慢；对纯虚特征值，显式 Euler 仍满足 $|1+ih\omega|>1$，所以条件稳定的本质不会改变。

中点法

显式中点法 Explicit Midpoint Method 先用 Euler 估计半步状态，再用中点斜率走完整一步：

$$\mathbf k_1 = \mathbf f(\mathbf y_n,t_n),$$

$$\mathbf y_{\mathrm{mid}} = \mathbf y_n+\frac h2\mathbf k_1,$$

$$\mathbf k_2 = \mathbf f \left( \mathbf y_{\mathrm{mid}}, t_n+\frac h2 \right),$$

$$\mathbf y_{n+1} = \mathbf y_n+h\mathbf k_2.$$

起点斜率只把状态推到时间区间中点，真正的整步更新使用中点处的斜率。相比只看起点，采样位置更接近区间平均，因此在相同步长下能明显降低轨迹弯曲带来的误差。

起点斜率只负责找到区间中点，在对称位置计算的 $\mathbf k_2$ 更接近整段平均斜率。Taylor 展开可验证其局部误差为 $O(h^3)$、全局误差为 $O(h^2)$。

改进欧拉

改进 Euler Modified Euler，也称 Heun Method，先预测终点，再平均起点和终点的斜率：

$$\mathbf k_1 = \mathbf f(\mathbf y_n,t_n), \qquad \tilde{\mathbf y} = \mathbf y_n+h\mathbf k_1,$$

$$\mathbf k_2 = \mathbf f(\tilde{\mathbf y},t_n+h),$$

$$\mathbf y_{n+1} = \mathbf y_n+\frac h2(\mathbf k_1+\mathbf k_2).$$

它近似梯形积分，同样是全局二阶方法。对恒定加速度，位置更新化为熟悉的

$$\mathbf x_{n+1} = \mathbf x_n +h\mathbf v_n +\frac{h^2}{2}\mathbf a_n.$$

中点法和 Heun 法提高了精度，但仍属于显式方法，稳定区域有限。方法阶数高并不等于可以对任意刚性系统使用大步长。

自适应步长

固定步长会在平缓区域浪费计算，在碰撞或高速变化附近又可能过粗。自适应步长 Adaptive Step Size 可用 步长加倍 Step Doubling：一次长度 $h$ 的更新得到 $\mathbf y_h$，再用两次长度 $h/2$ 的更新得到 $\mathbf y_{h/2}$。对全局 $p$ 阶方法，两者的差是局部误差指示量，较细结果的 Richardson 误差估计为

$$\mathbf e = \frac{\mathbf y_h-\mathbf y_{h/2}}{2^p-1}.$$

位置、速度和其他状态的量纲不同，不能直接用一个未缩放的欧氏范数。常用逐分量容差构造

$$\varepsilon = \sqrt{ \frac1m \sum_{j=1}^{m} \left( \frac{e_j} {\mathrm{atol}_j+\mathrm{rtol}_j|y_j|} \right)^2 }.$$

若 $\varepsilon>1$ 则拒绝该步、减小 $h$；误差远低于 1 则可在下一步增大 $h$。

对全局 $p$ 阶方法，局部误差主项通常为 $Ch^{p+1}$，可按

$$h_{\mathrm{new}} \approx \eta h \left( \frac1{\varepsilon} \right)^{1/(p+1)}$$

调整，其中 $\eta<1$ 是安全系数。这里的 $\varepsilon$ 已经除以绝对和相对容差；若使用未归一化误差，分子才应显式写成目标容差。实际还要限制单次放大或缩小倍数。自适应步长控制精度，却不能突破显式方法的稳定区域；刚性问题中即使解变化平缓，也可能被迫使用很小步长。

轨迹平缓处可以跨更大的时间步，碰撞或速度快速变化处则自动缩短步长。Step Doubling 控制的是精度，不会把不稳定的显式积分器变成无条件稳定。

隐式欧拉

隐式 Euler Implicit Euler 在时间步终点计算导数：

$$\mathbf y_{n+1} = \mathbf y_n +h\mathbf f(\mathbf y_{n+1},t_{n+1}).$$

$\mathbf y_{n+1}$ 同时出现在等式两侧，因此每一步都要解方程。定义残差

$$\mathbf R(\mathbf y) = \mathbf y-\mathbf y_n -h\mathbf f(\mathbf y,t_{n+1}),$$

目标是求 $\mathbf R(\mathbf y)=\mathbf0$。Newton 迭代在线性化后求

$$\left( I-h\frac{\partial\mathbf f}{\partial\mathbf y} \right) \Delta\mathbf y = -\mathbf R(\mathbf y),$$

这就是 Newton 迭代 Newton’s Method。再更新 $\mathbf y\leftarrow\mathbf y+\Delta\mathbf y$，直到残差足够小。大系统中的 Jacobian 通常稀疏，线性方程求解成本会主导每一步。

对测试方程 $\dot y=\lambda y$，

$$y_{n+1} = y_n+h\lambda y_{n+1},$$

所以

$$g_I = \frac1{1-h\lambda}.$$

当 $\operatorname{Re}\lambda\le0$ 时，隐式 Euler 对任意 $h>0$ 都不会放大该衰减模态。对 $\lambda=i\omega$，

$$|g_I| = \frac1{\sqrt{1+h^2\omega^2}}\le1,$$

$\omega\ne0$ 时严格小于 1，振荡会被人工衰减；$\omega=0$ 时等于 1。它的稳定性优于显式 Euler，精度仍只有全局一阶，并带有明显数值耗散。稳定意味着结果不爆炸，并不等于结果准确或守恒。

RK4

经典 四阶 Runge–Kutta RK4 在一步内查询四次斜率：

$$\mathbf k_1 = \mathbf f(\mathbf y_n,t_n),$$

$$\mathbf k_2 = \mathbf f \left( \mathbf y_n+\frac h2\mathbf k_1, t_n+\frac h2 \right),$$

$$\mathbf k_3 = \mathbf f \left( \mathbf y_n+\frac h2\mathbf k_2, t_n+\frac h2 \right),$$

$$\mathbf k_4 = \mathbf f \left( \mathbf y_n+h\mathbf k_3, t_n+h \right).$$

再按

$$\mathbf y_{n+1} = \mathbf y_n +\frac h6 (\mathbf k_1+2\mathbf k_2+2\mathbf k_3+\mathbf k_4)$$

组合。权重使 Taylor 展开直到四阶都与精确解匹配，因此局部误差为 $O(h^5)$、全局误差为 $O(h^4)$。

RK4 对平滑、非刚性问题很有效，但它仍是显式方法，没有无条件稳定性，也不保持 Hamiltonian 系统的能量或辛结构。长时间轨道模拟可能更适合 symplectic integrator。

Verlet

Verlet Integration 从位置的前向和后向 Taylor 展开开始：

$$\mathbf x(t+h) = \mathbf x(t)+h\dot{\mathbf x}(t) +\frac{h^2}{2}\ddot{\mathbf x}(t) +O(h^3),$$

$$\mathbf x(t-h) = \mathbf x(t)-h\dot{\mathbf x}(t) +\frac{h^2}{2}\ddot{\mathbf x}(t) +O(h^3).$$

两式相加消去速度，得到经典 Verlet：

$$\mathbf x_{n+1} = 2\mathbf x_n-\mathbf x_{n-1} +h^2\mathbf a_n +O(h^4).$$

速度可用中心差分近似：

$$\mathbf v_n \approx \frac{\mathbf x_{n+1}-\mathbf x_{n-1}}{2h}.$$

Velocity Verlet 显式保存速度：

$$\mathbf x_{n+1} = \mathbf x_n+h\mathbf v_n +\frac{h^2}{2}\mathbf a_n,$$

$$\mathbf v_{n+1} = \mathbf v_n +\frac h2 (\mathbf a_n+\mathbf a_{n+1}).$$

这一写法对只依赖位置和时间的加速度 $\mathbf a(\mathbf x,t)$ 最直接，因为先得到 $\mathbf x_{n+1}$ 后就能计算 $\mathbf a_{n+1}$。若加速度还依赖未知速度，例如黏性阻尼 $\mathbf a(\mathbf x,\mathbf v,t)$，第二式不再是直接显式更新，需要迭代、隐式处理或专用变体。

在适用条件下，Velocity Verlet 具有时间对称性，并在许多保守系统中比普通 Euler 更好地控制长期能量漂移。

位置约束

Position-Based Dynamics 直接在位置上处理碰撞、距离和不可压缩等约束：

1. 根据外力预测新位置 $\mathbf x_i^*$；
2. 反复投影位置，使约束 $C_j(\mathbf x)=0$ 或 $C_j(\mathbf x)\ge0$ 尽量满足；
3. 用修正后的位置更新速度：

$$\mathbf v_i^{n+1} = \frac{\mathbf x_i^{n+1}-\mathbf x_i^n}{h}.$$

同一质点常同时参与多个约束，一次遍历中修正一个约束可能破坏前面刚满足的另一个约束。因此需要多轮迭代：Gauss–Seidel 风格立即使用最新位置，收敛较快；Jacobi 风格先收集所有修正再批量应用，更易并行但通常需要更多轮次。

图中先由外力预测粒子位置，再在约束投影阶段逐轮消除穿透和密度误差，最后把修正前后的位移差回算为速度。投影轮数越少，约束残差越明显；轮数增加会提高等效刚度和计算成本。

对单个等式约束，先在预测位置处线性化：

$$C(\mathbf x+\Delta\mathbf x) \approx C(\mathbf x) +\sum_i \nabla_{\mathbf x_i}C\cdot\Delta\mathbf x_i =0.$$

为让质量较小的粒子承担更多修正，取质量加权的梯度方向

$$\Delta\mathbf x_i = -w_i\lambda\nabla_{\mathbf x_i}C,$$

其中 $w_i=1/m_i$。代回线性化方程可得

$$\lambda = \frac{C(\mathbf x)} {\displaystyle\sum_i w_i\lVert\nabla_{\mathbf x_i}C\rVert^2},$$

再用上式计算每个 $\Delta\mathbf x_i$。分母接近零时，约束梯度无法给出稳定修正，应跳过该约束或加入小的正则项。对不等式 $C(\mathbf x)\ge0$，只有 $C<0$、约束被违反时才执行边界投影；已经满足约束时保持 inactive。更完整的接触求解还要让约束间隙与接触乘子满足互补条件。

位置投影能快速消除明显穿透和密度误差，迭代次数、步长和约束顺序会影响等效刚度。它通常稳定、易控，却会引入数值耗散，不等同于精确求解原始力学方程。

刚体

刚体模拟 Rigid-Body Simulation 除了质心位置 $\mathbf x$ 和线速度 $\mathbf v$，还需要姿态 $R$ 或单位四元数、角速度 $\boldsymbol\omega$、质量 $m$ 与转动惯量 $I$。以下角动量 $\mathbf L$ 和力矩 $\boldsymbol\tau$ 都相对质心定义。线动量和角动量满足

$$\dot{\mathbf p}=\mathbf F, \qquad \dot{\mathbf L}=\boldsymbol\tau,$$

物体坐标中的惯量张量 $I_{\mathrm{body}}$ 随刚体保持不变。若 $R$ 把 body space 变换到 world space，则

$$I_{\mathrm{world}} = R I_{\mathrm{body}}R^T.$$

因此

$$\mathbf v=\frac{\mathbf p}{m}, \qquad \boldsymbol\omega=I_{\mathrm{world}}^{-1}\mathbf L.$$

若 $R$ 把 body space 向量变换到 world space，且 $\boldsymbol\omega$ 用 world space 表示，旋转矩阵的连续更新为

$$\dot R=[\boldsymbol\omega]_\times R,$$

$[\boldsymbol\omega]_\times$ 是叉积的反对称矩阵。若采用 world-to-body 矩阵或 body-space 角速度，乘法侧与符号会相应改变。直接对 $R$ 的元素积分会逐渐破坏正交性，需要重新正交化或使用归一化四元数。

碰撞会在很短时间内产生冲量，瞬间改变线动量和角动量。接触求解还要同时满足不穿透、恢复系数和摩擦约束，因此比单粒子积分多出姿态、力矩与接触系统。

位置流体

Position-Based Fluids 把液体表示为粒子，并用平滑核估计粒子 $i$ 附近的密度：

$$\rho_i = \sum_j m_jW(\lVert\mathbf x_i-\mathbf x_j\rVert,h_s),$$

$W$ 是支持半径为 $h_s$ 的核函数。不可压缩条件要求密度接近静止密度 $\rho_0$，可写成约束

$$C_i(\mathbf x) = \frac{\rho_i}{\rho_0}-1 =0.$$

记 $\nabla W_{ij}$ 为 $W(\mathbf x_i-\mathbf x_j,h_s)$ 对 $\mathbf x_i$ 的梯度。密度约束对粒子位置的梯度为

$$\nabla_{\mathbf x_k}C_i = \begin{cases} \dfrac1{\rho_0} \displaystyle\sum_jm_j\nabla W_{ij},&k=i,\\[6pt] -\dfrac{m_k}{\rho_0}\nabla W_{ik},&k\ne i\text{ 且 }k\in\mathcal N_i,\\[6pt] \mathbf0,&\text{其他}. \end{cases}$$

因此每个密度约束的乘子可取

$$\lambda_i = -\frac{C_i} {\displaystyle \sum_k w_k\lVert\nabla_{\mathbf x_k}C_i\rVert^2 +\varepsilon},$$

其中 $w_k=1/m_k$，$\varepsilon$ 防止邻域退化时分母过小。对常见的等质量粒子，把相邻约束对粒子 $i$ 的修正合并，可写成

$$\Delta\mathbf x_i = \frac1{\rho_0} \sum_{j\in\mathcal N_i} (\lambda_i+\lambda_j+s_{\mathrm{corr}}) \nabla W_{ij}.$$

$s_{\mathrm{corr}}$ 是抑制低密度区域粒子成团的人工压力项，常取

$$s_{\mathrm{corr}} = -k \left[ \frac{W(\lVert\mathbf x_i-\mathbf x_j\rVert,h_s)} {W(\Delta q,h_s)} \right]^n,$$

其中 $k>0$、$n>0$，参考距离 $\Delta q$ 位于核支持域内。算法预测位置后，反复计算邻域、密度、乘子和位置修正，直到密度误差足够小。

自由表面粒子天然缺少邻居，直接强迫 $\rho_i=\rho_0$ 可能把表面错误拉拢；常见处理是只强烈纠正过密状态、加入人工压力，并使用边界粒子、ghost sample 或 SDF 约束补足固体边界的核支持。PBF 能高效产生近似不可压缩流体，但不是对 Navier–Stokes 方程的完整离散；黏性、涡量和表面张力需要额外模型。

两种视角

Lagrangian 方法 跟随物质运动，粒子携带质量、速度、温度和材料状态。它容易追踪自由表面与历史信息，但邻域不断变化，大形变后采样可能不均匀。

Eulerian 方法 把速度、压力和密度存放在固定空间网格上，物质穿过网格。求空间导数和压力投影方便，处理自由表面时则要额外追踪界面，也会产生数值耗散。

两种名称描述观察方式：跟着物质走，或站在固定位置看物质流过。流体方法经常在二者之间转移信息，以同时利用粒子的材料追踪和网格的规则计算结构。

图中 Eulerian 方法让固定网格记录不同位置的速度和压力，粒子穿过网格；Lagrangian 方法让粒子携带质量和材料历史，网格只在需要时参与计算。前者便于求空间导数，后者便于追踪自由表面和大形变，实际算法会在两种表示之间交换信息。

物质点法

Material Point Method，简称 MPM，是粒子与背景网格结合的混合方法。粒子携带质量、速度、体积、形变梯度和材料参数；网格只在当前时间步承担邻域交互和数值更新。

粒子到网格阶段用形函数权重 $w_{ip}=N_i(\mathbf x_p)$ 转移质量和动量。后文的

$$\nabla w_{ip} = \left.\nabla_{\mathbf x}N_i(\mathbf x)\right|_{\mathbf x=\mathbf x_p}$$

表示网格形函数对空间坐标的梯度，在粒子位置求值。质量与动量转移为

$$m_i = \sum_p w_{ip}m_p,$$

$$m_i\mathbf v_i = \sum_p w_{ip}m_p\mathbf v_p.$$

这是最简的 PIC-style 动量传递；APIC 等方法还会传递粒子局部仿射速度，以减少旋转信息损失。材料本构关系根据形变梯度 $F_p$ 计算粒子的 Cauchy stress $\boldsymbol\sigma_p$。若 $V_p$ 是当前粒子体积，网格节点的内力为

$$\mathbf f_i^{\mathrm{int}} = -\sum_p V_p\boldsymbol\sigma_p \nabla w_{ip}.$$

负号来自内力对形变能的负梯度。加上外力后，只对 $m_i>0$ 的占用节点执行

$$\mathbf v_i^{n+1} = \mathbf v_i^n +\Delta t\, \frac{ \mathbf f_i^{\mathrm{int}} +\mathbf f_i^{\mathrm{ext}} }{m_i},$$

再处理固定边界、摩擦和碰撞。最简化的网格到粒子插值为

$$\mathbf v_p^{n+1} = \sum_i w_{ip}\mathbf v_i^{n+1},$$

网格速度在粒子处的梯度为

$$\nabla\mathbf v_p = \sum_i \mathbf v_i^{n+1} (\nabla w_{ip})^T.$$

它更新形变梯度和粒子位置：

$$F_p^{n+1} = (I+\Delta t\,\nabla\mathbf v_p)F_p^n,$$

$$\mathbf x_p^{n+1} = \mathbf x_p^n +\Delta t \sum_iw_{ip}\mathbf v_i^{n+1}.$$

$F_p$ 记录材料历史，$\boldsymbol\sigma_p(F_p)$ 把材料模型重新带回下一步网格内力，形成 P2G、网格更新、G2P 的闭环。实际 PIC、FLIP、APIC 等传输方式在数值耗散与噪声之间具有不同权衡。

粒子保存材料历史，网格避免直接维护不断扭曲的连接关系，因此 MPM 很适合雪、沙、泥浆和弹塑性体的大变形。它仍需解决网格接触、材料本构、时间步稳定性和高效稀疏存储等问题。