Loop problem

% I have got this simple program runnning a loop to satisfy the condition 'nvvni = 0' stop the program
PT= 50;
c = [0.20,0.10,0.10;0.20,0.20,0.20;0 0 0];
P= [190 72.2 51.6; 20.44 15.57 4.956;0 0 0];
K = P/PT;
A = c(1,1)*(K(1,1)-1)+ c(1,2)*(K(1,2)-1)+c(1,3)*(K(1,3)-1)+c(2,1)*(K(2,1)-1)+c(2,2)*(K(2,2)-1)+c(2,3)*(K(2,3)-1);
B = (c(1,1)*(K(1,1)-1)/K(1,1))+(c(1,2)*(K(1,2)-1)/K(1,2))+(c(1,3)*(K(1,3)-1)/K(1,3))+(c(2,1)*(K(2,1)-1)/K(2,1))+(c(2,2)*(K(2,2)-1)/K(2,2))+(c(2,3)*(K(2,3)-1)/K(2,3));
%disp(B)
for i=1:10
n=0;
nvni = (A )/(A-B);
n=n+1;
nvvni = (c(1,1)*(K(1,1)))/(nvni*(K(1,1)-1)+1)+(c(1,2)*(K(1,2)))/(nvni*(K(1,2)-1)+1)+(c(1,3)*(K(1,3)))/(nvni*(K(1,3)-1)+1)+(c(2,1)*(K(2,1)))/(nvni*(K(2,1)-1)+1)+(c(2,2)*(K(2,2)))/(nvni*(K(2,2)-1)+1)+(c(2,3)*(K(2,3)))/(nvni*(K(2,3)-1)+1);
if nvvni == 0;
disp(nvni);
break
else
nvvvi = (c(1,1)*(K(1,1)^2)/(nvni*(K(1,1) + 1)^2))+(c(1,2)*(K(1,2)^2)/(nvni*(K(1,2) + 1)^2))+(c(1,3)*(K(1,3)^2)/(nvni*(K(1,3) + 1)^2))+(c(2,1)*(K(2,1)^2)/(nvni*(K(2,1) + 1)^2))+(c(2,2)*(K(2,2)^2)/(nvni*(K(2,2) + 1)^2))+(c(2,3)*(K(2,3)^2)/(nvni*(K(2,3) + 1)^2));
nvvni = nvni - (nvvni/nvvvi);
nvni =nvvni;
end
end

1 Commento

Jan
Jan il 19 Feb 2012
Please format your code properly. I've done this for you this time.
What is your question?

Accedi per commentare.

Risposte (1)

Jan
Jan il 19 Feb 2012

0 voti

It is very unlikely that the result of such a complicated formula is exactly 0.0 due to rounding errors. Better check for a certain limit:
if abs(nvvni) < 1e-14
The size of the limit depends on the physical meaning of the formula and cannot be defined uniquely. An analysis of the sensitivity is recommended in addition: How sensitive is nvvni to tiny variations of the inputs, e.g. when using c+eps.

1 Commento

Walter Roberson
Walter Roberson il 19 Feb 2012
http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F

Accedi per commentare.

Categorie

Scopri di più su Loops and Conditional Statements in Centro assistenza e File Exchange

Richiesto:

il 19 Feb 2012

Modificato:

il 16 Ott 2013

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by