[SciPy-user] Jacobian matrix

Nils Wagner nwagner at mecha.uni-stuttgart.de
Mon Mar 21 01:29:22 CST 2005


John Gleeson wrote:

>
> On Mar 17, 2005, at 9:59 AM, Nils Wagner wrote:
>
>> Hi all,
>>
>> Is there a simpler way (built-in function) to compute the Jacobian 
>> of  a vector function ?
>> How about higher order derivatives of vector-valued functions ?
>>
>>
>> from scipy import *
>>
>> def f(x):
>> #
>>   tmp = zeros(5,Float)
>>   tmp[0] = x[0] + x[1] + x[2]**2 + x[3]**2 + x[4]**2 - 2
>>   tmp[1] = x[0] - x[1] + x[2]**2 + x[3]**2 + x[4]**2
>>   tmp[2] = -x[2]**2 + x[3]**2 + x[4]**2
>>   tmp[3] =  x[2]**2 - x[3]**2 + x[4]**2
>>   tmp[4] =  x[2]**2 + x[3]**2 - x[4]**2
>>   return tmp
>>
>> x0 = array(([1.02,1.02,0.02,0.02,0.02]))
>> eps = 1.e-5
>> J = zeros((5,5),Float)
>> for i in arange(0,5):
>>
>>  ei = zeros(5,Float)
>>  ei[i] = 1.0
>>  J[:,i] = (f(x0+eps*ei)-f(x0))/eps
>>
>>
>
> This approach doesn't uses the ScientificPython extension to scipy, 
> so  it doesn't satisfy your request for a built-in function, but it 
> is  fairly simple, and even provides higher order derivatives:

Afaik, ScientificPython is an extra package and not a part of scipy.
http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Scientific_8.html
However, it would be nice, if this functionality is directly available 
in scipy as well.

NIls

>
> >>> def f(x):
> ... #
> ...   tmp = zeros(5,PyObject)  # --- Notice change in typecode ---
> ...   tmp[0] = x[0] + x[1] + x[2]**2 + x[3]**2 + x[4]**2 - 2
> ...   tmp[1] = x[0] - x[1] + x[2]**2 + x[3]**2 + x[4]**2
> ...   tmp[2] = -x[2]**2 + x[3]**2 + x[4]**2
> ...   tmp[3] =  x[2]**2 - x[3]**2 + x[4]**2
> ...   tmp[4] =  x[2]**2 + x[3]**2 - x[4]**2
> ...   return tmp
>
> >>> x0 =  
> array(([DerivVar(1.02,0,2),DerivVar(1.02,1,2),DerivVar(0.02,2,2),DerivVa 
> r(0.02,3,2),DerivVar(0.02,4,2)]),PyObject)
> >>> D = f(x0)
> >>> D
> array([(0.041199999999999903, [1.0, 1.0, 0.040000000000000001,  
> 0.040000000000000001, 0.040000000000000001], [[0.0, 0.0, 0.0, 0.0,  
> 0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0, 0.0], [0.0, 
> 0.0,  0.0, 2.0, 0.0], [0.0, 0.0, 0.0, 0.0, 2.0]]) ,
>             (0.0012000000000000001, [1.0, -1.0, 0.040000000000000001,  
> 0.040000000000000001, 0.040000000000000001], [[0.0, 0.0, 0.0, 0.0,  
> 0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0, 0.0], [0.0, 
> 0.0,  0.0, 2.0, 0.0], [0.0, 0.0, 0.0, 0.0, 2.0]]) ,
>             (0.00040000000000000002, [0.0, 0.0, 
> -0.040000000000000001,  0.040000000000000001, 0.040000000000000001], 
> [[0.0, 0.0, 0.0, 0.0,  0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 
> -2.0, 0.0, 0.0], [0.0, 0.0,  0.0, 2.0, 0.0], [0.0, 0.0, 0.0, 0.0, 
> 2.0]]) ,
>             (0.00040000000000000002, [0.0, 0.0, 0.040000000000000001,  
> -0.040000000000000001, 0.040000000000000001], [[0.0, 0.0, 0.0, 0.0,  
> 0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0, 0.0], [0.0, 
> 0.0,  0.0, -2.0, 0.0], [0.0, 0.0, 0.0, 0.0, 2.0]]) ,
>             (0.00040000000000000002, [0.0, 0.0, 0.040000000000000001,  
> 0.040000000000000001, -0.040000000000000001], [[0.0, 0.0, 0.0, 0.0,  
> 0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0, 0.0], [0.0, 
> 0.0,  0.0, 2.0, 0.0], [0.0, 0.0, 0.0, 0.0, -2.0]]) ],'O')
>
> # To get them into a more natural form:
>
> >>> f0 = array(([x[0] for x in D]))
> >>> f0
> array([ 0.0412,  0.0012,  0.0004,  0.0004,  0.0004])
> >>> Df0 = array(([x[1] for x in D]))
> >>> Df0
> array([[ 1.  ,  1.  ,  0.04,  0.04,  0.04],
>        [ 1.  , -1.  ,  0.04,  0.04,  0.04],
>        [ 0.  ,  0.  , -0.04,  0.04,  0.04],
>        [ 0.  ,  0.  ,  0.04, -0.04,  0.04],
>        [ 0.  ,  0.  ,  0.04,  0.04, -0.04]])
> >>> D2f0 = array(([x[2] for x in D]))
> >>> D2f0
> array([[[ 0.,  0.,  0.,  0.,  0.],
>         [ 0.,  0.,  0.,  0.,  0.],
>         [ 0.,  0.,  2.,  0.,  0.],
>         [ 0.,  0.,  0.,  2.,  0.],
>         [ 0.,  0.,  0.,  0.,  2.]],
>        [[ 0.,  0.,  0.,  0.,  0.],
>         [ 0.,  0.,  0.,  0.,  0.],
>         [ 0.,  0.,  2.,  0.,  0.],
>         [ 0.,  0.,  0.,  2.,  0.],
>         [ 0.,  0.,  0.,  0.,  2.]],
>        [[ 0.,  0.,  0.,  0.,  0.],
>         [ 0.,  0.,  0.,  0.,  0.],
>         [ 0.,  0., -2.,  0.,  0.],
>         [ 0.,  0.,  0.,  2.,  0.],
>         [ 0.,  0.,  0.,  0.,  2.]],
>        [[ 0.,  0.,  0.,  0.,  0.],
>         [ 0.,  0.,  0.,  0.,  0.],
>         [ 0.,  0.,  2.,  0.,  0.],
>         [ 0.,  0.,  0., -2.,  0.],
>         [ 0.,  0.,  0.,  0.,  2.]],
>        [[ 0.,  0.,  0.,  0.,  0.],
>         [ 0.,  0.,  0.,  0.,  0.],
>         [ 0.,  0.,  2.,  0.,  0.],
>         [ 0.,  0.,  0.,  2.,  0.],
>         [ 0.,  0.,  0.,  0., -2.]]])
>
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user

 



More information about the SciPy-user mailing list