[Haskell-cafe] Re: (Chaos) [An interesting toy]

[Vincent Kraeutler]

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


-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 249 bytes
Desc: OpenPGP digital signature
Url : http://www.haskell.org/pipermail/haskell-cafe/attachments/20070507/25360bf7/signature.bin