Persistent output problems with a for loop block!

Hello everyone. I was trying to run a certain code block:

clear;
E1 = 70e9;
E2 = E1/2;
E3 = E1/4;
t  = 0.0125;
L1 = 0.1;
L2 = 0.3;
L3 = 0.2;
P1 = 50e3;
P2 = 20e3;
P3 = 30e3;
w1 = 0.01;
w2 = 0.025;
el_L=0.1;
%BAR - 1 
B1_LWN=(L2/el_L)+1;
%BAR - 2 
B2_LWN=(L2/el_L)+1;
%BAR - 3 
B3_LWN=(L1/el_L)+1;
%BAR - 4 
B4_LWN=(L3/el_L)+1;
totalNodes=B1_LWN+B2_LWN+B3_LWN+B4_LWN-3;
x=zeros(1,totalNodes);
y=zeros(1,totalNodes);
NWL=[B1_LWN,B2_LWN,B3_LWN,B4_LWN;   %number of length-wise nodes of beams 
     L2,L2,L1,L3;                   % length of beams in (m)
     0,0,0.3,0.4;                   % starting x coordinate of beam
     0.3,0.3,0.4,0.6;               % ending x coordinate of beam
     0.005,0.02,0.0125,0.0125; ];    % y coordinate of element cg
%BEAM NODE COORDINATES
node=0;
for k=1:4
    inc=NWL(2,k)/((NWL(1,k)-1));
    Xcoords=NWL(3,k):inc:NWL(4,k);
                   for i=1:(NWL(1,k)-1)
                   node = node + 1;
                             if(k==4&&i==NWL(1,4)-1)
                             x(1,node) = Xcoords(i); 
                             y(1,node) = NWL(5,k); 
                             node=node+1;
                             x(1,node) = Xcoords(i+1);
                               y(1,node) = NWL(5,k); 
                             else
                             x(1,node) = Xcoords(i); % nodal x-position
                             y(1,node) = NWL(5,k);   % nodal y-position
                             end
                     end 
  end

So when I run the code I have no syntax errors. However, if I replace the values of the 1st and 2nd row of matrix NWL with the SAME constant values (i.e: B1_LWN replaced with 4 as when you calculate it's value from the %BAR-1 commented section, you end up getting 4) for all of the 1st and 2nd rows, I get an entirely different output for vectors x and y!!! Can someone explain why this is happening? It makes no sense at all. I tried everything. Thanks.

[Guillaume]

Rather than adding a blank line between each line of code which does not make it any easier to read, just copy/paste your code then press the {} Code button (which will just put two spaces at the start of each line).

Answer 1

[Tim Berk]

That's crazy!!

The problem is caused by numerical precision. When Matlab calculates 0.3/0.1, this is not exactly 3. I'm sure someone else can explain why, but apparently Matlab estimates divisions. Try

isequal(0.3/0.1,3)

or even better

0.3/0.1-3

(The latter does not return zero).

A workaround to solve your problem is to round the values in B1_LWN, B2_LWN etc. Because

isequal(round(0.3/0.1),3)

does return 1.

[José-Luis]

I haven't gone through OP's post because it is awfully formatted. If you are right, that's not the way to solve the problem.

This is due to the limited precision of ieee 754 doubles:

https://en.wikipedia.org/wiki/Double-precision_floating-point_format

Only a subset of decimal numbers can be exactly represented that way. When you are using logical comparisons on doubles, you should use a tolerance. round is too blunt of a tool. You should look into eps() instead.

[Guillaume]

Matlab does not "estimate divisions".

The same way that you cannot write down exactly the numer 3/10 in base 10 (0.333333.... you'll have to stop writing the 3s somewhere), computers cannot store all numbers exactly in base 2. 0.1 is actually such a number. The nearest value is stored instead (off by about 1e-17).

That's a very well known issue that you'll face when doing computation on a computer in any language unless you use symbolic computation which is orders of magnitude slower than floating point maths.

In most cases the infinitesimal error is negligible and can be ignored but it does accumulate in computations. This is why you can never use == (or isequal and co.) to compare floating point numbers. They may differ at the last decimal due to accumulated rounding even though they theoretically should be equal).

The way to compare floating point numbers is to check that their difference is smaller than an arbitrary value (dependent on their magnitude):

abs((0.3/0.1) - 3) <= 1e-15

[Stephen23]

"but apparently Matlab estimates divisions"

No, actually MATLAB performs perfectly correct binary divisions of the binary data that all numeric variables are stored as. In just the same way that you cannot write 1/3 as a finite decimal, there is no way to store 0.1 as a finite binary. This is not specific to MATLAB, but applies to all floating-point variables on almost all computers in the world. That some beginners have not yet learned that binary floating-point numbers are NOT decimal numbers is not any fault of MATLAB's, and does not mean that MATLAB "estimates division".

"A workaround to solve your problem is to round the values..."

If the values contain decimal fractions then using round introduces artifacts into the data that do not actually exist, in which case the real solution is to compare the absolute difference against a tolerance value. If the values should be integer then round or conversion to integer is quite appropriate. Read these to know why round is not a general solution for decimal values:

https://www.mathworks.com/matlabcentral/answers/321709-strange-behaviour-of-matlab-find-command#comment_444258

https://www.mathworks.com/matlabcentral/answers/350303-why-0-6-is-not-equal-to-0-2-3#comment_472139

Before giving such advice you might like to actually read about the topic that you are discussing, such as in the MATLAB documentation: "..computations sometimes yield mathematically nonintuitive results. It is important to note that these results are not bugs in MATLAB."

source: https://www.mathworks.com/help/matlab/matlab_prog/floating-point-numbers.html

and also in this forum:

https://www.mathworks.com/matlabcentral/answers/57444-faq-why-is-0-3-0-2-0-1-not-equal-to-zero

https://www.mathworks.com/matlabcentral/answers/316889-why-two-equal-numbers-are-not-equal

http://www.mathworks.com/matlabcentral/answers/102419-how-do-i-determine-if-the-error-in-my-answer-is-the-result-of-round-off-error-or-a-bug

http://www.mathworks.com/matlabcentral/answers/140656-about-numerical-difficulty-of-equality-comparison-operation-between-two-double-real-numbers-in-matlab

This excellent paper is highly reccomended:

https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

http://faculty.tarleton.edu/agapie/documents/cs_343_arch/papers/1991_Goldberg_FloatingPoint.pdf

[Tim Berk]

Thanks for pointing that out. I thought it was best to (even with my limited knowledge) give a solution for OP. Using round in this particular case actually solves their particular problem.

OP is comparing the output of lets say 0.3/0.1 to integers. In simplified form:

for i = 1:N
  if i==0.3/0.1

@Jose, your version would be (as I understand it)

if i==(0.3/0.1+eps(0.3/0.1))

@Guillaume, your version (as I understand it)

if abs(i-0.3/0.1)<1e-15

Would you use this when simply comparing (what are supposed to be) integers? Isn't round the best way to go to the closest integer? Even after all your explanations (for this particular problem) I would still use round as it you can actually take this out of the loop, i.e.

B = round(0.3/0.1);
for i = 1:N
  if i==B

@Stephen, what would you use? Apparently you know a lot about this but I didn't see a solution.

[José-Luis]

I would use:

 abs(3 - 0.3/0.1)  < eps(3)

which is rather stringent.

Using round() is a bad idea. Please refer to the links Stephen posted. However much you say that they are supposed to be integers, they are not. They are doubles.

You can probably cross a country road without looking left and right and get away with it most of the time. Until one day when that becomes a one way ticket to meet your maker.

Answer 2

[Rahul Pillai]

Hello everyone. Thank you for all your suggestions. I worked this solution out last night because I figured it was a precision problem that goes undetected to the human eye! I simply used the round() function while performing the division. It worked :)

[José-Luis]

It might have worked this time but it will come back and bite you eventually. The better approach is to use a tolerance.

That being said, please accept the answer that best solves your problem.