Error with vpa solve

Hello everyone.

The program computes for the first Tc but then gives me an error. Can you help me?

Tb=112; %K
Z=0.27;
R=0.0821;
DH=8200; %J/mol


for Tc=120:0.01:200
%clapeyron  Pc
Pc=exp((DH/8.314)*((1/Tb)-(1/Tc)));
syms a b c Vc;
X=R*Tc/(Vc-b)-(a/Tc)/((Vc+c)^2)-Pc;
eqn1=R*Tc/(Vc-b)-(a/Tc)/((Vc+c)^2)-Pc==0;
eqn2=diff(X,Vc)==0;
eqn3=diff(X,Vc,2)==0;
eqn4=(Pc*Vc)/(R*Tc)-Z==0;
sol=vpasolve([eqn1,eqn2,eqn3,eqn4],[a,b,c,Vc]);
Vc=sol.Vc;
a=sol.a 
b=sol.b 
c=sol.c 

syms V
P_star=1;
Pb=R*Tb/(V-b)-(a/Tb)/((V+c)^2);
Z=R*Tb/(V-b)-(a/Tb)/((V+c)^2)-P_star==0;
V_soluzione=real(vpasolve(Z,V));
Vliq=min(V_soluzione);
Vvs=max(V_soluzione);
Y=Pb-P_star;
integralearee=real(vpaintegral(Y,Vliq,Vvs))
 if integralearee<=0
   break
 end
end
vol_liq_sat=Vliq
vol_vap_sat=Vvs

[Steven Lord]

Please show us the the full and exact text of the error and/or warning messages (all the text displayed in red and/or orange in the Command Window.)

Answer 1

[David Hill]

Tc=fzero(@I,126);  
[~,Vliq,Vvs]=I(Tc);

function [integralearee,Vliq,Vvs]=I(Tc)
  Tb=112; %K
  Z=0.27;
  R=0.0821;
  DH=8200; %J/mol
  syms a b c Vc V;
  Pc=exp((DH/8.314)*((1/Tb)-(1/Tc)));
  X=R*Tc/(Vc-b)-(a/Tc)/((Vc+c)^2)-Pc;
  eqn1=R*Tc/(Vc-b)-(a/Tc)/((Vc+c)^2)-Pc==0;
  eqn2=diff(X,Vc)==0;
  eqn3=diff(X,Vc,2)==0;
  eqn4=(Pc*Vc)/(R*Tc)-Z==0;
  sol=vpasolve([eqn1,eqn2,eqn3,eqn4],[a,b,c,Vc]);
  a=sol.a; 
  b=sol.b;
  c=sol.c; 
  P_star=1;
  Pb=R*Tb/(V-b)-(a/Tb)/((V+c)^2);
  Z=R*Tb/(V-b)-(a/Tb)/((V+c)^2)-P_star==0;
  V_soluzione=real(vpasolve(Z,V));
  Vliq=min(V_soluzione);
  Vvs=max(V_soluzione);
  Y=Pb-P_star;
  integralearee=real(vpaintegral(Y,Vliq,Vvs));
end

[Alessio Falcone]

Hi thank you very much for the reply.

How do I set that the program must stop when 'integralearee' is equal to 0 or very close to it? furthermore 'Vliq' must be approximately equal to 0.036. Then I have to show the Tc found on the display, thank you very much