Verlet integration

Newton's equation of motion for conservative physical systems is

${\displaystyle M{\ddot {\vec {x}}}(t)=F{\big (}{\vec {x}}(t){\big )}=-\nabla V{\big (}{\vec {x}}(t){\big )},}$

or individually

${\displaystyle m_{k}{\ddot {\vec {x}}}_{k}(t)=F_{k}{\big (}{\vec {x}}(t){\big )}=-\nabla _{{\vec {x}}_{k}}V{\big (}{\vec {x}}(t){\big )},}$

where

t is the time,

${\displaystyle {\vec {x}}(t)={\big (}{\vec {x}}_{1}(t),\ldots ,{\vec {x}}_{N}(t){\big )}}$ is the ensemble of the position vector of N objects,

V is the scalar potential function,

F is the negative gradient of the potential, giving the ensemble of forces on the particles,

M is the mass matrix, typically diagonal with blocks with mass ${\displaystyle m_{k}}$ for every particle.

This equation, for various choices of the potential function V, can be used to describe the evolution of diverse physical systems, from the motion of interacting molecules to the orbit of the planets.

After a transformation to bring the mass to the right side and forgetting the structure of multiple particles, the equation may be simplified to

${\displaystyle {\ddot {\vec {x}}}(t)=A{\big (}{\vec {x}}(t){\big )}}$

with some suitable vector-valued function A representing the position-dependent acceleration. Typically, an initial position ${\displaystyle {\vec {x}}(0)={\vec {x}}_{0}}$ and an initial velocity ${\displaystyle {\vec {v}}(0)={\dot {\vec {x}}}(0)={\vec {v}}_{0}}$ are also given.

To discretize and numerically solve this initial value problem, a time step ${\displaystyle \Delta t>0}$ is chosen, and the sampling-point sequence ${\displaystyle t_{n}=n\,\Delta t}$ considered. The task is to construct a sequence of points ${\displaystyle {\vec {x}}_{n}}$ that closely follow the points ${\displaystyle {\vec {x}}(t_{n})}$ on the trajectory of the exact solution.

Where Euler's method uses the forward difference approximation to the first derivative in differential equations of order one, Verlet integration can be seen as using the central difference approximation to the second derivative:

${\displaystyle {\frac {\Delta ^{2}{\vec {x}}_{n}}{\Delta t^{2}}}={\frac {{\frac {{\vec {x}}_{n+1}-{\vec {x}}_{n}}{\Delta t}}-{\frac {{\vec {x}}_{n}-{\vec {x}}_{n-1}}{\Delta t}}}{\Delta t}}={\frac {{\vec {x}}_{n+1}-2{\vec {x}}_{n}+{\vec {x}}_{n-1}}{\Delta t^{2}}}={\vec {a}}_{n}={\vec {A}}({\vec {x}}_{n}).}$

Verlet integration in the form used as the Störmer method[3] uses this equation to obtain the next position vector from the previous two without using the velocity as

${\displaystyle {\vec {x}}_{n+1}=2{\vec {x}}_{n}-{\vec {x}}_{n-1}+{\vec {a}}_{n}\,\Delta t^{2},\qquad {\vec {a}}_{n}={\vec {A}}({\vec {x}}_{n}).}$

The time symmetry inherent in the method reduces the level of local errors introduced into the integration by the discretization by removing all odd-degree terms, here the terms in ${\displaystyle \Delta t}$ of degree three. The local error is quantified by inserting the exact values ${\displaystyle {\vec {x}}(t_{n-1}),{\vec {x}}(t_{n}),{\vec {x}}(t_{n+1})}$ into the iteration and computing the Taylor expansions at time ${\displaystyle t=t_{n}}$ of the position vector ${\displaystyle {\vec {x}}(t\pm \Delta t)}$ in different time directions:

${\displaystyle {\begin{aligned}{\vec {x}}(t+\Delta t)&={\vec {x}}(t)+{\vec {v}}(t)\Delta t+{\frac {{\vec {a}}(t)\Delta t^{2}}{2}}+{\frac {{\vec {b}}(t)\Delta t^{3}}{6}}+{\mathcal {O}}(\Delta t^{4})\\{\vec {x}}(t-\Delta t)&={\vec {x}}(t)-{\vec {v}}(t)\Delta t+{\frac {{\vec {a}}(t)\Delta t^{2}}{2}}-{\frac {{\vec {b}}(t)\Delta t^{3}}{6}}+{\mathcal {O}}(\Delta t^{4}),\end{aligned}}}$

where ${\displaystyle {\vec {x}}}$ is the position, ${\displaystyle {\vec {v}}={\dot {\vec {x}}}}$ the velocity, ${\displaystyle {\vec {a}}={\ddot {\vec {x}}}}$ the acceleration, and ${\displaystyle {\vec {b}}}$ the jerk (third derivative of the position with respect to the time).

Adding these two expansions gives

${\displaystyle {\vec {x}}(t+\Delta t)=2{\vec {x}}(t)-{\vec {x}}(t-\Delta t)+{\vec {a}}(t)\Delta t^{2}+{\mathcal {O}}(\Delta t^{4}).}$

We can see that the first- and third-order terms from the Taylor expansion cancel out, thus making the Verlet integrator an order more accurate than integration by simple Taylor expansion alone.

Caution should be applied to the fact that the acceleration here is computed from the exact solution, ${\displaystyle {\vec {a}}(t)=A{\big (}{\vec {x}}(t){\big )}}$, whereas in the iteration it is computed at the central iteration point, ${\displaystyle {\vec {a}}_{n}=A({\vec {x}}_{n})}$. In computing the global error, that is the distance between exact solution and approximation sequence, those two terms do not cancel exactly, influencing the order of the global error.

To gain insight into the relation of local and global errors, it is helpful to examine simple examples where the exact solution, as well as the approximate solution, can be expressed in explicit formulas. The standard example for this task is the exponential function.

Consider the linear differential equation ${\displaystyle {\ddot {x}}(t)=w^{2}x(t)}$ with a constant w. Its exact basis solutions are ${\displaystyle e^{wt}}$ and ${\displaystyle e^{-wt}}$.

The Störmer method applied to this differential equation leads to a linear recurrence relation

${\displaystyle x_{n+1}-2x_{n}+x_{n-1}=h^{2}w^{2}x_{n},}$

or

${\displaystyle x_{n+1}-2\left(1+{\tfrac {1}{2}}(wh)^{2}\right)x_{n}+x_{n-1}=0.}$

It can be solved by finding the roots of its characteristic polynomial ${\displaystyle q^{2}-2\left(1+{\tfrac {1}{2}}(wh)^{2}\right)q+1=0}$. These are

${\displaystyle q_{\pm }=1+{\tfrac {1}{2}}(wh)^{2}\pm wh{\sqrt {1+{\tfrac {1}{4}}(wh)^{2}}}.}$

The basis solutions of the linear recurrence are ${\displaystyle x_{n}=q_{+}^{n}}$ and ${\displaystyle x_{n}=q_{-}^{n}}$. To compare them with the exact solutions, Taylor expansions are computed:

${\displaystyle {\begin{aligned}q_{+}&=1+{\tfrac {1}{2}}(wh)^{2}+wh\left(1+{\tfrac {1}{8}}(wh)^{2}-{\tfrac {3}{128}}(wh)^{4}+{\mathcal {O}}(h^{6})\right)\\&=1+(wh)+{\tfrac {1}{2}}(wh)^{2}+{\tfrac {1}{8}}(wh)^{3}-{\tfrac {3}{128}}(wh)^{5}+{\mathcal {O}}(h^{7}).\end{aligned}}}$

The quotient of this series with the one of the exponential ${\displaystyle e^{wh}}$ starts with ${\displaystyle 1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}(h^{5})}$, so

${\displaystyle {\begin{aligned}q_{+}&=\left(1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}(h^{5})\right)e^{wh}\\&=e^{-{\frac {1}{24}}(wh)^{3}+{\mathcal {O}}(h^{5})}\,e^{wh}.\end{aligned}}}$

From there it follows that for the first basis solution the error can be computed as

${\displaystyle {\begin{aligned}x_{n}=q_{+}^{\;n}&=e^{-{\frac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}(h^{4})}\,e^{wt_{n}}\\&=e^{wt_{n}}\left(1-{\tfrac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}(h^{4})\right)\\&=e^{wt_{n}}+{\mathcal {O}}(h^{2}t_{n}e^{wt_{n}}).\end{aligned}}}$

That is, although the local discretization error is of order 4, due to the second order of the differential equation the global error is of order 2, with a constant that grows exponentially in time.

Note that at the start of the Verlet iteration at step ${\displaystyle n=1}$, time ${\displaystyle t=t_{1}=\Delta t}$, computing ${\displaystyle {\vec {x}}_{2}}$, one already needs the position vector ${\displaystyle {\vec {x}}_{1}}$ at time ${\displaystyle t=t_{1}}$. At first sight, this could give problems, because the initial conditions are known only at the initial time ${\displaystyle t_{0}=0}$. However, from these the acceleration ${\displaystyle {\vec {a}}_{0}={\vec {A}}({\vec {x}}_{0})}$ is known, and a suitable approximation for the position at the first time step can be obtained using the Taylor polynomial of degree two:

${\displaystyle {\vec {x}}_{1}={\vec {x}}_{0}+{\vec {v}}_{0}\Delta t+{\tfrac {1}{2}}{\vec {a}}_{0}\Delta t^{2}\approx {\vec {x}}(\Delta t)+{\mathcal {O}}(\Delta t^{3}).}$

The error on the first time step then is of order ${\displaystyle {\mathcal {O}}(\Delta t^{3})}$. This is not considered a problem because on a simulation over a large number of time steps, the error on the first time step is only a negligibly small amount of the total error, which at time ${\displaystyle t_{n}}$ is of the order ${\displaystyle {\mathcal {O}}(e^{Lt_{n}}\Delta t^{2})}$, both for the distance of the position vectors ${\displaystyle {\vec {x}}_{n}}$ to ${\displaystyle {\vec {x}}(t_{n})}$ as for the distance of the divided differences ${\displaystyle {\tfrac {{\vec {x}}_{n+1}-{\vec {x}}_{n}}{\Delta t}}}$ to ${\displaystyle {\tfrac {{\vec {x}}(t_{n+1})-{\vec {x}}(t_{n})}{\Delta t}}}$. Moreover, to obtain this second-order global error, the initial error needs to be of at least third order.

A disadvantage of the Störmer–Verlet method is that if the time step (${\displaystyle \Delta t}$) changes, the method does not approximate the solution to the differential equation. This can be corrected using the formula[4]

${\displaystyle {\vec {x}}_{i+1}={\vec {x}}_{i}+({\vec {x}}_{i}-{\vec {x}}_{i-1})(\Delta t_{i}/\Delta t_{i-1})+{\vec {a}}_{i}\Delta t_{i}^{2}.}$

A more exact derivation uses the Taylor series (to second order) at ${\displaystyle t_{i}}$ for times ${\displaystyle t_{i+1}=t_{i}+\Delta t_{i}}$ and ${\displaystyle t_{i-1}=t_{i}-\Delta t_{i-1}}$ to obtain after elimination of ${\displaystyle {\vec {v}}_{i}}$

${\displaystyle {\frac {{\vec {x}}_{i+1}-{\vec {x}}_{i}}{\Delta t_{i}}}+{\frac {{\vec {x}}_{i-1}-{\vec {x}}_{i}}{\Delta t_{i-1}}}={\vec {a}}_{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}},}$

so that the iteration formula becomes

${\displaystyle {\vec {x}}_{i+1}={\vec {x}}_{i}+({\vec {x}}_{i}-{\vec {x}}_{i-1}){\frac {\Delta t_{i}}{\Delta t_{i-1}}}+{\vec {a}}_{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}}\,\Delta t_{i}.}$

The velocities are not explicitly given in the basic Störmer equation, but often they are necessary for the calculation of certain physical quantities like the kinetic energy. This can create technical challenges in molecular dynamics simulations, because kinetic energy and instantaneous temperatures at time ${\displaystyle t}$ cannot be calculated for a system until the positions are known at time ${\displaystyle t+\Delta t}$. This deficiency can either be dealt with using the velocity Verlet algorithm or by estimating the velocity using the position terms and the mean value theorem:

${\displaystyle {\vec {v}}(t)={\frac {{\vec {x}}(t+\Delta t)-{\vec {x}}(t-\Delta t)}{2\Delta t}}+{\mathcal {O}}(\Delta t^{2}).}$

Note that this velocity term is a step behind the position term, since this is for the velocity at time ${\displaystyle t}$, not ${\displaystyle t+\Delta t}$, meaning that ${\displaystyle {\vec {v}}_{n}={\tfrac {{\vec {x}}_{n+1}-{\vec {x}}_{n-1}}{2\Delta t}}}$ is a second-order approximation to ${\displaystyle {\vec {v}}(t_{n})}$. With the same argument, but halving the time step, ${\displaystyle {\vec {v}}_{n+1/2}={\tfrac {{\vec {x}}_{n+1}-{\vec {x}}_{n}}{\Delta t}}}$ is a second-order approximation to ${\displaystyle {\vec {v}}(t_{n+1/2})}$, with ${\displaystyle t_{n+1/2}=t_{n}+{\tfrac {1}{2}}\Delta t}$.

One can shorten the interval to approximate the velocity at time ${\displaystyle t+\Delta t}$ at the cost of accuracy:

${\displaystyle {\vec {v}}(t+\Delta t)={\frac {{\vec {x}}(t+\Delta t)-{\vec {x}}(t)}{\Delta t}}+{\mathcal {O}}(\Delta t).}$

Since velocity Verlet is a generally useful algorithm in 3D applications, a general solution written in C++ could look like below. A simplified drag force is used to demonstrate change in acceleration, however it is only needed if acceleration is not constant.

struct Body
{
    Vec3d pos { 0.0, 0.0, 0.0 };
    Vec3d vel { 2.0, 0.0, 0.0 }; // 2m/s along x-axis
    Vec3d acc { 0.0, 0.0, 0.0 }; // no acceleration at first
    double mass = 1.0; // 1kg
    double drag = 0.1; // rho*C*Area - simplified drag for this example

    /**
     * Update pos and vel using "Velocity Verlet" integration
     * @param dt DeltaTime / time step [eg: 0.01]
     */
    void update(double dt)
    {
        Vec3d new_pos = pos + vel*dt + acc*(dt*dt*0.5);
        Vec3d new_acc = apply_forces(); // only needed if acceleration is not constant
        Vec3d new_vel = vel + (acc+new_acc)*(dt*0.5);
        pos = new_pos;
        vel = new_vel;
        acc = new_acc;
    }

    Vec3d apply_forces() const
    {
        Vec3d grav_acc = Vec3d{0.0, 0.0, -9.81 }; // 9.81m/s^2 down in the Z-axis
        Vec3d drag_force = 0.5 * drag * (vel * abs(vel)); // D = 0.5 * (rho * C * Area * vel^2)
        Vec3d drag_acc = drag_force / mass; // a = F/m
        return grav_acc - drag_acc;
    }
};