11
$\begingroup$

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?

$\endgroup$
7
  • 2
    $\begingroup$ This is not a Mathematica question - it is a numerical analysis question: en.wikipedia.org/wiki/Round-off_error $\endgroup$ May 15, 2012 at 12:41
  • 7
    $\begingroup$ It's not the For[] loop, really. Try 0.008 + 0.001 // InputForm. Welcome to the fun of floating-point arithmetic! $\endgroup$ May 15, 2012 at 12:41
  • $\begingroup$ Good,Thank you! $\endgroup$
    – howard
    May 15, 2012 at 13:16
  • $\begingroup$ What's the point of Module[{}, ...] with empty variable list? As far as I understand, it should not have any effect at all. $\endgroup$
    – celtschk
    May 16, 2012 at 13:22
  • $\begingroup$ @celtschk I've used it a few times as a cool parentheses, reserved for future use... Then I stopped doing that $\endgroup$
    – Rojo
    May 16, 2012 at 18:34

2 Answers 2

17
$\begingroup$

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*)
$\endgroup$
4
  • 2
    $\begingroup$ The more traditional example is 1 + $MachineEpsilon == 1... ;) $\endgroup$ May 15, 2012 at 12:59
  • 1
    $\begingroup$ @J.M. Tradition is for wimps :) $\endgroup$
    – acl
    May 15, 2012 at 13:01
  • 7
    $\begingroup$ Bear 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). $\endgroup$ May 15, 2012 at 13:55
  • $\begingroup$ @OleksandrR. I didn't know that, thanks $\endgroup$
    – acl
    May 15, 2012 at 14:05
5
$\begingroup$

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

$\endgroup$
4
  • 1
    $\begingroup$ Replacing TableForm with InputForm reveals that one does not totally avoid floating-point error... $\endgroup$ May 15, 2012 at 13:55
  • $\begingroup$ @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? $\endgroup$ May 15, 2012 at 15:26
  • 3
    $\begingroup$ it's general, see e.g. floating-point-gui.de $\endgroup$
    – acl
    May 15, 2012 at 15:45
  • $\begingroup$ @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. $\endgroup$ May 15, 2012 at 16:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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