Message Boards Message Boards

2
|
7360 Views
|
26 Replies
|
18 Total Likes
View groups...
Share
Share this post:

Is it possible to use explicit quartics solution for solving biquadratics?

Posted 10 years ago

I have just discovered that explicit solution for quartics in form Mathematica gives it, is not suitable for solving biquadratic equation (i.e. when B=D=0). Say we have the following equation:

Roots[A x^4 + B x^3 + C x^2 + D x + E == 0, x]

After solving it explicitly and making certain denotations we get all four roots:

X1= -(B/(4 A)) - K4/2 - Sqrt[K5 - K6]/2
X2= -(B/(4 A)) - K4/2 + Sqrt[K5 - K6]/2
X3= -(B/(4 A)) + K4/2 - Sqrt[K5 + K6]/2
X4= -(B/(4 A)) + K4/2 + Sqrt[K5 + K6]/2

Where

K1 = 2 C^3 - 9 B C D + 27 A D^2 + 27 B^2 E - 72 A C E,
K2 = K1 + Sqrt[- 4 (C^2 - 3 B D + 12 A E)^3 + K1^2],
K3 = (2^(1/3) (C^2 - 3 B D + 12 A E))/(3 A K2^(1/3)) + K2^(1/3)/(3 A 2^(1/3)),
K4 = Sqrt[B^2/(4 A^2) - (2 C)/(3 A) + K3],
K5 = B^2/(2 A^2) - (4 C)/(3 A) - K3,
K6 = (-(B^3/A^3) + (4 B C)/A^2 - (8 D)/A)/(4 K4).

Then if we put B=D=0, and try to solve it with Roots and with explicit formula, answers will appear to be different. It happens only for some A,C,E combinations though, for some other solutions agree (it is probably based on discriminant sign and some other coefficients as mentioned in wikipedia). The problem comes from K2 coefficient when expression under square root becomes negative hence making all the further computations complex. Does anybody know how to deal with that? Is there any other stable explicit solution for the quartics? I have to calculate explicitly over a range of A,B,C,D,E coefficients including case when B=D=0. I am attaching a piece of Mathematica code where this problem is put in.

Attachments:
POSTED BY: Yuriy Ivanov
26 Replies

I think the conclusion is that the most general formula for the solution of an equation is generated by Reduce, since its result handles special cases. However, Reduce can only handle small problems.

POSTED BY: Frank Kampas

Hi, all thanks a lot for your brilliant contributions. But what are the conclusions for the the ' not so deeply involved' users? Do we need further instructions ( more than the documented ones) about using Solve, Roots, Reduce, Setprecision, ... or will there be modifications be done inside of mathematica?

Thanks again, Peter

POSTED BY: Peter Bischet
Posted 10 years ago

Peter and all,

Despite the fact this discussion raised questions much deeper than I expected, what I gained from it with regards to my current problem, is: In order to find the solution of quartic in most general form it should be preserved throughout a code as a Root object:

root_n = Root[#1^4 A1 + #1^3 A2 + #1^2 A3 + #1 A4 + A5 &, n]

Explicit solution given by Solve is valid only for certain combination of A1...A5. Reduce outputs the solution in form of a root object. However if equation is not expressed in polynomial form, Reduce will have troubles solving it and Solve is a way to go.

It helped me partially, now I am struggling with consequence of using Root object - ordering of the roots. Separate discussion is created.

Thanks all involved in discussion, appreciate your responses!

Regards, Yuriy

POSTED BY: Yuriy Ivanov

In fact, if I don't use Sort but look at a slightly different distance, e.g.

sols = Table[{Atemp, Btemp, Ctemp, Dtemp, Etemp} = results[[i]];
   {Sort[Table[
      Chop[N[Roots[
          Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
            Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, 
          x]]][[i, 2]], {i, 1, 
       Length[N[
         Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 
            0 /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
            Ee -> Temp}, x]]]}]], 
    Sort[Table[
      Chop[N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
           x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Etemp}]][[i, 2]], {i, 1, 
       Length[N[
         Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
           x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Etemp}]]}]]}, {i, 1, Length[results]}];

and then

Total @ Chop[(Re[#][[1]]^2 - Re[#][[2]]^2) + (Im[#][[1]]^2 - Im[#][[2]]^2)] & /@ sols

I get all zeros, which say that Roots gets it right nearly all the time.

Cheers, Marco

POSTED BY: Marco Thiel

Perhaps the "errors" you see are due to imprecision in machine precision calculations.

POSTED BY: Frank Kampas

Somehow, the next failure is

sol1: {-0.600937,0.0132837 -1.45082 I,0.0132837 +1.45082 I,1.04832}
sol2: {-0.600937,0.0132837 -1.45082 I,0.0132837 +1.45082 I,1.04832}

lets get this in the long form (accents ommited in courtesy of the WC (* Wolfram Community *) Editor)

In[184]:= {-0.6009373760877655, 0.013283715568729561 - 1.4508231140635948 I, 
            0.013283715568729561 + 1.4508231140635948 I, 1.0483249164991404} == 
          {-0.6009373760877655, 0.01328371556872926 - 1.450823114063595 I, 
            0.01328371556872926 + 1.450823114063595 I, 1.0483249164991404}
Out[184]= False

it's false in a strict sense, but: has been computed with precision 16 and is wrong beginning with the 15-th digit only.

So it seems that it could be made correct to any meaningful precision and I tend to believe Daniel and Frank.

POSTED BY: Udo Krause

This is the error message:

Reduce::ratnz: Reduce was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

So it did need (I think) rationalise the system.

Thanks,

Marco

POSTED BY: Marco Thiel

I think you need to do the calculation with infinite precision inputs and not numericize (sic) the results.

POSTED BY: Frank Kampas
In[196]:= MCNumber = 10000;
Tally[Monitor[
  Table[Atemp = Rationalize[RandomReal[] - 0.5]; 
   Btemp = Rationalize[RandomReal[] - 0.5];
   Ctemp = Rationalize[RandomReal[] - 0.5]; 
   Dtemp = Rationalize[RandomReal[] - 0.5];
   Etemp = Rationalize[RandomReal[] - 0.5];
   Sort[Table[SetPrecision[Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
         Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, x], 14][[i, 2]], {i, 1, 4}]] == 
    Sort[Table[SetPrecision[Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, x] /. {Aa -> Atemp, 
         Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, 14][[i, 2]], {i, 1, 4}]], {k, 1, MCNumber}], 
  ProgressIndicator[k, {0, MCNumber}]]]

Out[197]= {{True, 10000}}

No warnings, no error messages, no failures, just numerics.

In[198]:= $MachinePrecision
Out[198]= 15.9546

it makes sense to set the precision to 14 if machine precision is around 16.

POSTED BY: Udo Krause

Hi,

I agree that there is a numerical part involved, but I suppose that one has to be very careful with this.

If you run

dist = MCNumber = 1000; Table[
 Atemp = RandomReal[] - 0.5; Btemp = RandomReal[] - 0.5; 
 Ctemp = RandomReal[] - 0.5; Dtemp = RandomReal[] - 0.5; 
 Etemp = RandomReal[] - 0.5; 
 Norm[Sort[
    Table[N[Roots[
        Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
          Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, x]][[i,
        2]], {i, 1, 
      Length[N[
        Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 
           0 /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Temp}, x]]]}]] - 
   Sort[Table[
     N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
         x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
         Ee -> Etemp}][[i, 2]], {i, 1, 
      Length[N[
        Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
          x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
          Ee -> Etemp}]]}]]], {k, 1, MCNumber}];

you get a list of distanced between the two solutions. In fact, most of them are, as expected, quite close to zero, but there are some large deviations as well:

Sort[dist][[-150 ;;]]

{0.447706, 0.468712, 0.478374, 0.51479, 0.536065, 0.537396, 0.543305, \
0.660131, 0.689731, 0.807025, 0.809497, 0.873262, 0.915736, 0.942873, \
1.0216, 1.06628, 1.07255, 1.07922, 1.07954, 1.08418, 1.13938, \
1.14699, 1.16781, 1.17355, 1.18269, 1.1934, 1.19466, 1.21008, \
1.25215, 1.27071, 1.28718, 1.29177, 1.3206, 1.32214, 1.33523, \
1.35991, 1.36773, 1.37844, 1.40719, 1.40724, 1.40763, 1.41086, \
1.43494, 1.4613, 1.46464, 1.47223, 1.47275, 1.47728, 1.49559, \
1.52547, 1.53006, 1.53025, 1.53896, 1.55807, 1.58614, 1.59941, \
1.63932, 1.64731, 1.65441, 1.65704, 1.6784, 1.68514, 1.71468, \
1.73314, 1.74677, 1.76276, 1.77466, 1.84156, 1.85621, 1.87894, \
1.88623, 1.90572, 1.92754, 1.94024, 1.94376, 1.96933, 1.98945, \
2.03146, 2.05536, 2.09187, 2.09309, 2.14448, 2.18454, 2.18603, \
2.19656, 2.20124, 2.22276, 2.24433, 2.24465, 2.2963, 2.37969, \
2.41349, 2.42688, 2.47719, 2.4975, 2.51811, 2.53008, 2.5723, 2.59401, \
2.63813, 2.644, 2.65012, 2.65131, 2.69886, 2.7626, 2.78525, 2.80045, \
2.81197, 2.87692, 2.90851, 2.93173, 2.94927, 3.04928, 3.073, 3.08039, \
3.11489, 3.15659, 3.16508, 3.20288, 3.21572, 3.21919, 3.23016, \
3.24037, 3.24488, 3.30734, 3.3232, 3.32944, 3.33474, 3.40173, \
3.41504, 3.41817, 3.48765, 3.52086, 3.54337, 3.58852, 3.599, 3.61942, \
3.77933, 3.97308, 4.13652, 4.1701, 4.28506, 4.57054, 4.66001, \
4.68233, 4.72912, 6.57014, 6.70144, 7.60674, 15.8699}

Of course, these might be numerical problems, but it is important to note that in many cases the solutions are off by quite a bit. If we look at the rather complicated formula for one of the solutions of the quartic

  Solve[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, x][[3]]

    {x->1/2 \[Sqrt](Bb^2/(4 Aa^2)+(1/(3 Power[2, (3)^-1] Aa))((Sqrt[(-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^2-4 (12 Aa Ee-3 Bb Dd+Cc^2)^3]-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^(1/3))+(Power[2, (3)^-1] (12 Aa Ee-3 Bb Dd+Cc^2))/(3 Aa (Sqrt[(-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^2-4 (12 Aa Ee-3 Bb Dd+Cc^2)^3]-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^(1/3))-(2 Cc)/(3 Aa))-1/2 \[Sqrt](Bb^2/(2 Aa^2)+(-(Bb^3/Aa^3)+(4 Bb Cc)/Aa^2-(8 Dd)/Aa)/(4 \[Sqrt](Bb^2/(4 Aa^2)+(1/(3 Power[2, (3)^-1] Aa))((\[Sqrt]((-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^2-4 (12 Aa Ee-3 Bb Dd+Cc^2)^3)-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^(1/3))+(Power[2, (3)^-1] (12 Aa Ee-3 Bb Dd+Cc^2))/(3 Aa (\[Sqrt]((-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^2-4 (12 Aa Ee-3 Bb Dd+Cc^2)^3)-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^(1/3))-(2 Cc)/(3 Aa)))-(1/(3 Power[2, (3)^-1] Aa))((Sqrt[(-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^2-4 (12 Aa Ee-3 Bb Dd+Cc^2)^3]-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^(1/3))-(Power[2, (3)^-1] (12 Aa Ee-3 Bb Dd+Cc^2))/(3 Aa (Sqrt[(-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^2-4 (12 Aa Ee-3 Bb Dd+Cc^2)^3]-72 Aa Cc Ee+27 Aa Dd^2+27 Bb^2 Ee-9 Bb Cc Dd+2 Cc^3)^(1/3))-(4 Cc)/(3 Aa))-Bb/(4 Aa)}

it becomes quite clear that sometimes we divide by terms that might be very close to zero, which might give large errors. But I think there is something else - which is my fault. Playing with this I found ar problem in my approach. The sort command does not do what I wanted it to do...

MCNumber = 1000; results = {}; Do[
 Atemp = RandomReal[] - 0.5; Btemp = RandomReal[] - 0.5; 
 Ctemp = RandomReal[] - 0.5; Dtemp = RandomReal[] - 0.5; 
 Etemp = RandomReal[] - 0.5; 
 normpol = 
  Norm[Sort[
     Table[N[Roots[
         Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
           Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, 
         x]][[i, 2]], {i, 1, 
       Length[N[
         Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 
            0 /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
            Ee -> Temp}, x]]]}]] - 
    Sort[Table[
      N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
          x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
          Ee -> Etemp}][[i, 2]], {i, 1, 
       Length[N[
         Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
           x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Etemp}]]}]]]; 
 If[normpol > 1, 
  AppendTo[results, {Atemp, Btemp, Ctemp, Dtemp, Etemp}]], {k, 1, 
  MCNumber}];

Then you execute:

{Atemp, Btemp, Ctemp, Dtemp, Etemp} = results[[3]];
{Sort[Table[
   N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
        Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, x]][[i, 
     2]], {i, 1, 
    Length[N[
      Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp,
          Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Temp}, x]]]}]],
  Sort[Table[
   N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
       x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
       Ee -> Etemp}][[i, 2]], {i, 1, 
    Length[N[
      Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
        x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
        Ee -> Etemp}]]}]]}

That gives:

{{-0.75599 - 0.773082 I, -0.75599 + 0.773082 I, 0.208992, 0.882124}, 
{-0.75599 + 0.773082 I, -0.75599 - 0.773082 I, 0.208992 - 5.55112*10^-17 I, 0.882124 + 0. I}}

so the first and the second solutions are exchanged. That leads to a large deviation in the norm. I suppose that this is why there were so many large deviations.

Cheers,

Marco

POSTED BY: Marco Thiel

In this last statement Ctemp appears twice, Etemp is missing:

MCNumber = 10000; Tally[Monitor[Table[
   Atemp = RandomReal[] - 0.5; Btemp = RandomReal[] - 0.5; 
   Ctemp = RandomReal[] - 0.5; Ctemp = RandomReal[] - 0.5; 
   Dtemp = RandomReal[] - 0.5; 
   Sort[Table[
      N[Reduce[
         Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
           Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, 
         x]][[i, 2]], {i, 1, 4}]] == 
    Sort[Table[
      N[Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
          x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
          Ee -> Etemp}][[i, 2]], {i, 1, 4}]], {k, 1, MCNumber}], 
  ProgressIndicator[k, {0, MCNumber}]]]

then Etemp is used.

POSTED BY: Udo Krause

Ups, sorry. Here's the fix; the problem remains.

MCNumber = 1000; Tally[
 Monitor[Table[Atemp = RandomReal[] - 0.5; 
   Btemp = RandomReal[] - 0.5;
   Ctemp = RandomReal[] - 0.5; Dtemp = RandomReal[] - 0.5;
   Etemp = RandomReal[] - 0.5;
   Sort[Table[
      N[Reduce[
         Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
           Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, 
         x]][[i, 2]], {i, 1, 4}]] == 
    Sort[Table[
      N[Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
          x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
          Ee -> Etemp}][[i, 2]], {i, 1, 4}]], {k, 1, MCNumber}], 
  ProgressIndicator[k, {0, MCNumber}]]]

results in:

{{True, 996}, {False, 4}}

I'll fix that in the code above not to confuse readers. Thanks a lot for your corrections.

It does not change the general result though.

Cheers,

Marco

POSTED BY: Marco Thiel

I did the corrected MCNumber thing Rationalized and got a

In[166]:= MCNumber = 10000;
Tally[Monitor[Table[Atemp = Rationalize[RandomReal[] - 0.5]; 
   Btemp = Rationalize[RandomReal[] - 0.5];
   Ctemp = Rationalize[RandomReal[] - 0.5]; 
   Dtemp = Rationalize[RandomReal[] - 0.5];
   Etemp = Rationalize[RandomReal[] - 0.5];
   Sort[Table[N[Reduce[
         Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
           Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, 
         x]][[i, 2]], {i, 1, 4}]] == 
    Sort[Table[N[Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
          x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
          Ee -> Etemp}][[i, 2]], {i, 1, 4}]], {k, 1, MCNumber}], 
  ProgressIndicator[k, {0, MCNumber}]]]

Out[167]= {{True, 9953}, {False, 47}}

OMG! So it's interesting to run into the first failure

goAgain = True;
While[goAgain, 
 Atemp = Rationalize[RandomReal[] - 0.5];
 Btemp = Rationalize[RandomReal[] - 0.5];
 Ctemp = Rationalize[RandomReal[] - 0.5];
 Dtemp = Rationalize[RandomReal[] - 0.5];
 Etemp = Rationalize[RandomReal[] - 0.5];
 sol1 = Sort[Table[N[Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
         Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, x]][[i, 2]], {i, 1, 4}]];
 sol2 = Sort[Table[N[Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, x] /. {Aa -> Atemp, 
         Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp,  Ee -> Etemp}][[i, 2]], {i, 1, 4}]];
 goAgain = (sol1 == sol2);
 If[!goAgain, 
  Print["sol1: ", sol1];
  Print["sol2: ", sol2]
  ]
 ]

sol1: {-0.866444-0.043348 I,-0.866444+0.043348 I,0.595995 -1.03621 I,0.595995 +1.03621 I}
sol2: {-0.866444-0.043348 I,-0.866444+0.043348 I,0.595995 -1.03621 I,0.595995 +1.03621 I}

with standard precision you don't spot it. Must be behind that.

POSTED BY: Udo Krause

I have just checked Frank's observation. It is quite true that Reduce seems to do much better (I ignored some error messages), but it is not always right either.

MCNumber = 10000; Tally[Monitor[Table[
   Atemp = RandomReal[] - 0.5; Btemp = RandomReal[] - 0.5; 
   Ctemp = RandomReal[] - 0.5; Dtemp = RandomReal[] - 0.5; 
   Etemp = RandomReal[] - 0.5; 
   Sort[Table[
      N[Reduce[
         Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
           Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, 
         x]][[i, 2]], {i, 1, 
       Length[N[
         Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 
            0 /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
            Ee -> Temp}, x]]]}]] == 
    Sort[Table[
      N[Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
          x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
          Ee -> Etemp}][[i, 2]], {i, 1, 
       Length[N[
         Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
           x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Etemp}]]}]], {k, 1, MCNumber}], 
  ProgressIndicator[k, {0, MCNumber}]]]

gave

 {{True, 9981}, {False, 19}}

Another thing is that Solve is performing about the same as Roots:

MCNumber = 1000; Tally[Monitor[Table[
   Atemp = RandomReal[] - 0.5; Btemp = RandomReal[] - 0.5; 
   Ctemp = RandomReal[] - 0.5; Dtemp = RandomReal[] - 0.5; 
   Etemp = RandomReal[] - 0.5; 
   Sort[Table[
      x /. N[Solve[
          Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
            Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, 
          x]][[i]], {i, 1, 
       Length[N[
         Solve[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 
            0 /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
            Ee -> Temp}, x]]]}]
     ] == 
    Sort[Table[
      x /. N[Solve[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
           x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Etemp}][[i]], {i, 1, 
       Length[N[
         Solve[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
           x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Etemp}]]}]

     ], {k, 1, MCNumber}], ProgressIndicator[k, {0, MCNumber}]]]

gives

{{True, 621}, {False, 379}}

So 38% wrong results. Also, even if Reduce is better, it is curious that Solve and Roots have problems in more than one third of the cases without giving any warning.

Cheers,

Marco

PS: In fact I think that the program is very inefficient. I solve the equation once more for the upper limit of the Table command which is not really required. A 4 there does the trick.

MCNumber = 10000; Tally[Monitor[Table[
   Atemp = RandomReal[] - 0.5; Btemp = RandomReal[] - 0.5; 
   Ctemp = RandomReal[] - 0.5; Dtemp = RandomReal[] - 0.5; 
   Etemp = RandomReal[] - 0.5; 
   Sort[Table[
      N[Reduce[
         Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> Atemp, 
           Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, Ee -> Etemp}, 
         x]][[i, 2]], {i, 1, 4}]] == 
    Sort[Table[
      N[Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
          x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
          Ee -> Etemp}][[i, 2]], {i, 1, 4}]], {k, 1, MCNumber}], 
  ProgressIndicator[k, {0, MCNumber}]]]
POSTED BY: Marco Thiel

Seems to me that Reduce is a better way to attack this problem. As I understand it, Reduce handles all the cases, including measure zero cases.

r = Reduce[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, x];

r /. {Aa -> -30, Bb -> 0, Cc -> 13, Dd -> 0, Ee -> -1}

x == -(1/Sqrt[3]) || x == -(1/Sqrt[10]) || x == 1/Sqrt[10] || x == 1/Sqrt[3]

POSTED BY: Frank Kampas

The soluions are generically correct. For a measure zero set of parameter values they may give bad resuts on substitution. One can often recover correct values using Limit. In this example I did so by substituting for all but the cubic coefficient, then taking Limit[...] as that remaining coefficient approaches zero.

POSTED BY: Daniel Lichtblau

Dear Daniel,

are you sure about that measure zero statement? I just ran the following little program:

MCNumber = 10000; Tally[Monitor[Table[
   Atemp = RandomReal[]; Btemp = RandomReal[]; Ctemp = RandomReal[]; 
   Dtemp = RandomReal[]; Etemp = RandomReal[]; 
   Sort[Table[
      N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 
           0 /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Etemp}, x]][[i, 2]], {i, 1, 
       Length[N[
         Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 
            0 /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
            Ee -> Temp}, x]]]}]] == 
    Sort[Table[
      N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
          x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
          Ee -> Etemp}][[i, 2]], {i, 1, 
       Length[N[
         Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
           x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Etemp}]]}]], {k, 1, MCNumber}], 
  ProgressIndicator[k, {0, MCNumber}]]]

It generates different parameters and compares the solutions. The results are usually that about 22-23% of times we get it wrong.

{{True, 7734}, {False, 2266}}

Or am I just "lucky" at picking out values of the measure zero region?

Cheers,

Marco

PS: If the parameters are also negative, i.e. the uniform distribution is centred at zero it gets even worse:

MCNumber = 1000; Tally[Monitor[Table[
   Atemp = RandomReal[] - 0.5; Btemp = RandomReal[] - 0.5; 
   Ctemp = RandomReal[] - 0.5; Dtemp = RandomReal[] - 0.5; 
   Etemp = RandomReal[] - 0.5; 
   Sort[Table[
      N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 
           0 /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Etemp}, x]][[i, 2]], {i, 1, 
       Length[N[
         Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 
            0 /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
            Ee -> Temp}, x]]]}]] == 
    Sort[Table[
      N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
          x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
          Ee -> Etemp}][[i, 2]], {i, 1, 
       Length[N[
         Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, 
           x] /. {Aa -> Atemp, Bb -> Btemp, Cc -> Ctemp, Dd -> Dtemp, 
           Ee -> Etemp}]]}]], {k, 1, MCNumber}], 
  ProgressIndicator[k, {0, MCNumber}]]]

gives about 37% wrong results.

{{True, 626}, {False, 374}}
POSTED BY: Marco Thiel

Marco,

Yes, I'm sure they are generically correct. If I get a chance I'll look into your experiment to see what's going on.

Here are a few remarks that might be useful.

(1) Solve is invoking Roots so they will show similar behavior.

(2) Reduce is likely also mostly invoking Roots, with the addition of accounting for the possibility of leading coefficients vanishing.

(3) If one sets Quartics->False and works with parametrized Root objects the substitution problems should disappear. Root objects will not have issues because roots are isolated on the fly, so to speak (that is, after rather than before substitution).

(4) The "problem" with the quartic formula is actually not in the formula itself. It has to do with basic evaluation in Mathematica and in particular with detection of an indeterminate form. For the example that was initially posted we have this.

Aa = -30; Bb = 0; Cc = 13; Dd = 0; Ee = -1;
rts = Roots[aA x^4 + bB x^3 + cC x^2 + dD x + eE == 0, x];

Now substitute for four of the parameters.

rt2 = rts /. {aA -> Aa, bB -> Bb, cC -> Cc, eE -> Ee};
rt2[[1, 2]] // Simplify

(* Out[67]= -((1/(6*Sqrt[10]))*
      (Sqrt[26 - 529/(-11843 - 405*dD^2 + 
           9*Sqrt[5]*Sqrt[(-4 + 5*dD^2)*(4802 + 81*dD^2)])^(1/3) - (-11843 - 405*dD^2 + 
          9*Sqrt[5]*Sqrt[(-4 + 5*dD^2)*(4802 + 81*dD^2)])^(1/3)] + 
         Sqrt[52 + 529/(-11843 - 405*dD^2 + 
           9*Sqrt[5]*Sqrt[(-4 + 5*dD^2)*(4802 + 81*dD^2)])^(1/3) + (-11843 - 405*dD^2 + 
          9*Sqrt[5]*Sqrt[(-4 + 5*dD^2)*(4802 + 81*dD^2)])^(1/3) - (18*Sqrt[10]*dD)/ Sqrt[26 - 
          529/(-11843 - 405*dD^2 + 
              9*Sqrt[5]*Sqrt[(-4 + 5*dD^2)*(4802 + 
                    81*dD^2)])^(1/3) - (-11843 - 
             405*dD^2 + 9*Sqrt[5]*
                          Sqrt[(-4 + 5*dD^2)*(4802 + 81*dD^2)])^(1/3)]]))) *)

Here is a denominator that shows up. The fraction has dD as a factor in the numerator.

Sqrt[26 - 529/(-11843 - 405*dD^2 + 
       9*Sqrt[5]*Sqrt[(-4 + 5*dD^2)*(4802 + 81*dD^2)])^(1/3) - (-11843 - 405*dD^2 + (1/3)] /. dD -> 0

(* Sqrt[26 - 529/(-11843 + 882ISqrt[10])^(1/3) - (-11843 + 882ISqrt[10])^(1/3)] *)

This is what we call a hidden zero. A remarkably well hidden zero. The evaluation of substituting dD->0 into the full expression certainly sees a zero in the numerator, and does not notice the one in the denominator. So it evaluates the fraction to zero. This is all as I believe it should be. But the consequence is that an indeterminate form has instead become zero. Hence the numerically incorrect result after substitution.

(5) Limit will work harder to tease out the behavior of that indeterminate form and will end up with a numerically correct result corresponding to the substitution at hand.

POSTED BY: Daniel Lichtblau

Dear Daniel,

thank you very much for your reply. That is very useful indeed. I would have thought that all of Solve, Roots, Reduce should generally be correct. Solving these polynomials is certainly some basic functionality and I am convinced that Mathematica gets it right. I also do understand the problem of these hidden zeros; in fact, I have produced lots of them in my own code myself.

It was just a bit surprised at the results of my simulation. More than 30% error rate is a bit high for something that is generically ok; of course it might be due to my choice of parameters. I will also check whether something goes wrong when I compare the results of the two procedures; substitute first, or solve first.

I suppose that in most problems I would rather substitute the coefficients first anyway, but I could see myself committing the mistake of solving first, with the argument that the "heavy lifting" would only be done once, and then only substitutions would be needed.

I'll also try to investigate further.

Thank you, Marco

POSTED BY: Marco Thiel

Update, I have checked the

Quartics-> False 

option and as you say it does work for Solve and Roots- (not quite for Reduce). It does circumvent the idea that Yuriy had originally, but it does seem to work.

Thanks, Marco

POSTED BY: Marco Thiel
Posted 10 years ago

Marco,

I am not sure what is the trick that Root object does sine I am quite fresh to Mathematica and still have some blurry understanding of Mathematica features, but regarding Reduce from help we can get that:

By default Reduce does not use the general formulas for solving quartics in radicals Reduce[x^4 + 2 x^2 + 3 x + 4 == 0, x] // Last

x == Root[4 + 3 #1 + 2 #1^2 + #1^4 &, 4]

So that is why Reduce works better having Quartics -> False by default.

PS I am so glad the post got so many brilliant responses, many thanks to all contributors, I have a long way to go yet with Mathematica and this post is extremely useful in terms of learning.

Thanks, Yuriy

POSTED BY: Yuriy Ivanov

Dear Yuriy,

regarding Quartics, you might want to have a look at :

http://reference.wolfram.com/language/ref/Quartics.html

It is about how the roots are represented, implicitly or explicitly. The documentation page explains that quite nicely.

Cheers,

Marco

POSTED BY: Marco Thiel

Dear Daniel,

I think I found the problem in my approach. The sort command is a bit different from what I thought it was. I get many false positives like that - see post below.

Cheers and sorry, Marco

POSTED BY: Marco Thiel

I am surprised at that behaviour as well. I have filed a report.

Cheers, Marco

POSTED BY: Marco Thiel

So basically you are saying that you are surprised that

In[1]:= ClearAll["Global`*"]

In[2]:= N[
  Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, x]] /. {Aa -> -30,
  Bb -> 0,
  Cc -> 13,
  Dd -> 0,
  Ee -> -1}

Out[2]= x == -0.465475 - 3.2262*10^-9 I || 
 x == 0.465475 - 3.2262*10^-9 I || x == -0.465475 + 3.2262*10^-9 I || 
 x == 0.465475 + 3.2262*10^-9 I

In[3]:= Aa = -30;
Bb = 0;
Cc = 13;
Dd = 0;
Ee = -1;
N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0, x]]

Out[8]= x == 0.316228 || x == -0.316228 || x == 0.57735 || 
 x == -0.57735

So it looks as if you get different results whether you first declare the parameters and solve, or first solve and then substitute the parameters? The first does not appear to be correct, when substituted.

Also

Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> -30,Bb -> 0,Cc -> 13,Dd -> 0,Ee -> -1}

is

-1 + 13 x^2 - 30 x^4 == 0

The solution of which is:

Solve[ -1 + 13 x^2 - 30 x^4 == 0, x]
{{x -> -(1/Sqrt[3])}, {x -> 1/Sqrt[3]}, {x -> -(1/Sqrt[10])}, {x -> 1/Sqrt[10]}}

Cheers,

Marco

PS: This one

N[Roots[Aa x^4 + Bb x^3 + Cc x^2 + Dd x + Ee == 0 /. {Aa -> -30,Bb -> 0,Cc -> 13,Dd -> 0,Ee -> -1}, x]]

i.e. "first substitute, then solve" also gives the expected result.

POSTED BY: Marco Thiel
Posted 10 years ago

Marco, thanks for your contribution. Basically that is exactly what I am surprised at. If the algebraic solution (symbolic) is claimed to be general, then it should work no matter what the sequence of substituting of coefficients is, if I am not mistaken.

POSTED BY: Yuriy Ivanov
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract