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

2 Ansichten (letzte 30 Tage)
Nasir Qazi
Nasir Qazi am 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 Kommentare
Jiro Doke
Jiro Doke am 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 am 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

Melden Sie sich an, um zu kommentieren.

Antworten (3)

Grzegorz Knor
Grzegorz Knor am 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 am 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 Kommentare
Nasir Qazi
Nasir Qazi am 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 am 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.

Melden Sie sich an, um zu kommentieren.


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

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Startup and Shutdown finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by