Accuracy issues when reading data from an input file

Hi

A simple question. Please find the code below. I'm experiencing an accuracy issue. Tried %f, %e and %g in fscanf but none solves my problem.

Input file contains:
    0
    0.1
    0.2
    0.3
    0.4
    0.1
    0.2
    0.3
    0.4
    0.5
Code:

meshf = fopen('./test_accuracy.dat', 'r');
all = fscanf(meshf, '%f', [1,10]); 
all = all';
startN = all(1:5);
endN = all(6:10);
dl = endN - startN;
for j = 1:4
    dl_check(j,1) = dl(j) - dl(j+1);
end
fclose(meshf)

Output: dl =

   0.100000000000000
   0.100000000000000
   0.100000000000000
   0.100000000000000
   0.100000000000000

dl_check =

1.0e-016 *

                   0
   0.277555756156289
  -0.555111512312578
   0.555111512312578

My Question: Typically we would expect dl_check to be all zeros right!? What format flags should I use to prevent this problem.

Many Thanks

Amar

Answer 1

[Andrew Newell]

You can't do better than that because the decimal numbers are really represented by binary numbers, which can't represent a lot of these fractions any better. If you try rounding the numbers to the first decimal, you'll get the same answer:

all - round(all*10)/10

returns all zeros.

For most purposes that's accurate enough: all the numbers in dl_check are at most 2*eps in magnitude. If you need to do better than that , you'll need extended precision arithmetic (as in, for example, the Symbolic Toolbox).

[Walter Roberson]

http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F

[Amarpal]

'round' just rounds off to the nearest integer, so how do you specify, which decimal place to round off to? I found that the function 'roundoff' is the one to be used, but I can't get it running on the version of MATLAB I have (7.6.0 - R2008a). Any work-arounds please?

[Andrew Newell]

I multiplied by 10 before rounding and then divided by 10. That is the work-around.

[Walter Roberson]

You cannot _really_ round to a decimal point in any binary floating point representation. [0.1 decimal] is [0.0(0011)* binary] where the (0011) is repeated infinitely. No matter how you round that in the binary, the result will not equal [0.1 decimal] because [0.1 decimal] requires an infinite binary number. http://en.wikipedia.org/wiki/Binary_numeral_system

[Andrew Newell]

I guess I didn't really make that clear in my response.

[Amarpal]

"Andrew Newell 31 minutes ago

I multiplied by 10 before rounding and then divided by 10. That is the work-around."

Sorry, I didn't get that. Could you please expand? Even without multiplying and dividing by 10, the function 'round' results in zeroes in dl_check.

Also, I found an external function 'roundoff' on the file exchange server (http://www.mathworks.com/matlabcentral/fileexchange/1374-roundoff/all_files). It results in all zeroes for upto 16 decimal places of accuracy. 17th decimal place onwards the same problem of having non-zeroes as before comes up.

I'm still interested to know about your *10 /10 logic.

Your next comment came up as I was writing this. Yeah these things do come up to mind....had it in electronics in my pre university. Has been a bit of a disconnection now. I'm doing a PhD in mechanical engineering these days!

[Andrew Newell]

Suppose you want to round 0.11 to get 0.1. Multiply by 10 to get 1.1; round to get 1; divide by 10 to get 0.1.

[Andrew Newell]

As for the file in the File Exchange, 16 decimals is the same as you are getting in your original question.

[Amarpal]

Actually not. 17 decimal places gives non-zero entries, but 16 works just fine. Checked it!

[Walter Roberson]

You will NOT be able to get rid of the problem of that 17th decimal place. Not as long as you are using binary floating point arithmetic. The problem is inherent in using any number base other than one which is an exact multiple of the denominator of the fraction you are trying to represent. Try, for example, *exactly* representing 1/7 in base 10 with using only a finite number of positions. Or exactly representing 1/3 in base 2, 4, or 5 -- but 1/3 in base 6 is no problem as [1/3 decimal] = [0.2 senary]