I am puzzled by the behavior of Mathematica machine precision with numbers approaching zero. This manifests itself, e.g., with FixedPoint and the like.
In the examples below I will use the following "service" function to take only a few values of long lists:
ClearAll[take];
take[n_Integer][li_List] := If[Length[li] >= Abs[n], Take[li, n], li];
Consider finding the solution of two different equations:
FindRoot[Sin[2.0 x] == x, {x, 1}]
FindRoot[Sin[0.5 x] == x, {x, 1}]
They return respectively:
{x -> 0.9477471335169905}
{x -> 2.6312502524050707*^-15}
where the second one is actually zero to machine precision. edit: no: where the second can be considered as zero in many machine precision computations, even if it is several orders of magnitude bigger than the smallest floating point number.
Now consider doing the same with FixedPoint:
a = FixedPoint[Sin[2.0 #] &, 1.]
b = FixedPoint[Sin[0.5 #] &, 1.]
The behavior of the functions is such that the two should converge to the respective solutions above. Actually both never end.
An analysis of the values obtained (with FixedPointList) shows that the problem with the first function is an oscillation in the last bit(s):
take[-5]@ FixedPointList[Sin[2.0 #] &, 1., 1000] // Differences
{-2.220446*^-16, 2.220446*^-16, -2.220446*^-16, 2.220446*^-16}
so that a real fixed point is never achieved. In this case, SameTest -> Equal is of help:
FixedPoint[Sin[2.0 #] &, 1., SameTest -> Equal]
0.9477471335169858
On the other hand, using SameTest -> Equal does not help with the other example, since:
FixedPoint[Sin[0.5 #] &, 1., SameTest -> Equal]
also never ends.
The same analysis as before, though, shows this:
take[-5]@ FixedPointList[Sin[0.5 #] &, 1., 10000] // Differences
{-3.7962812653512530891065876917712356`15.477468515471298*^-3010,
-1.8981406326756265445532938458856178`15.477468515471298*^-3010,
-9.490703163378132722766469229428089`15.477468515471298*^-3011,
-4.745351581689066361383234614714044`15.477468515471298*^-3011}
Actually there is no bottom! The negative exponents become arbitrarily large and two successive numbers are never seen to be equal to zero.
Now the question: is it reasonable that, having given as input MachinePrecision numbers, arbitrarily large exponents are not confined into the expected range? This is what Chop would do, but it is not applied and must be forced with a user-given SameTest, which seems a little unexpected here.
Why can MachinePrecision numbers have an arbitrary negative exponent?
Comments
Post a Comment