deriv2.m

matlab算法集 matlab算法集

function y = deriv2 (a,h,p,m,e,f)
%-------------------------------------------------------------------- 
% Usage:       y = deriv2 (a,h,p,m,e,f);
%
% Description: Numerically estimate the second derivative of f(x) at
%              x = a.
%
% Inputs:      a = point of evaluation
%              h  = step size (h > 0)
%              p  = number of points to use (2 <= p <= 5)
%              m  = type of approximation:
%
%                       m      Approximation                 p
%                      ------------------------------------------
%                      -1   p-point backward difference   3,4,5
%                       0   p-point central difference    3,5
%                       1   p-point forward difference    3,4,5
%                      ------------------------------------------
%
%              e  = extrapolation parameter.  If e is nonzero,
%                   extrapoaltion 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)
%
% Output:      y = f''(a) 
%-------------------------------------------------------------------- 

% Extrapolation 

   p = min(2,min(5,p));
   if e 
      x = 2^(p-1);
      y = (x*deriv2(a,h,p,m,0,f) - deriv2(a,2*h,p,m,0,f))/(x-1);
   end

   if m < 0

% Backward difference 

      p = max(3,min(4,p));
      switch p 
         case 3,
            y = (f(a) - 2*f(a-h) + f(a-2*h));
            break;
         case 4,
            y = (2*f(a) - 5*f(a-h) + 4*f(a-2*h) - f(a-3*h));
            break; 
      end
 
   elseif m == 0
    
% Central difference /

      if p <= 3 
         p = 3;
      else
         p = 5;
      end
      switch p 
         case 3,
            y = f(a-h) - 2*f(a) + f(a+h);
            break;
         case 5,
            y = (-f(a-2*h) + 16*f(a-h) - 30*f(a) + 16*f(a+h)...
                 -f(a+2*h))/12;
            break; 
      end

   else 
       
% Forward difference 

      p = max(3,min(4,p));
      switch p 
         case 3,
            y = f(a) - 2*f(a+h) + f(a+2*h);
            break;
         case 4,
            y = 2*f(a) - 5*f(a+h) +4*f(a+2*h) - f(a+3*h);
            break; 
       end
    end
           
    y = y/(h*h);