Unsure of what is happening in for loop

Dear all,

I have the following for loop,

T92to20_state_d         = [2.50926000000000; 5.04629166666667; 3.81944833333333];
Rn92to20_state_d        = [24.1978271186779; 24.8685463692206; 24.5147611929163];
LvETPM92to20_state_d    = [11.5146886953329; 10.5486657087851; 10.7259878112342];
GE92to20_state_d        = [-12.3855735457507; -12.1467974925575; -12.2727450153218]; 
windspd92to20_state_d_2 = [1.64395657249513; 1.16191506903469; 1.62723837006298];
T4                      = [276.391481481482; 277.285000000000; 277.706296296296];
Pressure_atmos          = [99.189335562566010];

k                   = 0.35;                                               
z1                  = 0.5;   z2 = 1.14;             
z0                  = 0.21756875;                   
T                   = T92to20_state_d + 273.16;     
H_old               = Rn92to20_state_d - LvETPM92to20_state_d - GE92to20_state_d;   
H92to20_state_d     = H_old;
delta1              = ones(size(H_old));
L                   = ones(size(1:length(Rn92to20_state_d)+1));
zeta_old            = zeros(size(H_old));
l                   = 0;

for i = 1:3
     while delta1(i) > 0.01 
        u_s(i)  = (k.*windspd92to20_state_d_2(i)*0.277778*z1)/(log(z1/z0) - 4.7*zeta_old(i));   
        L(i)    = (T(i)*u_s(i)^3*1010*1.225)/(k*9.81*H92to20_state_d(i));                       
        a(i)    = 122.44*Pressure_atmos/T(i);                                                   
        b(i)    = u_s(i)*abs(T4(i) - T(i));                                                     
        c(i)    = 0.74*log(z2/z1) + (4.7/L(i))*(z2-z1);                                         
        H(i)    = a(i)*b(i)/c(i);                                                               
        zeta(i) = z2/L(i);
        delta1(i) = (L(i+1)-L(i))./L(i);
        H92to20_state_d(i)   = H(i);
        zeta_old(i) = zeta(i);
        l = l+1;
     end
end

Using manual calculation, I can get the values for i = 1 but not for i = 2 and 3.

% u_s(2) = (k.*windspd92to20_state_d_2(2)*0.277778*z1)/(log(z1/z0) - 4.7*zeta_old(2))
%        = 0.2444 (but answer from loop is 0.0679).
% delta1(1) = -0.6625 (but answer from loop is -0.7151)

I can't seem to find the problem. Thank you in advance.

Answer 1

[Mathieu NOE]

hello

I don't get why the while loop is for ?

The for loop by itself will generate this result

delta1 =

-0.7151

-0.1559

-0.6961

so the condition while delta1(i) > 0.01 is never met

I a have a feeling the code implementation is maybe not what is should be but I don't understand what is the target here

I beleive the code does even not need that for loop as it should be vectorized

also why have created the variable "l" which is incremented ?

so code sligthly modified for my test purpose

clc
clearvars

T92to20_state_d         = [2.50926000000000; 5.04629166666667; 3.81944833333333];
Rn92to20_state_d        = [24.1978271186779; 24.8685463692206; 24.5147611929163];
LvETPM92to20_state_d    = [11.5146886953329; 10.5486657087851; 10.7259878112342];
GE92to20_state_d        = [-12.3855735457507; -12.1467974925575; -12.2727450153218]; 
windspd92to20_state_d_2 = [1.64395657249513; 1.16191506903469; 1.62723837006298];
T4                      = [276.391481481482; 277.285000000000; 277.706296296296];
Pressure_atmos          = [99.189335562566010];
k                   = 0.35;                                               
z1                  = 0.5;   z2 = 1.14;             
z0                  = 0.21756875;                   
T                   = T92to20_state_d + 273.16;     
H_old               = Rn92to20_state_d - LvETPM92to20_state_d - GE92to20_state_d;   
H92to20_state_d     = H_old;
delta1              = ones(size(H_old));
L                   = ones(size(1:length(Rn92to20_state_d)+1));
zeta_old            = zeros(size(H_old));
l                   = 0;

% for i = 1:3
%      while delta1(ci) > 0.01 
for ci = 1:3
        u_s(ci)  = (k.*windspd92to20_state_d_2(ci)*0.277778*z1)/(log(z1/z0) - 4.7*zeta_old(ci));   
        L(ci)    = (T(ci)*u_s(ci)^3*1010*1.225)/(k*9.81*H92to20_state_d(ci));                       
        a(ci)    = 122.44*Pressure_atmos/T(ci);                                                   
        b(ci)    = u_s(ci)*abs(T4(ci) - T(ci));                                                     
        c(ci)    = 0.74*log(z2/z1) + (4.7/L(ci))*(z2-z1);                                         
        H(ci)    = a(ci)*b(ci)/c(ci);                                                               
        zeta(ci) = z2/L(ci);
        delta1(ci) = (L(ci+1)-L(ci))./L(ci);
        H92to20_state_d(ci)   = H(ci);
        zeta_old(ci) = zeta(ci);
        l = l+1; % what is the purpose of this ? 
end
% end

[EYKL]

Hi @Mathieu NOE,

The function of the code is to perform iteration for the variable L until L converges, such that the condition

% delta1(ci) = (L(ci+1)-L(ci))./L(ci) < 0.01

To summarise:

1. Calculate an initial H [H_old(1)] to start the iteration.
2. Calculate u_s(1), L(1), H(1) and zeta(1).
3. Repeat point 2 until L converges, i.e., delta1 < 0.01.

I am aware that perform iteration more than 3 steps are required but since I am getting weird values, I provided a short example here.

[Mathieu NOE]

hello again

what disturbs me is why are most of the variables (s, H, L ...) in this for loop not scalars but defined as vectors ? (the 2nd dimension has not to do with the iteration length as far as I understand)

this start already at the initilaisation phase

for example

T92to20_state_d = [2.50926000000000; 5.04629166666667; 3.81944833333333];

Rn92to20_state_d = [24.1978271186779; 24.8685463692206; 24.5147611929163];

LvETPM92to20_state_d = [11.5146886953329; 10.5486657087851; 10.7259878112342];

GE92to20_state_d = [-12.3855735457507; -12.1467974925575; -12.2727450153218];

windspd92to20_state_d_2 = [1.64395657249513; 1.16191506903469; 1.62723837006298];

T4 = [276.391481481482; 277.285000000000; 277.706296296296];

and

L = ones(size(1:length(Rn92to20_state_d)+1));

the for loop is only iteratin along those vector dimension but is not updating a full 3x1 sized vector

like zeta is defined as a 3 x 1 vector and the goal is update at each step one of the 3 components (one by one) of the vector, or the full vector at once ?

tx

[EYKL]

Hi @Mathieu NOE,

The goal is to update the variables one by one.

1. When i=1, update the first index of each variable.
2. When i=2, update the second index of each variable.
3. Then test if the delta1 < 0.01 condition is met. If not, continue iterating until it is.

Thank you!

[Mathieu NOE]

hello again

I still don't understand

If I follow what you say there would be only 3 iterations and this is not what you want ?

[EYKL]

Hi @Mathieu NOE,

Yes, the iterations should be > 3 but right now the problem is that after I manually calculate the values for i=2 using the previous i=1 values, I get different results from the for loop. So while the for loop purpose is for iterating using previous values, I don't think it is actually doing that.

Thanks

[Mathieu NOE]

hello again

ok I think I start to see what we want

try this and let me know

clc
clearvars

T92to20_state_d         = [2.50926000000000; 5.04629166666667; 3.81944833333333];
Rn92to20_state_d        = [24.1978271186779; 24.8685463692206; 24.5147611929163];
LvETPM92to20_state_d    = [11.5146886953329; 10.5486657087851; 10.7259878112342];
GE92to20_state_d        = [-12.3855735457507; -12.1467974925575; -12.2727450153218]; 
windspd92to20_state_d_2 = [1.64395657249513; 1.16191506903469; 1.62723837006298];
T4                      = [276.391481481482; 277.285000000000; 277.706296296296];
Pressure_atmos          = [99.189335562566010];
k                   = 0.35;                                               
z1                  = 0.5;   z2 = 1.14;             
z0                  = 0.21756875;                   
T                   = T92to20_state_d + 273.16;     
H_old               = Rn92to20_state_d - LvETPM92to20_state_d - GE92to20_state_d;   
H92to20_state_d     = H_old;
delta1              = ones(size(H_old));
L                   = ones(size(Rn92to20_state_d));
zeta_old            = 0;
l                   = 0;



% iteration 
for ci = 1:3
        u_s(ci)  = (k.*windspd92to20_state_d_2(ci)*0.277778*z1)/(log(z1/z0) - 4.7*zeta_old);   
        L(ci)    = (T(ci)*u_s(ci)^3*1010*1.225)/(k*9.81*H92to20_state_d(ci));                       
        a(ci)    = 122.44*Pressure_atmos/T(ci);                                                   
        b(ci)    = u_s(ci)*abs(T4(ci) - T(ci));                                                     
        c(ci)    = 0.74*log(z2/z1) + (4.7/L(ci))*(z2-z1);                                         
        H(ci)    = a(ci)*b(ci)/c(ci);                                                               
        zeta(ci) = z2/L(ci);
        if  ci == 1 
            delta1(ci) = delta1(1);
        else
            delta1(ci) = (L(ci)-L(ci-1))./L(ci-1);
        end
        H92to20_state_d(ci)   = H(ci);
        zeta_old = zeta(ci);
%         l = l+1; % not used 
end

[EYKL]

Hi @Mathieu NOE,

I have tried and it is want I want.

Thanks!

[Mathieu NOE]

My pleasure !