trapcomp.m

来自「Integraton routines in matlab」· M 代码 · 共 54 行

function [int, T, nfeval] = trapcomp(fun, a, b, tol)
%
%   function [int, nfeval] = trapcomp(f, a, b, tol)
%
%   Compute an approximation of the integral of f over
%   [a,b]  using composite Trazepzodal 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(T(n-1)-T(n)) <= tol*abs(T(n))
%   
%   Output
%   int    approximation of the integral  (=T(n))
%   T      T(n)  composite trapezoidal rule with step size h = (b-a)/n
%   nfeval number of function values required
%
%   
%   Version April 9, 2002
%   Matthias Heinkenschloss




n_max  = 100;   % maximum number of steps 

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

while( abs(T(n-1)-T(n)) > tol*abs(T(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)   = 0.5*T(n-1) + h*sum(fx);
end

int = T(n);

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