why why why why why ..........​..........​..........​..........​..

2 visualizaciones (últimos 30 días)
Nasir Qazi
Nasir Qazi el 27 de Feb. de 2012
I have an code and when I try to run it, the status say busy , and it wont show me the answer to my solution , what exactly is the problem ?
[EDITDED, Jan Simon]: Copied from the comment section:
This program does not converge:
% Calculate Bubble Point and Dew Point SRK-EOS
% P = RT / V - b - alpha*a / V(V-b)
%----------------------------------------------
% liquid phase composition
y = [0.86 0.05 0.05 0.02 0.01 0.0095 0.0005]';
% vapor phase composition
x = [0.45 0.05 0.05 0.03 0.01 0.01 0.4]';
% critical pressure
Pc =[4.5947 4.8711 4.2472 3.6397 3.3688 3.1233 1.965]';
% critical temperature
Tc =[190.7389 305.5111 370.0333 425.3444 469.8889 512.7778 644.4444]';
% Accentric Factor
w = xlsread('myfile.xlsx','Sheet1','E2:E8');
%---------------------------------------------
L = 1;
R = 8.314e-3; % MPa m3 / mole.K
T = 200; % K
Tr = T./Tc; % K
P(L) = 5.0; %MPa
m = 0.480+1.574.*w -0.176.*w.^2;
al=1+ m.*(1-sqrt(Tr)).^2;
a = 0.42748.*(((R.*Tc).^2)./Pc);
b = 0.08664.*(R.*Tc./Pc);
q=0;
while q==0
% Z3 - Z2 + (A-B-B^2)Z-AB=0
% for gas phase
SA = 0;
RR = zeros(7,1);
for i = 1:7
for j = 1:7
SA = SA + y(i).*y(j).*sqrt(a(i).*a(j).*al(i).*al(j));
RR(i) = RR(i) + y(j).*sqrt(a(i).*a(j).*al(i).*al(j));
end
end
SB = sum(y.*b);
% For Liquid Phase
SAA = 0;
MM = zeros(7,1);
for i = 1:7
for j = 1:7
SAA = SAA + x(i).*x(j)*sqrt(a(i).*a(j)*al(i).*al(j));
MM(i) = MM(i) +x(j).*sqrt(a(i).*a(j).*al(i).*al(j));
end
end
SBB = sum(x.*b);
%for Gas Phase
A = SA.*P(L)./(R*T).^2;
B = SB.*P(L)./(R*T);
ZV=roots([1 -1 A-B-B^2 -A*B]);
Zv=max(ZV);
Zvc=isreal(Zv);
if Zvc == false
Zv=abs(Zv);
disp('error 1')
end
%for liquid Phase
AA = (SAA.*P(L))./(R*T).^2;
BB = (SBB.*P(L)./(R*T));
Z=roots([1 -1 AA-BB-BB^2 -AA*BB]);
Zl=min(Z);
Zvl=isreal(Zl);
if Zvl == false
Zl=abs(Zl);
disp('error 2')
end
phiv =exp((b.*(Zv-1)/SB) - log(Zv-B)-(A/B).*((2.*RR/SA)-
(b/SB)).*log(1 + B/Zv));
phil =exp((b.*(Zl-1)./SBB) - log(Zl-BB)-(AA./BB).*((2.*MM/SAA)-
(b./SBB)).*log(1 + BB./Zl));
K = phil./phiv;
RA1= sum(y./K)-1;
if abs(RA1) < 1e-5 && abs((y./K) - x) < 1e-5
break
else
if L==1
L = L+1;
P(L) = 1.01.*P(L-1);
RA1;
else
if L==1
P(L)=P(L-1)-(RA1./(RA1-RA2))*(P(L-1)-P(L-2));
RA2 = RA1;
end
%adjustment
x2 = (1-0.5)*(y./K) + 0.5*x;
x = x2;
end
end
end
disp(P(L));
  2 comentarios
Jiro Doke
Jiro Doke el 27 de Feb. de 2012
Can you please rephrase the subject line so that it is specific to your question? Please read here:
http://www.mathworks.com/matlabcentral/about/answers/
and read the "Tips" under "Getting Started" -> "Ask a Question"
Walter Roberson
Walter Roberson el 27 de Feb. de 2012
Duplicate is at http://www.mathworks.com/matlabcentral/answers/30397-trial-and-error-help-is-not-giving-me-solution-to-that-still-waiting

Iniciar sesión para comentar.

Respuestas (3)

Grzegorz Knor
Grzegorz Knor el 27 de Feb. de 2012
It means, that program is still running. Run it in debug mode, and you will see step by step how it works.

Jan
Jan el 27 de Feb. de 2012
You can insert a waitbar or some output in the code to be informed about the status of process. Either the computations take a lot of time or you have created an infinite loop by accident.
  4 comentarios
Nasir Qazi
Nasir Qazi el 27 de Feb. de 2012
% Calculate Bubble Point and Dew Point SRK-EOS
% P = RT / V - b - alpha*a / V(V-b)
%----------------------------------------------
% liquid phase composition
y = [0.86 0.05 0.05 0.02 0.01 0.0095 0.0005]';
% vapor phase composition
x = [0.45 0.05 0.05 0.03 0.01 0.01 0.4]';
% critical pressure
Pc =[4.5947 4.8711 4.2472 3.6397 3.3688 3.1233 1.965]';
% critical temperature
Tc =[190.7389 305.5111 370.0333 425.3444 469.8889 512.7778 644.4444]';
% Accentric Factor
w = xlsread('myfile.xlsx','Sheet1','E2:E8');
%---------------------------------------------
L = 1;
R = 8.314e-3; % MPa m3 / mole.K
T = 200; % K
Tr = T./Tc; % K
P(L) = 5.0; %MPa
m = 0.480+1.574.*w -0.176.*w.^2;
al=1+ m.*(1-sqrt(Tr)).^2;
a = 0.42748.*(((R.*Tc).^2)./Pc);
b = 0.08664.*(R.*Tc./Pc);
q=0;
while q==0
% Z3 - Z2 + (A-B-B^2)Z-AB=0
% for gas phase
SA = 0;
RR = zeros(7,1);
for i = 1:7
for j = 1:7
SA = SA + y(i).*y(j).*sqrt(a(i).*a(j).*al(i).*al(j));
RR(i) = RR(i) + y(j).*sqrt(a(i).*a(j).*al(i).*al(j));
end
end
SB = sum(y.*b);
% For Liquid Phase
SAA = 0;
MM = zeros(7,1);
for i = 1:7
for j = 1:7
SAA = SAA + x(i).*x(j)*sqrt(a(i).*a(j)*al(i).*al(j));
MM(i) = MM(i) +x(j).*sqrt(a(i).*a(j).*al(i).*al(j));
end
end
SBB = sum(x.*b);
%for Gas Phase
A = SA.*P(L)./(R*T).^2;
B = SB.*P(L)./(R*T);
ZV=roots([1 -1 A-B-B^2 -A*B]);
Zv=max(ZV);
Zvc=isreal(Zv);
if Zvc == false
Zv=abs(Zv);
disp('error 1')
end
%for liquid Phase
AA = (SAA.*P(L))./(R*T).^2;
BB = (SBB.*P(L)./(R*T));
Z=roots([1 -1 AA-BB-BB^2 -AA*BB]);
Zl=min(Z);
Zvl=isreal(Zl);
if Zvl == false
Zl=abs(Zl);
disp('error 2')
end
phiv =exp((b.*(Zv-1)/SB) - log(Zv-B)-(A/B).*((2.*RR/SA)-
(b/SB)).*log(1 + B/Zv));
phil =exp((b.*(Zl-1)./SBB) - log(Zl-BB)-(AA./BB).*((2.*MM/SAA)-
(b./SBB)).*log(1 + BB./Zl));
K = phil./phiv;
RA1= sum(y./K)-1;
if abs(RA1) < 1e-5 && abs((y./K) - x) < 1e-5
break
else
if L==1
L = L+1;
P(L) = 1.01.*P(L-1);
RA1;
else
if L==1
P(L)=P(L-1)-(RA1./(RA1-RA2))*(P(L-1)-P(L-2));
RA2 = RA1;
end
%adjustment
x2 = (1-0.5)*(y./K) + 0.5*x;
x = x2;
end
end
end
disp(P(L));
% I am trying to using trial and error to get the results of P(L) , but my solution is not converging and keep going on. I need to know what can I do to resolve the issue .
Help Please
Jan
Jan el 27 de Feb. de 2012
Well, this is an important detail: As far as I understand your problem is the lack of convergence. This is another story.

Iniciar sesión para comentar.


Sean de Wolski
Sean de Wolski el 27 de Feb. de 2012
The bald and not excessively bald and not excessively smart hamster obeyed a terrified and not excessively terrified hamster.
  1 comentario
Jan
Jan el 27 de Feb. de 2012
@Nasir: If you are wondering, what the deeper meaning of this answer could be, type "why" in the command line - repeatedly.

Iniciar sesión para comentar.

Categorías

Más información sobre Startup and Shutdown en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by