代码之家  ›  专栏  ›  技术社区  ›  Azlof

符号特征值的辛计算

  •  6
  • Azlof  · 技术社区  · 7 年前

    我试图计算符号复矩阵的特征值 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()
    

    返回3个正确的复数符号特征值 M . 然而,当我设置 x=1. ,我得到以下错误:

    我还试图计算特征值如下:

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

    但它给了我一个 ConditionSet 特征值()

    有人知道如何正确地解决这个特征值问题吗 x ?

    1 回复  |  直到 7 年前
        1
  •  3
  •   user6655984 user6655984    7 年前

    你对M的定义让symphy的日子过得太难过了,因为它引入了浮点数。当需要符号解时,应避免浮点运算。这意味着:

    • 而不是 1./3. sp.Rational(1, 3) (辛有理数)或 sp.S(1)/3 它具有相同的效果,但更容易键入。
    • 而不是 1j (Python的虚单位)使用 sp.I (辛的虚单位)
    • x = 1. x = 1 (Python 2.7 habits和SymPy一起表现不佳)。

    随着这些变化 solveset 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))
    

    from sympy import * 然后键入所有这些sp.)


    我不太清楚为什么即使经过上述修改,SymPy的特征值方法也会报告失败。正如你所见 in the source 关于特征多项式。不同之处似乎在于该多项式的创建方式: 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' . 随后 {} ,未找到根。看起来是实施的不足。