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.

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

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

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

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.

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.

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.

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]

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

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

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

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.