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

10 visualizzazioni (ultimi 30 giorni)
Nasir Qazi
Nasir Qazi il 27 Feb 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 Commenti
Jiro Doke
Jiro Doke il 27 Feb 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 il 27 Feb 2012
Duplicate is at http://www.mathworks.com/matlabcentral/answers/30397-trial-and-error-help-is-not-giving-me-solution-to-that-still-waiting

Accedi per commentare.

Risposte (3)

Grzegorz Knor
Grzegorz Knor il 27 Feb 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 il 27 Feb 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 Commenti
Nasir Qazi
Nasir Qazi il 27 Feb 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 il 27 Feb 2012
Well, this is an important detail: As far as I understand your problem is the lack of convergence. This is another story.

Accedi per commentare.


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

Accedi per commentare.

Categorie

Scopri di più su Programming in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by