deriv.m

matlab算法集 matlab算法集

function dy = deriv (a,h,n,m,e,f)
%----------------------------------------------------------------------
% Usage:       dy = deriv (a,h,n,m,e,f)
%
% Description: Numerically estimate the first derivative of f(x) at
%              x = a.
%
% Inputs:      a = evaluation point
%              h = step size (h > 0)
%              n = number of points 
%              m = type of approximation:
%
%                    m      Approximation                n
%                   ----------------------------------------
%                   -1   n-point backward difference   2,3,4
%                    0   n-point central difference    3,5
%                    1   n-point forward difference    2,3,4
%                   ----------------------------------------
%
%              e  = extrapolation parameter.  If e is nonzero,
%                   extrapolation from step sizes 2h and h is used.
%              f  = string containing name of user-supplied function 
%                   to be differentiated.  The form of f is:
%
%                      function y = f(x)
%
%                   When f is called with scalar input x, it must 
%                   return the scalar y = f(x). 
%
% Outputs:     dy =  df(x)/dx at x = a.
%----------------------------------------------------------------------

% Initialize

   h = args (h,eps,h,2,'deriv');

% Extrapolation 

   if e 
      x = 2^(n-1);	
      dy = (x*deriv(a,h,n,m,0,f) - deriv(a,2*h,n,m,0,f))/(x-1);
      return
   elseif m < 0 

% Backward difference 

      if (n < 2) | (n > 4)
         n = args (n,2,4,3,'deriv');
      end
      switch n 
         case 2,
            dy = feval(f,a) - feval(f,a-h);
         case 3,
            dy = (3*feva(f,a) - 4*feval(f,a-h) + feval(f,a-2*h))/2;
         case 4,
            dy = (11*feval(f,a) - 18*feval(f,a-h) + 9*feval(f,a-2*h) - ...
                 2*feval(f,a-3*h))/6;
      end   

   elseif m == 0 
    
% Central difference 

      if (n ~= 3) & (n ~= 5)
         disp ('\nArgument 3 to function deriv was out of range.\n');
      end
      switch n 
         case 3,
            dy = (feval(f,a+h) - feval(f,a-h))/2;
         case 5,
            dy = (feval(f,a-2*h) - 8*feval(f,a-h) + 8*feval(f,a+h) - ...
                 feval(f,a+2*h))/12;
      end

   else
       
% Forward difference 

      if (n < 2) | (n > 4)
         n = args (n,2,4,3,'deriv');
      end
      switch n 
         case 2,
            dy = -feval(f,a) + feval(f,a+h);
         case 3,
            dy = (-3*feval(f,a) + 4*feval(f,a+h) - feval(f,a+2*h))/2;
         case 4,
            dy = (-11*feval(f,a) + 18*feval(f,a+h) - 9*feval(f,a+2*h) + ...
                 2*feval(f,a+3*h))/6;
       end
    end
    dy = dy/h;
%----------------------------------------------------------------------