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$
  • $\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 '14 at 3:01
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}
| improve this answer | |
$\endgroup$
  • $\begingroup$ ubpdqn Thank you very very much, it was my mistake in understanding any way. $\endgroup$ – Pipe Mar 22 '14 at 11:29
  • 1
    $\begingroup$ @Pipe we are all learning here. I am glad we could resolve.:) $\endgroup$ – ubpdqn Mar 22 '14 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 '14 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 '14 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 '14 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.