Beeman's algorithm

Beeman's algorithm is a method for numerically integrating ordinary differential equations, generally position and velocity, which is closely related to Verlet integration. It produces identical positions to Verlet, but is more accurate in velocities. It is most commonly seen in molecular dynamics simulations.

Equation

The formula used to compute the positions at time $t\; +\; Delta\; t$ is:

:$r(t+Delta\; t)\; =\; r(t)\; +\; v(t)\; Delta\; t\; +\; frac\{2\}\{3\}a(t)\; Delta\; t^2\; -\; frac\{1\}\{6\}\; a(t\; -\; Delta\; t)\; Delta\; t^2\; +\; O(\; Delta\; t^4)$

and this is the formula used to update the velocities:

:$v(t\; +\; Delta\; t)\; =\; v(t)\; +\; frac\{1\}\{3\}a(t\; +\; Delta\; t)\; Delta\; t\; +\; frac\{5\}\{6\}a(t)\; Delta\; t\; -\; frac\{1\}\{6\}a(t\; -\; Delta\; t)\; Delta\; t\; +\; O(Delta\; t^3)$

where

*$t$ is present time (ie: independent variable)
*$Delta\; t$ is the time step size
*$r(t)$ is the position at time t
*$v(t)$ is the velocity at time t
*$a(t)$ is the acceleration at time t
*the last term is the error term, using the big O notation

Predictor-Corrector Modifications

In systems where the forces are a function of velocity in addition to position, the above equations need to be modified into a predictor-corrector form whereby the velocities at time $t\; +\; Delta\; t$ are predicted and the forces calculated, before producing a corrected form of the velocities.

An example is:

:$r(t+Delta\; t)\; =\; r(t)\; +\; v(t)\; Delta\; t\; +\; frac\{2\}\{3\}a(t)\; Delta\; t^2\; -\; frac\{1\}\{6\}\; a(t\; -\; Delta\; t)\; Delta\; t^2\; +\; O(\; Delta\; t^4).$

The velocities at time t =$t\; +\; Delta\; t$ are then calculated from the positions.

:$v(t\; +\; Delta\; t)\; (predicted)\; =\; v(t)\; +\; frac\{3\}\{2\}a(t)\; Delta\; t\; -\; frac\{1\}\{2\}a(t\; -\; Delta\; t)\; Delta\; t$

The accelerations at time t =$t\; +\; Delta\; t$ are then calculated from the positions and predicted velocities.

:$v(t\; +\; Delta\; t)\; (corrected)\; =\; v(t)\; +\; frac\{1\}\{3\}a(t\; +\; Delta\; t)\; Delta\; t\; +\; frac\{5\}\{6\}a(t)\; Delta\; t\; -\; frac\{1\}\{6\}a(t\; -\; Delta\; t)\; Delta\; t$

Error term

As shown above, the error term is $O(Delta\; t^4)$ for position and $O(Delta\; t^3)$ velocity. In comparison, Verlet is $O(Delta\; t^4)$ for position and $O(Delta\; t^2)$ for velocity. In exchange for greater accuracy, Beeman's algorithm is moderately computationally more expensive.

Memory Requirements

The simulation must keep track of position, velocity, acceleration and previous acceleration vectors per particle (though some clever work-arounds for storing the previous acceleration vector are possible), keeping its memory requirements on par with velocity Verlet and slightly more expensive than the original Verlet method.

References

* [http://www.earth.ox.ac.uk/~keithr/moldy-manual/node9.html Integration Algorithms]
* [http://mse-jl1.eng.ohio-state.edu/Archive/Papers/05/Li05-2.8.pdf Basic Molecular Dynamics]
* [http://www.mse.ufl.edu/~ssinn/Backgrounds/back04_integ.html Integrator: Predictor-Corrector Algorithm]