Verlet integration (French pronunciation: [vɛʁˈlɛ]) is a numerical method used to integrate Newton's equations of motion. It is frequently used to calculate trajectories of particles in molecular dynamics simulations and computer graphics. The algorithm was first used in 1791 by Jean Baptiste Delambre and has been rediscovered many times since then, most recently by Loup Verlet in the 1960s for use in molecular dynamics. It was also used by P. H. Cowell and A. C. C. Crommelin in 1909 to compute the orbit of Halley's Comet, and by Carl Størmer in 1907 to study the trajectories of electrical particles in a magnetic field (hence it is also called Størmer's method). The Verlet integrator provides good numerical stability, as well as other properties that are important in physical systems such as time reversibility and preservation of the symplectic form on phase space, at no significant additional computational cost over the simple Euler method.
Basic Størmer–Verlet For a second-order differential equation of the type x ¨ ( t ) = A ( x ( t ) ) {\displaystyle {\ddot {\mathbf {x} }}(t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}} with initial conditions x ( t 0 ) = x 0 {\displaystyle \mathbf {x} (t_{0})=\mathbf {x} _{0}} and x ˙ ( t 0 ) = v 0 {\displaystyle {\dot {\mathbf {x} }}(t_{0})=\mathbf {v} _{0}} , an approximate numerical solution x n ≈ x ( t n ) {\displaystyle \mathbf {x} _{n}\approx \mathbf {x} (t_{n})} at the times t n = t 0 + n Δ t {\displaystyle t_{n}=t_{0}+n\,\Delta t} with step size Δ t > 0 {\displaystyle \Delta t>0} can be obtained by the following method:
set x 1 = x 0 + v 0 Δ t + 1 2 A ( x 0 ) Δ t 2 {\textstyle \mathbf {x} _{1}=\mathbf {x} _{0}+\mathbf {v} _{0}\,\Delta t+{\tfrac {1}{2}}\mathbf {A} (\mathbf {x} _{0})\,\Delta t^{2}} , for n = 1, 2, ... iterate x n + 1 = 2 x n − x n − 1 + A ( x n ) Δ t 2 . {\displaystyle \mathbf {x} _{n+1}=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {A} (\mathbf {x} _{n})\,\Delta t^{2}.}
Equations of motion Newton's equation of motion for conservative physical systems is
M x ¨ ( t ) = F ( x ( t ) ) = − ∇ V ( x ( t ) ) , {\displaystyle {\boldsymbol {M}}{\ddot {\mathbf {x} }}(t)=F{\bigl (}\mathbf {x} (t){\bigr )}=-\nabla V{\bigl (}\mathbf {x} (t){\bigr )},}
or individually
m k x ¨ k ( t ) = F k ( x ( t ) ) = − ∇ x k V ( x ( t ) ) , {\displaystyle m_{k}{\ddot {\mathbf {x} }}_{k}(t)=F_{k}{\bigl (}\mathbf {x} (t){\bigr )}=-\nabla _{\mathbf {x} _{k}}V{\left(\mathbf {x} (t)\right)},}
where
t {\displaystyle t} is the time,
x ( t ) = ( x 1 ( t ) , … , x N ( t ) ) {\displaystyle \mathbf {x} (t)={\bigl (}\mathbf {x} _{1}(t),\ldots ,\mathbf {x} _{N}(t){\bigr )}} is the ensemble of the position vector of N {\displaystyle N} objects,
V {\displaystyle V} is the scalar potential function,
F {\displaystyle F} is the negative gradient of the potential, giving the ensemble of forces on the particles,
M {\displaystyle {\boldsymbol {M}}} is the mass matrix, typically diagonal with blocks with mass m k {\displaystyle m_{k}} for every particle. This equation, for various choices of the potential function V {\displaystyle 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
x ¨ ( t ) = A ( x ( t ) ) {\displaystyle {\ddot {\mathbf {x} }}(t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}}
with some suitable vector-valued function A ( x ) {\displaystyle \mathbf {A} (\mathbf {x} )} representing the position-dependent acceleration. Typically, an initial position x ( 0 ) = x 0 {\displaystyle \mathbf {x} (0)=\mathbf {x} _{0}} and an initial velocity v ( 0 ) = x ˙ ( 0 ) = v 0 {\displaystyle \mathbf {v} (0)={\dot {\mathbf {x} }}(0)=\mathbf {v} _{0}} are also given.
Verlet integration (without velocities) To discretize and numerically solve this initial value problem, a time step Δ t > 0 {\displaystyle \Delta t>0} is chosen, and the sampling-point sequence t n = n Δ t {\displaystyle t_{n}=n\,\Delta t} considered. The task is to construct a sequence of points x n {\displaystyle \mathbf {x} _{n}} that closely follow the points x ( t n ) {\displaystyle \mathbf {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:
Δ 2 x n Δ t 2 = x n + 1 − x n Δ t − x n − x n − 1 Δ t Δ t = x n + 1 − 2 x n + x n − 1 Δ t 2 = a n = A ( x n ) . {\displaystyle {\begin{aligned}{\frac {\Delta ^{2}\mathbf {x} _{n}}{\Delta t^{2}}}&={\frac {{\frac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}-{\frac {\mathbf {x} _{n}-\mathbf {x} _{n-1}}{\Delta t}}}{\Delta t}}\\[6pt]&={\frac {\mathbf {x} _{n+1}-2\mathbf {x} _{n}+\mathbf {x} _{n-1}}{\Delta t^{2}}}\\&=\mathbf {a} _{n}=\mathbf {A} (\mathbf {x} _{n}).\end{aligned}}}
Verlet integration in the form used as the Størmer method uses this equation to obtain the next position vector from the previous two without using the velocity as
x n + 1 = 2 x n − x n − 1 + a n Δ t 2 , a n = A ( x n ) . {\displaystyle {\begin{aligned}\mathbf {x} _{n+1}&=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {a} _{n}\,\Delta t^{2},\\[6pt]\mathbf {a} _{n}&=\mathbf {A} (\mathbf {x} _{n}).\end{aligned}}}
Discretisation error 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 Δ t {\displaystyle \Delta t} of degree three. The local error is quantified by inserting the exact values x ( t n − 1 ) , x ( t n ) , x ( t n + 1 ) {\displaystyle \mathbf {x} (t_{n-1}),\mathbf {x} (t_{n}),\mathbf {x} (t_{n+1})} into the iteration and computing the Taylor expansions at time t = t n {\displaystyle t=t_{n}} of the position vector x ( t ± Δ t ) {\displaystyle \mathbf {x} (t\pm \Delta t)} in different time directions:
x ( t + Δ t ) = x ( t ) + v ( t ) Δ t + a ( t ) Δ t 2 2 + b ( t ) Δ t 3 6 + O ( Δ t 4 ) x ( t − Δ t ) = x ( t ) − v ( t ) Δ t + a ( t ) Δ t 2 2 − b ( t ) Δ t 3 6 + O ( Δ t 4 ) , {\displaystyle {\begin{aligned}\mathbf {x} (t{+}\Delta t)&=\mathbf {x} (t)+\mathbf {v} (t)\Delta t+{\frac {\mathbf {a} (t)\Delta t^{2}}{2}}+{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}{\left(\Delta t^{4}\right)}\\[1ex]\mathbf {x} (t{-}\Delta t)&=\mathbf {x} (t)-\mathbf {v} (t)\Delta t+{\frac {\mathbf {a} (t)\Delta t^{2}}{2}}-{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}{\left(\Delta t^{4}\right)},\end{aligned}}}
where x {\displaystyle \mathbf {x} } is the position, v = x ˙ {\displaystyle \mathbf {v} ={\dot {\mathbf {x} }}} the velocity, a = x ¨ {\displaystyle \mathbf {a} ={\ddot {\mathbf {x} }}} the acceleration, and b = a ˙ = x … {\displaystyle \mathbf {b} ={\dot {\mathbf {a} }}={\overset {\dots }{\mathbf {x} }}} the jerk (third derivative of the position with respect to the time). Adding these two expansions gives
x ( t + Δ t ) = 2 x ( t ) − x ( t − Δ t ) + a ( t ) Δ t 2 + O ( Δ t 4 ) . {\displaystyle \mathbf {x} (t{+}\Delta t)=2\mathbf {x} (t)-\mathbf {x} (t{-}\Delta t)+\mathbf {a} (t)\Delta t^{2}+{\mathcal {O}}{\left(\Delta t^{4}\right)}.}
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, a ( t ) = A ( x ( t ) ) {\displaystyle \mathbf {a} (t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}} , whereas in the iteration it is computed at the central iteration point, a n = A ( x n ) {\displaystyle \mathbf {a} _{n}=\mathbf {A} (\mathbf {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.
A simple example 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 x ¨ ( t ) = w 2 x ( t ) {\displaystyle {\ddot {x}}(t)=w^{2}x(t)} with a constant w {\displaystyle w} . Its exact basis solutions are e w t {\displaystyle e^{wt}} and e − w t {\displaystyle e^{-wt}} . The Størmer method applied to this differential equation leads to a linear recurrence relation
x n + 1 − 2 x n + x n − 1 = h 2 w 2 x n , {\displaystyle x_{n+1}-2x_{n}+x_{n-1}=h^{2}w^{2}x_{n},}
or
x n + 1 − 2 ( 1 + 1 2 ( w h ) 2 ) x n + x n − 1 = 0. {\displaystyle x_{n+1}-2\left(1+{\tfrac {1}{2}}\left(wh\right)^{2}\right)x_{n}+x_{n-1}=0.}
It can be solved by finding the roots of its characteristic polynomial
q 2 − 2 ( 1 + 1 2 ( w h ) 2 ) q + 1 = 0 {\displaystyle q^{2}-2\left(1+{\tfrac {1}{2}}(wh)^{2}\right)q+1=0} . These are
q ± = 1 + 1 2 ( w h ) 2 ± w h 1 + 1 4 ( w h ) 2 . {\displaystyle q_{\pm }=1+{\tfrac {1}{2}}\left(wh\right)^{2}\pm wh{\sqrt {1+{\tfrac {1}{4}}\left(wh\right)^{2}}}.}
The basis solutions of the linear recurrence are x n = q + n {\displaystyle x_{n}=q_{+}^{n}} and x n = q − n {\displaystyle x_{n}=q_{-}^{n}} . To compare them with the exact solutions, Taylor expansions are computed:
q + = 1 + 1 2 ( w h ) 2 + w h ( 1 + 1 8 ( w h ) 2 − 3 128 ( w h ) 4 + O ( h 6 ) ) = 1 + ( w h ) + 1 2 ( w h ) 2 + 1 8 ( w h ) 3 − 3 128 ( w h ) 5 + O ( h 7 ) . {\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}}{\left(h^{6}\right)}\right)\\&=1+(wh)+{\tfrac {1}{2}}(wh)^{2}+{\tfrac {1}{8}}(wh)^{3}-{\tfrac {3}{128}}(wh)^{5}+{\mathcal {O}}{\left(h^{7}\right)}.\end{aligned}}}
The quotient of this series with the exponential e w h {\displaystyle e^{wh}} starts with 1 − 1 24 ( w h ) 3 + O ( h 5 ) {\displaystyle 1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}{\left(h^{5}\right)}} , so
q + = ( 1 − 1 24 ( w h ) 3 + O ( h 5 ) ) e w h = e − 1 24 ( w h ) 3 + O ( h 5 ) e w h . {\displaystyle {\begin{aligned}q_{+}&=\left(1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}{\left(h^{5}\right)}\right)e^{wh}\\&=e^{-{\frac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)}\,e^{wh}.\end{aligned}}}
From there it follows that for the first basis solution the error can be computed as
x n