romberg.m

matlab算法集 matlab算法集

function [y,r] = romberg (a,b,tol,m,f,dm)
%-----------------------------------------------------------------------
% Usage:       [y,r] = romberg (a,b,tol,m,f,dm);
%
% Description: Use the Romberg integration mathod to numerically 
%              integrate the function f(x) over the interval [a,b].
%
% Inputs:      a   = lower limit of integration
%              b   = upper limit of integration
%              tol = error tolerance for truncation error (tol >= 0)
%              m   = maximum extrapolation level (2 to 10) 
%              f   = string containing name of user-supplied function 
%                    to be integrated. The form of f is:
%
%                    function y = f(x)
%
%                    When f is called it must return the value f(x).
%              dm  = optional display mode.  If present, intermediate
%                    results are displayed.
%
% Outputs:     y = estimate of integral
%              r = extrapolation level used.  If r < m, then the 
%                  following termination criterion was satisfied 
%                  where E is the estimated truncation error:
%
%                  |E| < tol
%-----------------------------------------------------------------------

% Initialize

   tol = args (tol,0,tol,3,'romberg');
   m   = args (m,2,10,4,'romberg');

   display = nargin > 5;
   k = 1;
   n = 1;
   R = zeros (m,m);
   h = b - a;
   E = tol + 1;
   R(1,1) = h*(feval(f,a) + feval(f,b))/2;
   if display
      fprintf ('\n  1%11.7f',R(1,1));
   end
      
% Compute estimates 
   
   if m > 1
      while (abs(E) >= tol) & (k < m) 
         k = k + 1;
         n = 2*n;
         h = h/2;
         R(k,1) = (feval(f,a) + feval(f,b))/2;
         for j = 1 : n-1
            R(k,1) = R(k,1) + feval(f,a+j*h);
         end 	
         R(k,1) = h*R(k,1);  	
         for j = 2 : k
            x =  4^(j-1);
            R(k,j) = (x*R(k,j-1) - R(k-1,j-1))/(x - 1);
         end
         E = (R(k,k-1) - R(k-1,k-1))/(4^(k-1) - 1);
         if display
            fprintf ('\n %2i',k);
            for j = 1 : k
               fprintf ('%11.7f',R(k,j));
            end
         end
      end
   end

% Finalize 

   r = k;
   y = R(k,k);
%-----------------------------------------------------------------------