adaptint.m

用matlab写的一些数值算法

function  [I, fcnt] = adaptint(f,a,b,tol)
% Integral of function f over interval [a,b]
% tol : desired absolute accuracy.
% fcnt: number of function evaluations

% Version 4.06.2004.  INCBOX

% Initialize.  End and mid point values of integrand
ff = feval(f,[a  (a+b)/2  b]);  fcnt = 3;
% Initial Simpson approximation
I1 = (b-a) * (ff(1) + 4*ff(2) + ff(3)) / 6;
% Recursive computation
[I,fcnt] = adpi(f,a,b,tol,ff,I1,fcnt);

% Auxiliary function
function  [I,fcnt] = adpi(f,a,b,tol,ff,I1,fcnt)
% Check contribution from (sub)interval [a,b]
h = (b-a)/2;  m = (a+b)/2;
% Mid point values in half intervals
fm = feval(f, [(a+m)/2  (m+b)/2]);  fcnt = fcnt + 2;
IL = h*(ff(1) + 4*fm(1) + ff(2))/6;   % Left half interval
IR = h*(ff(2) + 4*fm(2) + ff(3))/6;  % Right half interval
% Refined approximation with extrapolation
I2 = IL + IR;  
I = I2 + (I2 - I1)/15;
% Check accuracy
if  abs(I-I2) > tol    % Refine both subintervals 
  [IL,fcnt] = adpi(f,a,m,tol/2,[ff(1) fm(1) ff(2)],IL,fcnt);
  [IR,fcnt] = adpi(f,m,b,tol/2,[ff(2) fm(2) ff(3)],IR,fcnt);
  I = IL + IR;
end