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
1 个评论
Steven Lord
2021-11-9
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.)
回答（1 个）
David Hill
2021-11-9
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
类别
在 帮助中心 和 File Exchange 中查找有关 Assumptions 的更多信息
2021-11-9
2021-11-9
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
