Sunday 15 February 2015

python - Solving a cubic equation -


I am writing as part of a program, I need to fix a cubic equation (a numerical root Finder):

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

I'm trying However, consider the following code (this is Python but this is very common code):

  a = 1.0 b = 0.0c = 0.2 - 1.0 d = -0.7 * 0.2 q = (3 * a * c - b ** 2) / (9 * one ** 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 the delta is less than zero, so we use other equations from the article: rho = (-q ** 3) ** The hypothetical portion for ** 0.5 # x1 is insignificant because it is s_real = rho ** (1./3.) T_real = rho ** (1./3 print "s [real] = ", S_real print" t [real] = ", t_lr 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 output is:

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

Then the output is not zero, and therefore x1 is not really a solution in Wikipedia article Mistake?

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 it to work on multiplier Arrays of.

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

, I make these changes in 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 also cos .)


No comments:

Post a Comment