for loop does not save the values in the respective array

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 

[KSSV]

Two problems in your code.

1. The line Vf(i)= vol; should be inside the for loop.
2. The while loop is not working. You need to make few changes.

[Juliana Quintana]

Hey, thanks for answer. Could you help me with the changes I have to make ? Or do you know where i have to made them ?

[Ersad Mert Mutlu]

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]

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]

Most welcome 🤜

Answer 1

[Torsten]

%________________
%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 

[Juliana Quintana]

It works as expected now, thank you so much for the help.

[Torsten]

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

[Juliana Quintana]

Perfect, I will check it out.

Thank you