python-2.7 - 用 sympy 计算符号特征值

Question

我正在尝试计算符号复矩阵的特征值 M尺寸3x3 .在某些情况下，eigenvals()完美地工作。例如，以下代码:

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()

返回 M 的 3 个正确的复符号特征值.但是，当我设置 x=1. 时，我收到以下错误:

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

我还尝试按如下方式计算特征值:

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

但它返回给我一个 ConditionSet无论如何，即使 eigenvals()可以胜任。

有谁知道如何正确解决这个特征值问题，对于 x 的任何值?

Answer 1

你对 M 的定义让 SymPy 过得很艰难，因为它引入了浮点数。当您需要符号解决方案时，应避免 float 。这意味着:

而不是 1./3. (Python 的浮点数)使用 sp.Rational(1, 3) (SymPy 的有理数)或 sp.S(1)/3它具有相同的效果，但更容易输入。

而不是 1j (Python 的虚数单位)使用 sp.I (SymPy 的虚数单位)

而不是 x = 1. , 写 x = 1 (Python 2.7 习惯和 SymPy 在一起很糟糕)。

随着这些变化 solveset或 solve找到特征值，虽然 solve让他们更快。此外，您可以制作 Poly 对象并应用 roots到它，这可能是最有效的:

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))

(执行 from sympy import * 比输入所有这些 sp 更容易。)

我不太清楚为什么 SymPy 的 eigenvals 方法即使进行了上述修改也会报告失败。如您所见 in the source ，它并没有比上面的代码做的更多:调用 roots在特征多项式上。不同之处似乎在于创建此多项式的方式: M.charpoly(lam)返回

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')

与神秘(对我而言) domain='EX' .随后，申请 roots返回 {} ，没有找到根。看起来像实现的缺陷。