Questions tagged [verlet-integration]

The problem to solve

Verlet integration is a method for the numerical solution of differential equations

x''(t) = a( x(t) ),    x(t₀)=x₀, x'(t₀)=v₀

that simulate mechanical systems and other conservative or Hamiltonian systems. The objects of the system are assumed to be enumerated and their positions assembled into the vector x(t), their velocities assembled in the vector valued function v(t)=x'(t). a(t) consists of the accelerations of the single objects corresponding to the force acting on each object, divided by its mass.

Verlet integration is a symplectic integrator and will only work as expected on conservative systems, which, for example, excludes the consideration of friction. More precisely this means that the system is required do possess an energy function

E(x,v) = 0.5 * v^T*M*v + P(x)

where the first term is the kinetic energy, M is the (symmetric, and often diagonal) mass matrix and P(x) the potential energy. Then the force acting on the objects of the system is minus the gradient of the potential energy and thus the acceleration is

a(x) = - M^-1 * grad P(x).

Usually the dynamically system is encoded via evaluation of the function a(x), such a procedure is named eval_a(x) in this article.

The solution obtained

Verlet integration, like any numerical integration of ODE, produces sequences

(x_n), (v_n), _{(n in IN)}

of sample points in the phase space of positions and velocities corresponding to a prescribed sequence (t_n) of time instants, so that x_n and v_n approximate the exact solutions x(t_n) and v(t_n)=x'(t_n).

In the formulation of the method the velocities are of secondary importance, the basic version even leaves them out.

The numerical method(s)

Verlet integration comes in three main flavors, all three for a constant time step dt, so that t_n=t₀+n*dt. All three compute the same position sequence, basic Verlet computes only this. The other two also compute velocities, while the leapfrog variant has the least operations, for exact velocities the velocity Verlet method should be used.

• basic Verlet: before each step n->n+1, the variables t, x, x_old are t_n, x_n, x_n-1, initialization returns with the state for n=1.

  verlet_init:
      x_old = x0
      a = eval_a( x0 );
      x = x0 + v0*dt + 0.5*a*dt*dt;   t=t0+dt;
  
  verlet_step:
      a = eval_a(x);
      x_new = 2*x-x_old + a*dt*dt;
      x_old = x;
      x = x_new; t+=dt;
      do_collisions(t,x,x_old);
  

• velocity verlet: before each step n->n+1, the variables t, x, v are t_n, x_n, v_n, initialization returns with the state for n=0.

  verlet_init:
      a = eval_a( x0 );
      v = v0;
      x = x0; t = t0;
  
  verlet_step:
      v += a*0.5*dt;
      x += v*dt; t += dt;
      do_collisions(t,x,v,dt);
      a = eval_a(x);
      v += a*0.5*dt;
      do_statistics(t,x,v); 
  

• leapfrog verlet: before each step n->n+1, the variables t, x, v are t_n, x_n, v_n+0.5, where t_n+0.5=t_n+0.5*dt, initialization returns with the state for n=0.

  verlet_init:
      a = eval_a( x0 );
      v = v0 + 0.5*a*dt;
      x = x0;     t=t0;
  
  verlet_step:
      x += v*dt; t += dt;
      do_collisions(t,x,v,dt);
      a = eval_a(x);
      v += a*dt;
  

The sequence of the computations must be preserved, however the lines inside `verlet_step' may be rotated. Then the initialization needs to be adapted so that the unrolled loop stays the same.

Numerical errors

The approximation error at time t behaves as O(e^Lt*dt²), it is a second order method.

The Verlet integration method belongs to the family of symplectic integration methods. These preserve the energy E and other first integrals of the system almost perfectly. In the Verlet method the energy oscillates with time constant amplitude O(dt²) around the initial energy level.