Group Abstract Group Abstract

Message Boards Message Boards

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

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

Yuriy Ivanov, NTNU
Posted 11 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
Sort By:
POSTED BY: Daniel Lichtblau

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

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
POSTED BY: Marco Thiel
Yuriy Ivanov
Yuriy Ivanov, NTNU
Posted 11 years ago
POSTED BY: Yuriy Ivanov

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 BY: Udo Krause
POSTED BY: Marco Thiel
POSTED BY: Frank Kampas

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

Cheers, Marco
POSTED BY: Marco Thiel
Yuriy Ivanov
Yuriy Ivanov, NTNU
Posted 11 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
Yuriy Ivanov
Yuriy Ivanov, NTNU
Posted 11 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

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

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

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[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

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

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

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

Perhaps the "errors" you see are due to imprecision in machine precision calculations.
POSTED BY: Frank Kampas

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

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

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

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
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard