用 sympy 计算符号特征值
Computation of symbolic eigenvalues with sympy
我正在尝试计算大小为 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()
returns 我 3 个正确的复符号特征值 M
。但是,当我设置 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)
但它 returns 我 ConditionSet
无论如何,即使 eigenvals()
可以完成这项工作。
有谁知道如何正确解决这个特征值问题,对于 x
的任何值?
你对 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
可以更快地找到它们。此外,您可以制作一个 Poly 对象并对其应用 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))
(做 from sympy import *
比输入所有这些 sp 更容易。)
我不太清楚为什么即使进行了上述修改,SymPy 的特征值方法也会报告失败。如您所见 in the source,它所做的并不比上面的代码多多少:在特征多项式上调用 roots
。区别似乎在于创建此多项式的方式: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
returns{}
,没有找到根。看起来像是实施的缺陷。
我正在尝试计算大小为 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()
returns 我 3 个正确的复符号特征值 M
。但是,当我设置 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)
但它 returns 我 ConditionSet
无论如何,即使 eigenvals()
可以完成这项工作。
有谁知道如何正确解决这个特征值问题,对于 x
的任何值?
你对 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
可以更快地找到它们。此外,您可以制作一个 Poly 对象并对其应用 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))
(做 from sympy import *
比输入所有这些 sp 更容易。)
我不太清楚为什么即使进行了上述修改,SymPy 的特征值方法也会报告失败。如您所见 in the source,它所做的并不比上面的代码多多少:在特征多项式上调用 roots
。区别似乎在于创建此多项式的方式: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
returns{}
,没有找到根。看起来像是实施的缺陷。