Initial Value Problems - Numerical Methods in Engineering with MATLAB - page 288

x=x+h;
y2 = y0 + 2.0*h*feval(dEqs,x,y1);
y0=y1;
y1=y2;

end
y = 0.5*(y1 + y0 + h*feval(dEqs,x,y2));

Bulirsch-Stoer Algorithm

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-
rithm shines.

bulStoer

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)
% INPUT:
% dEqs = handle of function that returns the first-order
% differential equations F(x,y) = [dy1/dx,dy2/dx,...].

19

Stoer,J., and Bulirsch, R., Introduction to Numerical Analysis , Springer, 1980.