8
$\begingroup$

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:

A curve under low precision

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:

A curve under high presision

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?

Improve this question
$\endgroup$
17
  • 1
    $\begingroup$ 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.. $\endgroup$
    – halirutan
    Commented Sep 19, 2012 at 11:57
  • 1
    $\begingroup$ @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? $\endgroup$
    – xzczd ♦
    Commented Sep 19, 2012 at 12:21
  • 3
    $\begingroup$ And thanks to @J.M. who found this answer. $\endgroup$
    – halirutan
    Commented Sep 19, 2012 at 12:49
  • 2
    $\begingroup$ @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. $\endgroup$
    – celtschk
    Commented Sep 19, 2012 at 15:09
  • 2
    $\begingroup$ @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). $\endgroup$
    – celtschk
    Commented Sep 20, 2012 at 9:19

1 Answer 1

Reset to default
11
$\begingroup$

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

Mathematica graphics

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}

Mathematica graphics

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.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ 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... $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Sep 19, 2012 at 15:20
  • $\begingroup$ 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. $\endgroup$
    – xzczd ♦
    Commented Sep 20, 2012 at 5:43
  • $\begingroup$ 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. $\endgroup$
    – xzczd ♦
    Commented Sep 20, 2012 at 5:43

Not the answer you're looking for? Browse other questions tagged or ask your own question.