Beeman's algorithm is a method for numerically integrating ordinary differential equations of order 2, more specifically Newton's equations of motion
\ddotx=A(x)
The formula used to compute the positions at time
t+\Deltat
x(t+\Deltat)
tandt-\Deltat
x(t+\Deltat) =x(t)+v(t)\Deltat +
1 | |
6 |
l(4a(t)-a(t-\Deltat)r)\Deltat2 +O(\Deltat4)
t+\Deltat
tandt+\Deltat
a(t+\Deltat)
\begin{align} x(t+\Deltat) &=x(t)+v(t)\Deltat +
1 | |
6 |
l(a(t+\Deltat)+2a(t)r)\Deltat2+O(\Deltat4);\\ v(t+\Deltat)\Deltat &=x(t+\Deltat)-x(t) +
16 | |
l(2a(t+\Delta |
t)+a(t)r)\Deltat2 +O(\Deltat4); \end{align}
In tests it was found that this corrector step needs to be repeated at most twice. The values on the right are the old values of the last iterations, resulting in the new values on the left.
Using only the predictor formula and the corrector for the velocities one obtains a direct or explicit method which is a variant of the Verlet integration method:
\begin{align} x(t+\Deltat) &=x(t)+v(t)\Deltat +
1 | |
6 |
l(4a(t)-a(t-\Deltat)r)\Deltat2 +O(\Deltat4)\\ v(t+\Deltat) &=v(t) +
16 | |
l(2a(t+\Delta |
t)+5a(t)-a(t-\Deltat)r)\Deltat +O(\Deltat3); \end{align}
This is the variant that is usually understood as Beeman's method.
Beeman also proposed to alternatively replace the velocity update in the last equation by the second order Adams–Moulton method:
v(t+\Deltat) =v(t) +
1 | |
12 |
l(5a(t+\Deltat)+8a(t)-a(t-\Deltat)r)\Deltat +O(\Deltat3)
where
t
\Deltat
x(t)
v(t)
a(t)
x(t)
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+\Deltat
An example is:
x(t+\Deltat)=x(t)+v(t)\Deltat+
2 | |
3 |
a(t)\Deltat2-
1 | |
6 |
a(t-\Deltat)\Deltat2+O(\Deltat4).
The velocities at time
t=t+\Deltat
v(t+\Deltat)~(predicted)=v(t)+
3 | |
2 |
a(t)\Deltat-
1 | |
2 |
a(t-\Deltat)\Deltat+O(\Deltat3).
The accelerations
a(t+\Deltat)
t=t+\Deltat
v(t+\Deltat)~(corrected)=v(t)+
5 | |
12 |
a(t+\Deltat)\Deltat+
2 | |
3 |
a(t)\Deltat-
1 | |
12 |
a(t-\Deltat)\Deltat+O(\Deltat3).
As shown above, the local error term is
O(\Deltat4)
O(\Deltat3)
O(\Deltat3)
O(\Deltat2)
The simulation must keep track of position, velocity, acceleration and previous acceleration vectors per particle (though some clever workarounds 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.