Numerical Integration: Verlet integration

Position Verlet

 To a close approximation (actually performing a Taylor Series Expansion about  x(t±Δt) ), that might look like this:

\[x(t+\Delta t) = x(t) + v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 + \frac{1}{6}b(t) \Delta t^3 + \mathcal{O}(\Delta t^4)\]

This means that if we need to find the next  x, we need the current  x ,  v ,  a , etc. However, because few people calculate the jerk term, our error is typically  \(\mathcal{O}(\Delta t^3)\) . That said, we can calculate  x  with less knowledge and higher accuracy if we play a trick! Let's say we want to calculate  x  of the previous timestep. Again, to a close approximation, that might look like this:

\[x(t-\Delta t) = x(t) - v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 - \frac{1}{6}b(t) \Delta t^3 + \mathcal{O}(\Delta t^4)\]

Now, we have two equations to solve for two different timesteps in x, one of which we already have. If we add the two equations together and solve for x(t+Δt), we find that

\[x(t+ \Delta t) = 2x(t) - x(t-\Delta t) + a(t)\Delta t^2 + \mathcal{O}(\Delta t^4)\]

So, this means we can find our next  x simply by knowing our current  x , the  x  before that, and the acceleration! No velocity necessary! In addition, this drops the error to \(\mathcal{O}(\Delta t^4)\) , which is great! Here is what it looks like in code:

  
  XXXXXXXXXXXXXXXXX

Now, obviously this poses a problem; what if we want to calculate a term that requires velocity, like the kinetic energy, \(\frac{1}{2}mv^2\). In this case, we certainly cannot get rid of the velocity! Well, we can find the velocity to  O(Δt2)  accuracy by using the Stormer-Verlet method. We have the equations for  x(t+Δt)  and  x(t−Δt) above, so let's start there. If we subtract the latter from the former, we get the following:

\[x(t+\Delta t) - x(t - \Delta t) = 2v(t)\Delta t + \frac{1}{3}b(t)\Delta t^3\]

When we solve for  v(t), we get

\[\begin{align} v(t) &= \frac{x(t+\Delta t) - x(t-\Delta t)}{2\Delta t} + \frac{b(t) \Delta t^3}{3 \Delta t} \\  v(t) &= \frac{x(t+\Delta t) - x(t-\Delta t)}{2\Delta t} + \mathcal{O}(\Delta t^2). \end{align}\]

Note that the 2 in the denominator makes sense because we are going over 2 timesteps. It's essentially solving  \(v=\frac{Δx}{Δt}\). In addition, we can calculate the velocity of the next timestep like so

\[v(t+\Delta t) = \frac{x(t+\Delta t) - x(t)}{\Delta t} + \mathcal{O}(\Delta t)\]

Velocity Verlet

Now, let's say we actually need the velocity to calculate out next timestep. Well, in this case, we simply cannot use the above approximation and instead need to use the Velocity Verlet algorithm.

In some ways, this algorithm is even simpler than above. We can calculate everything like

\[\begin{align} x(t+\Delta t) &=x(t) + v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 \\ a(t+\Delta t) &= f(x(t+\Delta t)) \\ v(t+\Delta t) &= v(t) + \frac{1}{2}(a(t) + a(t+\Delta t))\Delta t \end{align}\]

which is literally the kinematic equation above, solving for  x ,  v , and  a  every timestep. You can also split up the equations like so

\[\begin{align} v(t+\frac{1}{2}\Delta t) &= v(t) + \frac{1}{2}a(t)\Delta t \\ x(t+\Delta t) &=x(t) + v(t+\frac{1}{2}\Delta t)\Delta t \\ a(t+\Delta t) &= f(x(t+\Delta t)) \\ v(t+\Delta t) &= v(t+\frac{1}{2}\Delta t) + \frac{1}{2}a(t+\Delta t)\Delta t \end{align}\]

Here is the velocity Verlet method in code:

  XXXXXXXXXXXXXXXXXX

Even though this method is more widely used than the simple Verlet method mentioned above, it unfortunately has an error term of \(\mathcal{O}(\Delta t^2)\), which is two orders of magnitude worse. 

Verlet Leapfrog Algorithm

Takes the form

\[\begin{aligned}x_{i+1}&=x_{i}+v_{i}\Delta t+{\frac {1}{2}}a_{i}\Delta t^{2},\\v_{i+1}&=v_{i}+{\frac {1}{2}}(a_{i}+a_{i+1})\Delta t.\end{aligned}\]

Useful References

https://en.wikipedia.org/wiki/Verlet_integration

https://julesjacobs.com/2019/03/15/leapfrog-verlet.html

https://pikuma.com/blog/verlet-integration-2d-cloth-physics-simulation

https://www.algorithm-archive.org/contents/verlet_integration/verlet_integration.html