All geek questions in one place

Calculating symbolic eigenvalues ​​with sympy

I am trying to calculate the eigenvalues ​​of a symbolic complex matrix M size 3x3 . In some cases, eigenvals() works fine. For example, the following code:

 import sympy as sp kx = sp.symbols('kx') x = 0. M = sp.Matrix([[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]]) M[0, 0] = 1. M[0, 1] = 2./3. M[0, 2] = 2./3. M[1, 0] = sp.exp(1j*kx) * 1./6. + x M[1, 1] = sp.exp(1j*kx) * 2./3. M[1, 2] = sp.exp(1j*kx) * -1./3. M[2, 0] = sp.exp(-1j*kx) * 1./6. M[2, 1] = sp.exp(-1j*kx) * -1./3. M[2, 2] = sp.exp(-1j*kx) * 2./3. dict_eig = M.eigenvals()

returns me 3 correct complex symbolic eigenvalues ​​of M However, when I set x=1. I get the following error:

raise MatrixError ("Could not calculate eigenvalues ​​for {}". format (self))

I also tried to calculate the eigenvalues ​​as follows:

 lam = sp.symbols('lambda') cp = sp.det(M - lam * sp.eye(3)) eigs = sp.solveset(cp, lam)

but it returns a ConditionSet to me anyway, even when eigenvals() can do the job.

Does anyone know how to properly solve this problem with an eigenvalue, for any value of x ?

+5
source share
1 answer

Your definition of M made life too complicated for SymPy because it introduced floating point numbers. When you need a symbolic solution, you should avoid floats. It means:

  • instead of 1./3. (Python floating-point number) use sp.Rational(1, 3) (Sympy is a rational number) or sp.S(1)/3 , which has the same effect but is easier to type.
  • instead of 1j (the imaginary unit of Python) use sp.I (the imaginary unit of SymPy)
  • instead of x = 1. , write x = 1 (Python 2.7 and SymPy habits go bad together).

With these changes, solveset or solve find eigenvalues, although solve gets them much faster. In addition, you can create a Poly object and apply roots to it, which is probably the most efficient:

 M = sp.Matrix([ [ 1, sp.Rational(2, 3), sp.Rational(2, 3), ], [ sp.exp(sp.I*kx) * sp.Rational(1, 6) + x, sp.exp(sp.I*kx) * sp.Rational(1, 6), sp.exp(sp.I*kx) * sp.Rational(-1, 3), ], [ sp.exp(-sp.I*kx) * sp.Rational(1, 6), sp.exp(-sp.I*kx) * sp.Rational(-1, 3), sp.exp(-sp.I*kx) * sp.Rational(2, 3), ] ]) lam = sp.symbols('lambda') cp = sp.det(M - lam * sp.eye(3)) eigs = sp.roots(sp.Poly(cp, lam))

(It would be easier to do from sympy import * than to enter all these sp.)

I don’t quite understand why the SymPy eigenvalue method reports an error even with the above modifications. As you can see in the source , this is not much more than what the above code does: call roots in the characteristic polynomial, The difference seems to be how this polynomial is created: M.charpoly(lam) returns

 PurePoly(lambda**3 + (I*sin(kx)/2 - 5*cos(kx)/6 - 1)*lambda**2 + (-I*sin(kx)/2 + 11*cos(kx)/18 - 2/3)*lambda + 1/6 + 2*exp(-I*kx)/3, lambda, domain='EX')

with the mysterious (to me) domain='EX' . Subsequently, the roots application returns {} , no roots were found. Looks like a lack of implementation.

+2
source

Source: https://habr.com/ru/post/1272035/

More articles:


All Articles