Editing Verlet integration (section)

==Basic Størmer–Verlet==

For a [[second-order differential equation]] of the type <math>\ddot{\mathbf x}(t) = \mathbf A\bigl(\mathbf x(t)\bigr)</math> with initial conditions <math>\mathbf x(t_0) = \mathbf x_0</math> and <math>\dot{\mathbf x}(t_0) = \mathbf v_0</math>, an approximate numerical solution <math>\mathbf x_n \approx \mathbf x(t_n)</math> at the times <math>t_n = t_0 + n\,\Delta t</math> with step size <math>\Delta t > 0</math> can be obtained by the following method:
# set <math display="inline">\mathbf x_1 = \mathbf x_0 + \mathbf v_0\,\Delta t + \tfrac 1 2 \mathbf A(\mathbf x_0)\,\Delta t^2</math>,
# for ''n'' = 1, 2, ... iterate <math display="block">
\mathbf x_{n+1} = 2 \mathbf x_n - \mathbf x_{n-1} + \mathbf A(\mathbf x_n)\,\Delta t^2.
</math>

===Equations of motion===
Newton's equation of motion for conservative physical systems is
:<math>\boldsymbol M \ddot{\mathbf x}(t) = F\bigl(\mathbf x(t)\bigr) = -\nabla V\bigl(\mathbf x(t)\bigr),</math>
or individually
:<math>m_k \ddot{\mathbf x}_k(t) = F_k\bigl(\mathbf x(t)\bigr) = -\nabla_{\mathbf x_k} V\left(\mathbf x(t)\right),</math>
where
* <math>t</math> is the time,
* <math>\mathbf x(t) = \bigl(\mathbf x_1(t), \ldots, \mathbf x_N(t)\bigr)</math> is the ensemble of the position vector of <math>N</math> objects,
* <math>V</math> is the scalar potential function,
* <math>F</math> is the negative [[Potential gradient|gradient of the potential]], giving the ensemble of forces on the particles,
* <math>\boldsymbol M</math> is the [[mass matrix]], typically diagonal with blocks with mass <math>m_k</math> for every particle.

This equation, for various choices of the potential function <math>V</math>, can be used to describe the evolution of diverse physical systems, from the motion of [[molecular dynamics|interacting molecules]] to the [[N-body problem|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
:<math>\ddot{\mathbf x}(t) = \mathbf A\bigl(\mathbf x(t)\bigr)</math>
with some suitable vector-valued function <math>\mathbf A(\mathbf x)</math> representing the position-dependent acceleration. Typically, an initial position <math>\mathbf x(0) = \mathbf x_0</math> and an initial velocity <math>\mathbf v(0) = \dot{\mathbf x}(0) = \mathbf v_0</math> are also given.

===Verlet integration (without velocities)===
To discretize and numerically solve this [[initial value problem]], a time step <math>\Delta t > 0</math> is chosen, and the sampling-point sequence <math>t_n = n\,\Delta t</math> considered. The task is to construct a sequence of points <math>\mathbf x_n</math> that closely follow the points <math>\mathbf x(t_n)</math> 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:
:<math>\begin{align}
\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{align}</math>

''Verlet integration'' in the form used as the ''Størmer method''<ref>[http://www.fisica.uniud.it/~ercolessi/md/md/node21.html webpage] {{Webarchive|url=https://web.archive.org/web/20040803212552/http://www.fisica.uniud.it/~ercolessi/md/md/node21.html |date=2004-08-03 }} with a description of the Størmer method.</ref> uses this equation to obtain the next position vector from the previous two without using the velocity as
:<math>\begin{align}
\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{align}</math>

===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 <math>\Delta t</math> of degree three. The local error is quantified by inserting the exact values <math>\mathbf x(t_{n-1}), \mathbf x(t_n), \mathbf x(t_{n+1})</math> into the iteration and computing the [[Taylor expansion]]s at time <math>t = t_n</math> of the position vector <math>\mathbf{x}(t \pm \Delta t)</math> in different time directions:
:<math>\begin{align}
\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)\\
\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{align}</math>
where <math>\mathbf{x}</math> is the position, <math>\mathbf{v} = \dot{\mathbf x}</math> the velocity, <math>\mathbf{a} = \ddot{\mathbf x}</math> the acceleration, and <math>\mathbf{b} = \dot{\mathbf a} = \overset{\dots}{\mathbf x}</math> the [[Jerk (physics)|jerk]] ([[third derivative]] of the position with respect to the time).

Adding these two expansions gives
:<math>\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).</math>
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, <math>\mathbf a(t) = \mathbf A\bigl(\mathbf x(t)\bigr)</math>, whereas in the iteration it is computed at the central iteration point, <math>\mathbf a_n = \mathbf A(\mathbf x_n)</math>. 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 <math>\ddot x(t) = w^2 x(t)</math> with a constant <math>w</math>. Its exact basis solutions are <math>e^{wt}</math> and <math>e^{-wt}</math>.

The Størmer method applied to this differential equation leads to a linear [[recurrence relation]]
:<math>x_{n+1} - 2x_n + x_{n-1} = h^2 w^2 x_n,</math>
or
:<math>x_{n+1} - 2\left(1 + \tfrac12(wh)^2\right) x_n + x_{n-1} = 0.</math>
It can be solved by finding the roots of its characteristic polynomial
<math>q^2 - 2\left(1 + \tfrac12(wh)^2\right)q + 1 = 0</math>. These are
:<math>q_\pm = 1 + \tfrac 1 2 (wh)^2 \pm wh \sqrt{1 + \tfrac 1 4 (wh)^2}.</math>
The basis solutions of the linear recurrence are <math>x_n = q_+^n</math> and <math>x_n = q_-^n</math>. To compare them with the exact solutions, Taylor expansions are computed:
:<math>\begin{align}
q_+ &= 1 + \tfrac12(wh)^2 + wh\left(1 + \tfrac18(wh)^2 - \tfrac{3}{128}(wh)^4 + \mathcal O\left(h^6\right)\right)\\
    &= 1 + (wh) + \tfrac12(wh)^2 + \tfrac18(wh)^3 - \tfrac{3}{128}(wh)^5 + \mathcal O\left(h^7\right).
\end{align}</math>
The quotient of this series with the exponential <math>e^{wh}</math> starts with <math>1 - \tfrac1{24}(wh)^3 + \mathcal O\left(h^5\right)</math>, so
:<math>\begin{align}
q_+ &= \left(1 - \tfrac1{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{align}</math>
From there it follows that for the first basis solution the error can be computed as
:<math>\begin{align}
x_n = q_+^{n}
 &= e^{-\frac{1}{24}(wh)^2\,wt_n + \mathcal O\left(h^4\right)}\,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\left(h^2 t_n e^{wt_n}\right).
\end{align}</math>
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.

===Starting the iteration===
Note that at the start of the Verlet iteration at step <math>n = 1</math>, time <math>t = t_1 = \Delta t</math>, computing <math>\mathbf x_2</math>, one already needs the position vector <math>\mathbf x_1</math> at time <math>t = t_1</math>. At first sight, this could give problems, because the initial conditions are known only at the initial time <math>t_0 = 0</math>. However, from these the acceleration <math>\mathbf a_0 = \mathbf A(\mathbf x_0)</math> is known, and a suitable approximation for the position at the first time step can be obtained using the [[Taylor polynomial]] of degree two:
:<math>\mathbf x_1 = \mathbf{x}_0 + \mathbf{v}_0 \Delta t + \tfrac12 \mathbf a_0 \Delta t^2
\approx \mathbf{x}(\Delta t) + \mathcal{O}\left(\Delta t^3\right).</math>

The error on the first time step then is of order <math>\mathcal O\left(\Delta t^3\right)</math>. 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 <math>t_n</math> is of the order <math>\mathcal O\left(e^{Lt_n} \Delta t^2\right)</math>, both for the distance of the position vectors <math>\mathbf x_n</math> to <math>\mathbf x(t_n)</math> as for the distance of the divided differences <math>\tfrac{\mathbf x_{n+1} - \mathbf x_n}{\Delta t}</math> to <math>\tfrac{\mathbf x(t_{n+1}) - \mathbf x(t_n)}{\Delta t}</math>. Moreover, to obtain this second-order global error, the initial error needs to be of at least third order.

===Non-constant time differences===
A disadvantage of the Størmer–Verlet method is that if the time step (<math>\Delta t</math>) changes, the method does not approximate the solution to the differential equation. This can be corrected using the formula<ref>{{cite web| last=Dummer|first=Jonathan|title=A Simple Time-Corrected Verlet Integration Method| url=http://lonesock.net/article/verlet.html}}</ref>
:<math>
\mathbf x_{i+1}
= \mathbf x_i + \left(\mathbf x_i - \mathbf x_{i-1}\right) \frac{\Delta t_i}{\Delta t_{i-1}} + \mathbf a_i \Delta t_i^2.
</math>

A more exact derivation uses the Taylor series (to second order) at <math>t_i</math> for times <math>t_{i+1} = t_i + \Delta t_i</math> and <math>t_{i-1} = t_i - \Delta t_{i-1}</math> to obtain after elimination of <math>\mathbf v_i</math>
:<math>
\frac{\mathbf x_{i+1} - \mathbf x_i}{\Delta t_i}
+ \frac{\mathbf x_{i-1} - \mathbf x_i}{\Delta t_{i-1}}
= \mathbf a_i\,\frac{\Delta t_{i} + \Delta t_{i-1}}2,
</math>
so that the iteration formula becomes
:<math>
\mathbf x_{i+1}
= \mathbf x_i
+ (\mathbf x_i - \mathbf x_{i-1}) \frac{\Delta t_i}{\Delta t_{i-1}}
+ \mathbf a_i\,\frac{\Delta t_{i} + \Delta t_{i-1}}2\,\Delta t_i.
</math>

===Computing velocities – Størmer–Verlet method===
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 <math>t</math> cannot be calculated for a system until the positions are known at time <math>t + \Delta t</math>. This deficiency can either be dealt with using the [[#Velocity Verlet|velocity Verlet]] algorithm or by estimating the velocity using the position terms and the [[mean value theorem]]:
:<math>
\mathbf{v}(t) =
\frac{\mathbf{x}(t + \Delta t) - \mathbf{x}(t - \Delta t)}{2\Delta t}
+ \mathcal{O}\left(\Delta t^2\right).
</math>

Note that this velocity term is a step behind the position term, since this is for the velocity at time <math>t</math>, not <math>t + \Delta t</math>, meaning that <math>\mathbf v_n = \tfrac{\mathbf x_{n+1} - \mathbf x_{n-1}}{2\Delta t}</math> is a second-order approximation to <math>\mathbf{v}(t_n)</math>. With the same argument, but halving the time step, <math>\mathbf v_{n+\frac12} = \tfrac{\mathbf x_{n+1} - \mathbf x_n}{\Delta t}</math> is a second-order approximation to <math>\mathbf{v}\left(t_{n+\frac12}\right)</math>, with <math>t_{n+\frac12} = t_n + \tfrac12 \Delta t</math>.

One can shorten the interval to approximate the velocity at time <math>t + \Delta t</math> at the cost of accuracy:
:<math>\mathbf{v}(t + \Delta t) = \frac{\mathbf{x}(t + \Delta t) - \mathbf{x}(t)}{\Delta t} + \mathcal{O}(\Delta t).</math>