simpson.m

matlab算法集 matlab算法集

function y = simpson (a,b,n,f)
%-----------------------------------------------------------------------
% Usage:       y = simpson (a,b,n,f);
%
% Description: Use Simpson's rule to numerically integrate the
%              function f(x) over the interval [a,b] using n panels.
%
% Inputs:      a = lower limit of integration
%              b = upper limit of integration
%              n = number of panels (n >= 2 and even)
%              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).
%
% Outputs:     y = estimate of integral
%
% Note:        The truncation error is proportional to n^(-4). 
%-----------------------------------------------------------------------

% Initialize

   if (n < 2) | mod(n,2)
      disp ('Argumant 4 of simpson must be even and >= 2.')
      return
   end
   h = (b - a)/n;
   y = feval(f,a) + feval(f,b);

% Compute integral

   for i = 1 : n-1
      if mod(i,2)
         y = y + 4*feval(f,a + i*h);
      else
         y = y + 2*feval(f,a + i*h); 
      end
   end  
   y = h*y/3;
%-----------------------------------------------------------------------