三次方程的 Sympy 解似乎是错误的并且包含虚数单位。怎么了？

Question

from sympy import *
from sympy.abc import x

f = x**3 - 3*x + 1
res = solve(f)
print(res)

但是，它给出的虚数单位答案不正确。检查 viasubs也不返回 0。

更新：

资源：

[-3/((-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)/3,
 -(-1/2 + sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)/3 - 3/((-1/2 + sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)),
 -(27/2 + 27*sqrt(3)*I/2)**(1/3)/3 - 3/(27/2 + 27*sqrt(3)*I/2)**(1/3)]

子检查：

>>> print(f.subs(x, res[0]))
1 + (-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3) + (-3/((-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)/3)**3 + 9/((-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3))

剧情：

>>> plot(f, xlim=(-2, 2), ylim=(-2, 2))

Answer 1

让我们尝试将 sympy 计算的符号解转换为浮点数：

from sympy import *
var("x")
f = x**3 - 3*x + 1
sol = solve(f)
print([s.n() for s in sol])
# [0.347296355333861 - 0.e-23*I, 1.53208888623796 + 0.e-20*I, -1.87938524157182 + 0.e-23*I]

在这里你可以看到非常小的虚部。这些是舍入误差。我们可以这样删除它们：

print([s.n(chop=True) for s in sol])
# [0.347296355333861, 1.53208888623796, -1.87938524157182]

现在，我将把第一个符号解代入方程：

print(f.subs(x, sol[0]))
# 1 + (-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3) + (-3/((-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)/3)**3 + 9/((-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3))

它看起来并不为零，但我们可以要求 sympy 简化它，如下所示：

print(f.subs(x, sol[0]).simplify())
# 0

正如您所看到的，sympy 计算的符号解是正确的。

Answer 2

另一个答案很好，表明 SymPy 给出的表达式是正确的，即使它们涉及复数。我只是想澄清一下为什么solve返回这个结果。

\n

给出这些表达式的原因是因为由于不可约原因，在不使用复数的情况下不可能用根式来表达给定多项式的根：

\n

https://en.wikipedia.org/wiki/Casus_irreducibilis

\n

SymPy 提供了两种简单的方法来避免在寻求不可约三次方根的精确表示时出现不可约原因。一种是要求涉及三角函数的表达式而不是纯根式表达式：

\n

In [2]: roots(x**3 - 3*x + 1, x, trig=True)\nOut[2]: \n\xe2\x8e\xa7      \xe2\x8e\x9b\xcf\x80\xe2\x8e\x9e          \xe2\x8e\x9b2\xe2\x8b\x85\xcf\x80\xe2\x8e\x9e          \xe2\x8e\x9b4\xe2\x8b\x85\xcf\x80\xe2\x8e\x9e   \xe2\x8e\xab\n\xe2\x8e\xa8-2\xe2\x8b\x85cos\xe2\x8e\x9c\xe2\x94\x80\xe2\x8e\x9f: 1, 2\xe2\x8b\x85cos\xe2\x8e\x9c\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x8e\x9f: 1, 2\xe2\x8b\x85cos\xe2\x8e\x9c\xe2\x94\x80\xe2\x94\x80\xe2\x94\x80\xe2\x8e\x9f: 1\xe2\x8e\xac\n\xe2\x8e\xa9      \xe2\x8e\x9d9\xe2\x8e\xa0          \xe2\x8e\x9d 9 \xe2\x8e\xa0          \xe2\x8e\x9d 9 \xe2\x8e\xa0   \xe2\x8e\xad\n

\n

这里的三角表达式-2*cos(pi/9)是纯实数，并且可以精确地表示解，而无需在表达式中使用中间复数。并不总是能够像这样表达任何多项式的根（取决于例如多项式的次数），就像并不总是能够表示根式的根一样（参见阿贝尔-鲁菲尼定理）。

\n

表示不可约多项式的根的更通用方法是使用 RootOf：

\n

In [3]: real_roots(x**3 - 3*x + 1)\nOut[3]: \n\xe2\x8e\xa1       \xe2\x8e\x9b 3             \xe2\x8e\x9e         \xe2\x8e\x9b 3             \xe2\x8e\x9e         \xe2\x8e\x9b 3             \xe2\x8e\x9e\xe2\x8e\xa4\n\xe2\x8e\xa3CRootOf\xe2\x8e\x9dx  - 3\xe2\x8b\x85x + 1, 0\xe2\x8e\xa0, CRootOf\xe2\x8e\x9dx  - 3\xe2\x8b\x85x + 1, 1\xe2\x8e\xa0, CRootOf\xe2\x8e\x9dx  - 3\xe2\x8b\x85x + 1, 2\xe2\x8e\xa0\xe2\x8e\xa6\n\nIn [4]: [r.n() for r in real_roots(x**3 - 3*x + 1)]\nOut[4]: [-1.87938524157182, 0.347296355333861, 1.53208888623796]\n

\n

有些人不喜欢 RootOf，因为它在某种程度上并不“明确”，但它是用符号表示多项式根的最佳方式。RootOf 表示结合了代数知识（该数字是给定多项式的精确根）、根是否实数的精确确定以及计算根到任何所需精度的数值能力：

\n

In [5]: for r in real_roots(x**3-3*x+1):\n   ...:     print(r, "=", r.evalf(50))\n   ...: \nCRootOf(x**3 - 3*x + 1, 0) = -1.8793852415718167681082185546494629398724162685289\nCRootOf(x**3 - 3*x + 1, 1) = 0.34729635533386069770343325353862959200075135436814\nCRootOf(x**3 - 3*x + 1, 2) = 1.5320888862379560704047853011108333478716649141608\n

\n

使用 RootOf，我们总是知道任何根是否是实数，因此当对 RootOf 进行数值计算时，不存在小的虚部。对于非实根，我们确切地知道哪些根是彼此共轭的等。代数数的代数闭域的每个元素都可以用 RootOf 表示形式表示，这意味着具有任何次数的有理系数的任何多项式的所有根。有了RootOf，我们可以同时避免不可约原因、阿贝尔-鲁菲尼定理和卡尔达诺公式的混乱的局限性。

\n

solve如果 SymPy\ 的函数在这里返回 RootOf 而不是非实部首表达式会更好，但对于许多 SymPy 用户来说，它可能比 casus irreducibilis 更令人困惑。当前solve将返回RootOf，但仅当roots未能找到所有根的根式表达式时才返回。

\n