标签：应该 性能 ips 稳定性 问题 思路 常用 进一步 离散化
介绍常用的时间积分方法，及最终的求解过程。
0 物理系统描述
在物理引擎中，借助牛顿第二运动定律对系统进行描述，即
\[\begin{aligned}
\boldsymbol{f} &= \boldsymbol{f}(\boldsymbol{x}) \\frac{\partial\boldsymbol{v}_i}{\partial t} &= \frac{\boldsymbol{f}_i}{m_i} \\frac{\partial\boldsymbol{x}_i}{\partial t} &= \boldsymbol{v}_i
\end{aligned}
\]
有时候也会用 \(\boldsymbol{q}\) 来表示模型中节点的位置，那么系统描述即为：
\[\begin{aligned}
\boldsymbol{f} &= \boldsymbol{f}(\boldsymbol{x}) \\ddot{\boldsymbol{q}} &= \frac{\boldsymbol{f}}{\boldsymbol{M}}
\end{aligned}
\]
上述方程就是物理引擎中的 ODE 部分。该方程组的解法主要有显式（Forward Euler、Semi-implicit Euler）、隐式（Backward Euler）等。
1、时间积分方法
在仿真计算过程中，已知模型中节点在 \(t\) 时刻的位置、速度等信息，进一步求解其在 \(t+1\) 时刻的速度、位置。
1.1 显式时间积分（Explicit / Forward Euler）
计算方式如下：
\[\begin{aligned}
\boldsymbol{v}_{t+1} &= \boldsymbol{v}_{t} + \Delta t \frac{\boldsymbol{f}_{t}}{m}\\boldsymbol{x}_{t+1} &= \boldsymbol{x}_{t} + \Delta t \boldsymbol{v}_{t}
\end{aligned}
\]
在该方法中，由模型中节点的位置 \(\boldsymbol{x}_{t}\) 直接计算得到受力 \(\boldsymbol{f}_{t}\) ，进而可直接计算得到 \(t+1\) 时刻模型中节点的速度、位置。
1.2 半隐式积分（Explicit / Semi-implicit Euler aka. Symplectic Euler）
计算方式如下：
\[\begin{aligned}
\boldsymbol{v}_{t+1} &= \boldsymbol{v}_{t} + \Delta t \frac{\boldsymbol{f}_{t}}{m}\\boldsymbol{x}_{t+1} &= \boldsymbol{x}_{t} + \Delta t \boldsymbol{v}_{t+1}
\end{aligned}
\]
在该方法中，同样由模型中节点的位置 \(\boldsymbol{x}_{t}\) 直接计算得到受力 \(\boldsymbol{f}_{t}\) ，进而可直接计算得到 \(t+1\) 时刻模型中节点的速度、位置。
Tips：Forward Euler 和 Semi-implicit Euler 略微不同，在计算 \(\boldsymbol{x}_{t+1}\) 的时候，一个是用了 \(\boldsymbol{v}_{t}\)，另一个是用了 \(\boldsymbol{v}_{t+1}\)。现在，Forward Euler 用的较少，Semi-implicit Euler 用的较多一些。虽然，Semi-implicit Euler 只差别了一点点，但是准确性上会有本质性的提升。
1.3 仿真流程（显式积分）
在使用显式（Forward Euler or Semi-implicit Euler）进行仿真的时候，仿真流程有如下几个步骤：
- 计算节点受力 \(\boldsymbol{f}_t = \boldsymbol{f}(\boldsymbol{x}_t)\)
- 计算新的速度 \(\boldsymbol{v}_{t+1} = \boldsymbol{v}_t + \Delta t\frac{\boldsymbol{f}_t}{m}\)
- 碰撞检测（此时会更正速度）
- 计算新的位置 \(\boldsymbol{x}_{t+1} = \boldsymbol{x}_t + \Delta t \boldsymbol{v}_{t+1}\) （Semi-implicit Euler）
显式时间积分器的性能缺陷： Easy to explore
\[\Delta t \le c\sqrt{\frac{m}{k}} \quad (c \sim 1)
\]
关于稳定性 Stability 和爆炸 Explode 问题：
（略）
1.4 隐式积分（implicit Euler）
计算方式如下：
\[\begin{aligned}
\boldsymbol{v}_{t+1} &= \boldsymbol{v}_{t} + \Delta t \frac{\boldsymbol{f}_{t + 1}}{m}\\boldsymbol{x}_{t+1} &= \boldsymbol{x}_{t} + \Delta t \boldsymbol{v}_{t + 1}
\end{aligned}
\]
亦或者记作如下的形式
\[\begin{aligned}
\boldsymbol{x}_{t+1} &= \boldsymbol{x}_{t} + \Delta t \boldsymbol{v}_{t + 1} \\boldsymbol{v}_{t+1} &= \boldsymbol{v}_{t} + \Delta t \mathbf{M}^{-1} \boldsymbol{f}(\boldsymbol{x}_{t + 1})
\end{aligned}
\]
上述系统方程，可以转化成一个非线性的偏微分方程（PDE），通常有两种思路：（1）化简消去 \(\boldsymbol{x}_{t+1}\)；（2）化简消去 \(\boldsymbol{v}_{t+1}\)
（1）隐式积分求解 - 消去 \(\boldsymbol{x}_{t+1}\)
化简消去 \(\boldsymbol{x}_{t+1}\) 后，得到系统方程，即
\[\boldsymbol{v}_{t+1} = \boldsymbol{v}_{t} + \Delta t \mathbf{M}^{-1} \boldsymbol{f}(\boldsymbol{x}_{t} + \Delta t \boldsymbol{v}_{t + 1})
\]
这是一个关于 \(\boldsymbol{v}_{t+1}\) 的偏微分方程 PDE 。
（2）隐式积分求解 - 消去 \(\boldsymbol{v}_{t+1}\)
化简消去 \(\boldsymbol{v}_{t+1}\) 后，得到系统方程，即
\[\boldsymbol{x}_{t+1} = \boldsymbol{x}_{t} + \Delta t \boldsymbol{v}_t + \Delta t^2 \mathbf{M}^{-1} \boldsymbol{f}(\boldsymbol{x}_{t +1} )
\]
这是一个关于 \(\boldsymbol{x}_{t+1}\) 的偏微分方程 PDE 。
Tips：在这里，消去 \(\boldsymbol{x}_{t+1}\) 而保留 \(\boldsymbol{v}_{t+1}\) 的用意，应该是便于进行碰撞处理。在一些物理引擎中，会先采用 \(\boldsymbol{x}_{t}\) 或 \(\boldsymbol{x}_{t} + \Delta t \boldsymbol{v}_t\) 作为模型中节点的位置，进行碰撞检测。得到碰撞结果后，进一步更新/限制节点的速度 \(\boldsymbol{v}_{t+1}\) ，这样的好处好象是稳定性会好一些，不会出现节点的位置穿过碰撞界限等。
2 物理引擎中的 PDE 求解
如第 1 节中看到，在物理仿真中，通过空间上的离散化，计算得到了模型中节点上的受力 \(\boldsymbol{f}\) ；通过运动方程，描述模型的运动规律，得到了一组 ODE；对模型的运动状态在时间上进行离散，并通过显式/隐式积分器，可以得到一组 PDE。那么，最终，就需要进行 PDE 的求解。（更为复杂的，比如带约束等情况，后面会再整理。这里只描述最基本的模型运动仿真）
物理引擎中的时间积分方法及求解
标签：应该 性能 ips 稳定性 问题 思路 常用 进一步 离散化
原文地址：https://www.cnblogs.com/wghou09/p/13245923.html
