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)