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

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?

-

## marked as duplicate by whuber, Oleksandr R., Brett Champion, Verbeia♦, Mr.Wizard♦Sep 23 '12 at 0:19

The last sentence in my previous comment sounds harsh. It is not meant like that! It's just a hint where you could start looking in your lecture notes or when searching the net. I don't know how to say this in a different way.. I'm German. Bratwurst and Rucksack, you know.. –  halirutan Sep 19 '12 at 11:57
@halirutan Haha, never mind, and in fact I really don't have numerical error analysis in my university…at least it's not a required course for my major and in my memory I never see this in the list of elective course, too… our lessons always talk about those analytical solutions and say little about numeric solutions, usually the books just say "the numeric solution should be done with computer". So, can you explain the details? –  xzczd Sep 19 '12 at 12:21
And thanks to @J.M. who found this answer. –  halirutan Sep 19 '12 at 12:49
@xzczd: The basic principle is this: If you only have a certain number of digits of precision, you'll lose something. Consider for example the expression $10^{100}(1+10^{-100}-1)$. If you have it as analytical formula, it's obvious that the result is $1$. But if you calculate e.g. to 5 digits, then 1.0000+1.0000e-100 = 1.0000 since the 1.0000e-100 is rounded away. Then subtracting 1.0000 again gives 0.0000 and multiplying with 1.0000e100 won't change that. –  celtschk Sep 19 '12 at 15:09
@xzczd: The problem is that exact arithmetcs is often orders of magnitude slower than machine arithmetics (and in many cases, an exact solution will not even be possible where numerics gives good results). I think the right thing to do is to identify key parts where precision may be needed (look out for those sums and differences) and only rewrite those to either use more precise (or, where possible, exact) methods, or to avoid lossy operations by rewriting the functions in question (for example, by using the series close to a critical point, as you did in your last block). –  celtschk Sep 20 '12 at 9:19

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 420*^-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.

-
To see why 2 - 2 E^x + 2 x + x^2 is so troublesome to evaluate numerically (due to subtractive cancellation), consider the Maclaurin series for the exponential function... –  Ｊ. Ｍ. Sep 19 '12 at 15:20
I've read the answer and the links given by J.M and @Mr.Wizard. Let me try to retell it to see if I've understand it: The error occurs because when doing numeric calculations Mma will keep the precision step by step i.e. several steps of approximation have been done before the final result comes out. –  xzczd Sep 20 '12 at 5:43
If so, well, it's a little different from my imagination, since I think the approximation will only be done in the last step and the precision before the last step will be as high as possible (in my mind it almost equals to infinity…), just like what I turn to do when calculating something manually…and now I notice that my imagination is unreasonable. –  xzczd Sep 20 '12 at 5:43