Re: Round-off error?
- To: mathgroup at smc.vnet.net
- Subject: [mg131932] Re: Round-off error?
- From: Itai Seggev <itais at wolfram.com>
- Date: Sat, 2 Nov 2013 02:25:43 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-outx@smc.vnet.net
- Delivered-to: mathgroup-newsendx@smc.vnet.net
- References: <20131031034715.CCF8C6A6B@smc.vnet.net>
> Hi everyone,
>
> The following lines of code are stumping me.
>
> ClearAll[y]
> dt = 0.15;
> Do[y[i*dt] = "hi there", {i, 0, 20}]
> Table[{i, y[(i + 1)*dt], y[i*dt + dt]}, {i, 0, 19}]
>
> Here is the output:
>
> 0 hi there hi there
> 1 hi there hi there
> 2 hi there hi there
> 3 hi there hi there
> 4 hi there hi there
> 5 hi there y(0.9)
> 6 hi there y(1.05)
> 7 hi there hi there
> 8 hi there hi there
> 9 hi there hi there
> 10 hi there hi there
> 11 hi there hi there
> 12 hi there y(1.95)
> 13 hi there hi there
> 14 hi there hi there
> 15 hi there hi there
> 16 hi there hi there
> 17 hi there hi there
> 18 hi there y(2.85)
> 19 hi there hi there
>
>
> It seems that y[i*dt+dt] does not match y[(i+1)*dt], perhaps due to
> numerical round-off? If you switch to dt=15/100, the problem goes away.
> You might also try different values for dt to see what happens.
>
>
> If it is numerical round-off, then why does this evaluate to true?
>
> 6*.15 == 5*.15 + .15
> True
The first thing to do whenever Mathematica does something unexpected is to look
at the FullForm.
In[1]:= 6*.15 //FullForm
Out[1]//FullForm= 0.8999999999999999`
In[2]:= 5*.15 + .15 //FullForm
Out[2]//FullForm= 0.9`
So clearly these are two different numbers. And the answer to your last
question is easy--Equal by default ignores the last seven bits in comparison,
pricesly to avoid too many False due to machine precision error. You can
control this using the variable Internal`$EqualTolerance:
In[3]:= Block[{Internal`$EqualTolerance = 0}, 6*.15 == 5*.15 + .15]
Out[3]= False
Note that this variable is specified using decimal digits, which is why its
default value is 2.10721 == N[Log[10, 2^7].
It's also easy to understand why it goes away with 15/100--you're now using
exact numbers and you end up with the exact result 9/10 in both cases.
However, you have an added twist, which is that you're using a function
definition / pattern. (And this is something I didn't fully understand until I
started poking at this, so thanks for the good question!) All patterns,
including literal patterns like 0.9`, are hashed for efficiency reasons. On
the other hand Equal is no longer strickly transitive (since two numbers which
different from a third by 6 bits need not differ from each other by less than
7). So machine numbers are effectively binned in the hash process to provide
something fast and transitive. You are suffereing because sometimes the two
computations are in the same bin and sometimes they are not.
In[31]:= Internal`HashSameQ[.15*6, .15*5 + .15]
Out[31]= False
In[32]:= Internal`HashSameQ[.15*5, .15*4 + .15]
Out[32]= True
This binning only affects patterns/definitions tied to specific approximate
numbers. If you want to avoid the binning, you could rewrite your definition
in the form
y[x_] /; x == num = "Hi there"
This will make sure that two numbers which are Equal trigger the same
definition, but at performance cost. How to balance those is dependent on your
needs and preferences.
Here I show that this change actually works. Notice that I'm using With in the
Do loop to make sure that num is evaluated before the definition of y is
created, which is necessary when generating definitions programmatically.
In[27]:= ClearAll[y]
dt=0.15;
Do[With[{num=dt*i}, y[x_] /; x==num="hi there"],{i,0,20}]
Table[{i,y[(i+1)*dt],y[i*dt+dt]},{i,0,19}]//TableForm
Out[30]//TableForm=
0 hi there hi there
1 hi there hi there
2 hi there hi there
3 hi there hi there
4 hi there hi there
5 hi there hi there
6 hi there hi there
7 hi there hi there
8 hi there hi there
9 hi there hi there
10 hi there hi there
11 hi there hi there
12 hi there hi there
13 hi there hi there
14 hi there hi there
15 hi there hi there
16 hi there hi there
17 hi there hi there
18 hi there hi there
19 hi there hi there
--
Itai
Itai Seggev
Mathematica Algorithms R&D