[Math] Fast and robust root of a cubic polynomial with constraints

cubicsnumerical methodspolynomialsroots

I'm looking for a fast and robust method for finding a root of a cubic polynomial

$$x^3 + px^2 + qx + r$$

To make the search more robust and faster, I'd like to leverage these properties:

  1. The polynomial has exactly one real positive root (two other roots are either complex or real and negative).
  2. Only the value of the positive root is needed.
  3. There's a decent initial guess on the value of the root.

So far my approach was to directly apply Newton's method to the function using the initial guess and that would give me a decent result in just a couple of iterations.

However in some cases the iteration would cause the current guess for the root to jump closer to one of the negative roots and the method would incorrectly start converging towards that root instead. While it's possible to detect this situation, it's hard to bring the iteration back on the right track.

There's an interesting article about solving quartics and cubics and a also an example implementation, but the methods are very generic, too slow for my needs and have robustness issues of their own.

It would be interesting to know if there's a way to make my iterative search more robust or possibly that there's a faster analytical method taking advantage of the extra properties.

Best Answer

I think you can avoid the Newton-Raphson altogether, since cubic is solvable in radicals.

Here is the complete algorithm (working under constraints outlined in the problem), in Mathematica:

Clear[PositiveCubicRoot];
PositiveCubicRoot[p_, q_, r_] := 
 Module[{po3 = p/3, a, b, det, abs, arg},
  b = ( po3^3 - po3 q/2 + r/2);
  a = (-po3^2 + q/3);
  det = a^3 + b^2;
  If[det >= 0,
   det = Power[Sqrt[det] - b, 1/3];
   -po3 - a/det + det
   ,
   (* evaluate real part, imaginary parts cancel anyway *)
   abs = Sqrt[-a^3];
   arg = ArcCos[-b/abs];
   abs = Power[abs, 1/3];
   abs = (abs - a/abs);
   arg = -po3 + abs*Cos[arg/3]
   ]
  ]

Then

In[222]:= PositiveCubicRoot[-2.52111798, -71.424692, -129.51520]

Out[222]= 10.499

However, if the Newton-Raphson method must be used, then a good initial guess is imperative. Binary division is a good method to isolate the root in the case at hand. We start with an arbitrary point $x$, I chose $x=1$, and double it while the polynomial at that point is negative. Then do binary division a certain number of times (the code below does it twice). Ultimately, polish it off with Newton-Raphson iterations:

In[283]:= 
NewtonRaphsonStartingPoint[{p_, q_, r_}] := Module[{x1=0, x2=1,f1,f2,xm,fm},
   f1 = r + x1 (q + (p + x1) x1);
   While[(f2 = r + x2 (q + (p + x2) x2)) <= 0, 
      x1 = x2; f1 = f2; x2 = 2 x2];
   Do[xm = (x1 + x2)/2; fm = r + xm (q + (p + xm) xm); 
    If[fm <= 0, f1 = fm; x1 = xm, f2 = fm; x2 = xm], {i, 2}];
   (f2 x2 - f1 x1)/(f2 - f1)
];
NewtonRaphsonIterate[{p_, q_, r_}, x0_Real] := 
 FixedPoint[
  Function[x, x - (r + x (q + (p + x) x))/(q + (2 p + 3 x) x)], x0]

In[285]:= 
NewtonRaphson[p_, q_, r_] := 
 NewtonRaphsonIterate[{p, q, r}, NewtonRaphsonStartingPoint[{p, q, r}]]

In[286]:= NewtonRaphson[-2.52111798, -71.424692, -129.51520]

Out[286]= 10.499
Related Question