使用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 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過於艱難，因爲它引入了浮點數。當你想要一個符號解決方案時，浮動應該避免。這意味着：

• 而不是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讓他們快得多。此外，你可以做一個多邊形對象並應用roots它，這可能是最有效的：（這將是更容易做from sympy import *比輸入所有這些SP）

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

---

我對於爲什麼SymPy的特徵值方法即使在進行上述修改的情況下報告失敗也不太清楚。正如你可以看到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的申請返回{}，找不到根。看起來像是實施的不足之處。