关于EI8703-计算机动画原理与技术课程
SJTU-EI8703-计算机动画原理与技术 课程笔记
1. ODE（常微分方程）的数值求解
1.1. 常微分方程（Ordinary Differential Equation，简称 ODE）
Ordinary 指单变量，即一个独立变量，在计算机动画中一般是时间 t 。
核心抽象定义：\(\dot{x}(t) = f(t, x(t))\)
在普适的数学语境下，一个常微分方程（ODE）就是一个等式，它表达了“某个未知函数的变化率，等于该函数当前状态的某个已知映射”。
- x(t)：未知的状态函数 (Unknown function that evaluates the state given time)
- 在最高抽象下，x 就是一个存在于 n 维状态空间（State Space）中的向量：\(x \in \mathbb{R}^n\)。
- 这个函数长什么样是要“解”出来的目标。
- f：已知的映射函数 (Known function / Vector Field)
- 这是一个我们确切知道的数学规则：\(f: \mathbb{R}^n \rightarrow \mathbb{R}^n\)。
- 只要输入当前的状态 x(t) ，函数会输出当前状态应该有的变化趋势。
- 在几何意义上，f 在状态空间中定义了一个向量场 (Vector Field)。
- \(\dot{x}(t)\)：未知函数关于 t 的导数 (Time derivative of the unknown function)
- 这就是状态 x 随 t 的变化率，即 \(\frac{dx}{dt}\)。
- 符号解法 (Symbolic/Analytical solutions)： 目标是运用数学技巧（如积分变换）推导出一个完美的、精确的数学公式。只要代入任意 t，就能得出绝对准确的结果。
- 数值解法 (Numerical solutions)： 计算机动画和物理引擎的核心，把变化规律 f 降级当成一个黑盒 (Black box)。不需要知道黑盒里是什么公式，只要向黑盒输入“现在的 t 和状态 x”，黑盒就会输出“当前的变化率（也就是导数/速度/加速度）”。然后利用这个变化率，一小步一小步地向前推算。
1.2. 初值问题（Initial Value Problem，IVP）
[定义 1] 含初值条件 \(x(t_0) = x_0\) 的常微分方程问题：
\[\begin{cases} \dot{x}(t) = f(t,x(t)), & t \ge t_0 \\ x(t_0) = x_0 \end{cases} \quad (1)\]
为初值问题 (Initial value problem, IVP)。
[定理 1] 若函数 f(t,x) 关于 x 满足 Lipschitz 条件，即存在常数 L > 0，使得对于任意的 \(t \ge t_0\)，任意的 x 和 \(\hat{x}\)，有：
\[|f(t,x) – f(t,\hat{x})| \le L|x – \hat{x}|\]
则常微分方程初值问题(1)存在唯一的解。
当我们知道一个常微分方程初值问题存在唯一解后，我们下一个要考虑的问题就是它的敏感性问题，即考虑初值发生扰动对结果产生的影响。
[定义 2] 对于常微分方程初值问题 (1)，考虑初值 \(x_0\) 的扰动使问题的解 x(t) 发生偏差的情形。当 \(t \to \infty\) 时，如果 x(t) 的偏差被控制在有界的范围内，则称该初值问题是稳定 (Stable) 的，否则该初值问题是不稳定 (Unstable) 的。进一步地，如果还有 \(t \to \infty\) 时，x(t) 的偏差收敛到 0，则称该初值问题是渐近稳定 (Asymptotically stable) 的。
- 注 1： 渐近稳定是比稳定更强的概念，即若一个问题是渐近稳定的，则必然是稳定的。
- 注 2： 对于不稳定的初值问题，初始数据的扰动将使 \(t \to \infty\) 时的误差结果无穷大。因此，为了保证数值求解的有效性，初值问题的稳定性至关重要。
1.3. 显式数值方法（Explicit Numeric Method）
在显式数值方法的过程中，未知数（未来的状态）已经被完全孤立在等式的左边。只需要把右边已知的数字代入公式，做简单的基础运算，就可以解得下一帧的结果。
假设我们要计算下一个时间步 \(t_{n+1}\) 的状态 \(x_{n+1}\)，步长为 \(\Delta t\)。显式方法的通用形式是：
\[x_{n+1} = x_n + \Delta t \cdot \Phi(t_n, x_n)\]
这里的 \(\Phi\) 是根据已知状态 \(x_n\) 计算出的，代表该时间步内的综合代表性斜率 / 变化率。
1.3.1. 龙格-库塔法（Runge-Kutta Methods, RK）总述
龙格-库塔法不是单一的算法，而是一整个显式算法家族。它的核心思想是：在 \(t_n\) 到 \(t_{n+1}\) 的这一个时间步 \(\Delta t\) 内，多次试探并采集不同点位的瞬时斜率（变化率 k），然后为这些斜率分配不同的权重进行加权平均，从而得到极其精准的综合斜率 \(\Phi\)。
其通用公式可以表示为：
\[x_{n+1} = x_n + \Delta t \cdot (a_1k_1 + a_2k_2 + \cdots + a_mk_m)\]
其中：
- \(a_i\) 是分配给各个斜率的权重，且所有权重相加等于 1（即 \(\sum a_i = 1\)）。
- \(k_1, k_2, \dots, k_m\) 是在区间内不同试探点通过调用“黑盒函数” f 算出的瞬时斜率。试探的次数越多、权重分配越合理，方法的阶数（准确度）就越高。
1.3.2. 一阶 RK（RK1）：显式/前向 欧拉（Explicit / Forward Euler）
显式欧拉法是最基础的RK 方法，它在区间内只试探了 1 次，即假设整个 \(\Delta t\) 时间段内的变化率都完全等于起点的变化率（权重 \(a_1 = 1\)）。
斜率采集：\[k_1 = f(t_n, x_n)\]
迭代公式：\[x_{n+1} = x_n + \Delta t \cdot k_1\]
由于只采样了一次，该方法仅具有 1 阶准确度，局部截断误差较大。
1.3.3. 二阶 RK（RK2）：中点法（Mid-Point Method）
为了提高精度，中点法采取了“先试探半步，再走整步”的策略，在区间内试探了 2 次。它完全舍弃了起点的斜率权重（\(a_1 = 0\)），将 100% 的权重交给了区间中点的斜率（\(a_2 = 1\)），因为数学上中点斜率更能代表整段区间的平均变化趋势。
斜率采集：\[\begin{aligned} k_1 &= f(t_n, x_n) \quad \text{(起点斜率)} \\ x_{mid} &= x_n + \frac{\Delta t}{2} \cdot k_1 \quad \text{(中点)} \\ k_2 &= f\left(t_n + \frac{\Delta t}{2}, x_{mid}\right) \quad \text{(中点斜率)} \end{aligned}\]
迭代公式：\[x_{n+1} = x_n + \Delta t \cdot k_2\]
它只需要多调用一次黑盒函数 f，就能将局部截断误差降低至 \(O(\Delta t^3)\)，具有 2 阶准确度。
1.3.4. 四阶 RK（RK4）：经典的龙格-库塔法
在单步内试探了 4 次（起点 1 次，中点 2 次，终点 1 次），并采用经典的 1:2:2:1 权重分配法则（即 \(\frac{1}{6}, \frac{2}{6}, \frac{2}{6}, \frac{1}{6}\)）。
斜率采集：\[\begin{aligned} k_1 &= f(t_n, x_n) \\ k_2 &= f\left(t_n + \frac{\Delta t}{2}, x_n + \frac{\Delta t}{2} k_1\right) \\ k_3 &= f\left(t_n + \frac{\Delta t}{2}, x_n + \frac{\Delta t}{2} k_2\right) \\ k_4 &= f(t_n + \Delta t, x_n + \Delta t \cdot k_3) \end{aligned}\]
迭代公式：\[x_{n+1} = x_n + \frac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4)\]
它的局部截断误差极小，达到了 4 阶准确度。
1.3.5. （特例）半隐式/辛 欧拉（Semi-Implicit / Symplectic Euler）
半隐式欧拉法不属于标准的龙格-库塔家族，但在游戏物理引擎中占据着极其核心的地位。
标准的 RK 家族严格遵循“完全用旧状态推导新状态”的纯显式套路。而半隐式欧拉法在更新物理状态（通常拆分为位置 x 和速度 v）时，打破了这一规则：它先算出下一帧的“新速度”，然后立刻利用这个“未来速度”来更新位置。
迭代公式：\[\begin{aligned} v_{n+1} &= v_n + \Delta t \cdot f_v(t_n, x_n, v_n) \\ x_{n+1} &= x_n + \Delta t \cdot v_{n+1} \end{aligned}\](注：式中的 \(f_v\) 代表计算速度变化率/加速度的黑盒函数。)
核心优势（能量守恒）： 虽然它的计算代价和 RK1 完全一样（1 阶准确度），但在模拟弹簧、单摆或轨道运动等保守力场时，RK 家族的方法（即使是 RK4）最终都会导致系统能量衰减或发散爆炸。而半隐式欧拉法属于辛积分器 (Symplectic Integrators)，它具有神奇的保结构特性，能够近乎完美地保持系统的能量守恒。
1.4. 数值方法的稳定性与误差分析
在使用数值方法求解初值问题时，除了需要考虑原问题本身的稳定性，还必须重点考虑数值方法本身的稳定性和误差。
1.4.1. 数值稳定性 (Numerical Stability)
[定义 3] 利用某个数值方法求解初值问题时，若在节点 \(t_n\) 上的函数近似值存在扰动 \(\delta_n\)，其引起的后续节点上的误差 \(\delta_m\)（其中 m > n）均不超过 \(\delta_n\)，即：
\[|\delta_m| \le |\delta_n|, \quad m > n\]
则称该方法是数值稳定 (Numerically stable) 的。
- 注 ： 此处的稳定性针对的是数值求解过程中的近似误差，这与 1.2 节中“原初值问题自身的稳定性”是两个独立的概念。
- 例： 单步前向欧拉法（显式欧拉）要保证数值稳定，必须满足 \(|1 + \lambda \Delta t| \le 1\)。这在物理仿真中的直观体现就是：采样步长 \(\Delta t\) 必须足够小，否则系统能量会发散爆炸。
1.4.2. 局部截断误差与准确度 (Local Truncation Error)
[定义 4] 设常微分方程初值问题的数值解法为 \(x_{n+1} = G(x_{n+1}, x_n, x_{n-1}, \cdots, x_{n-k})\)，假设其中之前的步骤 \(x_{n-i} \ (0 \le i \le k)\) 均等于完美的客观真实值 \(x(t_{n-i})\)，则我们称：
\[l_{n+1} = x(t_{n+1}) – x_{n+1}\]
为该方法的局部截断误差 (Local truncation error)。
局部截断误差的本质是“在假设之前的步骤完全准确的前提下，单次迭代计算所产生的新误差”，因此严格的定义必须是真实值 \(x(t_{n+1})\) 减去数值计算值 \(x_{n+1}\)。
[定义 5] 若一个求解常微分方程初值问题的数值解法，其局部截断误差 \(l_{n+1} = O(\Delta t^{p+1})\)，则称该方法具有 p 阶准确度。
注： \(O(\Delta t^{p+1})\) 表示和 \(\Delta t^{p+1}\) 同阶的小量，即存在常数 A，使得 \(\lim_{\Delta t \to 0} \frac{O(\Delta t^{p+1})}{\Delta t^{p+1}} = A\)。
核心结论： 只要方法的阶数 \(p \ge 1\)，随着仿真步长 \(\Delta t\) 的不断减小，最终的结果误差都将收敛到零。