时间校正的 Verlet 数值积分公式

The time-corrected Verlet numerical integration formula

网络上有一个常用的 verlet-integration 公式,由 Johnathan Dummer 编写,称为时间校正 Verlet。但是我读过几篇论坛帖子,人们在某些情况下会得到奇怪或意想不到的结果。

Formula by Johnathan Dummer:

x1 = x + (x – x0) * dt / dt0 + a * dt^2

还有一个Whosebug answer,其中指出 Dummer 的时间修正公式被破坏,张贴者提出他自己的推导是正确的。

Suggested correct formula by a Whosebug answer

x1 = x + (x – x0) * dt / dt0 + a * dt * (dt + dt0) / 2

嗯,Dummer的公式真的坏了吗?是的话楼主的推导是不是比较好?

PS:同样奇怪的是,Dummer 在他的网站上使用了 verlet 积分公式 x1 = x - x0 + a * dt^2 而不是正确的 x1 = 2x - x0 + a * dt^2.

维基百科页面Verlet integration - Non-constant time differences介绍了这两个公式,没有引用。我自己没有检查推导,但第二个改进公式的推理看起来很合理。

我已经下载了 Dummer 的电子表格并修改了其中一个公式以使用更正。结果好多了。

确切的结果是黄色的,我们看到只使用帧率波动的普通 Verlet 算法是不好的。 Dummer 的红色时间校正变体非常好,但有点偏差。改进校正后的深绿色版本要好得多。

对于具有二次解的重力弹丸,您可能会发现改进后的版本是精确的。当度数变得更高时,它会与真实路径不同,可能值得测试一下我们是否还能得到更好的近似值。

对 sin 曲线进行相同的计算表明改进后的方法要好得多。在这里,时间正确的 Verlet 漂移了很多。改进版更好,只差一点点准确答案。


为PS。请注意,如果您在 TCV 公式中设置 dt=dt0

x1 = x + (x – x0) * dt / dt0 + a * dt^2

你得到

x1 = x + x – x0 + a * dt^2
   = 2 x – x0 + a * dt^2

原始的 Verlet 公式。

我决定不再偷懒了,并展示了对原始 Verlet 方法在可变步长下的外观的某种推导。因为 Dummer 的这种错误改编似乎比我想象的更普遍,这令人难过。我还注意到,正如上面的答案所指出的,正确的版本现在与 Dummer 版本一起出现在维基百科上,尽管它是在我的 "suggested correct answer".

之后添加的

当我查看 Verlet 方法时,我发现它看起来很像 leapfrog、velocity Verlet、implicit Euler 等,它们看起来像修改中点的 second-order 版本,其中一些可能是相同的.在其中的每一个中,他们在某种程度上都有一个跨越式的想法,其中加速度的积分(到速度)和恒定速度的积分(到位置)每个都是交错的,以便它们重叠一半。这带来了 time-reversibility 和稳定性之类的东西,这对于模拟的 'realism' 比准确性更重要。而'realism',可信度,对于电子游戏来说更为重要。我们不在乎某物移动到的位置是否与它的精确质量真正造成的位置略有不同,只要它看起来并且感觉真实.我们没有计算将 high-powered 卫星望远镜指向何处以观察遥远 objects 或未来天体事件的特征。在这里,稳定性和效率优先于数学准确性。所以,越级法似乎是合适的。当您为可变时间步调整 leapfrog 时,它会失去一些优势,并且会失去一些对游戏物理的吸引力。 Stormer-Verlet 类似于 leapfrog,只不过它使用的是前一步的平均速度,而不是单独维护的速度。您可以采用与 leapfrog 相同的方式调整此 Stormer-Verlet。要将速度向前与固定加速度相结合,您可以使用前一步长度的一半和下一步长度的一半,因为它们是交错的。如果台阶像真正的蛙跳一样固定,它们的长度将相同,因此两个半长之和为一。我使用 h 表示步长,a/v/p 表示 acceleration/velocity/position,hl/pl 表示 'last',就像上一步一样。这些并不是真正的方程式,更像是赋值运算。

原创越级:

v = v + a*h
p = p + v*h

可变时间步长:

v = v + a*hl/2 + a*h/2
p = p + v*h

因素a/2:

v = v + a*(hl + h)/2
p = p + v*h

使用之前的位置(p - pl)/hl作为初始速度:

v = (p - pl)/hl + a*(hl + h)/2
p = p + v*h

代入,我们不需要v:

p = p + ( (p - pl)/hl + a*(hl + h)/2)*h

分发h:

p = p + (p - pl)*h/hl + a*h*(h + hl)/2

结果不像 Verlet 的原始 Stormer 形式那样简单或快速,2p - pl + a*h^2。我希望这是有道理的。您将省略实际代码中的最后一步,无需乘以 h 两次。

真正的推导是基于泰勒公式

x(t-h0) = x(t) - x'(t)*h0 + 0.5*x''(t)*h0^2 + O(h0^3)
x(t+h1) = x(t) + x'(t)*h1 + 0.5*x''(t)*h1^2 + O(h1^3)

现在从这两个公式中消除 x'(t) 以获得类似 Verlet 的公式

h0*x(t+h1) + h1*x(t-h0) = (h0+h1)*x(t) + 0.5*a(t)*h0*h1*(h0+h1) +O(h^3)

这使得传播公式

x(t+h1) = (1+h1/h0)*x(t) - h1/h0*x(t-h0) + 0.5*a(t)*h1*(h0+h1)
        = x(t) + (x(t)-x(t-h0))*h1/h0 + 0.5*a(t)*h1*(h0+h1)

所以修正后的公式确实是正确的。


请注意,如果您使用速度 Verlet 步骤

Verlet(dt) {
    v += a * 0.5*dt
    x += v*dt
    a = acceleration(x)
    v += a * 0.5*dt
}

那么每一步都是独立辛的,这样步与步之间步长的变化是绝对没有问题的。 然后