Calculating Jacobian matrix in sympy

2,409 views
Skip to first unread message

Kjetil brinchmann Halvorsen

unread,
Aug 23, 2012, 3:21:58 PM8/23/12
to sy...@googlegroups.com
I am using sympy to (try to) calculate a Jacobian matrix. It seems to
work, and I beleave correctly, but}is giving an error I dont
understand.
I don't know much python, so these surely can be simplified!
Code:

a11,a12,a13,a22,a23,a33= symbols('a11 a12 a13 a22 a23 a33', real=True)
x11,x12,x13,x22,x23,x33= symbols('x11 x12 x13 x22 x23 x33', real=True)
A = Matrix([[a11,a12,a13],[a12,a22,a23],[a13,a23,a33]])

In [7]: X = Matrix([[x11,x12,x13],[x12,x22,x23],[x13,x23,x33]])
Y = A*X+X*A
J = Matrix(6,6,lambda i,j:0) # This is for the Jacobian
vars = [x11,x12,x13,x22,x23,x33]

In [38]: for i in [0,1,2]:
....: for j in [0,1,2]:
....: k = 3*i+j
....: for l in [0,1,2,3,4,5]:
....: J[k,l] = diff(Y[i,j], vars[l])
....:
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
<ipython-input-38-46a5a73fca3d> in <module>()
3 k = 3*i+j
4 for l in [0,1,2,3,4,5]:
----> 5 J[k,l] = diff(Y[i,j], vars[l])
6

/home/kjetil/py/sympy/sympy/sympy/sympy/matrices/matrices.pyc in
__setitem__(self, key, value)
3537
3538
-> 3539 i, j = self.key2ij((i,j))
3540 if not (i>=0 and i<self.rows and j>=0 and j <
self.cols):
3541 raise IndexError("Index out of range:
a[%s]" % (key,))

/home/kjetil/py/sympy/sympy/sympy/sympy/matrices/matrices.pyc in
key2ij(self, key)
1102 i, j = [(k + n) if k < 0 else k for k, n in zip(key,
(self.rows, self.cols))]
1103 if not (i>=0 and i<self.rows and j>=0 and j < self.cols):
-> 1104 raise IndexError("Index out of range: a[%s]"%repr(key))
1105 return i,j
1106

IndexError: Index out of range: a[(6, 0)]

But J seems to be correct:

In [40]: J
Out[40]:
⎡2⋅a₁₁ 2⋅a₁₂ 2⋅a₁₃ 0 0 0 ⎤
⎢ ⎥
⎢ a₁₂ a₁₁ + a₂₂ a₂₃ a₁₂ a₁₃ 0 ⎥
⎢ ⎥
⎢ a₁₃ a₂₃ a₁₁ + a₃₃ 0 a₁₂ a₁₃⎥
⎢ ⎥
⎢ a₁₂ a₁₁ + a₂₂ a₂₃ a₁₂ a₁₃ 0 ⎥
⎢ ⎥
⎢ 0 2⋅a₁₂ 0 2⋅a₂₂ 2⋅a₂₃ 0 ⎥
⎢ ⎥
⎣ 0 a₁₃ a₁₂ a₂₃ a₂₂ + a₃₃ a₂₃⎦

In [41]:

???

Kjetil




--
"If you want a picture of the future - imagine a boot stamping on the
human face - forever."

George Orwell (1984)

Kjetil brinchmann Halvorsen

unread,
Aug 23, 2012, 3:30:58 PM8/23/12
to sy...@googlegroups.com
Followup. After the code posted in the previous post, I did:

In [43]: J.det()
Out[43]: 0

I don't beleave that! The jacobian of that transformation cannot be zero!

Kjetil

Matthew Rocklin

unread,
Aug 23, 2012, 3:34:01 PM8/23/12
to sy...@googlegroups.com
J = Matrix(6,6,lambda i,j:0)                 #  This is for the Jacobian
 vars = [x11,x12,x13,x22,x23,x33]

I'm not totally clear what you're doing here
 

In [38]: for i in [0,1,2]:
   ....: for j in [0,1,2]:
   ....:     k = 3*i+j
   ....:     for l in [0,1,2,3,4,5]:
   ....:         J[k,l] = diff(Y[i,j], vars[l])
   ....:
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-38-46a5a73fca3d> in <module>()
      3             k = 3*i+j
      4             for l in [0,1,2,3,4,5]:
----> 5                     J[k,l] = diff(Y[i,j], vars[l])
In this line Y is of shape 3, 3 but i goes from 0 to 5. So you're going outside the bounds of the Matrix. You're asking for the 3rd, 4th, and 5th, row of a 3x3 matrix. 

 

Matthew Rocklin

unread,
Aug 23, 2012, 3:44:21 PM8/23/12
to sy...@googlegroups.com
There is also a Matrix.jacobian method that might interest you. 
Reply all
Reply to author
Forward
0 new messages