For Loop Not Working??

Hi,

I am designing a programme which will operate using a for loop starting at a value of 4 and ending at 10 as shown in the code below:

for TSR=4:10;

I am unsure as to whether I should put the following line of code in the end of my program:

TSR=TSR+i;

I have also secified that i is a counting index of 1. If I dont put it in, the loop seems to keep interating for a TSR value of 4. But if I do put it in the for loop surpasses the final value of 10 and increases continuously. Could anybody please help with this problem?

Thanks, Kevin

Answer 1

[dpb]

Your code snippet is incomplete; you have no i defined as shown what more as "a counting index of 1" .

The for index variable in your case is TSR, there's no i involved and, "Yes, Virginia, you do not artificially index the for loop variable inside a for construct."

doc for % and all will be explained...read all the doc including the _Tips_ section at the bottom.

Then try the following exercise at the command line --

for i=4:10
  disp(i)
  i=i+1;
end

and explain the results.

NB: while often used as "traditional" loop indices, in Matlab i and j are predefined as the imaginary sqrt(-1). You can override them as above and it is often done even in Matlab, but be aware of what you are doing when do so particularly if there's any chance whatsoever that there are real (so to speak :) ) complex values lurking...

Answer 2

[dpb]

...I put an end if statement at the end which seems to have solved it...

Solved what???

for TSR=4:10 % TSR from 4 to 10
   while abs(Check)>=tolerance
     ...
     if TSR==10
     else
       TSR=TSR+i;
     end
   end
end

The solution is as documented to NOT try to update the loop variable inside the loop at all.

Remove the code snippet

     if TSR==10
     else
       TSR=TSR+i;
     end

entirely; it isn't needed and doesn't do anything, anyway. As the demo I suggested you run shows, the value of the loop index (TSR here in your case; i in mine) is overwritten anyway by the next loop iteration. Since after you updated TSR it isn't ever referenced until the beginning of the next loop iteration, that local copy that is changed is never referenced. To see this even more clearly, modify the sample I gave before to--

for i=4:10,disp(i),i=i+2;disp(i),end

and observe the results.

As the above conclusively demonstrates, whatever the actual problem you had before, it was NOT owing to the cause you've surmised; there's insufficient data provided to tell what, precisely, you did incorrectly before as you didn't post an example code that reproduced the problem to see...

I'll make some additional comments on the code in making it more "Mablaby"...

Rewrite

r_local=[(R./N)*1; (R./N)*2; ... (R./N)*8; (R./N)*9];

as

r_local=bsxfun(R./N),[1:9]); r_local=r_local(:);

Use inbuilt constant pi instead of crude approximation for pi...multiple places.

C=((8*pi.*r_local)./(B.*Cl_design)).*(1-cos(Phi));

for TSR=4:10
  while abs(Check)>=tolerance
    axial_induction_old=axial_induction;        
    TSR_local=TSR*(r_local./R);
    Phi=(2/3)*atan(1./TSR_local);
    axial_induction=1./(((4.*(sin(Phi).^2))./(sigma.*Cl_design.*cos(Phi)))+1);
    angular_induction=(1-(3*axial_induction))./((4.*axial_induction)-1);
    relative_wind=(1-axial_induction)./((1+axial_induction).*TSR);
    F=2/pi*acos(exp(-((B/2.*(1-r_over_R))./(r_over_R.*sin(relative_wind)))));
    C_T=(sigma.*((1-axial_induction).^2).*((Cl.*cos(relative_wind))+ ...
        (Cd.*sin(relative_wind))))./((sin(relative_wind)).^2);
    if C_T<0.96
      axial_induction=1./(1+(4.*F.*(sin(relative_wind).^2))./ ...
                            (sigma.*Cl.*cos(relative_wind)));  
    else
      axial_induction=1./(((4.*F.*cos(relative_wind))./(sigma.*Cl))-1); 
    end
    D=(8/(TSR.*N)).*(F.*(sin(Phi).^2).*(cos(Phi)- ...
        ((TSR_local).*(sin(Phi)))).*(sin(Phi)+ ...
        ((TSR_local).*(cos(Phi)))).*(1-(Cd./Cl).*atan(Phi)).*(TSR_local.^2));
    Cp=sum(D);
    Diff=axial_induction-axial_induction_old;
    Check=max(Diff(:));
    store_Phi(:,i)=Phi;
    store_TSR_local(:,i)=TSR_local;
    store_axial_induction(:,i)=axial_induction;
    store_angular_induction(:,i)=angular_induction;
    store_relative_wind(:,i)=relative_wind;
    store_Check(:,i)=Check;
    store_Diff(:,i)=Diff;
    store_Cp(:,i)=Cp;
    store_TSR(:,i)=TSR;
  end
end

Answer 3

[Ben11]

That is very strange; the default increment for for-loops is 1 so you should not have to add the TSR = TSR+i. By the way is it really +i or +1 ? Maybe you have something in your code which makes TSR smaller than 10 at every iteration so the loop keeps going. Can you post your whole code so we can better help?

[dpb]

... Maybe you have something in your code which makes TSR smaller than 10 at every iteration so the loop keeps going.

That's not how Matlab loops work...the demo I posted for Kevin illustrates this just as well with a slight modification...try

for i=4:10,disp(i),i=i-2;disp(i),end

at the command line.

Again I commend

doc for

to your attention--read the "Tips" section carefully...

Clearly if his code had an infinite loop he had "something" in his code that caused it but whatever it was, it was not the for loop index itself. As noted earlier, it's not possible to tell what his error was that caused the symptom; he didn't post any code that demonstrated the claimed behavior.

[Ben11]

Yep you're so right I need more coffee I guess thanks for correcting my non sense!

Answer 4

[Kevin]

Thanks for the help. I put an end if statement at the end which seems to have solved it. This is my entire code:

clc
clear

filename = 'book2.xlsx';
Input = xlsread(filename)
N = xlsread(filename,'C7:C15')
r = xlsread(filename,'D7:D15')
C = xlsread(filename,'E7:E15')
S = xlsread(filename,'F7:F15')
Non = xlsread(filename,'G7:G15')

% Inputs
R=0.4;                                      % Radius of Rotor
B=3;                                        % Number of blades
V=2;                                        % Fluid velocity
Rho=998;                                    % Fluid Density
N=9;                                        % Number of Blade Elements
Cp_estimate=0.5;                            % Estimate power coefficient
Alpha_design=4;                             % Design alpha
Cl_design=1.04;                             % Design lift coefficient

% Variables
TSR=4;
Cp=0;
i=1;                                        % Counter
alpha_new=0;
axial_induction=0;
tolerance=0.1;                              % Tolerance Value
Check=1;
axial_induction_old=0;                      % Initial value for old axial induction factor

r_local=[(R./N)*1; (R./N)*2; (R./N)*3; (R./N)*4; (R./N)*5; (R./N)*6; (R./N)*7; (R./N)*8; (R./N)*9];
r_over_R=r_local./R;
TSR_local=r_over_R.*TSR;

Phi=(2/3)*atan(1./TSR_local);
C=((8.*3.14.*r_local)./(B.*Cl_design)).*(1-cos(Phi));
sigma=(B*C)./(3.14.*r_local.*2);

Cl=[1.3; 1.1; 1; 0.9; 0.86; 0.83; 0.8; 0.75; 0.5];                          % Lift Coefficients
Cd=[0.027; 0.024; 0.02; 0.019; 0.018; 0.016; 0.013; 0.012; 0.01];           % Drag Coefficients

F=(2/3.14).*acos(exp(-(((B/2).*(1-(r_over_R)))./((r_over_R).*sin(Phi)))));

for TSR=4:10;                                                           % TSR from 4 to 10

     while abs(Check)>=tolerance

axial_induction_old=axial_induction;        
TSR_local=TSR*(r_local./R);                                                     % Local Tip Speed Ratio
Phi=(2/3)*atan(1./TSR_local);                                                   % Angle of Relative Fluid

axial_induction=1./(((4.*(sin(Phi).^2))./(sigma.*Cl_design.*cos(Phi)))+1);
angular_induction=(1-(3*axial_induction))./((4.*axial_induction)-1);

relative_wind=(1-axial_induction)./((1+axial_induction).*TSR);

F=(2/3.14).*acos(exp(-(((B/2).*(1-(r_over_R)))./((r_over_R).*sin(relative_wind)))));          % Tip Loss Factor

C_T=(sigma.*((1-axial_induction).^2).*((Cl.*cos(relative_wind))+(Cd.*sin(relative_wind))))./((sin(relative_wind)).^2);

if C_T<0.96
    axial_induction=1./(1+(4.*F.*(sin(relative_wind).^2))./(sigma.*Cl.*cos(relative_wind)));  
else
    axial_induction=1./(((4.*F.*cos(relative_wind))./(sigma.*Cl))-1); 
end

D=(8./(TSR.*N)).*(F.*(sin(Phi).^2).*(cos(Phi)-((TSR_local).*(sin(Phi)))).*(sin(Phi)+((TSR_local).*(cos(Phi)))).*(1-(Cd./Cl).*atan(Phi)).*(TSR_local.^2));
Cp=sum(D);

Diff=axial_induction-axial_induction_old;
Check=max(Diff(:));

store_Phi(:,i)=Phi;
store_TSR_local(:,i)=TSR_local;
store_axial_induction(:,i)=axial_induction;
store_angular_induction(:,i)=angular_induction;
store_relative_wind(:,i)=relative_wind;
store_Check(:,i)=Check;
store_Diff(:,i)=Diff;
store_Cp(:,i)=Cp;
store_TSR(:,i)=TSR;

if TSR==10
else TSR=TSR+i;
end

     end

end

figure(1)
plot(store_Cp,store_TSR)
hold all
title('Cp vs Tip Speed Ratio')
xlabel('TSR')
ylabel('Cp')

Answer 5

[dpb]

I seriously recast your code below...I may have ended up w/ a mismatched parens or two as I tried to remove many of the superfluous pairs of them and the unneeded dot operators to make it simpler to try to parse the actual code.

The biggest change, however, was to preset a vector of TSR values over which to iterate and then replaced the loop index to select from that vector. Then, the results are stored in a (preallocated) set of arrays by column from the first to last instead of off in the yonder somewhere as before.

Then, it appears you reversed the role of x- and y- axes in your plot command so I swapped them around.

filename = 'book2.xlsx';
data = xlsread(filename);   % only read the file once...
N = data(7:15,3);           % assumes numeric array is at origin;
r = data(7:15,4);           % adjust row/col offset to the actual
C = data(7:15,5);           % locations based on what text is also in 
S = data(7:15,6);           % the spreadsheet
Non = data(7:15,7);
%
% Inputs
R=0.4;                       % Radius of Rotor
B=3;                         % Number of blades
V=2;                         % Fluid velocity
Rho=998;                     % Fluid Density
N=9;                         % Number of Blade Elements
Cp_estimate=0.5;             % Estimate power coefficient
Alpha_design=4;              % Design alpha
Cl_design=1.04;              % Design lift coefficient
%
% Variables
TSR=4;                       % Initial tip speed ratio
Cp=0;                        % Initial power coefficient
i=1;                         % Counter
alpha_new=0;                 % Initial value for alpha new
axial_induction=0;           % Initial axial induction factor
tolerance=0.01;              % Tolerance Value
Check=1;                     % Initial check value
axial_induction_old=0;       % Initial value for old axial induction factor
r_local=bsxfun(R./N),[1:9]); r_local=r_local(:);
r_over_R=r_local./R;
TSR_local=r_over_R.*TSR;
Phi=(2/3)*atan(1./TSR_local);
C=(8*pi*r_local./(B*Cl_design)).*(1-cos(Phi));
sigma=B*C./(2*pi*r_local);
Cl=[1.3; 1.1; 1; 0.9; 0.86; 0.83; 0.8; 0.75; 0.5];  % Lift Coefficients
Cd=[0.027; 0.024; 0.02; 0.019; 0.018; 0.016; 0.013; 0.012; 0.01];
F=2/pi*acos(exp(-(B/2*(1-r_over_R)./(r_over_R.*sin(Phi)))));
axial_induction=1/(4*sin(Phi).^2./(sigma.*Cl_design.*cos(Phi))+1);
angular_induction=(1-3*axial_induction)./(4*axial_induction-1);
TSRvec=4:10;             % set the desired TSR values
[nr nc]=[length(r_local) length(TSRvec)];  % result sizes for preallocation
store_Phi=zeros(nr, nc);
store_TSR_local=zeros(nr, nc);
store_axial_induction=zeros(nr, nc);
store_angular_induction=zeros(nr, nc);
store_relative_wind=zeros(nr, nc);
store_Check=zeros(nr, nc);
store_Diff=zeros(nr, nc);
store_Cp=zeros(nr, nc);
store_TSR=zeros(nr, nc);
for i=1:length(TSRvec)   % and iterate over them
TSR=TSRvec(i);         % set local temporary from vector
while abs(Check)>=tolerance
  axial_induction_old=axial_induction;        
  TSR_local=TSR*(r_local./R);
  Phi=(2/3)*atan(1./TSR_local);
  relative_wind=(1-axial_induction)./((1+axial_induction)*TSR(i));
  F=2/pi*acos(exp(-(B/2*(1-r_over_R)./(r_over_R.*sin(relative_wind)))));
  C_T=(sigma.*((1-axial_induction).^2).*(Cl.*cos(relative_wind)+ ...
        Cd.*sin(relative_wind)))./((sin(relative_wind)).^2);
  if C_T(TSR,:)<0.96
    axial_induction=1./(1+(4.*F.*(sin(relative_wind).^2))./ ...
                               (sigma.*Cl.*cos(relative_wind)));  
  else
    axial_induction=1./(((4.*F.*cos(relative_wind))./(sigma.*Cl))-1); 
  end
  D=(8/(TSR*N))*(F.*sin(Phi).^2.*(cos(Phi)-(TSR_local.*sin(Phi))).* ...
    (sin(Phi)+(TSR_local.*cos(Phi))).*(1-Cd./Cl.*atan(Phi)).*TSR_local.^2);
  Cp=sum(D);
  Diff=axial_induction-axial_induction_old;
  Check=max(Diff(:));
  store_Phi(:,i)=Phi;
  store_TSR_local(:,i)=TSR_local;
  store_axial_induction(:,i)=axial_induction;
  store_angular_induction(:,i)=angular_induction;
  store_relative_wind(:,i)=relative_wind;
  store_Check(:,i)=Check;
  store_Diff(:,i)=Diff;
  store_Cp(:,i)=Cp;
  store_TSR(:,i)=TSR;
end
end
%
figure
plot(store_TSR,store_Cp)
title('Cp vs Tip Speed Ratio')
xlabel('TSR')
ylabel('Cp')

There is one serious problem in your code that undoubtedly is not doing what you intend--since I don't know what the intent is for certain I can only point it out for your attention--

if C_T(TSR,:)<0.96
  axial_induction=1./(1+(4.*F.*(sin(relative_wind).^2))./ ...
                               (sigma.*Cl.*cos(relative_wind)));  
else
  axial_induction=1./(((4.*F.*cos(relative_wind))./(sigma.*Cl))-1); 
end

In Matlab, the test above will be TRUE iff all values in the vector C_T are true so the above branch will only trigger when that occurs. If it is need for the case to be taken for any one value, then you would need

    if any(C_T(TSR,:))<0.96

As noted, you'll need to fix this logic up as you intend.

[Amin Mohammed]

Well explained, I'm facing similar issue, but still not able to work it out. I'll try to take the advantage of this answer. In answer provided by @dpb, in line 56:

relative_wind=(1-axial_induction)./((1+axial_induction)*TSR(i));

Should read:

relative_wind=(1-axial_induction)./((1+axial_induction)*TSR);

Answer 6

[Kevin]

Hi dpb,

Thank you so much for the help. I'm still relatively new to the world of Matlab and your advice really makes my life a bit easier.

Best regards, Kevin

[dpb]

This part of it would be easier yet if you would take and apply some of the advice given... :)

Answer 7

[Maximilian]

For the for loop below it works when I do not include the (e^(-n/r) part and I replace

it with (e*(-n/r). How can I make it work with the raised to the power part.

e = exp(1);
N= 1001;
Ic= zeros(1,N);
Vs = 1;
for n = 1:N
    Vc(n) = 1 - (e^(-n/r));
    Ic(n) = (Vs - Vc(n))/R;
end

[Walter Roberson]

exp(-n/r)

[Maximilian]

Okay I fixed that but it only returns the first four values and then the rest are zero.

N= 1001;
Ic= zeros(1,N);
Vs = 1;
for n = 1:N
    Ic(n) = (Vs - 1 - (exp(-n/r)))/R;
end

[dpb]

What is r?

exp(-z) rapidly --> 0 as z--> inf and n is 1:10001 so the expression

Vs-1-exp(-z) --> 1-(1-exp(-z)) --> 1-1+exp(-z) = exp(-z)

will become zero to machine precision.

NOTA BENE: The first version was (Vs-Vc)/R where Vc=1-exp(n/R). That expands to +exp(z) when the inner parentheses are removed...

And, depending upon how you're looking at the results, the answer may not be identically zero, but small enough to be rounded to zero.

For example if

 n=5;
 r=1;
exp(-n/r)

ans = 0.0067

r=r/10;
exp(-n/r)

ans = 1.9287e-22

the latter of which is close enough to zero to probably count... :)

[Walter Roberson]

Note that Vs is 1, and Vs - 1 would be 0, so your terms would be exp(-n/r)/R . When -n/r reaches -744.4 then the exp() would be too small to represent. And if you are using the default display representation of format short then that would try to express everything in a common power of 10 and would run out of precision very quickly. If you were to try to plot the exp(-n/r)/R then it only takes a few terms before the rest cannot be distinguished from 0 because pixels have a finite size.

[dpb]

Btw, @Walter Roberson got the -744 number from

log(eps(0))

ans = -744.4401

anything > will underflow explicitly.

I was going to note that "the MATLAB way" as I pointed out to the original poster of the Q to which this subthread is attached would be to use vectorized expressions rather than for...end loop

N=1000;
R=;
r=???;
VS=1;
n=0:N;
Ic=(Vs-(1-exp(-n/r)))/R;

[Walter Roberson]

I tend to forget about eps(0) so what I actually used was log(eps(realmin)) which is identical to log(eps(0))

[dpb]

Chuckle...