Re: why extending numbers by zeros instead of dropping precision is a good idea
- To: mathgroup at smc.vnet.net
- Subject: [mg117846] Re: why extending numbers by zeros instead of dropping precision is a good idea
- From: Bill Rowe <readnews at sbcglobal.net>
- Date: Fri, 1 Apr 2011 02:35:55 -0500 (EST)
On 3/31/11 at 4:06 AM, fateman at eecs.berkeley.edu (Richard Fateman)
wrote:
>It is occasionally stated that subtracting nearly equal
>quantities from each other is a bad idea and somehow unstable
>or results in noise. (JT Sardus said it on 3/29/2011, for example.)
>This is not always true; in fact it may be true hardly ever.
>It is, for example, the essence of Newton iteration.
>Here's simple example. If you have an approximation x for the
>square root of n, then
>r[x,n] gives you a better approximation where
>r[x_,n_]:= x- (x^2-n)/(2*x).
>Let's try it out.
>x=1.0
>r[x,2]
>r[%,2]
>r[%,2]
>r[%,2]
>gets us to 1.41421, which is pretty good.
>If we try
>Do[x = r[x, 2], {100}]
>we see that x is still 1.41421.
>That's using machine hardware floats.
>Now let's try it with significance arithmetic, using more precision
>x=1.0`20
>two=2.0`20
>r[x,two]
>.... this gives (a highly precision version of) 1.5
>Let's see if we can get a really good approximation for sqrt(2).
>And print out the results each iteration
>Do[x = r[x, two];Print[x], {100}]
>Oh dear, we get lower and lower accuracy and precision,
>eventually resulting in a division by zero. This loop
>converges to the square root of 2 being Indeterminate.
>Now this is no surprise to some people, but it may be to some
>others. Indeed one could say "there is no bug there,
>Mathematica is doing exactly what it is documented to do."
>Nevertheless, some other people might find it surprising. Maybe
>even distressing. If you figure out what is going on (it is
>not that complicated...)... you see:
>It is the consequence of not computing x^2-n, the difference
>of two nearly equal quantities, accurately.
>Now you can avoid this problem by setting $MinPrecision to (say) 20.
An alternative would be to use FixedPoint. That is
In[27]:= FixedPoint[r[#, 2] &, 1.0`20]
Out[27]= 1.414213562373095049
In[28]:= Precision[%]
Out[28]= 18.2351
In[29]:= FixedPoint[r[#, 2.0`20] &, 1.0`20]
Out[29]= 1.41421356237309505
In[30]:= Precision[%]
Out[30]= 18.0154
In[31]:= N[Sqrt[2], 20]
Out[31]= 1.4142135623730950488
It seems FixedPoint is doing something behind the scene that
avoids the problem you describe above.
Note, I am not arguing either Mathematica is doing the right
thing or the wrong thing when using a Do as you did. Simply, my
experience with Mathematica tells me I generally get better
behavior when I use functions intended to do exactly what I want
than I do when I attempt to code my own version.
It seems FixedPoint is doing something behind the scene that
avoids the problem you describe above.