python - Solve for eigenvectors "manually" using Sympy

Question

Setup

from sympy import Matrix, symbols, eye, factor, solve, simplify, solveset
from fractions import Fraction
_lambda = symbols("lambda")
x1, x2 = symbols("x_1, x_2")

Here's the matrices I'm working with

S = Matrix([[5,4],[4,5]])
M = Matrix([[1,0],[0,4]])
H = M**Fraction(-1,2)*S*M**Fraction(-1,2)

The next step I'd like to do is something like this

eigen = []
for eigenvalue in solve((S-_lambda*M).det()):
    eigen.append((
                   eigenvalue,
                   solve((S-_lambda*M).subs(_lambda,eigenvalue) * Matrix([x1,x2]))
                ))

But sympy can't find the answer... It probably just doesn't understand what I'm asking for but I'm not sure how to tell it more clearly.

print(eigen)
> [(25/8 - sqrt(481)/8, {x_1: nan, x_2: nan}),
 (sqrt(481)/8 + 25/8, {x_1: nan, x_2: nan})]

标签： pythonsympy

Answer 1

I'm not sure I quite understand what you are trying to do. The matrix H is unused. It looks like you want the eigenvalues and eigenvectors of the matrix S but I'm not sure what the M matrix is doing. I'll assume that you want to find the eigenvalues and eigenvectors of S.

We have the matrix S:

In [60]: S = Matrix([[5,4],[4,5]])                                                                                                             

In [61]: S                                                                                                                                     
Out[61]: 
⎡5  4⎤
⎢    ⎥
⎣4  5⎦

We first find its eigenvalues:

In [62]: l = Symbol('lambda')                                                                                                                  

In [63]: det(S - l*eye(2))                                                                                                                     
Out[63]: 
       2     
(5 - λ)  - 16

In [64]: l1, l2 = roots(_, l)                                                                                                                  

In [65]: l1                                                                                                                                    
Out[65]: 9

In [66]: l2                                                                                                                                    
Out[66]: 1

Now each of these eigenvalues has a corresponding eigenvector. Actually it's not a single eigenvector but rather a 1D family of vectors. If the vector is [x1, x2] then we can't uniquely solve for x1 and x2 because there is an infinite family of solutions. Using solve therefore gives a parametrised solution for x1 in terms of x2:

In [67]: x1, x2 = symbols('x1, x2')                                                                                                            

In [68]: solve((S - l1*eye(2)) * Matrix([x1, x2]), [x1, x2])                                                                                   
Out[68]: {x₁: x₂}

This tells us that a vector [x1, x2] will be a solution provided x1 == x2. Hence we have the family of eigenvectors [a, a] with eignevalue l1 (9):

In [69]: S*Matrix([1, 1])                                                                                                                      
Out[69]: 
⎡9⎤
⎢ ⎥
⎣9⎦

More directly what we actually want is the nullspace of S - l*I which we can compute as

In [70]: (S - l1*eye(2)).nullspace()                                                                                                           
Out[70]: 
⎡⎡1⎤⎤
⎢⎢ ⎥⎥
⎣⎣1⎦⎦

In [71]: (S - l2*eye(2)).nullspace()                                                                                                           
Out[71]: 
⎡⎡-1⎤⎤
⎢⎢  ⎥⎥
⎣⎣1 ⎦⎦

The nullspace method returns a list of vectors that are a basis for the nullspace. In this context that means that they are a basis for the eigenspace of the given eigenvalue. Since we have only one vector it is a 1D eigenspace consisting of all scalar multiples of [1, 1] (or [-1, 1] for l2). If there were repeated eigenvalues then it would be possible that the nullspace could have greater dimension.