物理引擎中的时间积分方法及求解

Posted 2020-12-14 wghou09

介绍常用的时间积分方法，及最终的求解过程。

---

0 物理系统描述

在物理引擎中，借助牛顿第二运动定律对系统进行描述，即

[egin{aligned} oldsymbol{f} &= oldsymbol{f}(oldsymbol{x}) \frac{partialoldsymbol{v}_i}{partial t} &= frac{oldsymbol{f}_i}{m_i} \frac{partialoldsymbol{x}_i}{partial t} &= oldsymbol{v}_i end{aligned} ]

有时候也会用 (oldsymbol{q}) 来表示模型中节点的位置，那么系统描述即为：

[egin{aligned} oldsymbol{f} &= oldsymbol{f}(oldsymbol{x}) \ddot{oldsymbol{q}} &= frac{oldsymbol{f}}{oldsymbol{M}} end{aligned} ]

上述方程就是物理引擎中的 ODE 部分。该方程组的解法主要有显式（Forward Euler、Semi-implicit Euler）、隐式（Backward Euler）等。

---

1、时间积分方法

在仿真计算过程中，已知模型中节点在 (t) 时刻的位置、速度等信息，进一步求解其在 (t+1) 时刻的速度、位置。

1.1 显式时间积分（Explicit / Forward Euler）

计算方式如下：

[egin{aligned} oldsymbol{v}_{t+1} &= oldsymbol{v}_{t} + Delta t frac{oldsymbol{f}_{t}}{m}\boldsymbol{x}_{t+1} &= oldsymbol{x}_{t} + Delta t oldsymbol{v}_{t} end{aligned} ]

在该方法中，由模型中节点的位置 (oldsymbol{x}_{t}) 直接计算得到受力 (oldsymbol{f}_{t}) ，进而可直接计算得到 (t+1) 时刻模型中节点的速度、位置。

1.2 半隐式积分（Explicit / Semi-implicit Euler aka. Symplectic Euler）

计算方式如下：

[egin{aligned} oldsymbol{v}_{t+1} &= oldsymbol{v}_{t} + Delta t frac{oldsymbol{f}_{t}}{m}\boldsymbol{x}_{t+1} &= oldsymbol{x}_{t} + Delta t oldsymbol{v}_{t+1} end{aligned} ]

在该方法中，同样由模型中节点的位置 (oldsymbol{x}_{t}) 直接计算得到受力 (oldsymbol{f}_{t}) ，进而可直接计算得到 (t+1) 时刻模型中节点的速度、位置。

Tips：Forward Euler 和 Semi-implicit Euler 略微不同，在计算 (oldsymbol{x}_{t+1}) 的时候，一个是用了 (oldsymbol{v}_{t})，另一个是用了 (oldsymbol{v}_{t+1})。现在，Forward Euler 用的较少，Semi-implicit Euler 用的较多一些。虽然，Semi-implicit Euler 只差别了一点点，但是准确性上会有本质性的提升。

1.3 仿真流程（显式积分）

在使用显式（Forward Euler or Semi-implicit Euler）进行仿真的时候，仿真流程有如下几个步骤：

• 计算节点受力 (oldsymbol{f}_t = oldsymbol{f}(oldsymbol{x}_t))
• 计算新的速度 (oldsymbol{v}_{t+1} = oldsymbol{v}_t + Delta tfrac{oldsymbol{f}_t}{m})
• 碰撞检测（此时会更正速度）
• 计算新的位置 (oldsymbol{x}_{t+1} = oldsymbol{x}_t + Delta t oldsymbol{v}_{t+1}) （Semi-implicit Euler）

显式时间积分器的性能缺陷： Easy to explore

[Delta t le csqrt{frac{m}{k}} quad (c sim 1) ]

关于稳定性 Stability 和爆炸 Explode 问题：

（略）

1.4 隐式积分（implicit Euler）

计算方式如下：

[egin{aligned} oldsymbol{v}_{t+1} &= oldsymbol{v}_{t} + Delta t frac{oldsymbol{f}_{t + 1}}{m}\boldsymbol{x}_{t+1} &= oldsymbol{x}_{t} + Delta t oldsymbol{v}_{t + 1} end{aligned} ]

亦或者记作如下的形式

[egin{aligned} oldsymbol{x}_{t+1} &= oldsymbol{x}_{t} + Delta t oldsymbol{v}_{t + 1} \boldsymbol{v}_{t+1} &= oldsymbol{v}_{t} + Delta t mathbf{M}^{-1} oldsymbol{f}(oldsymbol{x}_{t + 1}) end{aligned} ]

上述系统方程，可以转化成一个非线性的偏微分方程（PDE），通常有两种思路：（1）化简消去 (oldsymbol{x}_{t+1})；（2）化简消去 (oldsymbol{v}_{t+1})

（1）隐式积分求解 - 消去 (oldsymbol{x}_{t+1})

化简消去 (oldsymbol{x}_{t+1}) 后，得到系统方程，即

[oldsymbol{v}_{t+1} = oldsymbol{v}_{t} + Delta t mathbf{M}^{-1} oldsymbol{f}(oldsymbol{x}_{t} + Delta t oldsymbol{v}_{t + 1}) ]

这是一个关于 (oldsymbol{v}_{t+1}) 的偏微分方程 PDE 。

（2）隐式积分求解 - 消去 (oldsymbol{v}_{t+1})

化简消去 (oldsymbol{v}_{t+1}) 后，得到系统方程，即

[oldsymbol{x}_{t+1} = oldsymbol{x}_{t} + Delta t oldsymbol{v}_t + Delta t^2 mathbf{M}^{-1} oldsymbol{f}(oldsymbol{x}_{t +1} ) ]

这是一个关于 (oldsymbol{x}_{t+1}) 的偏微分方程 PDE 。

Tips：在这里，消去 (oldsymbol{x}_{t+1}) 而保留 (oldsymbol{v}_{t+1}) 的用意，应该是便于进行碰撞处理。在一些物理引擎中，会先采用 (oldsymbol{x}_{t}) 或 (oldsymbol{x}_{t} + Delta t oldsymbol{v}_t) 作为模型中节点的位置，进行碰撞检测。得到碰撞结果后，进一步更新/限制节点的速度 (oldsymbol{v}_{t+1}) ，这样的好处好象是稳定性会好一些，不会出现节点的位置穿过碰撞界限等。

---

2 物理引擎中的 PDE 求解

如第 1 节中看到，在物理仿真中，通过空间上的离散化，计算得到了模型中节点上的受力 (oldsymbol{f}) ；通过运动方程，描述模型的运动规律，得到了一组 ODE；对模型的运动状态在时间上进行离散，并通过显式/隐式积分器，可以得到一组 PDE。那么，最终，就需要进行 PDE 的求解。（更为复杂的，比如带约束等情况，后面会再整理。这里只描述最基本的模型运动仿真）