Introduction

Verlet integration is a nifty method for numerically integrating the equations of motion (typically linear, though you can use the same idea for rotational). It is a finite difference method that's popular with the Molecular Dynamics people. Actually, it comes in three flavors: the basic Position, the Leapfrog and the Velocity versions. We will be discussing the Position Verlet algorithm in this paper. It has the benefit of being quite stable, especially in the face of enforced boundary conditions (there is no explicit velocity term with which the position term can get out of sync). It is also very fast to compute (almost as fast as Euler integration), and under the right conditions it is 4^th order accurate (by comparison, the Euler method is only 1^st order accurate, and the second order Runge-Kutta method is only 2^nd order accurate [go figure]).

The disadvantages of the Verlet method are that it handles changing time steps badly, it is not a self-starter (it requires 2 steps to get going, so initial conditions are crucial), and it is unclear from the formulation how it handles changing accelerations. In this paper we will discuss all of these shortcomings, and see how to minimize their impact. The modified Verlet integrator is referred to as the Time-Corrected Verlet (TCV) and is shown below with its original counterpart. The computations used to generate the graphs for this paper are included in this Excel file.

Original Verlet:
x_i+1 = x_i + (x_i - x_i-1) ₊ a * dt * dt

Time-Corrected Verlet:
x_i+1 = x_i + (x_i - x_i-1) * (dt_i / dt_i-1) + a * dt_i * dt_i

To see why the TCV is an improvement, we need to see the math behind the original Verlet method.

Math

In this paper we will be talking exclusively about point masses, which are acted on by forces. Well, we all know about Newton's little equation: Force=d(Momentum)/dt. Momentum=mass*velocity and for non-relativistic speeds the mass is constant, removing our need for the chain rule, and yielding our familiar F=ma. So if all the forces acting on the point mass are summed (the vector F), then scaled by (1.0/m) we have the acceleration (a) of the point mass. Since we know how to go from force to acceleration, we will be starting from the acceleration term to keep the math less cluttered. In actual applications you will almost always be summing forces, then converting to accelerations.

Most of the graphs presented in this paper include the matching Euler simulation, just for reference. The Euler algorithm is extremely simple, and it will not be derived here. However, the equations are included below for your reference (order is important):

v = v + a * dt
x = x + v * dt

There is a more accurate version, but it is not strictly the Euler method, so it will not be used for this paper. However it does give 2^nd order accurate results, instead of merely 1^st order:

x = x + v * dt + 0.5 * a * dt * dt
v = v + a * dt

Please note that there is a faster way to derive the Verlet method's math* than what will be shown here, however it does not provide the insight needed to overcome the issues mentioned in the introduction. So we'll start from a few basic principles: a=dv/dt, and v=dx/dt. So for any given point mass, if we know the current position, velocity and acceleration (the acceleration must be constant over the given time step), we can compute exactly where it will be after the time step (dt) has elapsed.

(1)   x_i+1 = x_i + v_i * dt_i + 0.5 * a_i * dt_i * dt_i

or

(1a)   x_i+1 - x_i = v_i * dt_i + 0.5 * a_i * dt_i * dt_i

Now, if we wanted that equation formulated without the velocity term, we could replace it with some other known state variable, such as the position x.  But since we don't know x_i+1 yet, we shift the whole equation back a step:

(2)   x_i - x_i-1 = v_i-1 * dt_i-1 + 0.5 * a_i-1 * dt_i-1 * dt_i-1

and we can use simple integration to see that:

(3)   v_i = v_i-1 + a_i-1 * dt_i-1

or

(3a)   v_i-1 = v_i - a_i-1 * dt_i-1

and we substitute that into equation (2):

(4)   x_i - x_i-1 = (v_i - a_i-1 * dt_i-1) * dt_i-1 + 0.5 * a_i-1 * dt_i-1 * dt_i-1

or

(4a)   x_i - x_i-1 = v_i * dt_i-1 - 0.5 * a_i-1 * dt_i-1 * dt_i-1

For this next step we need to make a big assumption, the importance of which will be seen later: If we assume that neither the acceleration nor the time step vary between steps (i.e. that a_i-1 = a_i = a and that dt_i-1 = dt_i = dt) then we note that:

(5)   x_i - x_{i-1 +} a * dt * dt = v_i * dt - 0.5 * a * dt * dt + a * dt * dt

or

(5a)   x_i - x_{i-1 +} a * dt * dt = v_i * dt + 0.5 * a * dt * dt

You'll notice that the right hand side of equation (5a) is exactly the last half of equation (1), so we can work the modified equation (5a) back into equation (1):

(6)   x_i+1 = x_i + (x_i - x_i-1) ₊ a * dt * dt

and there you have the traditional Verlet Position integration method.

* Hint: remember the central difference 2^nd derivative approximation?
d²x / dt² = a_i²= (x_i+1 - 2*x_i + x_i+1) / dt²

Fundamental Problems