Solving a cubic equation

Question

As part of a program I'm writing, I need to solve a cubic equation exactly (rather than using a numerical root finder):

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

I'm trying to use the equations from here. However, consider the following code (this is Python but it's pretty generic code):

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

But the output is:

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

so the output is not zero, and so x1 isn't actually a solution. Is there a mistake in the Wikipedia article?

ps: I know that numpy.roots will solve this kind of equation but I need to do this for millions of equations and so I need to implement this to work on arrays of coefficients.

Answer 1

I've looked at the Wikipedia article and your program.

I also solved the equation using Wolfram Alpha and the results there don't match what you get.

I'd just go through your program at each step, use a lot of print statements, and get each intermediate result. Then go through with a calculator and do it yourself.

I can't find what's happening, but where your hand calculations and the program diverge is a good place to look.

Answer 2

Wikipedia's notation (rho^(1/3), theta/3) does not mean that rho^(1/3) is the real part and theta/3 is the imaginary part. Rather, this is in polar coordinates. Thus, if you want the real part, you would take rho^(1/3) * cos(theta/3).

I made these changes to your code and it worked for me:

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

(Of course, s_real = t_real here because cos is even.)

Answer 3

from cmath import *

solutions = set()


def cbrt(polynomial):
    solution = set()
    root1 = polynomial ** (1 / 3)
    root2 = (polynomial ** (1 / 3)) * (-1 / 2 + (sqrt(3) * 1j) / 2)
    root3 = (polynomial ** (1 / 3)) * (-1 / 2 - (sqrt(3) * 1j) / 2)
    solution.update({root1, root2, root3})
    return solution


def cardano(a, b, c, d):
    p = (3 * a * c - b ** 2) / (3 * a ** 2)
    q = (2 * b ** 3 - 9 * a * b * c + 27 * a ** 2 * d) / (27 * a ** 3)
    alpha = cbrt(-q / 2 + sqrt((q / 2) ** 2 + (p / 3) ** 3))
    beta = cbrt(-q / 2 - sqrt((q / 2) ** 2 + (p / 3) ** 3))
    for i in alpha:
        for j in beta:
            if abs((i * j) + p / 3) <= 0.00001:
                x = i + j - b / (3 * a)
                solutions.add(x)


def quadratic(a, b, c):
    D = b ** 2 - 4 * a * c
    x1 = (-b + sqrt(D)) / 2 * a
    x2 = (-b - sqrt(D)) / 2 * a
    solutions.update({x1, x2})


def linear(a, b):
    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)


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

if a != 0:
    cardano(a, b, c, d)

elif b != 0:
    quadratic(b, c, d)

else:
    linear(c, d)

print(solutions)