ReplaceAll[] and Limit[] don't give correct results for this expression under extreme variables

Question

Possible Duplicate:
Funny behaviour when plotting a polynomial of high degree and large coefficients

1/x^2 + (3 + x)/(6 (1 - Exp[x] + x))

——This is a expression that seems nothing special, so does its limit when x->0:

Limit[1/x^2 + (3 + x)/(6 (1 - Exp[x] + x)), x -> 0]
(* => 1/12 *)

But it becomes strange when trying to calculate a x near zero,

1/x^2 + (3 + x)/(6 (1 - Exp[x] + x)) /. x -> 0.00001
(* => -19403.7 *)
1/0.00001^2 + (3 + 0.00001)/(6 (1 - Exp[0.00001] + 0.00001))
(* => -19403.7 *)
Limit[1/x^2 + (3 + x)/(6 (1 - Exp[x] + x)), x -> 0.00001]
(* => 0. *)

Many people (including me) may check the curve of the expression as the first reaction for these results:

Plot[1/x^2 + (3 + x)/(6 (1 - Exp[x] + x)), {x, 0, 0.00001}]

And see:

Aha, it explains the matter! The curve is actually oscillating near zero!… but wait, try this:

Plot[1/x^2 + (3 + x)/(6 (1 - Exp[x] + x)), {x, 0, 0.001}, WorkingPrecision -> 30]

Now we see:

---

So all of these turn out to be a story about precision once again (why I always encounter this! ). I also found some resource for this error:

Series[1/x^2 + (3 + x)/(6 (1 - Exp[x] + x)), {x, 0, 10}];
Normal[%] /. x -> 0.00001
(* => 0.0833332 *)

N[1/x^2 + (3 + x)/(6 (1 - Exp[x] + x)) /. x -> 10^-5, 6]
(* => 0.0833332 *)

But, can this problem be solved only with ReplaceAll and Limit? Also, we see that Mma doesn't give warnings for these great error (while this sort of thing frequently occurs when I use NDSolve though the error seems not that big…it's another story), any way to make Mma give a warning message or something?

Answer 1

Lets summarize my comments into an answer.

When you inspect your expression a bit more you can see that you can drop some parts and still get an expression that's causing trouble:

expr = Numerator[
  Together[Most@ExpandAll[1/x^2 + (3 + x)/(6 (1 - Exp[x] + x))]]]

(*
  2 - 2 E^x + 2 x + x^2
*)

Plot[expr, {x, 0, 10^-6}]

Lets have a look what happens when we work with expr. First, we can just insert a numerical value which has machine precision and where the errors that happen are not traced.

expr /. x -> 4.*^-7
(*
  Out[22]= 1.0488*10^-16
*)

By using the . in the number we say please use the machine arithmetic. We don't know how reliable this result is but we can use Mathematicas arbitrary-precision arithmetics to get a n answer where we know how reliable it is. You have several options to input numbers for arbitrary-precision. If you give more digits than a machine-precision number has then Mathematica switches automatically

x1 = 4.00000000000000000*^-7;
Precision[x1]
expr /. x -> x1
Precision[%]

(*
Out[100]= 17.6021

Out[101]= -2.133*10^-20

Out[102]= 3.727
*)

See how the precision drops due to the evaluation of the expression! Usually, arbitrary-precision numbers are assembled with backticks. For instance 4`20*^-7 gives a number with precision 20. Please have a look at this tutorial and all the other pages about numerics in the documentation of Mathematica.

The most used function for numerical evaluation is maybe N. Here you can input an exact expression and evaluate it to whatever precision you like. Note that I use an exact number 4*^-7 with $\infty$ precision as input.

N[expr /. x -> 4*^-7, 20]
Precision[%]
(*    
Out[107]= -2.1333335466666837333*10^-20

Out[108]= 20.
*)

Therefore, N calculated this expression having a precision of 20 in the end. All this can now be applied to your Limit example

N[Limit[1/x^2 + (3 + x)/(6 (1 - Exp[x] + x)), x -> 10*^-6], 15]
Limit[1/x^2 + (3 + x)/(6 (1 - Exp[x] + x)), x -> 10.`35*^-6]
(*
Out[145]= 0.0833332222221759

Out[146]= 0.08333322222217592601410947237
*)

Or of course the Plot. Read the next example carefully. The first one says use machine numbers, the second one says use arbitrary-precision numbers with a precision of $MachinePrecision

Plot[1/x^2 + 3/(6 - 6*E^x + 6*x), {x, 0, 1.*^-6}, 
   WorkingPrecision -> #] & /@ {MachinePrecision, $MachinePrecision}

I'm not sure whether I really answered all of your questions and whether every detail is correct to the point. I kind of left my circle of competence here. I know numerics and how to watch out for errors, but I'm far from being an numerical analyst.