octave_packages/optim-1.2.0/deriv.m

   1 ## Copyright (C) 2000 Ben Sapp <bsapp@nua.lampf.lanl.gov>
   2 ## Copyright (C) 2011 Joaquín Ignacio Aramendía <samsagax@gmail.com>
   3 ## Copyright (C) 2011 Carnë Draug <carandraug+dev@gmail.com>
   4 ##
   5 ## This program is free software; you can redistribute it and/or modify it under
   6 ## the terms of the GNU General Public License as published by the Free Software
   7 ## Foundation; either version 3 of the License, or (at your option) any later
   8 ## version.
   9 ##
  10 ## This program is distributed in the hope that it will be useful, but WITHOUT
  11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  13 ## details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License along with
  16 ## this program; if not, see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## @deftypefn {Function File} {@var{dx} =} deriv (@var{f}, @var{x0})
  20 ## @deftypefnx {Function File} {@var{dx} =} deriv (@var{f}, @var{x0}, @var{h})
  21 ## @deftypefnx {Function File} {@var{dx} =} deriv (@var{f}, @var{x0}, @var{h}, @var{O})
  22 ## @deftypefnx {Function File} {@var{dx} =} deriv (@var{f}, @var{x0}, @var{h}, @var{O}, @var{N})
  23 ## Calculate derivate of function @var{f}.
  24 ##
  25 ## @var{f} must be a function handle or the name of a function that takes @var{x0}
  26 ## and returns a variable of equal length and orientation. @var{x0} must be a
  27 ## numeric vector or scalar.
  28 ##
  29 ## @var{h} defines the step taken for the derivative calculation. Defaults to 1e-7.
  30 ##
  31 ## @var{O} defines the order of the calculation. Supported values are 2 (h^2 order)
  32 ## or 4 (h^4 order). Defaults to 2.
  33 ##
  34 ## @var{N} defines the derivative order. Defaults to the 1st derivative of the
  35 ## function. Can be up to the 4th derivative.
  36 ##
  37 ## Reference: Numerical Methods for Mathematics, Science, and Engineering by
  38 ## John H. Mathews.
  39 ## @end deftypefn
  40
  41 function dx = deriv (f, x0, h = 0.0000001, O = 2, N = 1)
  42
  43   if (ischar(f))
  44     f = str2func(f); # let's also support a string with str2func
  45   endif
  46
  47   if (nargin < 2)
  48     error ("Not enough arguments.");
  49   elseif (!isa (f, 'function_handle'))
  50     error ("The first argument 'f' must be a function handle.");
  51   elseif (!isvector (x0) || !isnumeric (x0))
  52     ## a scalar is 1x1 therefore counts as a vector too
  53     error ("The second argument 'x0' must be a numeric vector.");
  54   elseif (!isscalar (h) || !isnumeric (h))
  55     error ("The third argument 'h' must be a scalar.");
  56   elseif (!isscalar (O) || !isnumeric (O))
  57     error ("The fourth argument 'O' must be a scalar.");
  58   elseif (O != 2 && O != 4)
  59     error ("Only order 2 or 4 is supported.");
  60   elseif (!isscalar (N) || !isnumeric (N))
  61     error ("The fifth argument 'N' must be a scalar.");
  62   elseif ((N > 4) || (N < 1))
  63     error("Only 1st,2nd,3rd or 4th order derivatives are acceptable.");
  64   elseif (nargin > 5)
  65     warning("Ignoring arguements beyond the 5th.");
  66   endif
  67
  68   switch O
  69     case (2)
  70       switch N
  71         case (1)
  72           dx = (feval(f,x0+h)-feval(f,x0-h))/(2*h);
  73         case (2)
  74           dx = (feval(f,x0+h)-2*feval(f,x0)+feval(f,x0-h))/(h^2);
  75         case (3)
  76           dx = (feval(f,x0+2*h)-2*feval(f,x0+h)+2*feval(f,x0-h)-feval(f,x0-2*h))/(2*h^3);
  77         case (4)
  78           dx = (feval(f,x0+2*h)-4*feval(f,x0+h)+6*feval(f,x0)-4*feval(f,x0-h)+feval(f,x0-2*h))/(h^4);
  79         otherwise
  80           error("Only 1st,2nd,3rd or 4th order derivatives are acceptable.");
  81       endswitch
  82     case (4)
  83       switch N
  84         case (1)
  85           dx = (-feval(f,x0+2*h)+8*feval(f,x0+h)-8*feval(f,x0-h)+feval(f,x0-2*h))/(12*h);
  86         case (2)
  87           dx = (-feval(f,x0+2*h)+16*feval(f,x0+h)-30*feval(f,x0)+16*feval(f,x0-h)-feval(f,x0-2*h))/(12*h^2);
  88         case (3)
  89           dx = (-feval(f,x0+3*h)+8*feval(f,x0+2*h)-13*feval(f,x0+h)+13*feval(f,x0-h)-8*feval(f,x0-2*h)+feval(f,x0-3*h))/(8*h^3);
  90         case (4)
  91           dx = (-feval(f,x0+3*h)+12*feval(f,x0+2*h)-39*feval(f,x0+h)+56*feval(f,x0)-39*feval(f,x0-h)+12*feval(f,x0-2*h)-feval(f,x0-3*h))/(6*h^4);
  92         otherwise
  93           error("Only 1st,2nd,3rd or 4th order derivatives are acceptable.");
  94       endswitch
  95     otherwise
  96       error ("Only order 2 or 4 is supported.");
  97   endswitch
  98 endfunction