Verlet Integration

Posted on March 02, 2023

Derivation

Störmer-Verlet

For systems governed by conservative forces, Newton's laws specify a second order differential equation.

$$ \begin{equation} f(x) = m \ddot{x} \end{equation} $$

To discretize this system, we can use the central difference for the second derivative.

$$\begin{equation} \ddot{x}_n \approx \frac{1}{\Delta t^2} \left( x_{n - 1} - 2 x_n + x_{n + 1} \right) \end{equation}$$

Geometrically, this estimation is the second derivative of the parabola that passes through $(-\Delta t, x_{n - 1})$, $(0, x_n)$, and $(\Delta t, x_{n + 1})$.

If we solve for $x_{n+1}$, we find that we can use this method to step forward in time, as long as we start with two samples $x_0$ and $x_1$.

$$\begin{equation} x_{n+1} = 2 x_n - x_{n - 1} + \Delta t^2 \frac{f(x_n)}{m} \end{equation}$$

Geometrically, this corresponds to constructing the parabola that has second derivative $f(x_n)/m$ and passes through $(-\Delta t, x_{n - 1})$ and $(0, x_n)$, then extrapolating it out to $t = \Delta t$ to get our estimate for $x_{n+1}$.

But for physics problems, we usually start with a position and a velocity, rather than two consecutive position samples. So how can we get this integrator started? We can fit our parabola a third way: construct the parabola that passes through $(0, x_n)$, has slope $v_0$ at $t = 0$, and has second derivative $f(x_n)/m$. Then extrapolate it out to $t = \Delta t$.

$$\begin{equation} x_{n+1} = x_n + \Delta t \: v_n + \frac{\Delta t^2}{2} \frac{f(x_n)}{m} \end{equation}$$

(You can also think of this as the position of a particle at time $\Delta t$ if it starts at position $x_0$ with velocity $v_0$ and experiences a constant acceleration $f(x_n)/m$. It's the same, since a particle under constant acceleration follows a parabolic trajectory.)

So now, starting with $x_0$ and $v_0$, equation 4 gives us a way to compute $x_1$. Then we can use $x_0$ and $x_1$ to compute $x_2$ via equation 3, and so on.

Velocity Verlet

Often we'll want to know the velocity at the other time steps as well. To compute it, we can solve for the slope of yet another parabola. (Specifically, the slope at $t = \Delta t$ for the parabola passing through $(0, x_n)$ and $(\Delta t, x_{n+1})$ with second derivative $f(x_{n+1})/m$.)

$$\begin{equation} v_{n+1} = \frac{x_{n+1} - x_n}{\Delta t} + \frac{\Delta t}{2} \frac{f(x_{n+1})}{m} \end{equation}$$

If we're going to do this for every time step, we can use equations 4 and 5 to integrate forward, without ever using equation 3 directly. It's a different way to fit the same parabolas, so the two methods are equivalent.

But, if we just computed $x_{n+1}$ via equation 4, we can substitute it in to save ourselves some work.

$$\begin{equation} \begin{aligned} v_{n+1} &= \frac{\Delta t \: v_n + \frac{\Delta t^2}{2} \frac{f(x_n)}{m}}{\Delta t} + \frac{\Delta t}{2} \frac{f(x_{n+1})}{m} \\ &= v_n + \frac{\Delta t}{2} \frac{f(x_n)}{m} + \frac{\Delta t}{2} \frac{f(x_{n+1})}{m} \\ \end{aligned} \end{equation}$$

Note that the first two terms in this equation appear in equation 4 (with a factor of $\Delta t$). Let's factor that out as $v_{n + 1/2}$.

$$\begin{equation} \begin{aligned} v_{n+1/2} &= v_n + \frac{\Delta t}{2} \frac{f(x_n)}{m} \\ x_{n+1} &= x_n + \Delta t \: v_{n + 1/2} \\ v_{n+1} &= v_{n + 1/2} + \frac{\Delta t}{2} \frac{f(x_{n+1})}{m} \\ \end{aligned} \end{equation}$$

Cycling through these equations for each time step is equivlant to Störmer-Verlet, but it lets us inspect the velocity as well as the position at each time step. It's also more numerically stable. This method of integration is the default in molecular dynamics and other classical particle based simulations (i.e. those governed by Newton's laws rather than domain-specific partial differential equations).

As a Variation of Semi-Implicit Euler

Recall that the Euler method for a first order ODE $\dot{x} = f(x)$ simply takes steps along the tangents: $x_{n+1} = x_n + \Delta x \: f(x_n)$. We can break our second order problem into a coupled first order problem.

$$\begin{equation} \begin{aligned} \dot{v} &= \frac{f(x)}{m} \\ \dot{x} &= v \end{aligned} \end{equation}$$

Which gives the follwing Euler update rules.

$$\begin{equation} \begin{aligned} v_{n+1} &= v_n + \Delta t \frac{f(x)}{m} \\ x_{n+1} &= x_n + \Delta t \: v_n \end{aligned} \end{equation}$$

The semi-implicit Euler method just uses the latest velocity value when updating the position.

$$\begin{equation} \begin{aligned} v_{n+1} &= v_n + \Delta t \frac{f(x)}{m} \\ x_{n+1} &= x_n + \Delta t \: v_{n+1} \end{aligned} \end{equation}$$

We can also do it the other way around, and use the newest position value when updating the velocity.

$$\begin{equation} \begin{aligned} x_{n+1} &= x_n + \Delta t \: v_n \\ v_{n+1} &= v_n + \Delta t \frac{f(x_{n+1})}{m} \end{aligned} \end{equation}$$

If we use half the timestep, and apply the update rules from equations 10 and 11 in sequence, we end up with equation 7. So velocity Verlet can be seen as a variation of semi-implicit Euler. Both integrators are symplectic (more on this later), but semi-implicit Euler is order 1, while Verlet is order 2. The additional computational cost is negligible (in practice, all our time will be spent computing $f(x_n)$), so Verlet integration is usually preferable.

Properties

Order

Velocity Verlet is locally accurate to second order, as can be seen by plugging $v_{n + 1/2}$ into the expression for $x_{n+1}$ in equation 7 (you should get the exact Taylor expansion up to and including $O(h^2)$). As a single step method, its global error is then $O(h^2)$ instead of $O(h^3)$.

Locally, position Verlet is a third order method (meaning the error is $O(h^4)$).

$$\begin{aligned} x(t + h) &= x(t) + h x'(t) + \frac{h^2}{2} x''(t) + \frac{h^3}{3!} x'''(t) + O(h^4) \\ x(t - h) &= x(t) - h x'(t) + \frac{h^2}{2} x''(t) - \frac{h^3}{3!} x'''(t) + O(h^4) \\ x(t + h) + x(t - h) &= 2 x(t) + h^2 x''(t) + O(h^4) \\ x(t + h) &= 2 x(t) - x(t - h) + h^2 x''(t) + O(h^4) \\ \end{aligned}$$

Globally, however, the error is still $O(h^2)$. This happens because position Verlet is a multistep method.

For comparison, Euler methods (both forward and semi-implicit) both have local errors of order $O(h^2)$ and global errors of order $O(h)$.

Symplectic Properties

Both position and velocity Verlet are time reversible. See this notebook.

It nearly preserves energy. Interestingly, it doesn't exactly preserve energy. But it should be easy to compare with other methods numerically.

In technical terms, the Verlet method is a symplectic integrator. This means that it preserves the exterior product $dq \wedge dp$. In words, it preserves areas in phase space.