使用 Sympy 求解埃尔米特矩阵的特征值探针

Question

我正在尝试使用 Sympy 找到 Hermitian 矩阵的特征值。从线性代数中我们知道特征值应该是实数，但 Sympy 的输出总是包含无穷小的虚部。我知道乍一看这可能不是一个严重的问题，但是当我们尝试象征性地求解特征值时，它会带来很多麻烦。

from sympy import *
from sympy import init_printing
init_printing(use_latex=True) 
t12=symbols('t12')
t13=symbols('t13')
t16=symbols('t16')
t34=symbols('t34')
x=symbols('x')
y=symbols('y')
kx=symbols('kx')
ky=symbols('ky')
H0=Matrix([[0 for i in range(6)] for j in range(6)])

H0[0,1]=-t12
H0[0,2]=-t13
H0[0,5]=-t16*exp(-4*pi*I/3)
H0[1,3]=-t13
H0[1,4]=-t16*exp(4*pi*I/3)
H0[2,3]=-t34*exp(-4*pi*I/3)
H0[2,4]=-t13
H0[3,5]=-t13
H0[4,5]=-t12

H0[1,0]=-t12
H0[2,0]=-t13
H0[3,1]=-t13
H0[3,2]=-t34*exp(4*pi*I/3)
H0[4,1]=-t16*exp(-4*pi*I/3)
H0[4,2]=-t13
H0[5,0]=-t16*exp(4*pi*I/3)
H0[5,3]=-t13
H0[5,4]=-t12


eig=H0.evalf(5).eigenvals()

任何帮助，将不胜感激。

Answer 1

这是你的矩阵：

\n

In [6]: H0\nOut[6]: \n\xe2\x8e\xa1                                                                             2\xe2\x8b\x85\xe2\x85\x88\xe2\x8b\x85\xcf\x80\xe2\x8e\xa4\n\xe2\x8e\xa2                                                                             \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x8e\xa5\n\xe2\x8e\xa2                                                                               3  \xe2\x8e\xa5\n\xe2\x8e\xa2      0           -t\xe2\x82\x81\xe2\x82\x82          -t\xe2\x82\x81\xe2\x82\x83            0             0        -t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85\xe2\x84\xaf     \xe2\x8e\xa5\n\xe2\x8e\xa2                                                                                  \xe2\x8e\xa5\n\xe2\x8e\xa2                                                              -2\xe2\x8b\x85\xe2\x85\x88\xe2\x8b\x85\xcf\x80              \xe2\x8e\xa5\n\xe2\x8e\xa2                                                              \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80             \xe2\x8e\xa5\n\xe2\x8e\xa2                                                                 3                \xe2\x8e\xa5\n\xe2\x8e\xa2    -t\xe2\x82\x81\xe2\x82\x82            0             0           -t\xe2\x82\x81\xe2\x82\x83      -t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85\xe2\x84\xaf              0     \xe2\x8e\xa5\n\xe2\x8e\xa2                                                                                  \xe2\x8e\xa5\n\xe2\x8e\xa2                                                 2\xe2\x8b\x85\xe2\x85\x88\xe2\x8b\x85\xcf\x80                            \xe2\x8e\xa5\n\xe2\x8e\xa2                                                 \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80                            \xe2\x8e\xa5\n\xe2\x8e\xa2                                                   3                              \xe2\x8e\xa5\n\xe2\x8e\xa2    -t\xe2\x82\x81\xe2\x82\x83            0             0        -t\xe2\x82\x83\xe2\x82\x84\xe2\x8b\x85\xe2\x84\xaf           -t\xe2\x82\x81\xe2\x82\x83            0     \xe2\x8e\xa5\n\xe2\x8e\xa2                                                                                  \xe2\x8e\xa5\n\xe2\x8e\xa2                                  -2\xe2\x8b\x85\xe2\x85\x88\xe2\x8b\x85\xcf\x80                                          \xe2\x8e\xa5\n\xe2\x8e\xa2                                  \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80                                         \xe2\x8e\xa5\n\xe2\x8e\xa2                                     3                                            \xe2\x8e\xa5\n\xe2\x8e\xa2      0           -t\xe2\x82\x81\xe2\x82\x83      -t\xe2\x82\x83\xe2\x82\x84\xe2\x8b\x85\xe2\x84\xaf              0             0           -t\xe2\x82\x81\xe2\x82\x83    \xe2\x8e\xa5\n\xe2\x8e\xa2                                                                                  \xe2\x8e\xa5\n\xe2\x8e\xa2                     2\xe2\x8b\x85\xe2\x85\x88\xe2\x8b\x85\xcf\x80                                                        \xe2\x8e\xa5\n\xe2\x8e\xa2                     \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80                                                        \xe2\x8e\xa5\n\xe2\x8e\xa2                       3                                                          \xe2\x8e\xa5\n\xe2\x8e\xa2      0        -t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85\xe2\x84\xaf           -t\xe2\x82\x81\xe2\x82\x83            0             0           -t\xe2\x82\x81\xe2\x82\x82    \xe2\x8e\xa5\n\xe2\x8e\xa2                                                                                  \xe2\x8e\xa5\n\xe2\x8e\xa2      -2\xe2\x8b\x85\xe2\x85\x88\xe2\x8b\x85\xcf\x80                                                                      \xe2\x8e\xa5\n\xe2\x8e\xa2      \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80                                                                     \xe2\x8e\xa5\n\xe2\x8e\xa2         3                                                                        \xe2\x8e\xa5\n\xe2\x8e\xa3-t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85\xe2\x84\xaf              0             0           -t\xe2\x82\x81\xe2\x82\x83          -t\xe2\x82\x81\xe2\x82\x82            0     \xe2\x8e\xa6\n

\n

矩阵的特征值是其特征多项式的根。多项式根的符号表达方式存在局限性。我无法用显示的代码重现您所说的内容，因为我得到：

\n

In [3]: H0.eigenvals()\n---------------------------------------------------------------------------\nMatrixError: It is not always possible to express the eigenvalues of a matrix of size 5x5 or higher in radicals. We have CRootOf, but domains other than the rationals are not currently supported. If there are no symbols in the matrix, it should still be possible to compute numeric approximations of the eigenvalues using M.evalf().eigenvals() or M.charpoly().nroots().\n

\n

该错误消息是

\n

\n

大小为 5x5 或更大的矩阵的特征值并不总是能够用根式表示。我们有 CRootOf，但目前不支持有理数以外的域。如果矩阵中没有符号，仍然可以使用 M.evalf().eigenvals() 或 M.charpoly().nroots() 计算特征值的数值近似值。

\n

\n

原因是您有一个 6x6 矩阵，其特征多项式为 6 次，因此特征值是 6 次多项式的根。由于阿贝尔 - 5 次或以上的多项式可能无法用根式求解鲁菲尼定理：\n https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem

\n

我的猜测是，您看到这一点是因为您正在用符号替换值，然后计算特征值，例如如下所示：

\n

In [18]: H0.subs({t12:1, t13:2, t16:3, t34:4}).evalf().eigenvals()\nOut[18]: \n{-6.21170948161553 + 2.04287290941429e-64\xe2\x8b\x85\xe2\x85\x88: 1, -3.60555127546399 - 1.30254464798347e-64\xe2\x8b\x85\xe2\x85\x88: 1, -0.6\n43945118785511 + 1.4910290487412e-64\xe2\x8b\x85\xe2\x85\x88: 1, 0.643945118785511 - 2.93808243220787e-64\xe2\x8b\x85\xe2\x85\x88: 1, 3.6055512\n7546399 + 4.04818363587083e-64\xe2\x8b\x85\xe2\x85\x88: 1, 6.21170948161553 - 3.2972673293568e-64\xe2\x8b\x85\xe2\x85\x88: 1}\n

\n

这是因为您的矩阵中有复数，因此根是使用复浮点计算的，然后会因舍入误差而产生小的虚部。

\n

有多种解决此问题的方法，这些方法可能适用于不同的示例，也可能不适用于不同的示例。一种方法是首先精确计算特征多项式（不带浮点数），然后才要求求根：

\n

In [22]: l = symbols(\'lambda\')\n\nIn [23]: p = H0.subs({t12:1, t13:2, t16:3, t34:4}).charpoly(l).as_expr()\n\nIn [24]: p\nOut[24]: \n 6       4    2 \xe2\x8e\x9b             2/3      3 ____\xe2\x8e\x9e                2/3      3 ____\n\xce\xbb  - 52\xe2\x8b\x85\xce\xbb  + \xce\xbb \xe2\x8b\x85\xe2\x8e\x9d434 - 89\xe2\x8b\x85(-1)    + 89\xe2\x8b\x85\xe2\x95\xb2\xe2\x95\xb1 -1 \xe2\x8e\xa0 - 224 - 16\xe2\x8b\x85(-1)    + 16\xe2\x8b\x85\xe2\x95\xb2\xe2\x95\xb1 -1 \n\nIn [25]: nroots(p)\nOut[25]: \n[-6.21170948161553, -3.60555127546399, -0.643945118785511, 0.643945118785511, 3.60555127546399, 6.2\n1170948161553]\n

\n

这里nroots可以看出，根是完全实数，因为多项式的系数是实数。需要精确的算术才能达到这一点，以便确定系数是实数，即使它们的表达式涉及复数。您还可以计算特征多项式，然后仅替换这些值，因为获取多项式是其中最慢的部分。

\n

给定一个具有（精确）有理系数的多项式，SymPy 始终可以使用 RootOf 表示其根。你的系数是代数系数而不是有理系数，但这意味着它仍然可以用正确的参数进行因式分解：

\n

In [26]: pf = p.factor(extension=True)\n\nIn [27]: pf\nOut[27]: \n\xe2\x8e\x9b 2     \xe2\x8e\x9e \xe2\x8e\x9b 4       2     \xe2\x8e\x9e\n\xe2\x8e\x9d\xce\xbb  - 13\xe2\x8e\xa0\xe2\x8b\x85\xe2\x8e\x9d\xce\xbb  - 39\xe2\x8b\x85\xce\xbb  + 16\xe2\x8e\xa0\n

\n

这里的extension=True意思是构造一个可包含多项式系数的扩展域，然后在该域上因式分解多项式环。在这种情况下，结果足以对多项式进行因式分解。现在我们只有具有有理系数的多项式，因此所有根都可以使用 RootOf 来表示：

\n

In [28]: real_roots(pf)\nOut[28]: \n\xe2\x8e\xa1       \xe2\x8e\x9b 4       2        \xe2\x8e\x9e               \xe2\x8e\x9b 4       2        \xe2\x8e\x9e         \xe2\x8e\x9b 4       2        \xe2\x8e\x9e       \n\xe2\x8e\xa3CRootOf\xe2\x8e\x9dx  - 39\xe2\x8b\x85x  + 16, 0\xe2\x8e\xa0, -\xe2\x88\x9a13, CRootOf\xe2\x8e\x9dx  - 39\xe2\x8b\x85x  + 16, 1\xe2\x8e\xa0, CRootOf\xe2\x8e\x9dx  - 39\xe2\x8b\x85x  + 16, 2\xe2\x8e\xa0, \xe2\x88\x9a13, \n\n       \xe2\x8e\x9b 4       2        \xe2\x8e\x9e\xe2\x8e\xa4\nCRootOf\xe2\x8e\x9dx  - 39\xe2\x8b\x85x  + 16, 3\xe2\x8e\xa0\xe2\x8e\xa6\n\nIn [29]: [r.n() for r in real_roots(pf)]\nOut[29]: \n[-6.21170948161553, -3.60555127546399, -0.643945118785511, 0.643945118785511, 3.60555127546399, 6.2\n1170948161553]\n

\n

上面的因式分解可能无需用以下符号替换值：

\n

In [33]: p = H0.charpoly(l).as_expr()\n\nIn [34]: p\nOut[34]: \n 6    4 \xe2\x8e\x9b       2        2        2      2\xe2\x8e\x9e    2 \xe2\x8e\x9b   4        2    2       2/3    2    2   3 ____  \n\xce\xbb  + \xce\xbb \xe2\x8b\x85\xe2\x8e\x9d- 2\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82  - 4\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83  - 2\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86  - t\xe2\x82\x83\xe2\x82\x84 \xe2\x8e\xa0 + \xce\xbb \xe2\x8b\x85\xe2\x8e\x9dt\xe2\x82\x81\xe2\x82\x82  + 4\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83  - (-1)   \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86  + \xe2\x95\xb2\xe2\x95\xb1 -1 \xe2\x8b\x85t\n\n  2    2        2    2         2/3        2         3 ____        2             2/3        2       \n\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86  + 2\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  - 4\xe2\x8b\x85(-1)   \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 + 4\xe2\x8b\x85\xe2\x95\xb2\xe2\x95\xb1 -1 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 - 2\xe2\x8b\x85(-1)   \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 + \n\n  3 ____        2            4        2    2        2              4        2    2\xe2\x8e\x9e      4    2    \n2\xe2\x8b\x85\xe2\x95\xb2\xe2\x95\xb1 -1 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 + 4\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83  + 4\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86  - 4\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 + t\xe2\x82\x81\xe2\x82\x86  + 2\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 \xe2\x8e\xa0 - t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  - 2\n\n 3 ____    3    2             2/3    3    2            2    4         2/3    2    2             3 _\n\xe2\x8b\x85\xe2\x95\xb2\xe2\x95\xb1 -1 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 + 2\xe2\x8b\x85(-1)   \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 - 4\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83  - 2\xe2\x8b\x85(-1)   \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 + 2\xe2\x8b\x85\xe2\x95\xb2\xe2\x95\xb1 \n\n___    2    2           3 ____    2    2    2       2/3    2    2    2     3 ____        4         \n-1 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 - \xe2\x95\xb2\xe2\x95\xb1 -1 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  + (-1)   \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  - 4\xe2\x8b\x85\xe2\x95\xb2\xe2\x95\xb1 -1 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 + 4\xe2\x8b\x85\n\n    2/3        4             2/3        2    2         3 ____        2    2            4    2      \n(-1)   \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 - 2\xe2\x8b\x85(-1)   \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 + 2\xe2\x8b\x85\xe2\x95\xb2\xe2\x95\xb1 -1 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 - 4\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86  + 4\xe2\x8b\x85t\n\n  2    3          4    2\n\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 - t\xe2\x82\x81\xe2\x82\x86 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 \n\nIn [35]: %time pf = p.factor(extension=True)\n...\n

\n

虽然最后一步很慢，所以我不知道需要多长时间，甚至不知道它是否会成功，因为不能保证多项式因式分解。如果它确实成功并给出了二次和四次的乘积，那么应该可以用根号来象征性地表达根，而无需用值替换符号。不过，由于不可约原因，这些公式也可能存在虚部较小的问题：\n https://en.wikipedia.org/wiki/Casus_irreducibilis

\n

编辑：我现在看到，H0.evalf().eigenvals()实际上确实返回了一些涉及符号和大量小浮点数的根的复杂表达式。无论如何，我所说的所有原因都会发生这种情况，因为在调用 evalf 后矩阵中有复数。

\n

这里还有一种获得特征值的精确表示的方法。首先在您的代码中使用符号a代替exp(4*pi*I/3)，同样1/a适用于exp(-4*pi*I/3)。然后你可以找到特征多项式并对其进行因式分解：

\n

In [78]: p = H0.charpoly(l).as_expr().factor().as_numer_denom()[0]\n\nIn [79]: p\nOut[79]: \n \xe2\x8e\x9b   2              2        2        2          \xe2\x8e\x9e \xe2\x8e\x9b   2  2              2        2        2       \n-\xe2\x8e\x9d- a \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 - a\xe2\x8b\x85\xce\xbb  + a\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82  + a\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86  - t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86\xe2\x8e\xa0\xe2\x8b\x85\xe2\x8e\x9d- a \xe2\x8b\x85\xce\xbb \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 - 2\xe2\x8b\x85a \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 + a \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\n\n     2      4      2    2        2    2      2    2      2    2        2    2          4          2\n\xe2\x82\x86\xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  + a\xe2\x8b\x85\xce\xbb  - a\xe2\x8b\x85\xce\xbb \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82  - 4\xe2\x8b\x85a\xe2\x8b\x85\xce\xbb \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83  - a\xe2\x8b\x85\xce\xbb \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86  - a\xe2\x8b\x85\xce\xbb \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  + a\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  + 4\xe2\x8b\x85a\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83  - 4\xe2\x8b\x85a\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \n\n                2    2    2                    2                  2\xe2\x8e\x9e\n\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 + a\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  - \xce\xbb \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 - 2\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 + t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 \xe2\x8e\xa0\n

\n

您可以使用它roots(p, l)来获得一些相当复杂的根表达式。不过我们想代入a回去，这样就可以进行因式分解：

\n

In [96]: pf = Mul(*(f.factor(extension=True).collect(l) for f in Mul.make_args(p.subs(a, exp(4*pi*I/3\n    ...: ))))).as_independent(l)[1]\n\nIn [97]: pf\nOut[97]: \n\xe2\x8e\x9b   2      2                2\xe2\x8e\x9e \xe2\x8e\x9b 4    2 \xe2\x8e\x9b     2                  2      2      2\xe2\x8e\x9e      2    2      \n\xe2\x8e\x9d- \xce\xbb  + t\xe2\x82\x81\xe2\x82\x82  + t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 + t\xe2\x82\x81\xe2\x82\x86 \xe2\x8e\xa0\xe2\x8b\x85\xe2\x8e\x9d\xce\xbb  + \xce\xbb \xe2\x8b\x85\xe2\x8e\x9d- t\xe2\x82\x81\xe2\x82\x82  + t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 - 4\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83  - t\xe2\x82\x81\xe2\x82\x86  - t\xe2\x82\x83\xe2\x82\x84 \xe2\x8e\xa0 + t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  + 2\xe2\x8b\x85t\n\n      2                  2        4        2              2    2\xe2\x8e\x9e\n\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 - t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  + 4\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83  - 4\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 + t\xe2\x82\x81\xe2\x82\x86 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 \xe2\x8e\xa0\n

\n

这些可以为根提供更好的表达式，因为我们已经消除了所有复数：

\n

In [98]: r1, r2, r3, r4, r5, r6 = roots(pf, l)\n\nIn [99]: r1\nOut[99]: \n    _______________________\n   \xe2\x95\xb1    2                2 \n-\xe2\x95\xb2\xe2\x95\xb1  t\xe2\x82\x81\xe2\x82\x82  + t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 + t\xe2\x82\x81\xe2\x82\x86  \n\nIn [100]: r2\nOut[100]: \n   _______________________\n  \xe2\x95\xb1    2                2 \n\xe2\x95\xb2\xe2\x95\xb1  t\xe2\x82\x81\xe2\x82\x82  + t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 + t\xe2\x82\x81\xe2\x82\x86  \n\nIn [101]: r3\nOut[101]: \n       ____________________________________________________________________________________________\n      \xe2\x95\xb1                                            ________________________________________________\n     \xe2\x95\xb1     2                         2      2     \xe2\x95\xb1    4        3            2    2        2    2  \n    \xe2\x95\xb1   t\xe2\x82\x81\xe2\x82\x82    t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86        2   t\xe2\x82\x81\xe2\x82\x86    t\xe2\x82\x83\xe2\x82\x84    \xe2\x95\xb2\xe2\x95\xb1  t\xe2\x82\x81\xe2\x82\x82  - 2\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 + 8\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83  + 3\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86  -\n-  \xe2\x95\xb1    \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80 - \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80 + 2\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83  + \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80 + \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80 - \xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\n \xe2\x95\xb2\xe2\x95\xb1      2        2                2      2                                                        \n\n___________________________________________________________________________________________________\n___________________________________________________________________________________________________\n      2    2            2                2                3                2        2    2         \n 2\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  - 8\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86 - 8\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84 - 2\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86  + 2\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x82\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86\xe2\x8b\x85t\xe2\x82\x83\xe2\x82\x84  + 8\xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x83 \xe2\x8b\x85t\xe2\x82\x81\xe2\x82\x86  + 16\xe2\x8b\x85t\xe2\x82\x8