simpcomp.m

Integraton routines in matlab

function [int, S, nfeval] = simpcomp(fun, a, b, tol)
%
%   function [int, nfeval] = simpcomp(f, a, b, tol)
%
%   Compute an approximation of the integral of f over
%   [a,b]  using composite Simpson rule.
%
%   Input
%   fun    name of the matlab function that evaluates the integrand
%   a      lower limit of integration
%   b      upper limit of integration
%   tol    tolerance. Composite Trazepzodal rule is stopped if
%          abs(S(n-1)-S(n)) <= tol*abs(S(n))
%   
%   Output
%   int    approximation of the integral  (=S(n))
%   S      S(n)  composite trapezoidal rule with step size h = (b-a)/2^(n-1)
%   nfeval number of function values required
%
%   
%   Version April 9, 2002
%   Matthias Heinkenschloss


%   The following notation is used:
%   S(n) = value of the composite Simpson rule with step size 
%          h = (b-a)/2^(n-1)
%   T    = value of the composite Trapezoidal rule with step size 
%          h = (b-a)/2^(n-1)
%   M    = h * ( f(a+h/2) + f(a+3h/2) + ... + f(b-3h/2) + f(b-h/2) )
%          weighted sum of funcation values at the mid-points.
% 
%   Note that S(n) = (1/3)*T + (2/3)*M 
%   and that value of the composite Trapezoidal rule with step size h/2
%   is      T = 0.5*(T + M)
%   
n_max  = 100;   % maximum number of steps 

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


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

int = S(n);

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