基于物理的建模入门（三）：能量函数与刚性
前面的话
本系列一共7篇，目前已经发布的前两篇得到了超出预期的关注，感谢各位知友！本来只是圈地自萌的行为，没想到收到了不错的反响。现在我更有动力坚持把这个系列更完啦，大家一起加油吧。总有一天我们也会有自己的可以媲美Unity和Unreal的商业游戏引擎哒~!
前两篇传送门：
基于物理的建模入门（一）：微分方程基础
基于物理的建模入门（二）：粒子系统动力学
这一篇提纲挈领地讲了在实践中如何处理能量函数不稳定性的问题。比较短，也没有很复杂的公式，主要是指导性的内容，比前两篇更下饭~。原文见：https://www.cs.cmu.edu/~baraff/pbm/energons.pdf
正文
在本讲义的前几节课中，我们指出，微分方程的数值解在实践中会出现两种主要类型的问题：精度和稳定性 (accuracy and stability)。精度是指在尝试对连续微分方程的解进行数值逼近/近似时引入的小误差。稳定性是指用于数值逼近/近似的方法是否收敛的问题。在某些工程环境中，高精度非常重要，但对于基于物理的计算机动画，准确性往往不如稳定性重要。模拟/仿真中的小错误在计算机图形学中基本上不会被注意到，但不稳定性通常会使仿真结果完全无法使用。此外，稳定性通常是仿真计算速度的限制因素，因为我们虽然通常可以增加微分方程求解器的时间步长以使仿真尽可能快地运行，但最多只能到稳定性所允许的极限。下面我们将检查由基于物理的建模导致的微分方程不稳定的一些原因，并讨论处理它们的方法。
为了便于说明，我们不考虑受 F = ma 定律支配的「正确的」 牛顿力学（牛顿：物理学不存在了！），而是分析更简单的受 F = mv 定律指导的物理模型。根据古希腊亚里士多德的理论，这一更简单的模型产生的「物理学」中，任何东西都不会移动，除非有力作用于它。这一模型的优点是我们可以直接从能量函数构造有趣的行为，而无需引入摩擦力和阻尼力。因此，该模型特别容易分析，但仍包含了我们了解不稳定性的来源所需的重要元素。
表现出导致不稳定性的关键问题的最简单的物理系统，是在二维中受势能 E(x, y) 影响的单位质量粒子。粒子上的力由向量 (\partial E/\partial x,\partial E/\partial y) 给出，因此使用亚里士多德物理学，粒子的运动方程就很简单： \partial X/\partial t = (\partial E/\partial x) 和 \partial Y/\partial t = (\partial E/\partial y) 。 当然，系统的行为取决于能量函数的精确性质。
假设我们虑上述系统在能量最小值附近的行为。如果我们考虑一个足够小的区域，我们只需要其泰勒展开中的前几项就可以合理地近似能量函数。泰勒展开式中的常数项不会影响粒子的行为，因为只有势能的相对变化才是显着的。梯度在最小值处消失，所以一阶导数项没有影响。因此，二阶导数项是对最小值附近的行为有重要影响的第一项。由于二阶导数项产生二次能量 (quadratic energy)，因此我们对观察系统在能量函数为二次时的行为十分感兴趣。
让我们从考虑二次能量函数 E = (x^2 + y^2)/2 开始。能量在0处具有最小值，其值为0。作用在粒子上的合力就是 F = \nabla E = (x,y) 。因此，作用在粒子上的力是从粒子位置回到原点的向量。在亚里士多德动力学下，粒子无论从那里开始，都会走一条沿直线路径走向原点。粒子越接近原点时会越减速，因为梯度变得越来越小。如果粒子被放置在原点，它的速度将为0，因为不存在梯度。
假设我们尝试使用欧拉方法计算粒子的运动。如果 P_0 表示粒子的起始位置，则其在第一个时间步长结束时的位置由 P_1 = P_0 + h\nabla E = (1 h)P_0 给出，其中 h 是时间步长。显然，如果 h 很小，则需要进行大量迭代才能使粒子接近(0, 0)处的最小值。随着 h 变大，粒子在一次迭代中移动得更远，因此仿真进行得更快。在 0
图1：对称能量函数F = (x^2+y^2)/2
如果能量函数是如上的对称二次方程，并且步长设置得当，欧拉方法可以很好地工作。我们怎么能对拥有「一步跳到最小值」能力的函数产生怨言呢？即使我们不知道要使用的正确步长，也应该能够在相对较小的步数内达到最小值。如果迭代对于所有步长都是稳定的，那当然更好，但现下的情形依旧可以被视为欧拉方法的成功。
当我们把最初的对称能量函数拿来，将其沿一个方向拉伸时，就会发生一些非常不同的事情了。假设我们有 E = (\epsilon x^2 + y^2 )/2 ，那么作用在粒子上的合力为 F = \nabla E = (\epsilon x,y) ，且不再指向最小值（除非 x = 0 、 y = 0 或 \epsilon = 1 ）。现在考虑如果我们使用欧拉方法会发生什么。我们从某个点 P_0 = (x_0, y_0) 开始。下一个点是 P_1 = ((1 h \epsilon)x_0, (1 h)y_0) 。请注意，在这种情况下，如果 h>2 ，欧拉迭代的 y 分量将发散。因此，即使在保持稳定性的最大可能步长下， x 的值每次迭代也仅减少 (1 2\epsilon) ，它描述的是时间常数与 1/\epsilon 成正比的指数衰减。
请注意，仅仅是对能量函数的简单缩放，就导致了这么深刻的变化。对称能量的情况下，如果选择正确的步长，我们单次迭代就可以达到最小值。 就算没有「正确」步长，只要我们选择一个稳定的步长，也将以与步长的倒数成正比的时间常数接近最小值。然而，随着能量的拉伸，即使我们选择可能的最大稳定步长，我们接近最小值的时间常数也只能与 1/\epsilon 成正比。这有多糟糕？要多糟有多糟。在实际问题中， 1/\epsilon 可能是 10^6 或 10^9 或更糟的数量级。显然这是一个严峻的难题。
为了理解这个难题的本质，让我们考虑一下拉伸的能量是什么样子的。随着 \epsilon 变得非常小，能量函数不再像一个对称的球，而越来越像一个陡峭的峡谷。在峡谷的底部，有一个非常小的坡度朝向最小值，但在其他任何地方，主要的坡度都朝向峡谷的底部。这样一来，使用欧拉方法的仿真将主要使粒子在峡谷两侧来回上下移动，同时非常非常缓慢地向最小值漂移，如图2所示。我们的问题在于，欧拉法基于能量函数的局部平面逼近，当能量函数开始看起来像一个高度弯曲的峡谷时，除非邻域非常小（因此步长也非常小），否则该逼近会变得非常差。结果就是，在这些情况下，欧拉方法变得极其缓慢。
图2：拉伸能量函数F = (\epsilon x^2+y^2)/2
请注意，当我们拉伸能量函数时，我们改变了微分方程真解 (true solutions) 的时间常数。在对称情况下，解是 x = x_0e^{t} ， y = y_0e^{t} 。在拉伸的情况下，解是 x = x_0e^{\epsilon t} ， y = y_0e^{t} 。对称情况下的时间常数相同，但在拉伸情况下，时间常数的比率为 e （注：此处是否应为 \epsilon ）。这是微分方程呈「刚性」的常见原因。如我们之前所见，最大步长由具有最快时间常数的过程确定，因此这会使较长时间常数的仿真慢得无法忍受。
为了应对刚性的难题，我们有两个主要的选择。可以使用基于更准确的能源「地貌 (landscape)」局部模型的方法，或者尝试制定 (formulate) 一组不同但相关的微分方程组来解决原始问题。如果可能，最好重新制定问题以消除刚性，但通常除了（硬着头皮）模拟刚性方程外别无选择。
在基于物理的建模的许多用途中，人们一直很想使用各种弹簧来强制约束。问题是，要用弹簧满足约束，必须使弹簧常数非常大，而这会导致刚性方程组。在许多情况下，拉格朗日乘数或变量的改变可以让使用刚性小得多的微分方程获得同样结果变得可能。
如果无法重新制定微分方程以消除刚性，则解决它们的唯一可行方法是尝试其他技术，用比欧拉方法更好的局部力模型来对抗刚性。通常，它们被称为隐式方法 (implicit methods)，我们将在后续「计算机图形学的连续介质动力学」中详细处理其中一种方法。
隐式微分方程求解器的细节在各种文本中都有处理。在这里，我们试图给出一些关于它们如何工作的直觉。请记住，对称能量的梯度指向最小值，而拉伸能量的梯度则不是。这是我们难题的很大一部分。如果我们继续沿梯度方向移动，我们到达最小值会非常缓慢。因此，为了迈出非常大的步长，微分方程求解器应该直接向最小值移动，即使面对拉伸的能量。微分方程求解器如何知道要走哪条路才能达到最小值？通常情况下，这有多困难都可能，但对于二次能量来说，则特别简单。二次能量具有线性梯度，而最小值出现在梯度消失时，因此可以通过求解一组线性方程来找到出现这种情况的点。对于拉伸能量的情况，梯度为 (\epsilon x, y) ，很明显它仅在点(0, 0)处消失。因此，通过求解适当的线性方程，微分方程求解器可以向二次能量函数的最小值移动。
通常，点 x_0 附近能量函数梯度的最佳线性估计由表达式 \nabla E(x) \approx \nabla E(x_0) + H \cdot (x x_0) 给出，其中 H 是二阶导数矩阵，即由 H = \partial^2E/\partial x_i \partial x_j 给出的黑塞矩阵 (Hessian matrix)。由于 H 是对称的，通过基的适当正交变化，总能将 H 变成对角矩阵。对于对称能量函数， H 已经是对角矩阵：它就是个2x2的单位矩阵。对于拉伸能量， H 也已经对元素 (\epsilon, 1) 成对角矩阵。在对角化H的基上， H 的条件数 (condition number) 定义为 H 的模长最大和最小的对角线元素的模长之比。如果条件数很小，那么我们有一个几乎对称的能量函数，方程不会很刚性。如果条件数很大，则意味着如果我们考虑给到H中最小和最大对角线元素的两个维度，则我们处于 1/\epsilon 对应于条件数的拉伸能量函数的情形中。显然，如果条件数很大，我们很可能陷入麻烦。
对称和拉伸能量函数指出求解微分方程的典型困难题以及典型解。如果我们查看任何平滑能量函数的合适邻域，我们可以将其行为近似为二次函数。如果我们选择正确的坐标，那么黑塞矩阵将是对角矩阵，因此我们可以通过黑塞矩阵的对角线元素来表征这些坐标中的局部行为。如果这些项的模长都近似相同，那么能量函数则是在合理范围内对称的（至少在局部如此），我们的数值问题也就相对不难。如果这些项模长差异很大，那么如果我们选择对应于模长最小项和最大项的两个坐标，我们就会得到一个局部看起来像一条陡峭峡谷的能量函数。为了在这种形状的能量函数中走出很大的步长，我们需要使用隐式方法来检查矩阵 H 以找到峡谷的轴，使它可以顺峡谷而下，直接向最小值移动。我们将在下一篇讲义「计算机图形学的连续介质动力学简介」中看到这样一种方法的示例。
系列传送门：
基于物理的建模入门（一）：微分方程基础
基于物理的建模入门（二）：粒子系统动力学
基于物理的建模入门（三）：能量函数与刚性