2

I already read the other discussion here on SO but I still cannot understand why the code does not work as expected.

In my code I check that k is approximately equal to 1, but still it does not work:

if approx(k, 1, eps)
    display('Approx equal');
else
    k
    k == 1
    approx(k, 1, eps)
    abs(k - 1) < eps
end

Instead of displaying the string, this is the result (I have enabled format long):

k = 1.000000000000000
ans = 0
ans = 0
ans = 0

This is baffling! I also tried to increase the error by k * 1e20 but the result is still 1e20... What to do here?

Note: although not really relevant to the question, here is the definition of approx:

function r = approx(a, b, tol)
    r = a <= b + tol && a >= b - tol;
end

EDIT: I changed the approx function to avoid catastrophic cancellation:

function r = approx(a, b, tol)
    r = a <= b + tol && a + tol >= b; % assumes tol is positive
end
Community
  • 1
  • 1
rubik
  • 8,814
  • 9
  • 58
  • 88
  • 1
    I might just point out that I cannot reproduce this (R2015a OS X 64-bit). What version of Matlab are you using? – dfrib Apr 10 '16 at 08:45
  • @dfri I'm using R2015b on Windows 32-bit. – rubik Apr 10 '16 at 08:47
  • Please change your code to print at least 17 significant digits, this will allow us to reproduce your problem. Code to print 17 significant digits can be found [in my answer here](http://stackoverflow.com/a/35626253/2732801) – Daniel Apr 10 '16 at 09:03
  • @Daniel Oh, indeed it is `k = 0.99999999999999967`. – rubik Apr 10 '16 at 09:08
  • @dfri I removed the subtraction from `approx` (now it's `r <= b + tol && a + tol >= b`) and found that `approx(k, 1, 1e-15) = 1` but `approx(k, 1, 1e-16) = 0`. And since `eps ~ 2e-16` I was considering using `approx(k, 1, 10*eps)` in my ifs. Does that sound right? – rubik Apr 10 '16 at 09:11
  • @rubik: We don't know you application and how large the errors can be, but setting `tol` to the maximum possible error which can happen is probably the solution. – Daniel Apr 10 '16 at 16:28

1 Answers1

1

Try displaying:

k - 1.000000000000000

I think you will find that there is a tiny but non-zero residual value.

Raymond Hettinger
  • 216,523
  • 63
  • 388
  • 485