更精确的模拟：verlet积分详解

[phenomLi]

阅读matter.js源码看到刚体update部分时发现了一种新的积分方法，赶紧学习总结下来。

注意，今天这篇比较硬核，需要一点点高数基础。

欧拉积分的缺陷

大多数物理引擎，包括box2d，在进行刚体位置更新时，通常使用以下方法：

// 更新速度
body.velocity += body.acceleration*dt;
// 更新位置 
body.position += body.velocity*dt;

这种我们经常在高中物理课本看到的速度和位置计算方法叫做欧拉积分法（Euler integration）。欧拉积分最大的好处就是简单，易于实现，性能好。

但是在中学，我们学习的都是匀加速运动，即加速度不变的运动，欧拉积分在恒定加速度的情况下表现良好。然而现代物理引擎常常包含复杂的碰撞，约束和摩擦力等，此时的加速度是不断变化的。在数学角度上，欧拉也可以求解变加速运动，因为加速度再怎么变化，位移s关于时间t的函数始终都是一个连续函数，而连续函数必定可积。但是在计算机中情况就不一样了，计算机模拟运动是离散的，物理引擎会设定一个步长dt（通常为1/60）作为每次模拟的时间间隔，也就是每隔dt时间，物理引擎便执行一次模拟，更新物理数据。

欧拉积分用当前时刻的值y(t) + y'(t)*dt来求y(t + 1)时刻的值。这种方法主要的误差来源于： 离散情况下dt区间内的y'(t)是变化的，而这里我们都用y(t)替代了，所以欧拉求得的y(t + 1)只能是近似值。 所以当模拟时间长了之后，欧拉积分就会不够精确，甚至不够稳定。

verlet积分法

Verlet积分法在wiki百科的定义：

Verlet积分法是经典力学（牛顿力学）中的一种最为普遍的积分方法，被广泛运用在分子运动模拟（Molecular Dynamics Simulation），行星运动以及织物变形模拟等领域。Verlet算法要解决的问题是，给定粒子t时刻的位置r和速度v，得到t+dt时刻的位置r(t+dt)和速度v(t+dt)。最简单的方法是前向计算（考虑当前和未来）的速度位移公式，也就是欧拉积分法，但精度不够，且不稳定。Verlet积分是一种综合过去、现在和未来的计算方法（居中计算），精度为O(4), 稳定度好，且计算复杂度不比欧拉积分高多少。

balabala。。。总结一下就是：

• verlet积分与欧拉积分最大的差别就是：欧拉积分考虑当前和下一时刻，verlet积分考虑上一时刻，当前，下一时刻

• verlet积分有着更好的精度，即4阶精度

• verlet积分和欧拉积分性能差不多

verlet积分其实很简单，只有两条公式，重点是要理解两条公式怎么来的，下面我尽量用简单易懂的表达来推出这两条verlet积分公式。

1

首先，上面说到的，verlet积分考虑上一时刻，当前，下一时刻，我们将上一时刻和下一时刻的位置关于时间的函数表示为x(t+Δt)和x(t-Δt)，其中Δt为一次步长。

由于多种力的存在，位置x关于时间t是一个复杂的函数，我们没办法用初等函数的形式将其表示出来。但是，我们有一个强有力的工具：泰勒公式。泰勒公式可以把函数近似成多项式和的形式，供我们研究函数的性质。关于泰勒公式的知识这里不多表述，高数不好的看这里。so，我们把 x(t+Δt)和x(t-Δt) 进行泰勒展开：





通常展开到4阶已经够了，后面的皮亚诺余项我们忽略掉。

2

将上面两条表达式进行相加，我们得到公式1，此公式用作计算下一时刻位置：

3

将上面两条表达式进行相减，我们得到公式2，此公式用作计算当前速度：

由公式1可以得出：要计算下一时刻的位置，只要知道当前位置，上一时刻位置和当前加速度即可，不需要知道当前时刻速度。
由公式2可以得出：只有在知道上一时刻和下一时刻的位置时，才有可能知道当前时刻的速度。也就是说，在当前时刻不能获得速度信息，必须要等到下一时刻的位置确定以后，才能返回来计算当前的速度。

就这些内容，其实还是很好理解的。

代码实现

这个贴代码其实没多大意义，因为都是套公式，没什么逻辑性的东西，但是为了撑撑排版还是贴一下吧。

/**
 * verlet积分更新物理
 * @param body 刚体对象
 * @param dt 步长
 */
function updatePhysicsByVerlet(body: any, dt: number) {
        // 上一时刻位置
    let prevPosition = body.prevPosition,
        // 当前时刻位置
        curPosition = body.curPosition,
        // 下一时刻位置
        nextPosition = null,
        // 下一时刻速度
        nextVelocity = null,
        // 当前时刻速度
        curVelocity = null,
        // 加速度
        acceleration = body.acceleration;

    // 套用verlet公式计算下一时刻位置和当前时刻速度
    nextPosition = 2*curPosition - prevPosition + acceleration*dt*dt;
    curVelocity = 0.5*(nextPosition - prevPosition)/dt;

    // 欧拉积分计算下一时刻速度：v' = v + a
    nextVelocity = curVelocity + acceleration;

    // 更新刚体的速度
    body.velocity = nextVelocity;

    // 更新刚体的位置
    body.prevPosition = curPosition;
    body.curPosition = body.position;
    body.position = nextPosition;
}

预告

下一篇issue我会说一下物理引擎是如何做到精确地每次隔dt时间更新一次的，这里面也大有学问。

[iceiceice123]

求速度那里两式相减的三次项是怎么消掉的啊？

[phenomLi]

求速度那里两式相减的三次项是怎么消掉的啊？

是笔误，余项应该是三阶，不是二阶

[vkensou]

学习了

[vkensou]

我觉得欧拉积分的写法应该是：
p+=vt+0.5at^2
v+=at
这样才是匀加速运动

[937318481]

你说的是解析的解，并不是欧拉方法