Graphics Programs Reference
y2 = y0 + 2.0*h*feval(dEqs,x,y1);
y = 0.5*(y1 + y0 + h*feval(dEqs,x,y2));
Whenused onits own, the module midpoint has a major shortcoming: the solution
at points between the initial and final values of x cannot be refinedbyRichardson
extrapolation,sothat y is usable onlyat the last point. This deficiencyis rectifiedin
the Bulirsch-Stoermethod. The fundamental ideabehind themethodissimple: apply
the midpoint methodinapiecewise fashion. That is, advance the solutioninstages of
length H , using the midpoint methodwith Richardson extrapolation to perform the
integrationineach stage. The valueof H can bequite large, since the precision of the
result is determined mainlybythe step length h in the midpoint method, not by H
The original Bulirsch and Stoer technique 19 is a complexprocedurethat incorpo-
rates many refinements missing in our algorithm. However, the function bulStoer
givenbelowretains the essential ideasofBulirsch and Stoer.
What are the relative merits of adaptive Runge-Kutta and Bulirsch-Stoermeth-
ods? The Runge-Kuttamethodismore robust, having higher tolerance fornonsmooth
functions andstiffproblems. Inmost applicationswherehighprecisionisnot required,
it also tendstobe moreefficient. However,this is not the case in the computation of
high-accuracy solutions involving smooth functions, where the Bulirsch-Stoeralgo-
Thisfunction contains a simplifiedalgorithm for the Bulirsch-Stoermethod.
function [xSol,ySol] = bulStoer(dEqs,x,y,xStop,H,tol)
% Simplified Bulirsch-Stoer method for integration of y' = F(x,y).
% USAGE: [xSol,ySol] = bulStoer(dEqs,x,y,xStop,H,tol)
% dEqs = handle of function that returns the first-order
differential equations F(x,y) = [dy1/dx,dy2/dx,...].
Stoer,J., and Bulirsch, R., Introduction to Numerical Analysis , Springer, 1980.