Dr. StrangeNumbers or: How I Learned to Stop Worrying and Love Floating Point Arithmetic

11

1

The following is the program.

test[t_, dt_] := 
 Module[{}, For[ti = dt, ti <= t, ti = ti + dt, Print[ti];]; 
  Print[MemberQ[{0.01, 0.02}, ti]]; Return[0];]    
test[0.01, 0.001]

(* 0.001
...
False
0 *)

Obviously the result is wrong. Copy the above result, you will be surprised to see:

(* ...
0.009000000000000001`
0.010000000000000002`
... *)

Why is this so?

howard

Posted 2012-05-15T12:36:45.020

Reputation: 353

2

This is not a Mathematica question - it is a numerical analysis question: http://en.wikipedia.org/wiki/Round-off_error

– Mark McClure – 2012-05-15T12:41:22.617

7It's not the For[] loop, really. Try 0.008 + 0.001 // InputForm. Welcome to the fun of floating-point arithmetic! – J. M.'s ennui – 2012-05-15T12:41:45.037

Good,Thank you! – howard – 2012-05-15T13:16:15.577

What's the point of Module[{}, ...] with empty variable list? As far as I understand, it should not have any effect at all. – celtschk – 2012-05-16T13:22:09.917

@celtschk I've used it a few times as a cool parentheses, reserved for future use... Then I stopped doing that – Rojo – 2012-05-16T18:34:05.380

@Guesswho I see that you added the faq tag. I've typically been using it only when there are a minimum of eight duplicates. This question presently has only two. Maybe you can find others that should be marked as duplicates? (Possibly questions that have been deleted that should have been marked duplicate instead?)

– Mr.Wizard – 2015-08-04T08:38:46.733

I'll try to look for some, @Mr. Wizard; I believe this has been asked here in one form or another (which also makes searching difficult). – J. M.'s ennui – 2015-08-04T08:44:46.390

Answers

17

While it is a bad idea to compare floating point numbers, in this case I think something simpler is going on: you have an off-by-1 problem in your loop. See what your code does:

test[t_, dt_] := Module[{},
  For[ti = dt, ti <= t, ti = ti + dt, Print[ti];];
  Print[MemberQ[{0.01, 0.02}, ti]];
  Return[0];
  ]

then, after

test[0.01, .001];

look at ti:

ti
(*0.011*)

If you use < rather than <= in the condition, the final value will indeed be what you expect it to be:

test2[t_, dt_] := Module[{},
  For[ti = dt, ti < t, ti = ti + dt, Print[ti];];
  Print[MemberQ[{0.01, 0.02}, ti]];
  Return[0];
  ]

test2[0.01, .001]
ti

prints a list up to 0.009 and returns True.

In general, though, don't do things like that with floating-point numbers. Even if it works here, it will eventually fail. Observe:

(1 + $MaxMachineNumber) == $MaxMachineNumber
(*True*)

or, as JM points out,

1 + $MachineEpsilon == 1
(*True*)

acl

Posted 2012-05-15T12:36:45.020

Reputation: 19 146

2The more traditional example is 1 + $MachineEpsilon == 1... ;) – J. M.'s ennui – 2012-05-15T12:59:18.197

1@J.M. Tradition is for wimps :) – acl – 2012-05-15T13:01:07.930

7Bear in mind Equal applies an extra tolerance in Mathematica. The proper comparison is Block[{Internal`$EqualTolerance = -Infinity}, 1 == 1 + $MachineEpsilon] (= False) vs. Block[{Internal`$EqualTolerance = -Infinity}, 1 == 1 + $MachineEpsilon/2] (= True). – Oleksandr R. – 2012-05-15T13:55:50.163

@OleksandrR. I didn't know that, thanks – acl – 2012-05-15T14:05:50.437

5

This seems to work fine and is more in the spirit of the Mathematica way and perhaps illustrates why a Functional approach tends to lead to fewer bugs than the procedural approach:

test[t_, dt_] := {#, MemberQ[{0.01, 0.02}, #]} & /@ Range[dt, t, dt];

Which when run,

test[0.01, 0.001] // TableForm

gives:

Mathematica graphics

image_doctor

Posted 2012-05-15T12:36:45.020

Reputation: 9 964

1Replacing TableForm with InputForm reveals that one does not totally avoid floating-point error... – J. M.'s ennui – 2012-05-15T13:55:16.800

@J.M. that's interesting, it makes me curious as to why that happens, the increment doesn't seem to be anywhere near the limit of precision of of a floating point number. Is this a general effect with FP numbers or more specific to MMA? – image_doctor – 2012-05-15T15:26:07.167

3

it's general, see e.g. http://floating-point-gui.de/

– acl – 2012-05-15T15:45:29.627

@acl, Thanks, surprising how in a great number of years that's never knowingly tripped me up, but here because MemberQ is being used as a criterion, rather than a < or > based test, it become obvious. – image_doctor – 2012-05-15T16:00:43.627