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個の正しい複雑なシンボリック固有値を返します。場合でもeigenvals()ができ、

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

を、それは私にどのような場合でもConditionSetを返します。私も次のように固有値を計算してみました

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

：私はx=1.設定した場合しかし、私は次のエラーを取得します仕事をしなさい。

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を適用し、おそらく最も効率的であることができます。

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

（すべてのこれらのSPを入力するよりも、from sympy import *を行うことが容易になります。）

---

私はSymPyの固有値法が上記の修正を行っても失敗を報告する理由についてはっきりしていません。 in the sourceが表示されているので、上記のコードよりもはるかに機能しません。特性多項式でrootsを呼び出します。 M.charpoly(lam)リターンdomain='EX'（私には）不思議なと

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

：違いは、この多項式を作成する方法であるように思われます。その後、rootsのアプリケーションは{}を返します。ルートは見つかりません。実装の不足のように見えます。