The Bulirsch-Stoer method is an adaptive method which uses Gragg's modified midpoint method to estimate the solution of an initial value problem for various step sizes. The estimates are fit to a "diagonal" rational function or a polynomial as a function of the step size and the limit as the step size tends to zero is taken as the final estimate. If the absolute difference of two successive scaled estimates is less than a user predefined positive number epsilon, then the procedure is assumed to have converged.
This technique can be extremely useful especially in the neighborhood of poles or other singularities.
- int Gragg_Bulirsch_Stoer( double (*f)(double, double), double y0, double *y1,
double x0, double h, double *h_new, double epsilon, double yscale, int rational_extrapolate )
This function uses the GBS method to return the estimate of the solution of the initial value problem, y' = f(x,y); y = y[0] when x = x0, at x0 + h. Upon returning h_new contains the estimated step size so that the final answer is within epsilon * yscale of the actual solution at x = x0 + h, this value of h_new can be used as the initial step size h in the subsequent call to this function. On input, y0 is the value of y(x) at x = x0 and on output y1 is the value of y(x) at x = x0 + h. This function returns a 0 if a solution was found, i.e. if yn and yn+1 are two successive estimates then the solution is assumed to be found if |yn / yscale - yn+1 / yscale| < epsilon. The function returns -1 if the method failed to converge, -2 if an attempt was made to divide by zero, and -3 if yscale is zero.
- The file, bulirsch_stoer.c, contains the version of Gragg_Bulirsch_Stoer( ) written in C.
