why not 10^-3*10 = 10^-2? or why (10^-3 *10 < 10^-2) returns true? How to deal with this?

Hi everyone

Please help me.

I start a small loop to find some maximum feasible alpha for some bigger task in my code. alpha is started at 10^-16,and here is the summary of the code

alphak = 10^-16;

while condition if alphak < 0.01 (1) alphak = alphak*10; else alphak = alphak + 0.01; if alphak >= 0.98 break; end end end

Unfortunately, even when alphak reaches 0.01 (shown when I debug), the next value alphak takes is 0.1 instead 0.02. This creates a fatal error for my bigger code.

Please let me know why condition (1) above fails to work? And how to fix it?

Thank you very much

Nhat

[Matt J]

Please highlight your code (only) and use the

toolbar button to format it in a readable font.

Answer 1

[James Tursa]

Floating point arithmetic is in general not exact. Your 10^-16, for instance, cannot be represented exactly in IEEE double precision arithmetic. You need to rework your conditions and tests to account for this (use tolerances or switch to integers etc.). E.g., your exact numbers you are using are:

>> x = 10^-16
  x =
    1.0000e-016
>> for k=1:15;disp(num2strexact(x));x=x*10;end
9.999999999999999790977867240346035618411149408467364363417573258630000054836273193359375e-17
1.00000000000000007770539987666107923830718560119501514549256171449087560176849365234375e-15
1.00000000000000015659149039876225977948004212947363811281320522539317607879638671875e-14
1.00000000000000015659149039876225977948004212947363811281320522539317607879638671875e-13
1.0000000000000001818350393658346375526553562185227974623558111488819122314453125e-12
1.00000000000000010105568267120302867849435113356548754381947219371795654296875e-11
1.00000000000000016567916802690831577782315520153133547864854335784912109375e-10
1.000000000000000269076744596036775136749241710276692174375057220458984375e-9
1.000000000000000351794805851339542623890110917272977530956268310546875e-8
1.000000000000000351794805851339542623890110917272977530956268310546875e-7
1.000000000000000378264585453036428219775189063511788845062255859375e-6
1.0000000000000004206162328157514451731913140974938869476318359375e-5
1.0000000000000004544975507059234587359242141246795654296875e-4
1.0000000000000004544975507059234587359242141246795654296875e-3
1.00000000000000054123372450476381345652043819427490234375e-2

You can find num2strexact here:

http://www.mathworks.com/matlabcentral/fileexchange/22239-num2strexact-exact-version-of-num2str

Answer 2

[Jan]

Instead of the floating point operations, you can use more stable integer operations:

for k = -16:-2
   alphak = 10^k;
   if condition
     break;
   end

     if k == -2
       for k2 = 0.01:0.01:0.98
         alphak = k2;
         if condition
           break;
         end
       end
     end
  end

This has drawbacks: 1. it looks ugly, 2. alphak=0.01 is checked twice, 3. the k2 loop relies on floating point arithmetics again, but in a more stable way.