romberg.m

用matlab写的一些数值算法

function  [I, info, R] = romberg(f,a,b,tol,p)
% Romberg's method for the integral of  f  over  [a,b].
% Assume that  T(h) = T(0) + a1*h^p(1) + ... + aq*h^p(q) ,
% where  q = length(p) .
% Default:  p = [2 4 ... 14]
% Stop when two adjacent values in the same column differ 
% less than  tol .  
% info = [estimated error, m, r], where m+1 and r-1 are the 
%    number of function evaluations and extrapolations.

% Version 11.12.2003.  INCBOX

% Initialize the scheme
if  nargin < 5
  p = 2 : 2 : 14;
end
q = length(p);  R = zeros(q+2,q+1);
h = b - a;  m = 1;
R(1,1) = h*(feval(f,a) + feval(f,b))/2;
mdif = inf;   % initialize min difference
for  i = 2 : q+2
  % First element in next row
  h = h/2;  m = 2*m;
  s = 0;
  for  j = 1 : 2 : m
    s = s + feval(f, a+j*h);
  end
  R(i,1) = R(i-1,1)/2 + h*s;
  % Extrapolate
  jmax = min(i-1,q);
  for  j = 1 : jmax
    R(i,j+1) = R(i,j) + (R(i,j) - R(i-1,j))/(2^p(j) - 1);
  end
  % Check accuracy
  [md, jb] = min(abs(R(i,1:jmax) - R(i-1,1:jmax)));
  if  md < mdif  % better result
    info = [md m jb];  I = R(i,jb);
    if  md <= tol
      R = R(1:i,1:jmax+1);        % return active part of R
      return
    end
  end
end