Verlet integration

From Free net encyclopedia

(Redirected from Verlet)

Verlet integration is a method for calculating the trajectories of particles in molecular dynamics simulations. The verlet integrator offers greater stability than the much simpler Euler integration methods, as well as other properties that are important in physical systems such as time-reversibility and area preserving properties. Verlet integrators are also sometimes used in the physics engines in computer games. The method was developed by French physicist Loup Verlet in 1967.

1 The algorithm
2 Derivation
3 Error Terms
4 Constraints
5 Collision reactions
6 See also
7 External links

[edit]

The algorithm

At first it may seem natural to simply calculate trajectories using Euler integration, which is defined by

<math>x(t_0 + \Delta t) = x(t_0) + v(t_0) \Delta t</math>

<math>v(t_0 + \Delta t) = v(t_0) + a(t_0) \Delta t</math>

where <math>t_0</math> is the current time and <math>\Delta t</math> is the time step. However, this kind of integration suffers from many problems, as discussed at Euler integration. The Verlet algorithm reduces the level of errors introduced into the integration by calculating the position at the next time step from the positions at the previous and current time steps, without using the velocity.

<math>x(t_0 + \Delta t) = 2x(t_0) - x(t_0 - \Delta t) + a\Delta t^2</math>

The velocity at each time step is then not calculated until the next time step.

<math>v(t_0) = \frac{x(t_0 + \Delta t) - x(t_0 - \Delta t)}{2\Delta t}</math>

This can create technical challenges in molecular dynamics simulations, because kinetic energy and instantaneous temperatures at time <math>t_0</math> cannot be calculated for a system until the positions are known at time <math>t_0 + \Delta t</math>.

The Verlet equations can also be modified to create a very simple damping effect (for instance, to emulate air friction in computer games):

<math>x(t_0 + \Delta t) = (2-f) x(t_0) -(1-f) x(t_0 - \Delta t) + at^2.</math>

Where f is a number representing the percentage of velocity per update that is lost to friction (0-1).

[edit]

Derivation

The Verlet integrator can be derived from the Taylor expansion of a trajectory. Let <math>x(t)</math> be the trajectory of a particle at time <math>t</math>. The Taylor expansion around time <math>t_0</math> then gives:

<math>x(t_0 + \Delta t) = x(t_0) + \Delta t x'(t_0) + \frac{1}{2} \Delta t^2 x(t_0) + \frac{1}{6} \Delta t^3 x'(t_0) + O(\Delta t^4) </math>

and

<math>x(t_0 - \Delta t) = x(t_0) - \Delta t x'(t_0) + \frac{1}{2} \Delta t^2 x(t_0) - \frac{1}{6} \Delta t^3 x'(t_0) + O(\Delta t^4) </math>

Adding these together gives:

<math>x(t_0 + \Delta t) + x(t_0 - \Delta t) = 2x(t_0) + \Delta t^2 x(t_0) + O(\Delta t^4) \, </math>

Arranging into the more conventional format:

<math>x(t_0 + \Delta t) = 2x(t_0) - x(t_0 - \Delta t) + \Delta t^2 x(t_0) + O(\Delta t^4) \, </math>

where the term <math>O(\Delta t^4)</math> represents fourth-order and higher terms in the Taylor expansion.

This offers the clear advantage that the third-order term from the Taylor expansion cancels out, thus making the Verlet integrator an order more accurate than integration by simple Taylor expansion alone.

Furthermore, by subtracting the two original Taylors series, we obtain the velocity as:

<math>v(t_0) = \frac{x(t_0 + \Delta t) - x(t_0 - \Delta t)}{2\Delta t} + O(\Delta t^2). </math>

The velocity calculation is thus two orders less accurate than the trajectory calculation. However, this error does not accumulate because the velocity is recalculated from the more accurate trajectories at every timestep.

[edit]

Error Terms

The local error in position of the Verlet integrator is <math>O(\Delta t^4)</math> as described above, and the local error in velocity is <math>O(\Delta t^2)</math>.

The global error in position, in contrast, is <math>O(\Delta t^2)</math> and the global error in velocity is <math>O(\Delta t^2)</math>. These can be derived by noting the following:

<math>\mathit{error}(x(t_0 + \Delta t)) = O(\Delta t^4) \, </math>

and

<math>x(t_0 + 2\Delta t) = 2x(t_0 + \Delta t) - x(t_0) + \Delta t^2 x(t_0 + \Delta t) + O(\Delta t^4) \, </math>

Therefore:

<math>\mathit{error}(x(t_0 + 2\Delta t)) = 2\mathit{error}(x(t_0)) + O(\Delta t^4) = 3O(\Delta t^4) \, </math>

Similarly:

<math>\mathit{error}(x(t_0 + 3\Delta t)) = 6O(\Delta t^4) \, </math>

<math>\mathit{error}(x(t_0 + 4\Delta t)) = 10O(\Delta t^4) \, </math>

<math>\mathit{error}(x(t_0 + 5\Delta t)) = 15O(\Delta t^4) \, </math>

Which can be generalized to (it can be shown by induction, but it is given here without proof):

<math>\mathit{error}(x(t_0 + n\Delta t)) = \frac{n(n+1)}{2}O(\Delta t^4) \, </math>

If we consider the global error in position between x(t) and x(t + T), where <math>T = n\Delta t</math>, it is clear that:

<math>\mathit{error}(x(t_0 + T)) = (\frac{T^2}{2\Delta t^2} + \frac{T}{2\Delta t}) O(\Delta t^4) \, </math>

And therefore, the global (cumulative) error over a constant interval of time is given by:

<math>\mathit{error}(x(t_0 + T)) = O(\Delta t^2) \, </math>

Because the velocity is determined in a non-cumulative way from the positions in the Verlet integrator, the global error in velocity is also <math>O(\Delta t^2)</math>.

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.

[edit]

Constraints

The most notable thing that is now easier due to using Verlet integration rather than Newtonian 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 times 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 tension between two points usually gives the best results coupled with the verlet algorithm. Here's how:

The x variables are the positions of the points, 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 3. 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.

Here's what an optimized 2D verlet constraint would look like in C:

dx = x2-x1;
dy = y2-y1;
d1 = sqrt(dx*dx+dy*dy);
d2 = 0.5*(d1-r)/d1;
dx = dx*d2;
dy = dy*d2;
x1 += dx;
x2 -= dx;
y1 += dy;
y2 -= dy;

3D:

dx = x2-x1;
dy = y2-y1;
dz = z2-z1;
d1 = sqrt(dx*dx+dy*dy+dz*dz);
d2 = 0.5*(d1-r)/d1;
dx = dx*d2;
dy = dy*d2;
dz = dz*d2;
x1 += dx;
x2 -= dx;
y1 += dy;
y2 -= dy;
z1 += dz;
z2 -= dz;

Problems, however, arise when multiple constraints position a point. It is extremely hard not to violate at least one. One way to solve this is to loop through all the constraints multiple times, eventually coming to a set of positions that satisfies the constraints enough to be accurate, or have it simply loop a set number of times.

[edit]

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. Simple as that. The verlet integration takes care of the velocity imparted because of the collision. For more information on how to do this, click Advanced Character Physics by Thomas Jakobsen.

[edit]