Types of N-body ODE Integration Methods

June 9, 1996

39

The large step size H can be so large that the results from using even the largest

number of steps (8 in our example) does not have to be very accurate. When all the trial

movements are combined to form the polynomial approximation, the end result will be an

extremely accurate method. So, even though each movement of step size H requires a

large number of evaluations of the force function

*f()*

, it is still much more efficient that

most other methods. It can also be made to be much more accurate than almost any other

method because the final result does not suffer from the rounding errors that would be

created by using incredibly tiny steps in other methods. For this reason, the Gragg-Bulir-

sch-Stoer method is clearly the best method when extremely accurate calculations must be

done.

There several pieces to making the Gragg-Bulirsch-Stoer method work:

First, Richardson's idea that the final answer can be thought of as a function of the

step size and this function can be used to extrapolate to what the answer would be if step

size was zero.

Secondly, Gragg figured out that the mid-point method mentioned above, slightly

reworked, has error terms that has only even powers of

*h*

. By careful manipulation when

creating the polynomial approximation, it is possible to cancel out two orders of the error

term at a time.

Thirdly, Bulirsch and Stoer figured out that using rational approximation instead of

Gragg's polynomial approximation allowed you to break the limits of convergence of the

taylor series around singularities and poles in the complex plain. They also figured out a

good method of decreasing the step size (the Bulirsch sequence).

Lastly, Marciniak's method also implements a method of "Discrete Mechanics".

That is, this particular implementation tries to keep the constants of motion constant at the

expense of keeping the positions and velocities close to the exact answer. Since it is often

impossible to tell what the exact answer should be, having the system at least behave well

in other respects can be very useful.

Combining all these things together yields a very accurate, very fast and

*very*

com-

plicated method. It is so complicated in fact, that this is the only method that I didn't write

the code for myself. Instead I used the example code out of Marciniak's book, fixed a few

bugs and made a few tweaks. The only problem is that Marciniak only implemented

Gragg's method of polynomial extrapolation instead of using rational extrapolation. When

XStar uses this method and two stars pass close together, the combined use of polynomial

approximation and the discrete mechanics can give spectacularly wrong results. Try `xstar

-m gpemce8 -a .25' sometime...

Sources: (14:582-8,12:121-4,129-30)

This document is best viewed as n-body.pdf because the translation to html was buggy.