The time-corrected Verlet numerical integration fo

问题:

There is a commonly used verlet-integration formula on the web by Johnathan Dummer, called Time-Corrected Verlet. However I've read several forum posts, that people get weird or unexpected results with it in certain conditions.

Formula by Johnathan Dummer:

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

There is also a stackoverflow answer, which states that Dummer's time-corrected formula is broken and the poster presents his own derivation as the correct one.

Suggested correct formula by a stackoverflow answer

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

Well, is Dummer's formula really broken? If yes, is the derivation of the poster better?

PS: It is also weird that Dummer uses the verlet integration formula x1 = x - x0 + a * dt^2 on his website instead of the correct x1 = 2x - x0 + a * dt^2.

回答1:

The wikipedia page Verlet integration - Non-constant time differences presents the two formula, without referenced. I've not checked the derivation myself but the reasoning for the second improved formula looks sound.

I've downloaded Dummer's spreadsheet and modified one of the formula to use the correction. The results are much better.

The exact results are in yellow, we see that just using the normal Verlet algorithm with fluctuating frame-rate is bad. Dummer's time-correct varient in red is pretty good, but a little off. The dark green version with the improved correction is much better.

For projectiles under gravity which has a quadratic solution you may find that the improved version is exact. When the degree gets a bit higher it will vary from the true path and it might be worth testing to see if we still get a better approximation.

Doing the same calculation for a sin curve shows the improved method is considerably better. Here Time-correct Verlet is drifting quite a bit. The improved version is better, only a little bit off the exact answer.

---

For the PS. Note that if you set dt=dt0 in the TCV formula

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

you get

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

the original Verlet formula.

回答2:

I decided to quit being lazy and show some kind of derivation of how the original Verlet method looks with a variable step size. Because it seems like this faulty adaptation by Dummer is more pervasive than I thought, which is saddening. I also noticed that, as the above answer points out, the correct version is now on wikipedia alongside the Dummer one, though it was added after my "suggested correct answer".

When I look at Verlet method, I see that it looks a lot like leapfrog, velocity Verlet, implicit Euler etc., which look like second-order versions of modified midpoint, and some of them may be identical. In each of these, to some degree they have a leapfrog idea where integration of acceleration (into velocity) and integration of constant velocity (into position) are each staggered so that they overlap by half. This brings things like time-reversibility and stability, which are more important for the 'realism' of a simulation than accuracy is. And the 'realism', the believability, is more important for video games. We don't care if something moves to a slightly different position than it's exact mass would have truly caused, so long as it looks and feels realistic. We're not calculating where to point our high-powered satellite telescopes to look at features on distant objects, or future celestial events. Here, stability and efficiency takes priority here over mathematical accuracy. So, it seems like leapfrog method is appropriate. When you adapt leapfrog for variable time step, it loses some of this advantage, and it loses some of it's appeal for game physics. Stormer-Verlet is like leapfrog, except it uses the average velocity of the previous step instead of a separately maintained velocity. You can adapt this Stormer-Verlet in the same way as leapfrog. To integrate velocity forward with a fixed acceleration, you use half the length of the previous step and half the length of the next step, because they are staggered. If the steps were fixed like true leapfrog, they would be the same length and so the two half lengths sum to one. I use h for step size, a/v/p for acceleration/velocity/position, and hl/pl for 'last' as in previous step. These aren't really equations, more like assignment operations.

Original leapfrog:

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

With variable time step:

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

Factor a/2:

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

Use previous position (p - pl)/hl for initial velocity:

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

Substitute, we don't need v:

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

Distribute h:

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

The result is not as simple or fast as the original Stormer form of Verlet, 2p - pl + a*h^2. I hope this makes some sense. You would omit the last step in actual code, no need to multiply h twice.

回答3:

The true derivation is based on Taylor formulas

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)

and now eliminate x'(t) from these two formulas to get a Verlet-like formula

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

which makes a propagation formula

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)

so that indeed the corrected formula is the correct one.

---

Note that if you use velocity Verlet steps

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

then each step is indepentently symplectic, so that the change of step size between steps is absolutely unproblematic. then the