Computation of symbolic eigenvalues with sympy

Question

I'm trying to compute eigenvalues of a symbolic complex matrix Mof size 3x3. In some cases, eigenvals() works perfectly. 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 compute eigenvalues for {}".format(self))

I also tried to compute eigenvalues as follows:

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

but it returns me a ConditionSet in any case, even when eigenvals() can do the job.

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

Solly

Answer 1

Your definition of M made life too hard for SymPy because it introduced floating point numbers. When you want a symbolic solution, floats are to be avoided. That means:

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

With these changes either solveset or solve find the eigenvalues, although solve gets them much faster. Also, you can make a Poly object and apply roots to it, which is probably 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 type all these sp.)

---

I'm not quite clear on why SymPy's eigenvals method reports failure even with the above modifications. As you can see in the source, it doesn't do much more than what the above code does: call roots on the characteristic polynomial. The difference appears to be in the way 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 mysterious (to me) domain='EX'. Subsequently, an application of roots returns {}, no roots found. Looks like a deficiency of the implementation.

Hemicycle

Thanks a lot for your help. It seems that my problem came rather from the use of 1j instead of sp.I, but the use of Rational with surely help! Problem solved for me, but still there is an issue with SymPy's eigenvals... – Solly

I simplified your example and posted as a SymPy issue – Hemicycle

Issue solved on github. For those who are interested, the fix has been pushed in branch master of SymPy. Thanks Michelle! – Solly