用sympy计算符号特征值 [英] Computation of symbolic eigenvalues with sympy

问题描述

我正在尝试计算大小为3x3的符号复数矩阵M的特征值.在某些情况下，eigenvals()可以完美运行.例如，以下代码:

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

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

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

raise MatrixError(无法为{}".format(self)计算特征值

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)

，但是即使eigenvals()可以完成任务，在任何情况下都返回ConditionSet.

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

对于x的任何值，有人知道如何正确解决此特征值问题吗?

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

推荐答案

您对M的定义使SymPy的生活变得困难，因为它引入了浮点数.当您需要符号解决方案时，应避免使用浮点数.这意味着:

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:

• 代替1./3.(Python的浮点数)，使用sp.Rational(1, 3)(SymPy的有理数)或sp.S(1)/3具有相同的效果，但更易于键入.
• 使用sp.I(SymPy的虚构单元)代替1j(Python的虚构单元)
• 代替x = 1.编写x = 1(Python 2.7习惯和SymPy在一起使用效果不佳).

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

通过这些更改，solveset或solve可以找到特征值，尽管solve可以更快地获取它们.另外，您可以制作一个Poly对象并将roots应用于它，这可能是最有效的:

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

(执行from sympy import *比​​键入所有这些sp.要容易得多)

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

我不清楚为什么SymPy的特征值方法即使进行了上述修改也报告失败.您可以在源中看到 ，它的作用不只是上面的代码:在特征多项式上调用roots.区别似乎在于创建多项式的方式:M.charpoly(lam)返回

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

与神秘的(对我来说)domain='EX'.随后，roots的应用程序返回{}，未找到任何根.看起来是实现的不足.

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

这篇关于用sympy计算符号特征值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！