Third degree equation resolution

Question

On page 17 of the Handbook of Mathematical Functions With Formulas, Graphs and Mathematical Tables by Abramowitz and Stegun, it is written that given the third degree equation:

$$z^3 + a_2\,z^2 + a_1\,z + a_0 = 0$$

and calculating the following values first:

$$q = \frac{1}{3}\,a_1 - \frac{1}{9}\,a_2^2\,, \quad r = \frac{1}{6}\left(a_1\,a_2 - 3\,a_0\right) - \frac{1}{27}\,a_2^3$$

and then these other values:

$$s_1 = \left(r + \left(q^3 + r^2\right)^{\frac{1}{2}}\right)^{\frac{1}{3}} \,, \quad s_2 = \left(r - \left(q^3 + r^2\right)^{\frac{1}{2}}\right)^{\frac{1}{3}}\,,$$

the three desired roots of the cubic equation are:

$$\begin{align*} & z_1 = (s_1 + s_2) - \frac{a_2}{3}\,, \\ & z_2 = - \frac{1}{2}\,(s_1 + s_2) - \frac{a_2}{3} + \frac{\text{i}\,\sqrt{3}}{2}\,(s_1 - s_2)\,, \\ & z_3 = - \frac{1}{2}\,(s_1 + s_2) - \frac{a_2}{3} - \frac{\text{i}\,\sqrt{3}}{2}\,(s_1 - s_2)\,.\end{align*}$$

Unfortunately, writing the following code:

a2 = 0;
a1 = 1;
a0 = 1;

NSolve[z^3 + a2 z^2 + a1 z + a0 == 0, z]

q = 1/3 a1 - 1/9 a2^2;
r = 1/6 (a1 a2 - 3 a0) - 1/27 a2^3;

s1 = (r + (q^3 + r^2)^(1/2))^(1/3);
s2 = (r - (q^3 + r^2)^(1/2))^(1/3);

z1 = (s1 + s2) - a2/3 // N
z2 = -1/2 (s1 + s2) - a2/3 + I Sqrt[3]/2 (s1 - s2) // N
z3 = -1/2 (s1 + s2) - a2/3 - I Sqrt[3]/2 (s1 - s2) // N

I get:

{{z -> -0.682328}, {z -> 0.341164 - 1.16154 I}, {z -> 0.341164 + 1.16154 I}}

0.835342 + 0.876227 I

0.341164 - 0.590913 I

-1.17651 - 0.285314 I

so there is definitely something that does not work, but I do not know what! Thank you!

Answer 1

As already noted by other answers, the version of Cardano in Abramowitz and Stegun uses a branch cut for the cube root that is not the same as the branch cut used by Mathematica.

So, here's how to implement Cardano manually for this case, as adapted from Numerical Recipes:

myCubic[ac_, bc_, cc_, dc_] := 
  Module[{a, b, c, da, db, ds, dt, du, q, r},
         {a, b, c} = {bc, cc, dc}/ac;
         q = (a^2 - 3 b)/9; r = (2 a^3 + 9 (3 c - a b))/54;
         ds = Sqrt[r^2 - q^3];
         If[Re[Conjugate[r] ds] < 0, ds = -ds];
         da = -(r + ds)^(1/3); db = If[da == 0, 0, q/da];
         dt = da + db; du = I Sqrt[3] (da - db)/2;
         {-dt/2 - du, dt, -dt/2 + du} - a/3]

Test:

With[{a = 1, b = 2, c = 3, d = 4}, 
     FromDigits[{a, b, c, d}, #] & /@ N[myCubic[a, b, c, d]] // Chop]
   {0, 0, 0}

With[{a = I - 3, b = 5 + 4 I, c = 2 - 2 I, d = 6}, 
     FromDigits[{a, b, c, d}, #] & /@ N[myCubic[a, b, c, d]] // Chop]
   {0, 0, 0}

Answer 2

The cube root (1/3 power) in A&S is CubeRoot in Mathematica:

a2 = 0;
a1 = 1;
a0 = 1;

NSolve[z^3 + a2 z^2 + a1 z + a0 == 0, z]

q = 1/3 a1 - 1/9 a2^2;
r = 1/6 (a1 a2 - 3 a0) - 1/27 a2^3;

s1 = (r + (q^3 + r^2)^(1/2)) // CubeRoot;
s2 = (r - (q^3 + r^2)^(1/2)) // CubeRoot;

z1 = (s1 + s2) - a2/3 // N
z2 = -1/2 (s1 + s2) - a2/3 + I Sqrt[3]/2 (s1 - s2) // N
z3 = -1/2 (s1 + s2) - a2/3 - I Sqrt[3]/2 (s1 - s2) // N
(*
  {{z -> -0.682328}, {z -> 0.341164 - 1.16154 I}, {z -> 0.341164 + 1.16154 I}}

  -0.682328
  0.341164 + 1.16154 I
  0.341164 - 1.16154 I
*)

Answer 3

Because of the cube roots, there are three, in general complex, solutions for s1 and s2. You have a positive discriminant, $q^3+r^2=31/108$, corresponding to one real root and two complex conjugate roots.

Test all possible solutions for s1 and s2 with something like the following:

Block[{a2 = 0, a1 = 1, a0 = 1, q, r, s1, s2, x},
      q = 1/3 a1 - 1/9 a2^2;
      r = 1/6 (a1 a2 - 3 a0) - 1/27 a2^3;
      Print["Discriminant: ", q^3 + r^2, " > 0"];
      s1 = Chop[x /. {ToRules[Reduce[x^3 == r + (q^3 + r^2)^(1/2), x]]}];
      s2 = Chop[x /. {ToRules[Reduce[x^3 == r - (q^3 + r^2)^(1/2), x]]}];
      Flatten[
         Table[{
            Chop[N[(s1[[i]] + s2[[j]]) - a2/3]],
            Chop[N[-1/2 (s1[[i]] + s2[[j]]) - a2/3 + I Sqrt[3]/2 (s1[[i]]-s2[[j]])]], 
            Chop[N[-1/2 (s1[[i]] + s2[[j]]) - a2/3 - I Sqrt[3]/2 (s1[[i]]-s2[[j]])]]},
            {i, 1, 3}, {j, 1, 3}], 1]
    ] // Column

Pick the solution with one real root and two complex conjugate roots. If the discriminant is less than zero, then there are three real roots to pick. If the discriminant equals zero, then there are three real roots with at least two of them equal.