Fast and robust root of a cubic polynomial with constraints

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.

• It's possible to modify Newton-Raphson such that the iteration is forced back on track if it, say, starts returning a negative root. It is also possible to leverage the cubic formula so that only the positive root is returned (but you'll want to check the discriminant first). Maybe you can post the actual cubic equation you have? – J. M. is a poor mathematician Aug 28 '11 at 19:30
• This is an example of a problematic case: $0.69357258x^3-1.7485783x^2-49.538208x-89.828194$, but of course I'm interested in a solution working for any cubic in a form described above. Could you say something more about the modified N-R or the cubic root formula? – robert Aug 28 '11 at 19:54
• I was sort of thinking that the $p,q,r$ are in fact expressions you plug the numbers you have into... e.g. was the 0.69357258 actual input or derived from something else? – J. M. is a poor mathematician Aug 28 '11 at 20:24
• Ah, sorry for the confusion: I wrote the original cubic in the form $x^3+px^2+qx+r$ but in the example in the comment as $ax^3+bx^2+cx+d$. It's meant to be the same thing, though: $x^3 -2.5211179x^2-71.424690x-129.51520$. So $p,q,r$ are scalars, inputs to the problem. I know the input polynomial with these coefficients has the described properties, because I know the geometrical configuration it has been derived from (described in this question) – robert Aug 28 '11 at 20:39
• 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... – J. M. is a poor mathematician Aug 29 '11 at 1:21

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:= PositiveCubicRoot[-2.52111798, -71.424692, -129.51520]

Out= 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:=
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:=
NewtonRaphson[p_, q_, r_] :=
NewtonRaphsonIterate[{p, q, r}, NewtonRaphsonStartingPoint[{p, q, r}]]

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

Out= 10.499
• It's hard to make a pronouncement on whether the analytical or iterative approach is faster/more convenient, at least not without testing. – J. M. is a poor mathematician Aug 29 '11 at 1:23
• I am sure iterative method will be faster in this case. – Sasha Aug 29 '11 at 1:27
• For the current setup, I suppose yes; you are after all using two transcendentals in addition to the radical. – J. M. is a poor mathematician Aug 29 '11 at 1:43
• I run some quick tests with a C application compiled with msvc running on an i7 comparing 10 iterations of my current version of N-R (without any improvements) to the closed-form solution you posted. The latter is somewhat faster and seems robust enough - I couldn't break it with any tests I tried so far, which is very cool :) It's speed advantage would be even bigger compared to a patched N-R solution using binary search (slow convergence) or better initial guess (larger root of the derivative + margin?). I'll polish both methods and we'll see the final verdict. – robert Aug 29 '11 at 22:41
• Could you post an explanation/link on how the closed-form solution was derived? Again, thanks for help everyone :) – robert Aug 29 '11 at 22:44

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.

• This is of course true in the mathematical sense, but doesn't always hold when performing operations with limited precision. But yes, it might be helpful to use a closed form formula to find the larger root of the polynomial's derivative, which is just a quadratic, increase it by a margin and use as a better starting point if it's larger than the original guess. Btw, $x^3 + px^2 + qx + r$ for any $p, q, r$ is convex and increasing around the largest root. – robert Aug 29 '11 at 15:37
• Try $x^3 - x + 6$, which has one real root (namely $-2$) and is concave there. – Robert Israel Aug 29 '11 at 18:58
• Ah, good point! Narrowed my thinking too much to my problem's domain, in which it couldn't be the case (the inflection point is on the negative side of the root). – robert Aug 29 '11 at 22:55