# Solutions of cubic equation D>0 or D<0

I have problem to obtain three roots of cubic equation

 N[Solve[a x^3 + b*x^2 + c*x + d == 0, x], 15]


using analytical procedure http://en.wikipedia.org/wiki/Casus_irreducibilis

Discriminant of the system is positive and solutions which I obtained with the formula x1, x2 and x3 and Mathematica Solve are not the same?

 a = 1;
b = 0.00100025;
c = 0.100000250125;
d = 0.0000500125;

p = (3 a c - b^2)/(3 a^2);
q = (2 b^3 - 9 a b c + 27 a^2 d)/(27 a^3);

N[q^2/4 + p^3/27]

t1 = \[Psi]1 Power[-(q/2) + Sqrt[q^2/4 + p^3/27], (
3)^-1] + \[Psi]1^2 Power[-(q/2) - Sqrt[q^2/4 + p^3/27], (3)^-1];
t2 = \[Psi]2 Power[-(q/2) + Sqrt[q^2/4 + p^3/27], (
3)^-1] + \[Psi]2^2 Power[-(q/2) - Sqrt[q^2/4 + p^3/27], (3)^-1];
t3 = \[Psi]3 Power[-(q/2) + Sqrt[q^2/4 + p^3/27], (
3)^-1] + \[Psi]3^2 Power[-(q/2) - Sqrt[q^2/4 + p^3/27], (3)^-1];

\[Psi]1 = 1;
\[Psi]2 = -(1/2) + Sqrt[3]/2 I;
\[Psi]3 = -(1/2) - Sqrt[3]/2 I;

x1 = t1 - b/(3 a);
x2 = t2 - b/(3 a);
x3 = t3 - b/(3 a);

N[x1]
N[x2]
N[x3]
N[Solve[a x^3 + b*x^2 + c*x + d == 0, x], 15]


And which formula to use in the case if the discriminant is negative?

For example if the equations is in the following form

 NSolve[0.5 + 25.1 x + 10. x^2 + x^3==0,x]

• This does not appear to be a Mathematica question. Since Mathematica is giving you correct results, the problem lies elsewhere. Perhaps MATH.stackexchange.com?
– ciao
Mar 22, 2014 at 3:01

The problem is your implementation of the cubic formula. In particular the use of the third roots of unit. Let $r=\sqrt{q^2+4p^3/27}$, then solution of $z^3=r$ can be expressed as $r^{1/3},r^{1/3}e^{2\pi i/3},r^{1/3}e^{4\pi i/3}$ or using the nomenclature in the post $r^{1/3}\psi_j$ where $j\in\{1,2,3\}$. This is not how the code is implemented. Hence the difference in result.

Further, note the discriminant of $a x^3+b x^2+ c x+ d$<0, hence one real solution and two complex soultions are expected. The discriminant of 0.5 + 25.1 x + 10. x^2 + x^3 is greater than zero (-> three real solutions). The formula

Here is an implementation of the cubic formula (there are almost certainly better ones but I post it to illustrate):

cubsoln[pol_, x_] := Module[{coeff, a0, a1, a2, a3, p, q, w3, w, ps},
coeff = CoefficientList[pol, x];
{a0, a1, a2, a3} = (coeff/Last@coeff);
p = (3 a1 - a2^2)/3;
q = (9 a1 a2  - 27 a0 - 2 a2^3)/27;
w3 = (q/2 + {1/2, -1/2} Sqrt[q^2 + 4 p^3/27]);
w = (Flatten@(Solve[z^3 == #, z] & /@ w3))[[All, 2]];
ps = Union[# - p/(3 #) & /@ w,
SameTest -> (Abs[#1 - #2] < 10^-15 &)];
Simplify[Chop[ps - a2/3]]
]


Now applying:

f = 0.0000500125 + 0.100000250125 x + 0.00100025 x^2 + x^3;
cubsoln[f, x]


yields:

{-0.000500125, -0.000250062 - 0.316228 I, -0.000250062 + 0.316228 I}


compared with:

Solve[f == 0, x]


yields:

{{x -> -0.000500125}, {x -> -0.000250062 -
0.316228 I}, {x -> -0.000250062 + 0.316228 I}}


and for 0.5 + 25.1 x + 10. x^2 + x^3:

cubsoln[0.5 + 25.1 x + 10. x^2 + x^3, x]


yields:

{-5., -4.97992, -0.0200806}


as expected from the discriminant.

The implementation fails when p=0, as this is $x^3=q$ or $(z-a_2/3)=q$.

Here is an implementation using roots of unity:

cubru[pol_, x_] := Module[{coeff, a0, a1, a2, a3, p, q, w3, w, ps},
coeff = CoefficientList[pol, x];
{a0, a1, a2, a3} = (coeff/Last@coeff);
p = (3 a1 - a2^2)/3;
q = (9 a1 a2  - 27 a0 - 2 a2^3)/27;
w3 = (q/2 + {1/2, -1/2} Sqrt[q^2 + 4 p^3/27]);
w = Flatten[#^(1/3) Table[Exp[ I j], {j, 0, 4 Pi/3, 2 Pi/3}] & /@
w3];
ps = Union[# - p/(3 #) & /@ w,
SameTest -> (Abs[#1 - #2] < 10^-15 &)];
Simplify[Chop[ps - a2/3]]
]


Testing:

cubru[f, x]


yields:

{-0.000500125, -0.000250062 - 0.316228 I, -0.000250062 + 0.316228 I}

• ubpdqn Thank you very very much, it was my mistake in understanding any way.
– Pipe
Mar 22, 2014 at 11:29
• @Pipe we are all learning here. I am glad we could resolve.:) Mar 22, 2014 at 11:30
• @ubpdqna I have one thing which is not clear, why the last step Simplify[Chop[ps - a2/3]] is not giving solutions the same with NSolve[0 == 0.0000500125 + 0.100000250125 x + 0.00100025 x^2 + x^3,x] if I set parameters immediately on the beginning?
– Pipe
Mar 22, 2014 at 13:31
• @Pipe I will look at when I can...it is midnight in my tz...could you explain exactly what you did?, eg in edit to Q. Mar 22, 2014 at 13:50
• @ubpdqna code is great and everything is working. I just followed the code and in the last step I thought that if I substitute parameters I obtain in the last step {-1.03412, -1.03412, 0.467062 - 0.865364 I, 0.467062 + 0.865364 I, 0.467062 - 0.865364 I, 0.467062 + 0.865364 I} and solutions which code gave is {{x -> -0.000500125}, {x -> -0.000250062 - 0.316228 I}, {x -> -0.000250062 + 0.316228 I}}. My question is how to see roots in analytic form which you used in the code?
– Pipe
Mar 22, 2014 at 15:04