7
$\begingroup$

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.

  • 0
    Hmm, apparently you have an eigenvalue equation as your original problem; I would think it might be more convenient to use some canned routine for finding eigenvalues/eigenvectors and pick appropriately...2011-08-29

2 Answers 2

5

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 
  • 0
    @robert See Wikipedia article on finding roots of [cubic](http://en.wikipedia.org/wiki/Cubic_function) for variety of derivations of closed form solutions.2011-08-30
5

If your function $f$ is convex and increasing, or concave and decreasing, on an interval $[a,b]$ that contains a solution of $f(x) = 0$, Newton's method starting at $b$ will converge to a solution. Similarly if it is concave and increasing, or convex and decreasing, and you start at $a$. For a cubic it's easy to identify the critical points and inflection and thus find such an interval.

  • 0
    Ah, good point! Narrowed my thinking too much to my problem's domain, in which it couldn't be the case (the in$f$lection point is on the negative side o$f$ the root).2011-08-29