romberg.m

Integraton routines in matlab

function [int, T, nfeval] = romberg(fun, a, b, tol)
%
%   function [int, T, nfeval] = romberg(f, a, b, tol)
%
%   Compute an approximation of the integral of f over
%   [a,b]  using Romberg integration.
%
%   Input
%   fun    name of the matlab function that evaluates the integrand
%   a      lower limit of integration
%   b      upper limit of integration
%   tol    tolerance. Romberg integration is stopped if
%          abs(T(n,n-1)-T(n,n)) <= tol*abs(T(n,n))
%   
%   Output
%   int    approximation of the integral (=T(n,n))
%   T      Romberg table
%   nfeval number of function values required
%
%   
%   Version April 9, 2002
%   Matthias Heinkenschloss




n_max  = 100;   % maximum number of steps 

h      = b-a;
x      = (a:h/2:b);
fx     = feval( fun, x);
nfeval = 3;
T(1,1) = 0.5*h*(fx(1)+fx(3));
h      = h/2;
T(2,1) = 0.5*T(1,1) + h*fx(2);
T(2,2) = T(2,1) + (T(2,1)-T(1,1))/3;
n      = 2;

while( abs(T(n,n-1)-T(n,n)) > tol*abs(T(n,n)) & n < n_max )
    n      = n+1;
    h      = h/2;
    x      = (a+h:2*h:b-h);
    fx     = feval( fun, x);
    nfeval = nfeval + length(fx);
    T(n,1) = 0.5*T(n-1,1) + h*sum(fx);
    for j = 2:n
        T(n,j) = T(n,j-1) + (T(n,j-1)-T(n-1,j-1))/(4^(j-1)-1);
    end
end

int = T(n,n);

if( n == n_max )
    disp(' WARNING: Romberg:')
    disp(' Maximum number of steps in Romberg integration reached.')
    disp(' int is not within tolerance')
end