Calcul des valeurs propres symboliques avec sympy

Question

J'essaie de calculer des valeurs propres d'une matrice complexe symbolique M de taille 3x3. Dans certains cas, eigenvals() fonctionne parfaitement. Par exemple, le code suivant:

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

me retourne 3 valeurs propres symboliques correctes complexes de M. Cependant, quand je mets x=1., je reçois l'erreur suivante:

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

J'ai essayé aussi de calculer les valeurs propres comme suit:

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

mais il me renvoie un ConditionSet dans tous les cas, même si eigenvals() peut fait le travail.

Est-ce que quelqu'un sait comment résoudre correctement ce problème de valeur propre, pour une valeur de x?

Answer 1

Votre définition de M a rendu la vie trop difficile pour SymPy car elle a introduit des nombres à virgule flottante. Quand vous voulez une solution symbolique, les flotteurs doivent être évités. Cela signifie:

• au lieu de 1./3. (nombre à virgule flottante de Python) utilisent sp.Rational(1, 3) (nombre rationnel de sympy) ou sp.S(1)/3 qui a le même effet mais il est plus facile de taper.
• au lieu de 1j (unité imaginaire de Python) utilisent sp.I (unité imaginaire de sympy)
• au lieu de x = 1., écrire x = 1 (Python 2.7 habitudes et sympy vont mal ensemble).

Avec ces changements soit solveset ou solve trouver les valeurs propres, bien que les solve devient beaucoup plus rapide. En outre, vous pouvez faire un objet Poly et appliquer roots à elle, ce qui est probablement la plus efficace: (. Il serait plus facile de faire from sympy import * que le type tous ces 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)) 

---

I Je ne comprends pas très bien pourquoi la méthode eigenvals de SymPy signale un échec même avec les modifications ci-dessus. Comme vous pouvez le voir in the source, il ne fait pas beaucoup plus que ce que fait le code ci-dessus: appelez roots sur le polynôme caractéristique. La différence semble être dans la façon dont ce polynôme est créé: M.charpoly(lam) retours

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

avec mystérieux (pour moi) domain='EX'. Par la suite, une application de roots renvoie {}, aucune racine trouvée. On dirait une déficience de la mise en œuvre.