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.

Solution 1

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.)

Solution 2

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.

Comments

[astrofrog]

I did go through very carefully, and also tried using a calculator, and I'm pretty confident I'm not doing anything wrong. Could this be a problem with the wikipedia article?

[John]

There it is. Can't believe I read through that several times.

[John]

See A. Rex's answer. As they say, "Well, THERE'S your problem!"

[astrofrog]

Thanks so much for spotting this! I was going crazy trying to understand why it wouldn't work.

[P i]

Seeing as you use Math.pow a lot, you might consider this extension: Number.prototype.pow = function(a) { return Math.pow(this, a); }; 2..pow(3) // 8

[Toby Speight]

Welcome to Stack Overflow! While this code snippet may solve the question, including an explanation really helps to improve the quality of your post. Remember that you are answering the question for readers in the future, not just the person asking now! Please edit your answer to add explanation, and give an indication of what limitations and assumptions apply.

[Cletrix]

No! any computer language can answer with correct algorithmic.

[SergeyT]

solve_deg3 for a=3.0, b=3.0, c=1.0, d=1.0 gives an incorrect result = -nan (it should be -1.0)