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

转载于物理引擎中的时间积分方法及求解 - wghou09 - 博客园

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

---

0 物理系统描述

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

f∂vi∂t∂xi∂t=f(x)=fimi=vif=f(x)∂vi∂t=fimi∂xi∂t=vi

有时候也会用 qq 来表示模型中节点的位置，那么系统描述即为：

fq¨=f(x)=fMf=f(x)q¨=fM

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

---

1、时间积分方法

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

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

计算方式如下：

vt+1xt+1=vt+Δtftm=xt+Δtvtvt+1=vt+Δtftmxt+1=xt+Δtvt

在该方法中，由模型中节点的位置 xtxt 直接计算得到受力 ftft ，进而可直接计算得到 t+1t+1 时刻模型中节点的速度、位置。

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

计算方式如下：

vt+1xt+1=vt+Δtftm=xt+Δtvt+1vt+1=vt+Δtftmxt+1=xt+Δtvt+1

在该方法中，同样由模型中节点的位置 xtxt 直接计算得到受力 ftft ，进而可直接计算得到 t+1t+1 时刻模型中节点的速度、位置。

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

1.3 仿真流程（显式积分）

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

• 计算节点受力 ft=f(xt)ft=f(xt)
• 计算新的速度 vt+1=vt+Δtftmvt+1=vt+Δtftm
• 碰撞检测（此时会更正速度）
• 计算新的位置 xt+1=xt+Δtvt+1xt+1=xt+Δtvt+1 （Semi-implicit Euler）

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

Δt≤cmk−−−√(c∼1)Δt≤cmk(c∼1)

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

（略）

1.4 隐式积分（implicit Euler）

计算方式如下：

vt+1xt+1=vt+Δtft+1m=xt+Δtvt+1vt+1=vt+Δtft+1mxt+1=xt+Δtvt+1

亦或者记作如下的形式

xt+1vt+1=xt+Δtvt+1=vt+ΔtM−1f(xt+1)xt+1=xt+Δtvt+1vt+1=vt+ΔtM−1f(xt+1)

上述系统方程，可以转化成一个非线性的偏微分方程（PDE），通常有两种思路：（1）化简消去 xt+1xt+1；（2）化简消去 vt+1vt+1

（1）隐式积分求解 - 消去 xt+1xt+1

化简消去 xt+1xt+1 后，得到系统方程，即

vt+1=vt+ΔtM−1f(xt+Δtvt+1)vt+1=vt+ΔtM−1f(xt+Δtvt+1)

这是一个关于 vt+1vt+1 的偏微分方程 PDE 。

（2）隐式积分求解 - 消去 vt+1vt+1

化简消去 vt+1vt+1 后，得到系统方程，即

xt+1=xt+Δtvt+Δt2M−1f(xt+1)xt+1=xt+Δtvt+Δt2M−1f(xt+1)

这是一个关于 xt+1xt+1 的偏微分方程 PDE 。

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

---

2 物理引擎中 PDE 的求解

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

（在物理引擎中，会得到各种 ODE、PDE、线性系统、非线性系统，名称上比较混乱。）其中涉及的求解方法也非常多，比如，线性化、牛顿法、Jacobin、CG（Conjugate Gradient）、优化隐式方法等等，非常多样

2.1 线性化（one step of Newton's method）

简单地求解，可以实现一步牛顿法，将 ff 在 xtxt 处一阶泰勒展开，得到如下：

（1）速度上的求解

vt+1=vt+ΔtM−1[f(xt)+∂f∂x(xt)Δtvt+1]vt+1=vt+ΔtM−1[f(xt)+∂f∂x(xt)Δtvt+1]

操作之后，系统变成了线性系统。整理之后，为

[I−Δt2M−1∂f∂x(xt)]vt+1=vt+ΔtM−1f(xt)[I−Δt2M−1∂f∂x(xt)]vt+1=vt+ΔtM−1f(xt)

（2）位置上的求解

xt+1=xt+Δtvt+Δt2M−1[f(xt)+∂f∂x(xt)Δtvt]xt+1=xt+Δtvt+Δt2M−1[f(xt)+∂f∂x(xt)Δtvt]

这里不是很确定，还需要再对照资料确认一下。

如上所述，线性化之后会得到如下的线性系统：

AbAvt+1=[I−Δt2M−1∂f∂x(xt)]=vt+ΔtM−1f(xt)=bA=[I−Δt2M−1∂f∂x(xt)]b=vt+ΔtM−1f(xt)Avt+1=b

对该线性系统的求解，可以采用 Jacobin / Gauss-Seidel iterations，或者 Conjugate gradients 等方法进行求解。

Jacobin iteration 求解线性系统的方法

单步的 Jacobin 迭代过程如下所示：

A = []
x = []
new_x = []
b = []

@ti.kernel
def iterate():
  for i in range(n):
    r = b[i]
    for j in range(n):
      if i != j:
        r -= A[i,j] * x[j]

    new_x[i] = r / A[i,i]

  for i in range(n):
    x[i] = new_x[i]

Jacobin 迭代的思想就是，每次只让一行（一个点）满足该方程，依次计算完成，就能让所有的点都（曾经）满足方程。多迭代几次，就能接近线性方程组的解了。（大概是这么个意思吧）

2.2 线性化（one step of Newton's method） - 进阶

在上节所述的线性系统中，加入一个系数 ββ ，那么，就会得到如下线性系统：

[I−βΔt2M−1∂f∂x(xt)]vt+1=vt+ΔtM−1f(xt)[I−βΔt2M−1∂f∂x(xt)]vt+1=vt+ΔtM−1f(xt)

其中，当 β=0β=0 时，为 forward / semi-implicit Euler (explicit integrator)
当 β=1/2β=1/2 时，为 middle-point (implicit integrator)
当 β=1β=1 时，为 backward Euler (implicit integrator)