Bulirsch Stoer integrator

For our comparative evaluation of SAAM s Chu Berman (Qm and Berman, 1974) integrator (also referred to as model code 10 in SAAM and CONSAM) we will be considering as alternate integrators the Runge Kutta 4/2 integrator (Press et ai, 1987), the Bulirsch Stoer integrator (Press et al., 1987), and Petzold s DASSL integrator (Petzold, 1983). 

The Bulirsch Stoer integrator performed with consistent accuracy but again was too slow to compete with Petzold s scheme as an alternative to the CBCCDS integrator. 

When the quadrature of eq 2 cannot be performed analytically the integration should be carried out numerically by robust routines such as the Runge-Kutta, Adams-Moulton predictor-corrector or Bulirsch-Stoer methods with step size and error control [53, 55, 56], These routines can also be found in computer codings at Netlib and in standard books on computer codes [53]. 

The quasiclassical trajectory method was used to study this system, and the variable step size modified Bulirsch-Stoer algorithm was specially developed for recombination problems such as this one. Comparisons were made with the fourth order Adams-Bashforth-Moulton predictor-corrector algorithm, and the modified Bulirsch-Stoer method was always more efficient, with the relative efficiency of the Bulirsch-Stoer method increasing as the desired accuracy increased. We measure the accuracy by computing the rms relative difference between the initial coordinates and momenta and their back-integrated values. For example, for a rms relative difference of 0.01, the ratio of the CPU times for the two methods was 1.6, for a rms relative difference of 0.001 it was 2.0, and for a rms relative difference of 10 it was 3.3. Another advantage of the variable step size method is that the errors in individual trajectories are more similar, e.g. a test run of ten trajectories yielded rms errors which differed by a factor of 53 when using the modified Bulirsch-Stoer... 

In spite of the modifications introduced by Bulirsch and Stoer, the method derived from the extended trapezoid rule suffers from several shortcomings that make it less efficient than other methods described later, when implemented in a general integration program. 

---