Numerical Integration - Numerical Methods in Engineering with MATLAB - page 212

romberg

The algorithm forRomberg integrationis implementedinthe function romberg . It
returns the value of the integral and the required number of function evaluations.
Richardson'sextrapolationis performedbythe subfunction richardson .

function [I,numEval] = romberg(func,a,b,tol,kMax)
% Romberg integration.
% USAGE: [I,numEval] = romberg(func,a,b,tol,kMax)
% INPUT:
% func = handle of function being integrated.

% a,b = limits of integration.

% tol = error tolerance (default is 1.0e-8).

% kMax = limit on the number of panel doublings

% (default is 20).

% OUTPUT:
% I = value of the integral.
% numEval = number of function evaluations.

if nargin < 5; kMax = 20; end
if nargin < 4; tol = 1.0e-8; end
r = zeros(kMax);
r(1) = trapezoid(func,a,b,0,1);
rOld = r(1);
fork=2:kMax
r(k) = trapezoid(func,a,b,r(k-1),k);
r = richardson(r,k);
if abs(r(1) - rOld) < tol
numEval=2ˆ(k-1)+1;I=r(1);
return

end
rOld = r(1);

end
error('Failed to converge')

functionr=richardson(r,k)
% Richardson's extrapolation in Eq. (6.14).
forj=k-1:-1:1
c = 4ˆ(k-j); r(j) = (c*r(j+1) - r(j))/(c-1);

end