0
$\begingroup$

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]
$\endgroup$
1
  • $\begingroup$ This does not appear to be a Mathematica question. Since Mathematica is giving you correct results, the problem lies elsewhere. Perhaps MATH.stackexchange.com? $\endgroup$
    – ciao
    Mar 22, 2014 at 3:01

1 Answer 1

4
$\begingroup$

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}
$\endgroup$
7
  • $\begingroup$ ubpdqn Thank you very very much, it was my mistake in understanding any way. $\endgroup$
    – Pipe
    Mar 22, 2014 at 11:29
  • 1
    $\begingroup$ @Pipe we are all learning here. I am glad we could resolve.:) $\endgroup$
    – ubpdqn
    Mar 22, 2014 at 11:30
  • $\begingroup$ @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? $\endgroup$
    – Pipe
    Mar 22, 2014 at 13:31
  • $\begingroup$ @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. $\endgroup$
    – ubpdqn
    Mar 22, 2014 at 13:50
  • $\begingroup$ @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? $\endgroup$
    – Pipe
    Mar 22, 2014 at 15:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.