Unsure of what is happening in for loop

Question

EYKL 2021-12-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1605950-unsure-of-what-is-happening-in-for-loop

评论： Mathieu NOE 2021-12-13

采纳的回答： Mathieu NOE

在 MATLAB Online 中打开

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2021-12-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1605950-unsure-of-what-is-happening-in-for-loop#answer_850390

在 MATLAB Online 中打开

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

8 个评论
显示 6更早的评论隐藏 6更早的评论

Mathieu NOE 2021-12-8

编辑：Mathieu NOE 2021-12-8

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

Mathieu NOE 2021-12-9

在 MATLAB Online 中打开

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