求解三次方程

Question

作为我正在编写的程序的一部分，我需要精确求解三次方程（而不是使用数值求根器）：

 a*x**3 + b*x**2 + c*x + d = 0.

我正在尝试使用 此处 的方程式。但是，请考虑以下代码（这是 Python，但它是非常通用的代码）：

 a =  1.0
b =  0.0
c =  0.2 - 1.0
d = -0.7 * 0.2

q = (3*a*c - b**2) / (9 * a**2)
r = (9*a*b*c - 27*a**2*d - 2*b**3) / (54*a**3)

print "q = ",q
print "r = ",r

delta = q**3 + r**2

print "delta = ",delta

# here delta is less than zero so we use the second set of equations from the article:

rho = (-q**3)**0.5

# For x1 the imaginary part is unimportant since it cancels out
s_real = rho**(1./3.)
t_real = rho**(1./3.)

print "s [real] = ",s_real
print "t [real] = ",t_real

x1 = s_real + t_real - b / (3. * a)

print "x1 = ", x1

print "should be zero: ",a*x1**3+b*x1**2+c*x1+d

但输出是：

 q =  -0.266666666667
r =  0.07
delta =  -0.014062962963
s [real] =  0.516397779494
t [real] =  0.516397779494
x1 =  1.03279555899
should be zero:  0.135412149064

所以输出不为零，所以 x1 实际上不是解。维基百科文章有错误吗？

ps：我知道 numpy.roots 会解决这种方程，但我需要为数百万个方程做这个，所以我需要实现它来处理系数数组。

原文由 astrofrog 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 1

维基百科的符号 (rho^(1/3), theta/3) 并不意味着 rho^(1/3) 是实部， theta/3 是虚部。相反，这是在极坐标中。因此，如果您想要实部，您将采用 rho^(1/3) * cos(theta/3) 。

我对您的代码进行了这些更改并且对我有用：

 theta = arccos(r/rho)
s_real = rho**(1./3.) * cos( theta/3)
t_real = rho**(1./3.) * cos(-theta/3)

（当然， s_real = t_real 在这里是因为 cos 是偶数。）

原文由 A. Rex 发布，翻译遵循 CC BY-SA 2.5 许可协议

Answer 2

我看过维基百科文章和您的程序。

我还使用 Wolfram Alpha 求解了方程式，结果与你得到的不匹配。

我只是在每一步都通过你的程序，使用大量打印语句，并获得每个中间结果。然后通过计算器自己完成。

我找不到发生了什么，但你的手算和程序的分歧是一个很好的观察点。

原文由 John 发布，翻译遵循 CC BY-SA 2.5 许可协议