The time-corrected Verlet numerical integration formula

Question

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.

Feather

Answer 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.

Interstate

Answer 2

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.

---

Notice that the main advantages of Verlet and other similar symplectic schemes over Runge-Kutta methods etc. crucially depends on using a fixed step. In detail, the modified energy function (and other more than quadratic constants of motion) that is largely constant under the numerical method is a modification of the exact energy where the difference scales with the step size. So when changing the step size a different modification gives a constant energy. Frequent changes of the step size allow thus, possibly, arbitrary changes of the energy level.

Costive

This is a more exact derivation, since Taylor series approximation is explained on the internet. I was just giving something simple. – Crutch

@LutzLehmann "Note that if you use velocity Verlet steps then each step is indepentently symplectic, so that the change of step size between steps is absolutely unproblematic." - Then what is the advantage of the time-corrected verlet method? Improved accuracy? Also, as I know velocity Verlet method assumes that the acceleration only depends on the position and does not depend on velocity. – Feather

@Feather : That the acceleration is only a function of x, indeed derives from the gradient of a potential function, is still necessary in the velocity Verlet variant to get a symplectic method. The method in question just looks more simple, I'm not sure if that still holds after the second look. Also, I'm not sure if it is still symplectic, I highly doubt it (the error in that should only depend on the variance of the step sizes, so still be small). – Costive

@LutzLehmann So when the acceleration depends on the velocity (like in case of friction), then the system won't be symplectic no matter you use a symplectic integration method. And this is true for all symplectic methods. Also, a system conserveres energy iff it is symmetric (invariant to time reversal) and symplectic. Are these assumptions right? I'm also interested whether the time-corrected Verlet remains symplectic, I think I create an another question about it. – Feather

@Feather : A system is conservative if it is Hamiltonian, x'=H_p, p'=-H_x. Equations like y''+y'^2+y=0 are also symmetric under time reversal, but what is the energy? This new question would belong to scicomp.SE or math.SE, perhaps physics.SE, but not here. – Costive

@LutzLehmann Thanks, it is more clear now: a symplectic integrator does conserve energy iff the underlying system does. – Feather

@LutzLehmann I also found a reference on Wikipedia (en.wikipedia.org/wiki/Leapfrog_integration#Algorithm) for this "kick-drift-kick" form of velocity Verlet you recommend. However it is not clear to me how is it different from the standard velocity Verlet, since mathematically the two expressions are equivalent, it is only rearranged for an explicit half-step velocity term. I guess it changes something numerically when the result is computed in floating-point precision where roundoff error applies. Is my assumption correct? – Feather

For fixed-step Verlet you can translate all the variants into each other, each defines the same position sequence (up to floating point error accumulation). If you make the time step variable, you have to somehow keep this pattern. One should at least consistency, that is, the new values should be at least second order accurate extrapolations of the current values. – Costive

@LutzLehmann I made a separate question about the time step adaptivity of velocity Verlet on scicomp: scicomp.stackexchange.com/q/34491/23946 – Feather

Answer 3

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.

Crutch

You can not start with Leapfrog, since the times for v and x differ. – Costive

But you can start with the formula for a leapfrog step, which is all I did. – Crutch

As long as you consider that the values v=v[n+1/2] are for the times (t[n]+t[n+1])/2. You implicitly do that by the rather unmotivated replacement of h by (h+h1)/2. So in fact you split Leapfrog at the boundaries of velocity Verlet and jump back-and-forth between these two interpretations. – Costive

The motivation is explicitly stated as "variable time step", which is indeed exactly what it does (edit: it does so by splitting the two time steps into explicit pieces). And there's no "jump back-and-forth". The single jump to Stormer Verlet -like is given as "use previous position (to approximate) initial velocity". I'm certain my original adaptation had reasoning that didn't resemble this post, but it was a long time ago, and the point was to give something simple, something appropriate for people who failed to correctly derive it themselves... like Dummer. – Crutch