【质点弹簧】质点位置积分

在物理模拟中，物理系统每帧都需要根据粒子当前位置，速度，加速度等信息计算粒子下个阶段相关信息，而这种一点一点计算位置的操作被称作对粒子位置的积分。

具体而言实现这种积分效果的有两种方法：

欧拉积分法

https://zhuanlan.zhihu.com/p/355170943

欧拉积分法是根据现实中的物理公式来求解下一帧位置的，具体而言分三种：

显式欧拉方法（完全根据当前帧的位置、速度、加速度数据计算）
半隐式欧拉方法（位置和加速度采用当前帧数据，速度采用下一帧数据）
隐式欧拉方法（加速度和速度都采用下一帧数据）

这三种方法的计算难度和结果误差各不相同，隐式欧拉最难，但精度最高，显式最简单，但精度最差。

半隐式欧拉法的计算方法最符合直觉，因为完全就是按照现实的物理法则来计算位置：

\[p(t+\Delta t)= v(t)\Delta t + a(t)\Delta t^2 + p(t) \]

但计算机与现实不同，计算机的每帧之间存在时间差，这些时间空白就导致了计算精度的损失。而且隐式欧拉法需要保存每帧都带有误差的速度数据，这就导致误差会随时间越来越大，最终导致物理效果完全崩溃。

以质点弹簧模型为例，误差可能导致位置的计算跑过头，可一旦跑过头，下次就要跑的更多，弹簧晃荡越来越大，最终浮点计算崩溃。因此相应的一种解决办法就是，给质点添加阻力并减少弹力，使质点每次都是离目标差点而不是跑过头。随时间流逝，最终质点肯定也是能抵达目的地的，但也因此物理效果明显有迟钝感，像东西都在水里一样。

欧拉积分法真心不好用，还得是下面的 Verlet 积分法。

Verlet 积分法

https://zhuanlan.zhihu.com/p/64278344

Verlet 积分法是一种，常用作物理模拟中的计算方法。相比欧拉积分法，这种方法误差更少，更稳定（依赖的物理数据少，随时间累计的误差也少）。

其计算公式如下（t 为当前时间）：

\[p(t+\Delta t) = p(t)+(p(t)-p(t-\Delta t))+a(t)\Delta t^2 \]

\[v(t) = \frac{p(t+\Delta t)-p(t-\Delta t)}{2\Delta t} \]

推导

粒子最重要的属性是位移，速度、加速度的最终目的都是为了计算位移，而 Verlet 就是根据位移函数的泰勒展开式推导而来。

将位置、速度、加速度看成是对时间 $t$ 的函数，则可将它们以及它们的关系写成如下形式：

当前位置 $= p(t)$
当前速度 $= v(t) =p'(t) $
当前加速度 $= a(t) = v'(t)$

那么位置函数的泰勒展开式，可以写成如下形式

$ \begin{aligned} p(t) &=p(t_0)+\frac{p'(t_0)}{1!}(t-t_0)+\frac{p''(t_0)}{2!}(t-t_0)^2+\frac{p'''(t_0)}{3!}(t-t_0)^3+\dots\\ &=p(t_0)+v(t_0)(t-t_0)+\frac{a(t_0)}{2}(t-t_0)^2+\frac{p'''(t_0)}{6}(t-t_0)^3+\dots\\ \end{aligned} $

若 $t_0$ 为当前时间，则其下一阶段（$t=t_0+\Delta t$）和上一阶段（$t=t_0-\Delta t$）的位置函数展开式如下：

$ \begin{aligned} p(t_0+\Delta t) &=p(t_0)+v(t_0)(t_0+\Delta t-t_0)+\frac{a(t_0)}{2}(t_0+\Delta t-t_0)^2+\frac{p'''(t_0)}{6}(t_0+\Delta t-t_0)^3+\dots\\ &=p(t_0)+v(t_0)\Delta t+\frac{a(t_0)}{2}\Delta t^2+\frac{p'''(t_0)}{6}\Delta t^3+O(\Delta t^4)\\ \end{aligned} $

$ \begin{aligned} p(t_0-\Delta t) &=p(t_0)+v(t_0)(t_0-\Delta t-t_0)+\frac{a(t_0)}{2}(t_0-\Delta t-t_0)^2+\frac{p'''(t_0)}{6}(t_0-\Delta t-t_0)^3+\dots\\ &=p(t_0)-v(t_0)\Delta t+\frac{a(t_0)}{2}\Delta t^2-\frac{p'''(t_0)}{6}\Delta t^3+O(\Delta t^4)\\ \end{aligned} $

将上述两式相加化简可得：

$ \begin{aligned} p(t_0+\Delta t)+p(t_0-\Delta t) &=2p(t_0)+a(t_0)\Delta t^2+2O(\Delta t^4)\\ p(t_0+\Delta t) &= 2p(t_0)+a(t_0)\Delta t^2-p(t_0-\Delta t)+O(\Delta t^4) \end{aligned} $

因此忽略误差后并稍加整理可以得下一阶段的粒子位置公式如下（t 表示当前时间）：

\[p(t+\Delta t) = p(t)+(p(t)-p(t-\Delta t))+a(t)\Delta t^2 \]

该公式只需粒子当前的位置信息和加速度，以及上一帧的位置即可计算下一帧的位置。当前帧的加速度通过牛顿第二定律即可获得，完全不依赖速度，因此抗干扰能力强。此外误差只有$O(\Delta t^4)$，也比欧拉积分少。

计算速度

利用位移除时间即可得到速度，而恰好该方法中我们需要记录上一帧的位置信息，因此该速度公式可用，并且我们还能计算出误差：

$ \begin{aligned} \frac{p(t_0+\Delta t)-p(t_0-\Delta t)}{2\Delta t}=\frac{2v(t_0)\Delta t+2O(\Delta t^3)}{2\Delta t}=v(t_0)+O(\Delta t^2) \end{aligned} $

故速度公式如下，其误差为 $O(\Delta t^2)$：

\[v(t) = \frac{p(t+\Delta t)-p(t-\Delta t)}{2\Delta t} \]