Editing Verlet integration (section)

== Velocity Verlet ==
A related, and more commonly used algorithm is the velocity Verlet algorithm,<ref>{{cite journal|last=Swope|first=William C. |author2=H. C. Andersen |author3=P. H. Berens |author4=K. R. Wilson|title=A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters|journal=The Journal of Chemical Physics|date=1 January 1982|volume=76|issue=1|pages=648 (Appendix)|doi=10.1063/1.442716|bibcode=1982JChPh..76..637S }}</ref> similar to the [[leapfrog method]], except that the velocity and position are calculated at the same value of the time variable (leapfrog does not, as the name suggests). This uses a similar approach, but explicitly incorporates velocity, solving the problem of the first time step in the basic Verlet algorithm:

:<math>\begin{align}
\mathbf{x}(t + \Delta t) &= \mathbf{x}(t) + \mathbf{v}(t)\, \Delta t + \tfrac{1}{2} \,\mathbf{a}(t) \Delta t^2, \\[6pt]
\mathbf{v}(t + \Delta t) &= \mathbf{v}(t) + \frac{\mathbf{a}(t) + \mathbf{a}(t + \Delta t)}{2} \Delta t.
\end{align}</math>

It can be shown that the error in the velocity Verlet is of the same order as in the basic Verlet. Note that the velocity algorithm is not necessarily more memory-consuming, because, in basic Verlet, we keep track of two vectors of position, while in velocity Verlet, we keep track of one vector of position and one vector of velocity. The standard implementation scheme of this algorithm is:
# Calculate <math>\mathbf{v}\left(t + \tfrac12\,\Delta t\right) = \mathbf{v}(t) + \tfrac12\,\mathbf{a}(t)\,\Delta t</math>.
# Calculate <math>\mathbf{x}(t + \Delta t) = \mathbf{x}(t) + \mathbf{v}\left(t + \tfrac12\,\Delta t\right)\, \Delta t</math>.
# Derive <math>\mathbf{a}(t + \Delta t)</math> from the interaction potential using <math>\mathbf{x}(t + \Delta t)</math>.
# Calculate <math>\mathbf{v}(t + \Delta t) = \mathbf{v}\left(t + \tfrac12\,\Delta t\right) + \tfrac12\,\mathbf{a}(t + \Delta t)\Delta t</math>.
This algorithm also works with variable time steps, and is identical to the 'kick-drift-kick' form of [[leapfrog method]] integration.

Eliminating the half-step velocity, this algorithm may be shortened to
# Calculate <math>\mathbf{x}(t + \Delta t) = \mathbf{x}(t) + \mathbf{v}(t)\,\Delta t + \tfrac12 \,\mathbf{a}(t)\,\Delta t^2</math>.
# Derive <math>\mathbf{a}(t + \Delta t)</math> from the interaction potential using <math>\mathbf{x}(t + \Delta t)</math>.
# Calculate <math>\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \tfrac12\,\bigl(\mathbf{a}(t) + \mathbf{a}(t + \Delta t)\bigr)\Delta t</math>.

Note, however, that this algorithm assumes that acceleration <math>\mathbf{a}(t + \Delta t)</math> only depends on position <math>\mathbf{x}(t + \Delta t)</math> and does not depend on velocity <math>\mathbf{v}(t + \Delta t)</math>.

One might note that the long-term results of velocity Verlet, and similarly of leapfrog are one order better than the [[semi-implicit Euler method]]. The algorithms are almost identical up to a shift by half a time step in the velocity.  This can be proven by rotating the above loop to start at step 3 and then noticing that the acceleration term in step 1 could be eliminated by combining steps 2 and 4.  The only difference is that the midpoint velocity in velocity Verlet is considered the final velocity in semi-implicit Euler method.

The global error of all Euler methods is of order one, whereas the global error of this method is, similar to the [[midpoint method]], of order two. Additionally, if the acceleration indeed results from the forces in a conservative mechanical or [[Hamiltonian system]], the energy of the approximation essentially oscillates around the constant energy of the exactly solved system, with a global error bound again of order one for semi-explicit Euler and order two for Verlet-leapfrog. The same goes for all other conserved quantities of the system like linear or angular momentum, that are always preserved or nearly preserved in a [[symplectic integrator]].<ref name="Hairer2003" />

The velocity Verlet method is a special case of the [[Newmark-beta method]] with <math>\beta = 0</math> and <math>\gamma = \tfrac12</math>.

=== Algorithmic representation ===
Since '''velocity Verlet''' is a generally useful algorithm in 3D applications, a solution written in C++ could look like below. This type of position integration will significantly increase accuracy in 3D simulations and games when compared with the regular Euler method.<syntaxhighlight lang="c++" line="1">
struct Body
{
    Vec3d pos { 0.0, 0.0, 0.0 };
    Vec3d vel { 2.0, 0.0, 0.0 }; // 2 m/s along x-axis
    Vec3d acc { 0.0, 0.0, 0.0 }; // no acceleration at first
    double mass = 1.0; // 1kg

    /**
     * Updates 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();
        Vec3d new_vel = vel + (acc+new_acc)*(dt*0.5);
        pos = new_pos;
        vel = new_vel;
        acc = new_acc;
    }

    /**
     * To apply velocity to your objects, calculate the required Force vector instead
     * and apply the accumulated forces here.
     */
    Vec3d apply_forces() const
    {
        Vec3d new_acc = Vec3d{0.0, 0.0, -9.81 }; // 9.81 m/s² down in the z-axis
        // Apply any other forces here...
        // NOTE: Avoid depending on `vel` because Velocity Verlet assumes acceleration depends on position.
        return new_acc;
    }
};
</syntaxhighlight>