Changing the value of only one parameter causes a problem (possibly with while loop)

Hello everyone. Firstly, as I do not have access to my other computer right now, I am using the online compiler, if that makes any difference (but I do not think so). In the attaced code, in the first section, when I set the d(1:10) value to 0.04, it just can't proceed further after the while loop and keeps calculating forever, as seen in screenshot (1). At least I assume so; the thing is, I haven't used matlab much and when I open the debugger (screenshot(2)) and step in and continue, it keeps doing the while loop over and over again. Overall, I do not get any crucial error that seemingly affect the working mechanism of the code: Only "consider preallocating for speed suggestions", so I cannot understand what is wrong with it. By the way, this happens as the value gets smaller: When d is smaller than 0.07. When I set it to something like 150, it still works. What might be the reasoon?

Thanks in advance!

%% Constants and Assumptions
g=9.81;
hgelgit=0.5;
he=0.3;
Vmin=0.6;
Vmax=0.9;
f=0.03;
s=0.25;
Cq(1)=0.61;
rho0=1008.301;
rhoa=998.83;
H=20;
Lana=1500;
d(1:10)=0.04;
hf(1)=0;
%% Initial Hydraulic Head
h=(rhoa*(H+hgelgit+he))/(rho0);
E(1)=h-H;
%% First Port
a(1)=pi*(d(1))^2./4;
q(1)=Cq(1)*a(1)*sqrt(2*g*E(1));
Q(1)=q(1);
D(1)=floor(20*(sqrt(4*Q(1)/(pi*Vmin))))/20;
A(1)=pi*D(1)^2./4;
V(1)=Q(1)/(pi*(D(1)^2./4)); %[m/s]
hizyuk(1)=0;
U(1)=q(1)/(pi*d(1)^2./4);
Dmax(1)=(sqrt(4*Q(1)/(pi*Vmin)));
Dmin(1)=(sqrt(4*Q(1)/(pi*Vmax)));

%% Other Ports
for i=2:10
    hf(i)=f*(s/D(i-1))*(V(i-1)^2./(2*g));
    E(i)=E(i-1)+hf(i);
    Cq(i)=0.63-0.58*(((V(i-1))^2)/(2*g*E(i)));
    a(i)=pi*(d(i))^2./4; 
    x=cumsum(a);
    q(i)=Cq(i)*a(i)*sqrt(2*g*E(i));
    Q(i)=Q(i-1)+q(i); 
    Dmax(i)=(sqrt(4*Q(i)/(pi*Vmin)));
    Dmin(i)=(sqrt(4*Q(i)/(pi*Vmax)));
    D(i)=D(i-1); 
    A(i)=pi*D(i)^2./4;
    V(i)=Q(i)/(pi*(D(i)^2./4)); 
    hizyuk(i)=V(i-1)^2./(2*g*E(i));
    while (V(i) >= Vmax)    
    D(i)=floor(20*(sqrt(4*Q(i)/(pi*Vmin))))/20; 
    V(i)=Q(i)/(pi*(D(i)^2./4)); 
    end
  U(i)=q(i)/(pi*d(i)^2./4);
end
y=x./(pi*D.^2./4);
%% Results
results=zeros(10,16); 
results(:,1)=transpose(hf);
results(:,2)=transpose(E); 
results(:,3)=transpose(hizyuk);
results(:,4)=transpose(Cq);
results(:,5)=transpose(a);
results(:,6)=transpose(x);
results(:,7)=transpose(q); 
results(:,8)=transpose(U); 
results(:,9)=transpose(Q); 
results(:,10)=transpose(Dmin);
results(:,11)=transpose(Dmax);
results(:,12)=transpose(D); 
results(:,13)=transpose(A);
results(:,14)=transpose(y);
results(:,15)=transpose(V); 
results(:,16)=transpose(Q)*1000; 

max(q)/min(q)
results

[Cris LaPierre]

The issue I see is that nothing inside the while loop changes loop to loop. The calculation of D(i) will always be the same, as will the resulting V(i). So if the recalculated value of V(i) is still >= Vmax, you will have an infinite loop.

Typically, you would update an index each loop, or one of the input values so the equation output changes each loop.

[Mahmut Cenk]

@Cris LaPierre, thanks a lot for your response. But can you elaborate a little, I feel like I am a little slow right now. Because for larger values of d, D(i) and V(i) values change within the loop. I mean, the code runs perfectly for values like 0.07 or greater. And I can confirm that the code correctly calculates the values by taking a simple calculation as reference and calculating using a calculator. Does it have anything to do with the "d" value being smaller than a certain threshold? That when the value is too small, the output values approach to zero or something?

[Cris LaPierre]

The value of d is used to calculate a, which is used to calculate q, which is turned into Q, which is used to compute D.

It would appear, then, when d is small, V(i) is less the Vmax, and when it is large, V(i) is greater than or equal to Vmax. There is no difference is how the values of D or V are calculated inside or outside your loop, so I suspect all that is happening is your loop runs infinitely, or not at all.

[Mahmut Cenk]

Cheers mate.

Answer 1

[Steven Lord]

Let's look at the while loop.

    while (V(i) >= Vmax)    
    D(i)=floor(20*(sqrt(4*Q(i)/(pi*Vmin))))/20; 
    V(i)=Q(i)/(pi*(D(i)^2./4)); 
    end

In order for the truth of the condition to change, either V(i) and/or Vmax must change. Nothing inside the loop assigns a value to Vmax so that means V(i) must change. Since we assign a value to V(i) in the loop that means this can't be an infinite loop, right?

Well, in order for V(i) to change one of Q(i) or D(i) must change. [pi is not going to change :)] You assign a value to D(i) in the loop but don't assign to Q(i) so if V(i) changes it must be because D(i) did.

D(i) depends on Q(i), pi, and Vmin. None of those change in the loop. So you're always going to get the same value of D(i). That means that the next line will always give you the same value for V(i). That means the while condition's truth cannot change unless the value of D(i) and/or the value of V(i) changes at the first iteration through the loop.

The code in the loop matches the code that was used to define V(i) two lines before the loop, so it's not going to change in that first iteration. That means If you enter the while loop, unless the value of D(i) changes at that first iteration, you're never going to exit the loop.

From looking at your code I suspect that the value of D(i) doesn't change either, but you have enough variables involved in the calculations that it would take a little while to prove that it doesn't and I'll leave that for you to check when you determine how to modify your code.

[Mahmut Cenk]

Cheers mate.

Answer 2

[Alan Stevens]

You can see from the following that when i = 3, 4 and 7, V(i) stays at a constant value greater than Vmax

%% Constants and Assumptions
g=9.81;
hgelgit=0.5;
he=0.3;
Vmin=0.6;
Vmax=0.9;
f=0.03;
s=0.25;
Cq(1)=0.61;
rho0=1008.301;
rhoa=998.83;
H=20;
Lana=1500;
d=0.04;
hf(1)=0;
%% Initial Hydraulic Head
h=(rhoa*(H+hgelgit+he))/(rho0);
E(1)=h-H;
%% First Port
a(1)=pi*(d)^2./4;
q(1)=Cq(1)*a(1)*sqrt(2*g*E(1));
Q(1)=q(1);
D(1)=floor(20*(sqrt(4*Q(1)/(pi*Vmin))))/20;
A(1)=pi*D(1)^2./4;
V(1)=Q(1)/(pi*(D(1)^2./4)); %[m/s]
hizyuk(1)=0;
U(1)=q(1)/(pi*d^2./4);
Dmax(1)=(sqrt(4*Q(1)/(pi*Vmin)));
Dmin(1)=(sqrt(4*Q(1)/(pi*Vmax)));

%% Other Ports
for i=2:10
    hf(i)=f*(s/D(i-1))*(V(i-1)^2./(2*g));
    E(i)=E(i-1)+hf(i);
    Cq(i)=0.63-0.58*(((V(i-1))^2)/(2*g*E(i)));
    a(i)=pi*(d)^2./4; 
    x=cumsum(a);
    q(i)=Cq(i)*a(i)*sqrt(2*g*E(i));
    Q(i)=Q(i-1)+q(i); 
    Dmax(i)=(sqrt(4*Q(i)/(pi*Vmin)));
    Dmin(i)=(sqrt(4*Q(i)/(pi*Vmax)));
    D(i)=D(i-1); 
    A(i)=pi*D(i)^2./4;
    V(i)=Q(i)/(pi*(D(i)^2./4)); 
    hizyuk(i)=V(i-1)^2./(2*g*E(i));
    steps = 0;
    while (V(i) >= Vmax) & (steps<4) 
        steps= steps+1;
        D(i)=floor(20*(sqrt(4*Q(i)/(pi*Vmin))))/20; 
        disp(['i = ', int2str(i),'  step = ', int2str(steps)])
        disp(['V(i)in= ',num2str(V(i))])
        V(i)=Q(i)/(pi*(D(i)^2./4)); 
        disp(['V(i)out= ',num2str(V(i))])
    end
    
    
  U(i)=q(i)/(pi*d^2./4);
end

i = 2  step = 1

V(i)in= 2.5565

V(i)out= 0.63911

i = 3  step = 1

V(i)in= 0.97981

V(i)out= 0.97981

i = 3  step = 2

V(i)in= 0.97981

V(i)out= 0.97981

i = 3  step = 3

V(i)in= 0.97981

V(i)out= 0.97981

i = 3  step = 4

V(i)in= 0.97981

V(i)out= 0.97981

i = 4  step = 1

V(i)in= 1.307

V(i)out= 1.307

i = 4  step = 2

V(i)in= 1.307

V(i)out= 1.307

i = 4  step = 3

V(i)in= 1.307

V(i)out= 1.307

i = 4  step = 4

V(i)in= 1.307

V(i)out= 1.307

i = 5  step = 1

V(i)in= 1.6163

V(i)out= 0.71836

i = 7  step = 1

V(i)in= 1.0191

V(i)out= 1.0191

i = 7  step = 2

V(i)in= 1.0191

V(i)out= 1.0191

i = 7  step = 3

V(i)in= 1.0191

V(i)out= 1.0191

i = 7  step = 4

V(i)in= 1.0191

V(i)out= 1.0191

i = 8  step = 1

V(i)in= 1.1652

V(i)out= 0.65545

[Mahmut Cenk]

This also clarified things, cheers mate.