[Scipy-tickets] [SciPy] #564: Change derivative() to return the gradient of multidimensional functions

SciPy scipy-tickets@scipy....
Sun Dec 16 19:13:23 CST 2007


#564: Change derivative() to return the gradient of multidimensional functions
-------------------------+--------------------------------------------------
 Reporter:  robfalck     |       Owner:  somebody
     Type:  enhancement  |      Status:  new     
 Priority:  normal       |   Milestone:  0.7     
Component:  scipy.misc   |     Version:          
 Severity:  minor        |    Keywords:          
-------------------------+--------------------------------------------------
 The following version of derivative in scipy.misc will return an array
 representing the gradient of a multidimensional function.  Currently
 derivative() only works for functions with a single independent variable.
 The difference in computational speed for one-dimensional functions has
 not been assessed.

 {{{
 def derivative(func,x0,dx=1.0,n=1,args=(),order=3):
     """Given a function, use a central difference formula with spacing dx
 to
        compute the nth derivative at x0.

        order is the number of points to use and must be odd.

        Warning: Decreasing the step size too small can result in
        round-off error.
     """
     assert (order >= n+1), "Number of points must be at least the
 derivative order + 1."
     assert (order % 2 == 1), "Odd number of points only."
     # pre-computed for n=1 and 2 and low-order for speed.
     if n==1:
         if order == 3:
             weights = array([-1,0,1])/2.0
         elif order == 5:
             weights = array([1,-8,0,8,-1])/12.0
         elif order == 7:
             weights = array([-1,9,-45,0,45,-9,1])/60.0
         elif order == 9:
             weights = array([3,-32,168,-672,0,672,-168,32,-3])/840.0
         else:
             weights = central_diff_weights(order,1)
     elif n==2:
         if order == 3:
             weights = array([1,-2.0,1])
         elif order == 5:
             weights = array([-1,16,-30,16,-1])/12.0
         elif order == 7:
             weights = array([2,-27,270,-490,270,-27,2])/180.0
         elif order == 9:
             weights =
 array([-9,128,-1008,8064,-14350,8064,-1008,128,-9])/5040.0
         else:
             weights = central_diff_weights(order,2)
     else:
         weights = central_diff_weights(order, n)

     ho = order >> 1
     ax0 = asfarray(x0).flatten()
     lx0 = len(ax0)
     derivs = zeros(lx0)
     for i in range(lx0):
         val = 0.0
         adx = zeros(lx0)
         adx[i] = dx
         for k in range(order):
             val += weights[k]*func(x0+(k-ho)*adx,*args)
         derivs[i] = val / product((dx,)*n,axis=0)
     return derivs
 }}}

-- 
Ticket URL: <http://scipy.org/scipy/scipy/ticket/564>
SciPy <http://www.scipy.org/>
SciPy is open-source software for mathematics, science, and engineering.


More information about the Scipy-tickets mailing list