2017-09-22 14 views
5

Tôi đang cố gắng tính giá trị riêng của ma trận phức hợp tượng trưng M của kích thước 3x3. Trong một số trường hợp, eigenvals() hoạt động hoàn hảo. Ví dụ: mã sau:Tính toán các giá trị riêng biệt tượng trưng với sympy

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

trả lại cho tôi 3 giá trị riêng biệt biểu tượng phức tạp của M. Tuy nhiên, khi tôi đặt x=1., tôi nhận được lỗi sau:

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

Tôi cũng đã cố gắng để tính toán giá trị riêng như sau:

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

nhưng nó sẽ trả về cho tôi một ConditionSet trong mọi trường hợp, ngay cả khi eigenvals() thể thực hiện công việc.

Có ai biết cách giải quyết đúng vấn đề về giá trị riêng này không, với bất kỳ giá trị nào của x?

Trả lời

2

Định nghĩa của bạn về M khiến cuộc sống quá khó đối với SymPy vì nó giới thiệu các số dấu phẩy động. Khi bạn muốn có một giải pháp tượng trưng, ​​phao nổi phải tránh. Điều đó có nghĩa:

  • thay vì 1./3. (số điểm nổi Python) sử dụng sp.Rational(1, 3) (số hợp lý SymPy) hoặc sp.S(1)/3 trong đó có tác dụng tương tự nhưng là dễ dàng hơn để gõ.
  • thay vì 1j (đơn vị ảo của Python) sử dụng sp.I (đơn vị ảo SymPy của)
  • thay vì x = 1., viết x = 1 (Python 2,7 thói quen và SymPy đi kém với nhau).

Với những thay đổi hoặc solveset hay solve tìm ra giá trị riêng, mặc dù solve được chúng nhanh hơn nhiều. Ngoài ra, bạn có thể làm cho một đối tượng Poly và áp dụng roots với nó, mà có lẽ hiệu quả nhất: (. Nó sẽ được dễ dàng hơn để làm from sympy import * hơn gõ tất cả các 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)) 


tôi 'không hoàn toàn rõ ràng về lý do tại sao phương pháp riêng của SymPy báo cáo thất bại ngay cả với những sửa đổi ở trên. Như bạn có thể thấy in the source, nó không hoạt động nhiều hơn những gì mã trên làm: gọi roots trên đa thức đặc trưng. Sự khác biệt dường như là một cách đa thức này được tạo ra: M.charpoly(lam) lợi nhuận

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

với bí ẩn (với tôi) domain='EX'. Sau đó, một ứng dụng của roots trả về {}, không tìm thấy rễ. Có vẻ như thiếu sự thực thi.

+0

Cảm ơn rất nhiều sự giúp đỡ của bạn. Có vẻ như vấn đề của tôi xuất phát từ việc sử dụng 1j thay vì sp.I, nhưng việc sử dụng Rational chắc chắn sẽ giúp bạn! Vấn đề được giải quyết cho tôi, nhưng vẫn còn có một vấn đề với bản địa của SymPy ... – Azlof

+0

Tôi đã đơn giản hóa ví dụ của bạn và đăng [dưới dạng vấn đề SymPy] (https://github.com/sympy/sympy/issues/13340) – FTP

+0

Vấn đề được giải quyết trên github.Đối với những người quan tâm, bản sửa lỗi đã được đẩy vào nhánh chủ của SymPy. Cảm ơn Michelle! – Azlof

Các vấn đề liên quan