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 am trying to solve an equation using the following piece of code:

f[t_] = Aa*t^4 + Bb*t^3 + Cc*t^2 + Dd;

s = Solve[f[t] == 0, t];

t1 = t /. s[[1]];
t2 = t /. s[[2]];
t3 = t /. s[[3]];
t4 = t /. s[[4]];  

And it gives me four symbolic formulas.

Then, I try to give numbers to the Aa, Bb, Cc, Dd parameters:

Aa = 1.2347931365926803*10^30 ;
Bb = 4.6356322739291924*10^23 ;
Cc  = 4.350806821541392*10^16 ;
Dd = -1.2562823055999998*10^6 ;

and the formulas provide the results:

-7.18016*10^-7 - 6.17065*10^-7 I
5.30307*10^-7 + 6.17065*10^-7 I
-7.18016*10^-7 + 6.17065*10^-7 I
5.30307*10^-7 - 6.17065*10^-7 I

But, if I solve using the numerical values of Aa, Bb, Cc, Dd from the beginning, I get:

-1.10255*10^-6
-9.38545*10^-8 - 9.99928*10^-7 I
-9.38545*10^-8 + 9.99928*10^-7 I
 9.14844*10^-7

I have also noticed that this behavior is sensitive to the values I give to the parameters. Using these instead,

Aa = 2.8262007068794462*10^25;
Bb = -5.389762192271189*10^21;
Cc = 2.5696632603084144*10^17;
Dd = -1.2213607220736102*10^6;

the problem seems to still be there, but not as serious. Any ideas on what might be causing this and maybe a way to avoid it?

share|improve this question
add comment

2 Answers

We encounter here numerical incompatibility of reducing the exact solutions written in terms of radicals by some arbitrary numerical parameters and merely numerically solving of corresponding algebraic equations.
In fact it is not really surprising taking into account values of given parameters. We might guess this effect appears more significant when the parameters get higher absolute values and become more asymmetric.

The documentation of ToRadicals says (see e.g. Details and Options ):

If Root objects in expr contain parameters, ToRadicals[expr] may yield a result that 
is not equal to expr for all values of the parameters. 

At first we define:

f[t_] := a t^4 + b t^3 + c t^2 + d

ncoeff = { a :> 1.2347931365926803*10^30, b :> 4.6356322739291924*10^23, 
           c :> 4.350806821541392*10^16, d :> -1.2562823055999998*10^6  };

fn[t_] := f[t] /. ncoeff

If we prevent getting solutions in terms of radicals there'll be no problems anymore. We can use e.g. the Quartics -> False option of Solve:

{t1e, t2e, t3e, t4e} = t /. Solve[ f[t] == 0, t, Quartics -> False]
{Root[d + c #1^2 + b #1^3 + a #1^4 &, 1], Root[d + c #1^2 + b #1^3 + a #1^4 &, 2], 
 Root[d + c #1^2 + b #1^3 + a #1^4 &, 3], Root[d + c #1^2 + b #1^3 + a #1^4 &, 4]}

now we have:

{t1e, t2e, t3e, t4e} /. ncoeff
 {-1.10255*10^-6, 9.14844*10^-7, 
  -9.38545*10^-8 - 9.99928*10^-7 I, -9.38545*10^-8 + 9.99928*10^-7 I}

and this result is numerically equal to numerical solutions.

t /. Solve[ fn[t] == 0, t]
{ -1.10255*10^-6, -9.38545*10^-8 - 9.99928*10^-7 I, 
  -9.38545*10^-8 + 9.99928*10^-7 I, 9.14844*10^-7}   

as well as

t /. FindRoot[ fn[t], {t, #}] & /@ {-10^-6, 10^-6, 10^-7 (-1 - 10 I), 10^-7 (-1 + 10 I)}
{ -1.10255*10^-6, 9.14844*10^-7, 
  -9.38545*10^-8 - 9.99928*10^-7 I, -9.38545*10^-8 + 9.99928*10^-7 I}

Now there is no ambiguity of the results.

share|improve this answer
    
Thank you for your answer Artes! The thing is that I was going for the analytical solution, because I want to use it in a fortran code. The way Solve provided it mislead me...I'm guessing there is no way around these Root objects, right? I've been reading this post so... –  nahn Jan 1 at 23:11
    
Also, I've found this post on general formula for solving quartic (4th degree polynomial) equations. They give a general solution based on what solve gives, but in another answer someone goes the old fashion way. I'll see what I can find maybe some constraints or something. –  nahn Jan 1 at 23:21
    
@nahn In general the class of equations solvable in terms of radicals is much narrower than those solvable in terms of root objects (the former is a subset of the latter). The idea of putting symbolic solutions found with Mathematica to transform them in Fortran sounds as a total misconception. Nonetheless there are no fundamental barriers. –  Artes Jan 1 at 23:25
    
I think you are right Artes. The key word being "solutions". Taking, for example a Curl(V), where V a 3D vector, expanding it through mathematica into its components and transforming it to fortran is a safe procedure. Solving an equation (unless you really know how mathematica works) turned out to be not a good idea. I've implemented the actual analytical solution and still could not get the numbers right. Thank you! –  nahn Jan 2 at 1:43
    
Your comments are not quite relevant to the essence of your question. You haven't clarified whether my post answers your question or not. If this is the case you should update your question providing clear idea what you would like to get. –  Artes Jan 5 at 0:02
add comment

The correct result can be obtained by switching to arbitrary precision calculations:

f[t_] := a t^4 + b t^3 + c t^2 + d
{t1e, t2e, t3e, t4e} = t /. Solve[f[t] == 0, t];
ncoeff = Rationalize[#, 0] &@{a :> 1.2347931365926803*10^30, 
    b :> 4.6356322739291924*10^23, c :> 4.350806821541392*10^16, 
    d :> -1.2562823055999998*10^6};
N[{t1e, t2e, t3e, t4e} /. ncoeff, 6]
{-9.38545*10^-8 - 9.99928*10^-7 I, -9.38545*10^-8 + 9.99928*10^-7 I, 
-1.10255*10^-6, 9.14844*10^-7}
share|improve this answer
add comment

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.