Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]
share|improve this question
    
This does not appear to be a Mathematica question. Since Mathematica is giving you correct results, the problem lies elsewhere. Perhaps MATH.stackexchange.com? –  rasher Mar 22 at 3:01

1 Answer 1

up vote 3 down vote accepted

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}
share|improve this answer
    
ubpdqn Thank you very very much, it was my mistake in understanding any way. –  Pipe Mar 22 at 11:29
    
@Pipe we are all learning here. I am glad we could resolve.:) –  ubpdqn Mar 22 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 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. –  ubpdqn Mar 22 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 at 15:04

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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