求解三次方程

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并不是一个真正的解。难道维基百科的文章有错误吗？

补充一下：我知道numpy.roots可以解决这种方程，但我需要处理数百万个方程，所以我需要实现一个可以在系数数组上工作的方案。

Answer 1

在编程中，有时候我们需要让程序在特定的条件下执行某些操作。这就像给程序设定了一些规则，只有当这些规则被满足时，程序才会继续进行。

比如说，你可以设定一个条件，如果用户输入的密码是正确的，程序才会允许他们登录。这种方式可以帮助我们控制程序的行为，确保它按照我们预期的方式运行。

有些编程语言提供了简单的方式来检查这些条件，比如使用“如果”语句（if statement）。这就像是在问：“如果这个条件成立，我就做这个事情。”

通过这种方式，我们可以让程序变得更加智能，能够根据不同的情况做出不同的反应。

from cmath import *


def linear(a, b):
    solutions = set()
    if a == 0 and b == 0:
        solutions.add(True)

    if a == 0 and b != 0:
        solutions.add(False)

    if a != 0:
        solutions.add(-b / a)
    return solutions


def quadratic(a, b, c):
    solutions = set()
    if a != 0:
        D = b ** 2 - 4 * a * c
        x1 = (-b + sqrt(D)) / (2 * a)
        x2 = (-b - sqrt(D)) / (2 * a)
        solutions.update({x1, x2})
    else:
        solutions.update(linear(b, c))
    return solutions


def cubic(a, b, c, d):
    solutions = set()
    if a != 0:
        p = (3 * a * c - b ** 2) / (3 * a ** 2)
        q = (2 * b ** 3 - 9 * a * b * c + 27 * a ** 2 * d) / (27 * a ** 3)
        for n in range(3):
            solutions.add((2 * sqrt(-p / 3) * cos(acos((-3 * q) * sqrt(-3 * p) / (2 * p ** 2)) / 3 + 2 * pi * n / 3))
                          - (b / (3 * a)))
    else:
        solutions.update(quadratic(b, c, d))
    return solutions


print('ax^3+bx^2+cx+d=0')
a, b, c, d = map(complex, input().split())

print(cubic(a, b, c, d))

Answer 2

这是A. Rex在JavaScript中的解决方案：

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

q = (3*a*c - Math.pow(b, 2)) / (9 * Math.pow(a, 2));
r = (9*a*b*c - 27*Math.pow(a, 2)*d - 2*Math.pow(b, 3)) / (54*Math.pow(a, 3));
console.log("q = "+q);
console.log("r = "+r);

delta = Math.pow(q, 3) + Math.pow(r, 2);
console.log("delta = "+delta);

// here delta is less than zero so we use the second set of equations from the article:
rho = Math.pow((-Math.pow(q, 3)), 0.5);
theta = Math.acos(r/rho);

// For x1 the imaginary part is unimportant since it cancels out
s_real = Math.pow(rho, (1./3.)) * Math.cos( theta/3);
t_real = Math.pow(rho, (1./3.)) * Math.cos(-theta/3);

console.log("s [real] = "+s_real);
console.log("t [real] = "+t_real);

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

console.log("x1 = "+x1);
console.log("should be zero: "+(a*Math.pow(x1, 3)+b*Math.pow(x1, 2)+c*x1+d));

Answer 3

维基百科的表示法 (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 是偶函数。）