Andrew Coppin wrote: [snip]
Fact #1: I don't *know* of any other numerical integration algorithm. (I've heard of RK4, but it's too complicated for me to understand.)
higher-order runge-kutta algorithms typically have high integration accuracy, but may require multiple energy/force evaluations per integration step. also, they are generally not symplectic [1, 2], which is believed to be rather undesirable if you're looking for long-term stability in your simulation. stoermer/verlet/leap-frog is symplectic, and requires only one force evaluation per time-step. it is therefore _the_ default choice for numerical integration of the equations of motion. it is really quite simple [3, 4]. if you're _really_ looking for speed, you'll want to look at r-RESPA, which is a generalization of leapfrog that allows for multiple-timestep integration [5]. if the force evaluation is indeed the limiting step in your computation (which i'm not sure is the case), you'll certainly want to look into that. it basically boils down to splitting your interaction function into a (rugged) short-range part (which is zero beyond some distance) and which is evaluated every timestep, and a (smooth) long-range part, which is evaluated only every so many steps.
Fact #2: I have tried running the simulation with several different, non-comensurate time step values, and it always seems to produce the same output, so I'm reasonably confident there are no integration errors.
getting a good feeling for the maximum allowed timestep size is usually the first step to getting performance out of your program. as a rule of thumb, you'll want more than 10 steps per "period" of your fastest degree of freedom. you can get a feeling for that by finding the maximum second derivative of the potential energy function (which is going to be close to one of your magnets, i suppose. speaking of which -- does the energy diverge near the magnets? i think it does.). to validate the above, you'll want to look at energy conservation -- compute the total (potential and kinetic) energy of your system at the beginning and the end (for a given simulation length in terms of time, not in terms of timesteps) of your simulation, and compare to the starting energy. you'll find that there's a critical step size below which your energy stays constant, and above which it starts to diverge rather quickly. you might want to take that step size (and divide it by half if you're cautious). kind regards, v. [1] http://srv.chim.unifi.it/orac/MAN/node5.html [2] http://mitpress.mit.edu/SICM/book-Z-H-57.html#%_chap_5 [3] http://einstein.drexel.edu/courses/CompPhys/Integrators/leapfrog/ [4] http://shootout.alioth.debian.org/gp4/benchmark.php?test=nbody&lang=all [5] http://cat.inist.fr/?aModele=afficheN&cpsidt=15255925