用 sympy 计算符号特征值

Question

我正在尝试计算大小为3x3 的符号复矩阵M 的特征值。在某些情况下，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 calculate 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 来说太难了，因为它引入了浮点数。当你想要一个符号解决方案时，要避免浮动。这意味着：

• 而不是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 的应用程序返回 {}，没有找到根。看起来像实施的缺陷。