Verlet Method

One of the most common drift-free higher-order algorithms is commonly attributed to Verlet [L. Verlet, Computer experiments on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules, Physical Review 159, 98 (1967); in this paper he discusses the method and its application to molecular dynamics simulations].We write the Taylor series expansion for $x_{n-1}$ in a form similar to Eq. (3)

\begin{displaymath} x_{n-1}=x_n - v_n \Delta t + \frac{1}{2} a_n (\Delta t)^2. \end{displaymath} (19)

If we add the forward and reverse forms, Eq. (3) and Eq. (19) respectively, we obtain
\begin{displaymath} x_{n+1} + x_{n-1}=2x_n + a_n (\Delta t)^2 + \mathcal{O}[(\Delta t)^4], \end{displaymath} (20)

or
\begin{displaymath} x_{n+1}=2x_n - x_{n-1} + a_n (\Delta t)^2. \end{displaymath} (21)

Similarly, the subtraction of the Taylor series for $x_{n+1}$ and $x_{n-1}$ yields
\begin{displaymath} v_n = \frac{x_{n+1} - x_{n-1}}{2 \Delta t} \ {\rm (original Verlet)}. \end{displaymath} (22)

Thus, the global error associated with the Verlet algorithm is third order for the position and second-order for the velocity. However, the velocity plays no part in the integration of the equations of motions. In the numerical analysis literature, the Verlet method is also knows as the ``explicit central difference method''.

Because the Verlet algorithm is not self-starting, another algorithm must be used to obtain the first few terms. An additional problem is that the new velocity Eq. (22) found by computing the difference between two quantities of the same order of magnitude. When using computers which always operate with finite numerical precision, such an operation results in a loss of numerical precision and may give rise to substantial roundoff error.

A mathematically equivalent version of the original Verlet algorithm is given by

$\displaystyle x_{n+1}$ $\textstyle =$ $\displaystyle x_n + v_n \Delta t + \frac{1}{2} a_n (\Delta t)^2,$ (23)
$\displaystyle v_{n+1}$ $\textstyle =$ $\displaystyle v_n + \frac{1}{2}(a_{n+1} + a_n) \Delta t, \ {\rm (velocity Verlet)}.$ (24)

We see that these equations, known as the velocity form of the Verlet algorithm, is self-starting and minimizes roundoff errors. We can derive Eqs. (23), (24) from Eqs. (21), (22) by the following considerations. We first solve Eq. (22) for $x_{n-1}$ and write $x_{n-1}=x_{n+1} - 2 v_n (\Delta t)$. If we substitute this expression for $x_{n-1}$ into Eq. (21) and solve for $x_{n+1}$, we find the form Eq. (23). Then we use Eq. (22) to write $v_{n+1}$ as
\begin{displaymath} v_{n+1}=\frac{x_{n+2} - x_n}{2 \Delta t}, \end{displaymath} (25)

and use Eq. (21) to obtain $x_{n+2}=2x_{n+1}- x_n + a_{n+1} (\Delta t)^2$. If we substitute this form for $x_{n+2}$ into Eq. (25), we obtain
\begin{displaymath} v_{n+1} = \frac{x_{n+1} - x_n}{\Delta t} + \frac{1}{2} a_{n+1} \Delta t. \end{displaymath} (26)

Finally, we use Eq. (23) for $x_{n+1}$ to eliminate $x_{n+1} - x_n$ from Eq. (26); after some algebra we obtain the desired result Eq. (24).

Another useful algorithm that avoids the roundoff error of the original Verlet algorithm is due to Beeman and Schofield. We write the Beeman algorithm in the form

$\displaystyle x_{n+1}$ $\textstyle =$ $\displaystyle x_n + v_n \Delta t + \frac{1}{6} (4 a_n - a_{n-1})(\Delta t)^2,$ (27)
$\displaystyle v_{n+1}$ $\textstyle =$ $\displaystyle v_n + \frac{1}{6} (2a_{n+1} + 5 a_n - a_{n-1})\Delta t.$ (28)

Note that this algorithm does not calculate particle trajectories more accurately than the Verlet algorithm. Its advantage is that in general does a better job of maintaining energy conservation. However, the Beeman algorithm is not self-starting.