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.
(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)
s_real = t_real here because
cos is even.)
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.