使用sympy计算符号特征值

Question

使用sympy计算符号特征值

7

我正在尝试计算一个符号复矩阵 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的值吗？

- Azlof

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user6655984 · Accepted Answer

您对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 方法即使进行了上述修改也无法成功报告。正如您在源代码中所看到的，它并没有比上面的代码更多地做一些事情：在特征多项式上调用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函数返回{}，表示没有找到根。看起来这个实现可能存在缺陷。