for loop does not save the values in the respective array

Question

Juliana Quintana 2022-2-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1653945-for-loop-does-not-save-the-values-in-the-respective-array

评论： Ersad Mert Mutlu 2022-2-20

Hello everyone, I'm looking for help about a for loop. I think it works and makes the loop, but it does not save the values in an array I made. It just show the last value taht seems to be the same value as the inizialazation I made. I would really appreciate your help.

clc 
clear all 
close all 
%________________
%parámetros 
% Pc y Tc tomado de  https://www.linde-gas.es/es/images/n-Butano_tcm316-612756.pdf
% w tomado de https://termoapunefm.files.wordpress.com/2011/10/41-1-_constantes_y_coeficientes_de_propiedades_fisicas.pdf
n = 100; 
T = 800; %K 
P = linspace(0.005,0.02,n) ; %atm 
Pc = 37.46; %atm 
Tc =425.16 ; % K
w = 0.1954 ; %
R = 0.08206 ; %atm
%_______________________________
%valores de while 
h=1e-8;
tol=1e-6;
iterMax=1000;
error = 10; 
iter = 0;  
%_______vector de respuesta_______
Vf = zeros(n,1);
for i = 1:n 
    P = linspace(0.005,0.02,n) ; %atm 
    V0 = 13129;
    
    while (error>tol && iter <iterMax)
        fun = patel(P(i),V0);
        derivada = (patel(P(i),(V0+h))-patel(P(i),V0))./h ;
        vol = V0 - fun./derivada;
        error = abs((vol-V0)./V0); 
        iter = iter + 1; 
        V0 = vol;
    end
    
end 
Vf(i)= vol; 
plot(P, Vf)
function resp = patel(P,V)
T = 800; %K 
Pc = 37.46; %atm 
Tc =425.16 ; % K
w = 0.1954 ; %
R = 0.08206 ; %atm
Tr = T/Tc; 
%__________Funciones del factor acentrico_______- 
F = 0.42413 + 1.30982*w - 0.295937*(w^2);
zetac = 0.329032 - 0.076799*w + 0.0211947*(w^2);
%_____________omega b ______________
poli = [1 ,(2-(3.*zetac)), 3.*(zetac.^2), (-zetac.^3)]; 
sol = roots(poli); 
sol = double(sol);
[omega_b, other] = min(sol(imag(sol)==0 & real(sol)>0));
%_____________omega a _________
omega_a =3*zetac^2+3*(1-2*zetac)*omega_b +omega_b^2+(1-3*zetac^3); 
%___________ omega c _________
omega_c = 1-3*zetac ; 
%_________constantes 
a = omega_a*(R^2*Tc^2/Pc)*(1+F*(1-sqrt(Tr)))^2; 
b = omega_b*(R*Tc/Pc);
c = omega_c*(R*Tc/Pc);
%__________residual_______
der = R.*T./(V-b);
izq = a./(V.*(V+b)+c.*(V-b));
resp = der - izq - P; 
end 

5 个评论
显示 3更早的评论隐藏 3更早的评论

Ersad Mert Mutlu 2022-2-20

编辑：Ersad Mert Mutlu 2022-2-20

在 MATLAB Online 中打开

Hi Juliana,

clc 
clear all 
close all 
%________________
%parámetros 
% Pc y Tc tomado de  https://www.linde-gas.es/es/images/n-Butano_tcm316-612756.pdf
% w tomado de https://termoapunefm.files.wordpress.com/2011/10/41-1-_constantes_y_coeficientes_de_propiedades_fisicas.pdf
n = 100; 
T = 800; %K 
P = linspace(0.005,0.02,n) ; %atm 
Pc = 37.46; %atm 
Tc =425.16 ; % K
w = 0.1954 ; %
R = 0.08206 ; %atm
%_______________________________
%valores de while 
h=1e-8;
tol=1e-6;
iterMax=1000;
error = 10; 
iter = 0;  
V0 = 13129;
%_______vector de respuesta_______
Vf = zeros(n,1);
for i = 1:n 
    P = linspace(0.005,0.02,n) ; %atm  
    derivada = (patel(P(i),(V0+h))-patel(P(i),V0))./h ;
    fun = patel(P(i),V0);
    vol = V0 - fun./derivada;
    while (error>tol && iter<iterMax)
        error = abs((vol-V0)./V0); 
        iter = iter + 1; 
        V0 = vol;
    end
    Vf(i)= vol; 
end 
plot(P, Vf)
function resp = patel(P,V)
T = 800; %K 
Pc = 37.46; %atm 
Tc =425.16 ; % K
w = 0.1954 ; %
R = 0.08206 ; %atm
Tr = T/Tc; 
%__________Funciones del factor acentrico_______- 
F = 0.42413 + 1.30982*w - 0.295937*(w^2);
zetac = 0.329032 - 0.076799*w + 0.0211947*(w^2);
%_____________omega b ______________
poli = [1 ,(2-(3.*zetac)), 3.*(zetac.^2), (-zetac.^3)]; 
sol = roots(poli); 
sol = double(sol);
[omega_b, other] = min(sol(imag(sol)==0 & real(sol)>0));
%_____________omega a _________
omega_a =3*zetac^2+3*(1-2*zetac)*omega_b +omega_b^2+(1-3*zetac^3); 
%___________ omega c _________
omega_c = 1-3*zetac ; 
%_________constantes 
a = omega_a*(R^2*Tc^2/Pc)*(1+F*(1-sqrt(Tr)))^2; 
b = omega_b*(R*Tc/Pc);
c = omega_c*(R*Tc/Pc);
%__________residual_______
der = R.*T./(V-b);
izq = a./(V.*(V+b)+c.*(V-b));
resp = der - izq - P; 
end 

Could you check this script please?

Juliana Quintana 2022-2-20

I just checked the script and it works, but it does not give the correct answer. I will check again , maybe there is something wrong with my equations. Thank you soo much for the help.

Ersad Mert Mutlu 2022-2-20

Most welcome 🤜

请先登录，再进行评论。

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

Answer 1

Torsten 2022-2-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1653945-for-loop-does-not-save-the-values-in-the-respective-array#answer_900390

编辑：Torsten 2022-2-20

在 MATLAB Online 中打开

%________________
%parámetros 
% Pc y Tc tomado de  https://www.linde-gas.es/es/images/n-Butano_tcm316-612756.pdf
% w tomado de https://termoapunefm.files.wordpress.com/2011/10/41-1-_constantes_y_coeficientes_de_propiedades_fisicas.pdf
n = 100; 
T = 800; %K 
P = linspace(0.005,0.02,n) ; %atm 
Pc = 37.46; %atm 
Tc =425.16 ; % K
w = 0.1954 ; %
R = 0.08206 ; %atm
%_______________________________
%valores de while 
h=1e-8;
tol=1e-6;
iterMax=1000;
error = 10; 
iter = 0;  
V0 = 13129;
%_______vector de respuesta_______
Vf = zeros(n,1);
volold = V0;
for i = 1:n 
    while  iter<iterMax 
        derivada = (patel(P(i),(volold+h))-patel(P(i),volold))./h ;
        fun = patel(P(i),volold);
        volnew = volold - fun/derivada;
        error = abs((volold-volnew)/volold);
        volold = volnew;
        if error>tol
          iter = iter + 1;
        else
          break
        end
    end
    Vf(i) = volnew;
end 
plot(P, Vf)
end
function resp = patel(P,V)
T = 800; %K 
Pc = 37.46; %atm 
Tc =425.16 ; % K
w = 0.1954 ; %
R = 0.08206 ; %atm
Tr = T/Tc; 
%__________Funciones del factor acentrico_______- 
F = 0.42413 + 1.30982*w - 0.295937*(w^2);
zetac = 0.329032 - 0.076799*w + 0.0211947*(w^2);
%_____________omega b ______________
poli = [1 ,(2-(3.*zetac)), 3.*(zetac.^2), (-zetac.^3)]; 
sol = roots(poli); 
sol = double(sol);
[omega_b, other] = min(sol(imag(sol)==0 & real(sol)>0));
%_____________omega a _________
omega_a =3*zetac^2+3*(1-2*zetac)*omega_b +omega_b^2+(1-3*zetac^3); 
%___________ omega c _________
omega_c = 1-3*zetac ; 
%_________constantes 
a = omega_a*(R^2*Tc^2/Pc)*(1+F*(1-sqrt(Tr)))^2; 
b = omega_b*(R*Tc/Pc);
c = omega_c*(R*Tc/Pc);
%__________residual_______
der = R.*T./(V-b);
izq = a./(V.*(V+b)+c.*(V-b));
resp = der - izq - P; 
end 

3 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten 2022-2-20

I made a small change to the program. It runs properly now without producing Inf values.

Juliana Quintana 2022-2-20

Perfect, I will check it out.

Thank you

请先登录，再进行评论。

for loop does not save the values in the respective array

5 个评论
显示 3更早的评论隐藏 3更早的评论

采纳的回答

3 个评论
显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

for loop does not save the values in the respective array

5 个评论 显示 3更早的评论隐藏 3更早的评论

采纳的回答

3 个评论 显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

5 个评论
显示 3更早的评论隐藏 3更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论