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$
  • 2
    $\begingroup$ This is not a Mathematica question - it is a numerical analysis question: en.wikipedia.org/wiki/Round-off_error $\endgroup$ – Mark McClure May 15 '12 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$ – J. M. will be back soon May 15 '12 at 12:41
  • $\begingroup$ Good,Thank you! $\endgroup$ – howard May 15 '12 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 '12 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 '12 at 18:34
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$
  • 2
    $\begingroup$ The more traditional example is 1 + $MachineEpsilon == 1... ;) $\endgroup$ – J. M. will be back soon May 15 '12 at 12:59
  • 1
    $\begingroup$ @J.M. Tradition is for wimps :) $\endgroup$ – acl May 15 '12 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$ – Oleksandr R. May 15 '12 at 13:55
  • $\begingroup$ @OleksandrR. I didn't know that, thanks $\endgroup$ – acl May 15 '12 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$
  • 1
    $\begingroup$ Replacing TableForm with InputForm reveals that one does not totally avoid floating-point error... $\endgroup$ – J. M. will be back soon May 15 '12 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$ – image_doctor May 15 '12 at 15:26
  • 3
    $\begingroup$ it's general, see e.g. floating-point-gui.de $\endgroup$ – acl May 15 '12 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$ – image_doctor May 15 '12 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.