5

I need to diagonalize a symbolic matrix with python. In Mathematica it can be done easily, but when using the module numpy.linalg I get problems.

For concreteness, consider the matrix

[[2, x], [x, 3]]

where x is a symbolic variable. I guess I get problems because the numpy package is provided for numerical computations, not symbolic, but I cannot find how to do it with sympy.

user
  • 5,370
  • 8
  • 47
  • 75
dapias
  • 2,512
  • 3
  • 15
  • 23
  • What if you create `x = sympy.Symbol('x')`, then initialize a numpy matrix as `A = np.array([[2, x], [x, 3]])`? – wflynny Sep 09 '13 at 16:02
  • @Bill That fails since numpy can't safely coerce the type to something it can handle. – Hooked Sep 09 '13 at 19:40
  • 1
    Matrix([[2,x],[x,3]]).diagonalize() should be enough. Do NOT use numpy for this, the numerical algorithms are completely inappropriate for symbolic calculation even if you use the `dtype=object`. – Krastanov Sep 09 '13 at 21:38
  • @Krastanov, can you explain me why I get two matrices using Matrix([[2,x],[x,3]]).diagonalize()? The second one is the correct but no idea about the first one. Thank you – dapias Sep 09 '13 at 22:04

2 Answers2

7

You can compute it from the eigenvalues, but there is actually a method that will do it for you, diagonalize

In [13]: M.diagonalize()
Out[13]:
⎛                                        ⎡     __________                       ⎤⎞
⎜                                        ⎢    ╱    2                            ⎥⎟
⎜⎡      -2⋅x                2⋅x       ⎤  ⎢  ╲╱  4⋅x  + 1    5                   ⎥⎟
⎜⎢─────────────────  ─────────────────⎥, ⎢- ───────────── + ─          0        ⎥⎟
⎜⎢   __________         __________    ⎥  ⎢        2         2                   ⎥⎟
⎜⎢  ╱    2             ╱    2         ⎥  ⎢                                      ⎥⎟
⎜⎢╲╱  4⋅x  + 1  - 1  ╲╱  4⋅x  + 1  + 1⎥  ⎢                        __________    ⎥⎟
⎜⎢                                    ⎥  ⎢                       ╱    2         ⎥⎟
⎜⎣        1                  1        ⎦  ⎢                     ╲╱  4⋅x  + 1    5⎥⎟
⎜                                        ⎢         0           ───────────── + ─⎥⎟
⎝                                        ⎣                           2         2⎦⎠

M.diagonalize() returns a pair of matrices (P, D) such that M = P*D*P**-1. If it can't compute enough eigenvalues, either because the matrix is not diagonalizable or because solve() can't find all the roots of the characteristic polynomial, it will raise MatrixError.

See also this section of the SymPy tutorial.

user
  • 5,370
  • 8
  • 47
  • 75
asmeurer
  • 86,894
  • 26
  • 169
  • 240
4

Assuming the matrix is diagonalizable, you can get the eigenvectors and eigenvalues by

from sympy import *
x = Symbol('x')
M = Matrix([[2,x],[x,3]])
print M.eigenvects()
print M.eigenvals()

Giving:

[(-sqrt(4*x**2 + 1)/2 + 5/2, 1, [[-x/(sqrt(4*x**2 + 1)/2 - 1/2)]
[                            1]]), (sqrt(4*x**2 + 1)/2 + 5/2, 1, [[-x/(-sqrt(4*x**2 + 1)/2 - 1/2)]
[                             1]])]
{sqrt(4*x**2 + 1)/2 + 5/2: 1, -sqrt(4*x**2 + 1)/2 + 5/2: 1}

You should check out the documentation, there are many other decompositions listed there.

Note that not every matrix is diagonalizable, but you can put every matrix into Jordan Normal Form using the sympy command .jordan_form.

Hooked
  • 84,485
  • 43
  • 192
  • 261
  • Thank you Hooked, but now I am confused, why do we got three components of the eigenvector if the matrix is 2x2. It's supposed that the eigenvectors cannot have greater dimension than the matrix. What do you thinK? – dapias Sep 09 '13 at 21:37
  • You've linked to an out-of-date version of the documentation. Replace the version with "latest" in the url to get the latest version. – asmeurer Sep 09 '13 at 22:10
  • @asmeurer Fixed and thanks, I didn't know that. Is that specific to the backend it's running, i.e. does this work for other sites? – Hooked Sep 10 '13 at 01:44
  • @dapias From the docs: three values are returned "A list of triples (eigenval, multiplicity, basis)". If all you are interested in is the diagonalization, asmeurer's answer is better (but do take note of the Jordan Normal form!). – Hooked Sep 10 '13 at 01:48
  • Hooked no this is something the SymPy docs implemented. The backend is just files hosted on github pages. – asmeurer Sep 10 '13 at 14:09