I'm using a simple version of Verlet integration for a particle–particle interaction system with collisions. At the end of each iteration, I integrate like this:

verletX = ( 1.92*x[i] - 0.92*x_old[i] ) + fx[i]*0.002;
verletY = ( 1.92*y[i] - 0.92*y_old[i] ) + fy[i]*0.002;
x_old[i] = x[i];
y_old[i] = y[i];
x[i] = verletX;
y[i] = verletY;

Can I increase its order of error (to have more precision for decreasing timesteps) using even more old data of X and Y? For now it should be less precise than a RK4, but I’m sure it’s much better than a simple Euler integration. For example can I use the data from the last eight iterations to reach RK4 levels?

The reason I’m not trying RK4 is that the calculation of forces takes too much time. However, having a copy of the old data and updating it does not. Also I don’t need exactness; I just need the system stability such that particles stay in fixed points rather than having explosive close range leaps.


Use multistep methods in this case. See the Adams-Bashforth method if it's nonstiff, or Adams-Bashforth-Moulton methods, or if its stiff BDF methods. These use past timepoints like you want in order to increase the order of accuracy.

I would highly recommend checking out something like Sundials which has variable timestep plus variable order methods. This helps a lot with the stability of these methods and lets them get a lot more efficiency out.

  • $\begingroup$ what does h mean in all these equations? step size in time or step size in number like 1,2,3,4? if it is time, is it sum of steps or just one step? $\endgroup$ – huseyin tugrul buyukisik Mar 19 '17 at 9:40
  • $\begingroup$ $h$ is the stepsize in time. It's assumed constant for the entire integration in the way its usually written. Since you're using past values which are assumed equally spaced, the derivation is much easier. That's why adaptive methods are orders of magnitude more difficult for these kinds of algorithms. $\endgroup$ – Chris Rackauckas Mar 19 '17 at 15:19
  • $\begingroup$ I tried 5 step version with trivial h values smaller than 1 and it worked $\endgroup$ – huseyin tugrul buyukisik Mar 19 '17 at 15:21
  • $\begingroup$ Note that Verlet is symplectic (preserves energy, avoids long-term drift) whereas I don't believe the same holds for ABM. Sympleticity is why they use Verlet for billion-particle simulations. $\endgroup$ – Richard Zhang Mar 25 '17 at 19:29
  • $\begingroup$ Indeed ABM is symplectic. But there are multi-step conjugate-symplectic methods. They are very similar to what I suggested, but use different coefficients. If the OP needs to preserve symplectic-ness of integration, then that's probably the way to go. $\endgroup$ – Chris Rackauckas Mar 25 '17 at 19:36

Your Answer

