Verlet integration
Encyclopedia
Verlet integration 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 Verlet integrator offers greater 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 cost over the simple Euler method. Verlet integration was used by Carl Størmer
to compute the trajectories of particles moving in a magnetic field (hence it is also called Störmer's method) and was popularized in molecular dynamics by French physicist Loup Verlet
in 1967.
or individually
where
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
with some suitable vector valued function A representing the position dependent acceleration. Typically, an initial position
and an initial velocity 
are also given.
, a time step
is chosen and the sampling point sequence 
considered. The task is to construct a sequence of points 
that closely follow the points 
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:
The Störmer method uses this equation to obtain the next position vector from the previous two as
without using the velocity. 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 h of degree three. The local error is quantified by inserting the exact values
into the iteration and computing the Taylor expansions at time 
of the position vector 
in different time directions.
where
is the position, 
the velocity, 
the acceleration and 
the jerk (third derivative of the position with respect to the time) 
.
Adding these two expansions gives
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,
, whereas in the iteration it is computed at the central iteration point, 
. 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.
.
Consider the linear differential equation
with a constant w. Its exact basis solutions are 
and 
.
The Störmer method applied to this differential equation leads to a linear recurrence relation
It can be solved by finding the roots of its characteristic polynomial
. These are
.
The basis solutions of the linear recurrence are
and 
. To compare them with the exact solutions, Taylor expansions are computed.
The quotient of this series with the one of the exponential
starts with 
, so
From there it follows that for the first basis solution the error can be computed as
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.
, time 
, computing 
, one already needs the position vector 
at time 
. At first sight this could give problems, because the initial conditions are known only at the initial time 
. However, from these the acceleration 
is known, and a suitable approximation for the first time step position can be obtained using the Taylor polynomial of degree two:
The error on the first time step calculation then is of order
. This is not considered a problem because on a simulation of over a large amount of timesteps, the error on the first timestep is only a negligible small amount of the total error, which at time 
is of the order 
, both for the distance of the position vectors 
to 
as for the distance of the divided differences 
to 
. Moreover, to obtain this second order global error, the initial error needs to be of at least third order.
simulations, because kinetic energy and instantaneous temperatures at time
cannot be calculated for a system until the positions are known at time 
. This deficiency can either be dealt with using the Velocity Verlet algorithm, or estimating the velocity using the position terms and the mean value theorem
:
Note that this velocity term is a step behind the position term, since this is for the velocity at time
, not 
, meaning that 
is an order two approximation to 
. With the same argument, but halving the time step, 
is an order two approximation to 
, with 
.
One can shorten the interval to approximate the velocity at time
at the cost of accuracy:
It can be shown that the error on the Velocity Verlet is of the same order as the Basic Verlet. Note that the Velocity algorithm is not necessarily more memory consuming, because it's not necessary to keep track of the velocity at every timestep during the simulation. The standard implementation scheme of this algorithm is:
Eliminating the half-step velocity, this algorithm may be shortened to
Note, however, that this algorithm assumes that acceleration
only depends on position 
, and does not depend on velocity 
.
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 shifted by half of a timestep in the velocity. This is easily 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 conservered quantities of the system like linear or angular momentum, that are always preserved or nearly preserved in a symplectic integrator
.
as described above, and the local error in velocity is 
.
The global error in position, in contrast, is
and the global error in velocity is 
. These can be derived by noting the following:
and
Therefore:
Similarly:
Which can be generalized to (it can be shown by induction, but it is given here without proof):
If we consider the global error in position between
and 
, where 
, it is clear that:
And therefore, the global (cumulative) error over a constant interval of time is given by:
Because the velocity is determined in a non-cumulative way from the positions in the Verlet integrator, the global error in velocity is also
.
In molecular dynamics simulations, the global error is typically far more important than the local error, and the Verlet integrator is therefore known as a second-order integrator.
The
variables are the positions of the points i at time t, the 
are the unconstrained positions (i.e. the point positions before applying the constraints) of the points i at time t, the d variables are temporary (they are added for optimization as the results of their expressions are needed multiple times), and r is the distance that is supposed to be between the two points. Currently this is in one dimension; however, it is easily expanded to two or three. Simply find the delta (first equation) of each dimension, and then add the deltas squared to the inside of the square root of the second equation (Pythagorean theorem
). Then, duplicate the last two equations for the number of dimensions there are. This is where verlet makes constraints simple - instead of say, applying a velocity to the points that would eventually satisfy the constraint, you can simply position the point where it should be and the verlet integrator takes care of the rest.
Problems, however, arise when multiple constraints position a vertex. One way to solve this is to loop through all the vertices in a simulation in a criss cross manner, so that at every vertex the constraint relaxation of the last vertex is already used to speed up the spread of the information. Either use fine time steps for the simulation, use a fixed number of constraint solving steps per time step, or solve constrains until they are met by a specific deviation.
When approximating the constraints locally to first order this is the same as the Gauss–Seidel method. For small matrices
it is known that LU decomposition
is faster. Large systems can be divided into clusters (for example: each ragdoll
=cluster). Inside clusters the LU method is used, between clusters the Gauss–Seidel method is used. The matrix code can be reused: The dependency of the forces on the positions can be approximated locally to first order, and the verlet integration can be made more implicit.
For big matrices sophisticated solvers (look especially for "The sizes of these
small dense matrices can be tuned to match the sweet spot" in http://crd.lbl.gov/~xiaoye/SuperLU/superlu_ug.pdf) for sparse matrices exist, any self made Verlet integration has to compete with these. The usage of (clusters of) matrices is not generally more precise or stable, but addresses the specific problem, that a force on one vertex of a sheet of cloth should reach any other vertex in a low number of time steps even if a fine grid is used for the cloth http://www.cs.cmu.edu/~baraff/papers/index.html (link needs refinement) and not form a sound wave.
Another way to solve Holonomic constraints
is to use constraint algorithm
s.
The Verlet integration would automatically handle the velocity imparted via the collision in the latter case, however note that this is not guaranteed to do so in a way that is consistent with collision physics
(that is, changes in momentum are not guaranteed to be realistic). Instead of implicitly changing the velocity term, you would need to explicitly control the final velocities of the objects colliding (by changing the recorded position from the previous time step).
The two simplest methods for deciding on a new velocity are perfectly elastic collision
s and inelastic collision
s. A slightly more complicated strategy that offers more control would involve using the coefficient of restitution
.
Where f is a number representing the fraction of the velocity per update that is lost to friction (0-1).
Isaac Newton
Sir Isaac Newton PRS was an English physicist, mathematician, astronomer, natural philosopher, alchemist, and theologian, who has been "considered by many to be the greatest and most influential scientist who ever lived."...
equations of motion
. It is frequently used to calculate trajectories
Trajectory
A trajectory is the path that a moving object follows through space as a function of time. The object might be a projectile or a satellite, for example. It thus includes the meaning of orbit—the path of a planet, an asteroid or a comet as it travels around a central mass...
of particles in molecular dynamics
Molecular dynamics
Molecular dynamics is a computer simulation of physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a period of time, giving a view of the motion of the atoms...
simulations and computer graphics
Computer graphics
Computer graphics are graphics created using computers and, more generally, the representation and manipulation of image data by a computer with help from specialized software and hardware....
. The Verlet integrator offers greater stability
Numerical stability
In the mathematical subfield of numerical analysis, numerical stability is a desirable property of numerical algorithms. The precise definition of stability depends on the context, but it is related to the accuracy of the algorithm....
, as well as other properties that are important in physical systems such as time-reversibility and preservation of the symplectic form on phase space
Symplectic integrator
In mathematics, a symplectic integrator is a numerical integration scheme for a specific group of differential equations related to classical mechanics and symplectic geometry. Symplectic integrators form the subclass of geometric integrators which, by definition, are canonical transformations...
, at no significant additional cost over the simple Euler method. Verlet integration was used by Carl Størmer
Carl Størmer
Fredrik Carl Mülertz Størmer was a Norwegian mathematician and physicist, known both for his work in number theory and for studying the movement of charged particles in the magnetosphere and the formation of aurorae....
to compute the trajectories of particles moving in a magnetic field (hence it is also called Störmer's method) and was popularized in molecular dynamics by French physicist Loup Verlet
Loup Verlet
Loup Verlet is a French physicist who pioneered the computer simulation of molecular dynamics models. In a famous 1967 paper he used what is now known as Verlet integration and the Verlet list Loup Verlet is a French physicist who pioneered the computer simulation of molecular dynamics models. ...
in 1967.
Equations of motion
Newton's equation of motion for conservative physical systems isor individually
where
- t is the time,
-
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
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
Molecular dynamics
Molecular dynamics is a computer simulation of physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a period of time, giving a view of the motion of the atoms...
to the orbit of the planets
N-body problem
The n-body problem is the problem of predicting the motion of a group of celestial objects that interact with each other gravitationally. Solving this problem has been motivated by the need to understand the motion of the Sun, planets and the visible stars...
.
After a transformation to bring the mass to the right side and forgetting the structure of multiple particles, the equation may be simplified to
with some suitable vector valued function A representing the position dependent acceleration. Typically, an initial position
and an initial velocity are also given.Störmer's method
To discretize and numerically solve this initial value problemInitial value problem
In mathematics, in the field of differential equations, an initial value problem is an ordinary differential equation together with a specified value, called the initial condition, of the unknown function at a given point in the domain of the solution...
, a time step
is chosen and the sampling point sequence considered. The task is to construct a sequence of points that closely follow the points 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:
The Störmer method uses this equation to obtain the next position vector from the previous two as
without using the velocity. 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 h of degree three. The local error is quantified by inserting the exact values
into the iteration and computing the Taylor expansions at time of the position vector in different time directions.where
is the position, the velocity, the acceleration and the jerk (third derivative of the position with respect to the time) .Adding these two expansions gives
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,
, whereas in the iteration it is computed at the central iteration point, . 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 as well as the approximative solution can be expressed in explicit formulas. The standard example for this task is the exponential functionExponential function
In mathematics, the exponential function is the function ex, where e is the number such that the function ex is its own derivative. The exponential function is used to model a relationship in which a constant change in the independent variable gives the same proportional change In mathematics,...
.
Consider the linear differential equation
with a constant w. Its exact basis solutions are and .The Störmer method applied to this differential equation leads to a linear recurrence relation
Recurrence relation
In mathematics, a recurrence relation is an equation that recursively defines a sequence, once one or more initial terms are given: each further term of the sequence is defined as a function of the preceding terms....
It can be solved by finding the roots of its characteristic polynomial
. These are.The basis solutions of the linear recurrence are
and . To compare them with the exact solutions, Taylor expansions are computed.The quotient of this series with the one of the exponential
starts with , soFrom there it follows that for the first basis solution the error can be computed as
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, time , computing , one already needs the position vector at time . At first sight this could give problems, because the initial conditions are known only at the initial time . However, from these the acceleration is known, and a suitable approximation for the first time step position can be obtained using the Taylor polynomial of degree two:The error on the first time step calculation then is of order
. This is not considered a problem because on a simulation of over a large amount of timesteps, the error on the first timestep is only a negligible small amount of the total error, which at time is of the order , both for the distance of the position vectors to as for the distance of the divided differences to . Moreover, to obtain this second order global error, the initial error needs to be of at least third order.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 dynamicsMolecular dynamics
Molecular dynamics is a computer simulation of physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a period of time, giving a view of the motion of the atoms...
simulations, because kinetic energy and instantaneous temperatures at time
cannot be calculated for a system until the positions are known at time . This deficiency can either be dealt with using the Velocity Verlet algorithm, or estimating the velocity using the position terms and the mean value theoremMean value theorem
In calculus, the mean value theorem states, roughly, that given an arc of a differentiable curve, there is at least one point on that arc at which the derivative of the curve is equal to the "average" derivative of the arc. Briefly, a suitable infinitesimal element of the arc is parallel to the...
:
Note that this velocity term is a step behind the position term, since this is for the velocity at time
, not , meaning that is an order two approximation to . With the same argument, but halving the time step, is an order two approximation to , with .One can shorten the interval to approximate the velocity at time
at the cost of accuracy:Velocity Verlet
A related, and more commonly used, algorithm is the Velocity Verlet algorithm , 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 first-timestep problem in the Basic Verlet algorithm:It can be shown that the error on the Velocity Verlet is of the same order as the Basic Verlet. Note that the Velocity algorithm is not necessarily more memory consuming, because it's not necessary to keep track of the velocity at every timestep during the simulation. The standard implementation scheme of this algorithm is:
- Calculate:
- Calculate:
- Derive
from the interaction potential using 
- Calculate:
.
Eliminating the half-step velocity, this algorithm may be shortened to
- Calculate:
- Derive
from the interaction potential using 
- Calculate:
.
Note, however, that this algorithm assumes that acceleration
only depends on position , and does not depend on velocity .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
Semi-implicit Euler method
In mathematics, the semi-implicit Euler method, also called symplectic Euler, semi-explicit Euler, Euler–Cromer, and Newton–Størmer–Verlet , is a modification of the Euler method for solving Hamilton's equations, a system of ordinary differential equations that arises in classical mechanics...
. The algorithms are almost identical up to a shifted by half of a timestep in the velocity. This is easily 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
Midpoint method
In numerical analysis, a branch of applied mathematics, the midpoint method is a one-step method for solving the differential equation y' = f, \quad y = y_0...
, of order two. Additionally, if the acceleration indeed results from the forces in a conservative mechanical or Hamiltonian system
Hamiltonian system
In physics and classical mechanics, a Hamiltonian system is a physical system in which forces are momentum invariant. Hamiltonian systems are studied in Hamiltonian mechanics....
, 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 conservered quantities of the system like linear or angular momentum, that are always preserved or nearly preserved in a symplectic integrator
Symplectic integrator
In mathematics, a symplectic integrator is a numerical integration scheme for a specific group of differential equations related to classical mechanics and symplectic geometry. Symplectic integrators form the subclass of geometric integrators which, by definition, are canonical transformations...
.
Error terms
The local error in position of the Verlet integrator is as described above, and the local error in velocity is .The global error in position, in contrast, is
and the global error in velocity is . These can be derived by noting the following:and
Therefore:
Similarly:
Which can be generalized to (it can be shown by induction, but it is given here without proof):
If we consider the global error in position between
and , where , it is clear that:And therefore, the global (cumulative) error over a constant interval of time is given by:
Because the velocity is determined in a non-cumulative way from the positions in the Verlet integrator, the global error in velocity is also
.In molecular dynamics simulations, the global error is typically far more important than the local error, and the Verlet integrator is therefore known as a second-order integrator.
Constraints
The most notable thing that is now easier due to using Verlet integration rather than Eulerian is that constraints between particles are very easy to do. A constraint is a connection between multiple points that limits them in some way, perhaps setting them at a specific distance or keeping them apart, or making sure they are closer than a specific distance. Often physics systems use springs between the points in order to keep them in the locations they are supposed to be. However, using springs of infinite stiffness between two points usually gives the best results coupled with the verlet algorithm. Here's how:The
variables are the positions of the points i at time t, the are the unconstrained positions (i.e. the point positions before applying the constraints) of the points i at time t, the d variables are temporary (they are added for optimization as the results of their expressions are needed multiple times), and r is the distance that is supposed to be between the two points. Currently this is in one dimension; however, it is easily expanded to two or three. Simply find the delta (first equation) of each dimension, and then add the deltas squared to the inside of the square root of the second equation (Pythagorean theoremPythagorean theorem
In mathematics, the Pythagorean theorem or Pythagoras' theorem is a relation in Euclidean geometry among the three sides of a right triangle...
). Then, duplicate the last two equations for the number of dimensions there are. This is where verlet makes constraints simple - instead of say, applying a velocity to the points that would eventually satisfy the constraint, you can simply position the point where it should be and the verlet integrator takes care of the rest.
Problems, however, arise when multiple constraints position a vertex. One way to solve this is to loop through all the vertices in a simulation in a criss cross manner, so that at every vertex the constraint relaxation of the last vertex is already used to speed up the spread of the information. Either use fine time steps for the simulation, use a fixed number of constraint solving steps per time step, or solve constrains until they are met by a specific deviation.
When approximating the constraints locally to first order this is the same as the Gauss–Seidel method. For small matrices
Matrix (mathematics)
In mathematics, a matrix is a rectangular array of numbers, symbols, or expressions. The individual items in a matrix are called its elements or entries. An example of a matrix with six elements isMatrices of the same size can be added or subtracted element by element...
it is known that LU decomposition
LU decomposition
In linear algebra, LU decomposition is a matrix decomposition which writes a matrix as the product of a lower triangular matrix and an upper triangular matrix. The product sometimes includes a permutation matrix as well. This decomposition is used in numerical analysis to solve systems of linear...
is faster. Large systems can be divided into clusters (for example: each ragdoll
Ragdoll physics
In computer physics engines, ragdoll physics is a type of procedural animation that is often used as a replacement for traditional static death animations.-Introduction:Early video games used manually-created animations for characters' death sequences...
=cluster). Inside clusters the LU method is used, between clusters the Gauss–Seidel method is used. The matrix code can be reused: The dependency of the forces on the positions can be approximated locally to first order, and the verlet integration can be made more implicit.
For big matrices sophisticated solvers (look especially for "The sizes of these
small dense matrices can be tuned to match the sweet spot" in http://crd.lbl.gov/~xiaoye/SuperLU/superlu_ug.pdf) for sparse matrices exist, any self made Verlet integration has to compete with these. The usage of (clusters of) matrices is not generally more precise or stable, but addresses the specific problem, that a force on one vertex of a sheet of cloth should reach any other vertex in a low number of time steps even if a fine grid is used for the cloth http://www.cs.cmu.edu/~baraff/papers/index.html (link needs refinement) and not form a sound wave.
Another way to solve Holonomic constraints
Holonomic constraints
In a system of point particles, holonomic constraints are relations between the coordinates and time which can be expressed in the following form:...
is to use constraint algorithm
Constraint algorithm
In mechanics, a constraint algorithm is a method for satisfying constraints for bodies that obey Newton's equations of motion. There are three basic approaches to satisfying such constraints: choosing novel unconstrained coordinates , introducing explicit constraint forces, and minimizing...
s.
Collision reactions
One way of reacting to collisions is to use a penalty-based system which basically applies a set force to a point upon contact. The problem with this is that it is very difficult to choose the force imparted. Use too strong a force and objects will become unstable, too weak and the objects will penetrate each other. Another way is to use projection collision reactions which takes the offending point and attempts to move it the shortest distance possible to move it out of the other object.The Verlet integration would automatically handle the velocity imparted via the collision in the latter case, however note that this is not guaranteed to do so in a way that is consistent with collision physics
Collision
A collision is an isolated event which two or more moving bodies exert forces on each other for a relatively short time.Although the most common colloquial use of the word "collision" refers to accidents in which two or more objects collide, the scientific use of the word "collision" implies...
(that is, changes in momentum are not guaranteed to be realistic). Instead of implicitly changing the velocity term, you would need to explicitly control the final velocities of the objects colliding (by changing the recorded position from the previous time step).
The two simplest methods for deciding on a new velocity are perfectly elastic collision
Elastic collision
An elastic collision is an encounter between two bodies in which the total kinetic energy of the two bodies after the encounter is equal to their total kinetic energy before the encounter...
s and inelastic collision
Inelastic collision
An inelastic collision, in contrast to an elastic collision, is a collision in which kinetic energy is not conserved.In collisions of macroscopic bodies, some kinetic energy is turned into vibrational energy of the atoms, causing a heating effect, and the bodies are deformed.The molecules of a gas...
s. A slightly more complicated strategy that offers more control would involve using the coefficient of restitution
Coefficient of restitution
The coefficient of restitution of two colliding objects is a fractional value representing the ratio of speeds after and before an impact, taken along the line of the impact...
.
Applications
The Verlet equations can also be modified to create a very simple damping effect (for instance, to emulate air friction in computer games):Where f is a number representing the fraction of the velocity per update that is lost to friction (0-1).