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
Best Answer
I'm going to post the method involving trig functions. As an example, we will be using $0=x^3-3x^2+3$.
To use this method, you must know that
$$\cos(3\theta)=4\cos^3(\theta)-3\cos(\theta)$$
which is easily derived from the sum of angles formulas and the Pythagorean theorem.
We manipulate this to give us
$$\cos(3\arccos(r))=4r^3-3r\tag0$$
Anyways, we start with
$$0=ax^3+bx^2+cx+d$$
$$0=x^3-3x^2+3$$
We always start with the substitution $x=y-\frac{b}{3a}$, which for our example is $x=y+1$
$$0=(y+1)^3-3(y+1)^2+3\\=y^3-3y+1$$
In general, we get something along the lines of
$$0=a\left(y-\frac{b}{3a}\right)^3+b\left(y-\frac{b}{3a}\right)^2+c\left(y-\frac{b}{3a}\right)+d\\=ay^3+py+q$$
Where $p$ and $q$ are constants.
Make the substitution $y=uz$ and multiply both sides by $v$.
$$0=v(uz)^3-3v(uz)+3v\\=vu^3z^3-3vuz+v$$
And have $vu^3=4$ and $-3vu=-3$ so that it comes in our $(0)$ form. Solving this system of equations gives $u=2$ and $v=\frac12$.
$$0=4z^3-3z+\frac12$$
Solving for $z$ and recalling the period of cosine, we get
$$z=\cos\left(\frac{2\pi k}3+\frac13\arccos\left(\frac{-1}2\right)\right)$$
$$y=2z=2\cos\left(\frac{2\pi k}3+\frac13\arccos\left(\frac{-1}2\right)\right)$$
$$x=1+y=1+2\cos\left(\frac{2\pi k}3+\frac13\arccos\left(\frac{-1}2\right)\right)$$
This is not always possible, but under those cases, we can just switch to different trig functions and use their corresponding formulas. It's pretty easy to see which trig function you should use at the "solve for $u$ and $v$ step".
An advantage to this method is that it avoids casus irreducibilis, which is a fancy way of saying you can't factor a cubic polynomial using only real numbers, radicals, and basic arithmetic operations. This does not include trig functions, however, which can be used to obtain nice forms of the solution.
To compare, the radical solution for the above polynomial is given as
$$x_1=1-\frac12(1-i\sqrt3)\sqrt[3]{\frac12(-1+i\sqrt3)}-\frac{1+i\sqrt3}{\sqrt[3]{4(-1+i\sqrt3}}$$
$$x_2=1-\frac{1-i\sqrt3}{\sqrt[3]{4(-1+i\sqrt3)}}-\frac12(1+i\sqrt3)\sqrt[3]{\frac12(-1+i\sqrt3)}$$
$$x_3=1+\frac1{\sqrt[3]{\frac12(-1+i\sqrt3)}}+\sqrt[3]{\frac12(-1+i\sqrt3)}$$
Compared side by side, you might find one form much nicer.