Matlab help regarding code

% Calculate the vapor Pressure of n-heptane using SRK-EOS
% P = RT / V - b  -  alpha*a / V(V-b)
R = 8.134e-3; % MPa m3 / mole.K
n = 0;
n = n+1;
for i=1:5
Pn = 1.6714;
Tr = 0.9;
Tc = 540.13; % K
Pc = 2.74; % Mpa;
T = Tr*Tc;
w = 0.350;
m = 0.480+1.574*w -0.176*w^2;
al=(1+m*(1-Tr^0.5))^2;
a = 0.42748*(((R*Tc)^2) /Pc);
b = 0.08664*(R*Tc / Pc);
% Z3 - Z2 + (A-B-B^2)Z-AB=0
A = (a*Pn*al)/(R*T)^2;
B = (b*Pn)/(R*T);
Zv=max(roots([1 -1  A-B-B^2 -A*B]));
Zl=min(roots([1 -1  A-B-B^2 -A*B]));
disp(Zv)
disp(Zl)
phiv(i)=exp((Zv -1)-log(Zv -B)-(A/B) *(log(Zv+B/Zv )));
phil(i)=exp((Zl -1)-log(Zl -B)-(A/B) *(log(Zl+B/Zl )));
Fv=Pn*phiv;
Fl=Pn*phil;
disp(phiv)
disp(phil);
disp(Fv);
disp(Fl);
    A=Fl-Fv;
if (A < 10^-4);
    Pnew=Pn;
       else 
       Pn+1 =Pn * (phil/phiv);       
end
end

Basically in this code I am trying to assume a pressure 'Pn' to calculate the results which satisfy this relation (A<10^-4) if this satisfy's then our assume pressure is right else if not then we use the relation Pn+1 = Pn*(Phil/Phiv) to caluclate another pressure and update the above pressure (previously assume) with this one. I have error with this and need advice please help.

regards

Answer 1

[Walter Roberson]

The syntax

Pn+1 =Pn * (phil/phiv);

is not valid. The left side of an assignment must be a valid variable name or array reference. For example,

Pn_plus_1 =Pn * (phil/phiv);

Or

P(n+1) =Pn * (phil/phiv);

In your "while" loop, you do not recalculate "A", so if the loop is entered at all, it is going to be repeated indefinitely. When you assign a variable value, the assignment does not affect any values that were already calculated before that point, including not any assignments that used that same variable name in the formula.

foo = 5;
bar = foo;
foo = 7;
bar

At the end of this, foo will be 7 and bar will still be what it was at the time the assignment for it was executed, namely 5. The update to the variable "foo" that was used in initializing "bar" does not cause an update to any variable already assigned to.

[Nasir Qazi]

Note:- The correct codes :- please have a look

% Calculate the vapor Pressure of n-heptane using SRK-EOS

% P = RT / V - b - alpha*a / V(V-b)

R = 8.134e-3; % MPa m3 / mole.K

n = 0;

n = n+1;

for i=1:100

Pn = 1.6714;

Tr = 0.9;

Tc = 540.13; % K

Pc = 2.74; % Mpa;

T = Tr*Tc;

w = 0.350;

m = 0.480+1.574*w -0.176*w^2;

al=(1+m*(1-Tr^0.5))^2;

a = 0.42748*(((R*Tc)^2) /Pc);

b = 0.08664*(R*Tc / Pc);

% Z3 - Z2 + (A-B-B^2)Z-AB=0

A = (a*Pn*al)/(R*T)^2;

B = (b*Pn)/(R*T);

Zv=max(roots([1 -1 A-B-B^2 -A*B]));

Zl=min(roots([1 -1 A-B-B^2 -A*B]));

disp(Zv)

disp(Zl)

phiv(i)=exp((Zv -1)-log(Zv -B)-(A/B) *(log(Zv+B/Zv )));

phil(i)=exp((Zl -1)-log(Zl -B)-(A/B) *(log(Zl+B/Zl )));

Fv=Pn*phiv;

Fl=Pn*phil;

disp(phiv)

disp(phil);

disp(Fv);

disp(Fl);

error=Fl-Fv;

if (error < 10^-4);

Pnew=Pn;

else

P(n+1) =Pn * (phil/phiv);

Pn = P(n+1) % this is suppose to update above Pn , but it doesn't

end

end

--------------------------------

This loop should run 100 time as a whole but when it satisfy the Condition Error < 10-4 , it stops there, how this is possible. This is not a homework question and I am trying to learn kindly Answer by pointing the code where i did mistake ..

thanks

[Image Analyst]

I don't understand. I copied, pasted, and ran it. I added "i" in the loop to print out the value of the loop index. It ran 100 times and never satisfied the condition "error < 10^-4. So I can't reproduce what you said happened.

[Nasir Qazi]

Assume pressure Pn will lead to a error which must be less then 10^-4 , if not it will Assumed another pressure from P(n+1)=Pn(previously assumes) x ( phil (previously calculate) / (phiv(previously calculated) to give us a pressure which will update the first assume pressure . now you suggest where is the wrong with the loop .

[Image Analyst]

I still don't understand. Maybe you can give a Pn where it does get into that part of the if statement and then bails out of the loop. But I see no reason why it should exit the loop and stop the program. There is no break statement anywhere so why should it exit the loop or "stop there" whatever that means. You mean like stop with an error because an exception got thrown? And how/where is Pnew ever used, even if it did get assigned?

[Nasir Qazi]

% these are the correct codes, have a look again

% Calculate the vapor Pressure of n-heptane using SRK-EOS

% P = RT / V - b - alpha*a / V(V-b)

R = 8.134e-3; % MPa m3 / mole.K

n = 0;

n = n+1;

for i=1:5

Pn = 1.6714; % Assume Pressure

Tr = 0.9;

Tc = 540.13; % K

Pc = 2.74; % Mpa;

T = Tr*Tc;

w = 0.350;

m = 0.480+1.574*w -0.176*w^2;

al=(1+m*(1-Tr^0.5))^2;

a = 0.42748*(((R*Tc)^2) /Pc);

b = 0.08664*(R*Tc / Pc);

% Z3 - Z2 + (A-B-B^2)Z-AB=0

A = (a*Pn*al)/(R*T)^2;

B = (b*Pn)/(R*T);

Zv=max(roots([1 -1 A-B-B^2 -A*B]));

Zl=min(roots([1 -1 A-B-B^2 -A*B]));

disp(Zv)

disp(Zl)

phiv(i)=exp((Zv -1)-log(Zv -B)-(A/B) *(log(Zv+B/Zv )));

phil(i)=exp((Zl -1)-log(Zl -B)-(A/B) *(log(Zl+B/Zl )));

Fv=Pn*phiv;

Fl=Pn*phil;

disp(phiv)

disp(phil);

disp(Fv);

disp(Fl);

error=(Fl-Fv);

if (error < 10^-4); % this statement has to be verify if not

Pn; % |

else % |

Pn+1 =Pn * (phil/phiv); % V

P(n+1) = Pn ( this pressure should update the previous assume pressure )

end

end