MATLAB Answers

MATLAB code never stops running (infinite loop)

31 views (last 30 days)
Marcos Da Conceicao
Marcos Da Conceicao on 28 Sep 2020 at 20:40
Answered: Walter Roberson on 28 Sep 2020 at 22:17
My Matlab code stuck in infinite loop, error comes from line 157 to 242 (fourth while loop). I really don't know what is wrong. I would appreciate any help on this!
clear all
clear variables
global nos noc nof CompMW CompSigma FeedTemp FeedPres Ftotal ...
length ir functiontemp
errormax=1e-6;
maxiterations=500;
nos=100;
noc=2;
z(1)=0.21;
z(2)=0.79;
Q(2)=1.04e-6;
selectivity=5.85;
Q(1)=Q(2)*selectivity;
FeedPres=150*76/14.7; %psia to cmHg
PermPres=14.7*76/14.7;
FeedTemp=25+273.15;
CompMW(1)=32; %O2 (g/mol)
CompMW(2)=28.02; %N2
CompMW(3)=44.01; %CO2
CompSigma(1)=3.433; %Angstroms
CompSigma(2)=3.667;
CompSigma(3)=3.996;
Ptotal=zeros(1,nos); %Pre-allocations for faster code execution
Rtotal=zeros(1,nos);
R=zeros(nos,noc);
P=zeros(nos,noc);
Sweep=zeros(1,noc);
R_temp=zeros(nos,noc);
P_temp=zeros(nos,noc);
nof=1000; %number of fibers
length=50; %cm
or=0.04/2; %cm
now=6;
fmin=nof*2*pi*or*length*Q(1)*FeedPres;
for n=1:now
Ftotal(n,1)=fmin;
end
PresDropSet=0.04; %psia
PresDropMax=0.2;
x(1)=-2.3506;
x(2)=-1.33584;
x(3)=-0.43607;
x(4)=0.43607;
x(5)=1.33584;
x(6)=2.35060;
w(1)=0.00453;
w(2)=0.15706;
w(3)=0.724629;
w(4)=0.724629;
w(5)=0.157067;
w(6)=0.00453;
id_avg=0.021;
id_std=0.068; % percent
pcounter=1;
while PresDropSet<=PresDropMax
error_f=0.1;
f_iterations=0;
while error_f>1e-4 && f_iterations<100
error_f=0;
for n=1:now
id=id_avg*(1.4142*id_std*x(n)+1);
ir=id/2; %cm
TotalArea=2*pi*or*length*nof;
A=TotalArea/nos;
for m=1:2
for i=1:noc
F(i)=z(i)*Ftotal(n,m);
end
for j=1:nos %Do crossflow calculations
Ptotal(j)=0;
Rtotal(j)=0;
for i=1:noc %Initialize with permeate vacuum
if j==1
P(j,i)=Q(i)*A*FeedPres*F(i)/Ftotal(n,m);
R(j,i)=F(i)-P(j,i);
else
P(j,i)=Q(i)*A*FeedPres*R(j-1,i)/Rtotal(j-1);
R(j,i)=R(j-1,i)-P(j,i);
end
Ptotal(j)=Ptotal(j)+P(j,i);
Rtotal(j)=Rtotal(j)+R(j,i);
end
c_error=0.1;
iterations=0;
while c_error>errormax && iterations<maxiterations %Iterate crossflow solution
c_error=0;
for i=1:noc
if j==1
P_new=(Q(i)*A*FeedPres*F(i)/Rtotal(j))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*FeedPres/Ftotal(n,m));
R_new=F(i)-P_new;
else
P_new=(Q(i)*A*FeedPres*R(j-1,i)/Rtotal(j))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*FeedPres/Rtotal(j));
R_new=R(j-1,i)-P_new;
end
c_error=c_error+(abs(R_new-R(j,i))+abs(P_new-P(j,i)))/Ftotal(n,m);
iterations=iterations+1;
P(j,i)=P_new;
R(j,i)=R_new;
end
Ptotal(j)=0;
Rtotal(j)=0;
for i=1:noc
Ptotal(j)=Ptotal(j)+P(j,i);
Rtotal(j)=Rtotal(j)+R(j,i);
end
end
end
Sweeptotal=0; %Sum permeates from last stage to first
for i=1:noc
Sweep(i)=0;
end
for j=nos:-1:1
if j==nos
Ptotal(j)=Ptotal(j)+Sweeptotal;
else
Ptotal(j)=Ptotal(j)+Ptotal(j+1);
end
for i=1:noc
if j==nos
P(j,i)=P(j,i)+Sweep(i);
else
P(j,i)=P(j,i)+P(j+1,i);
end
end
end
%Permeate flow
if m ==1
for j = 1:nos
for n = 1:now
for i=1:noc
Permf(n,i)=P(j,i);
end
end
end
end
for j=1:nos
RetPres(j)=FeedPres;
end
d_maxiterations=500;
c_error=0.1; %Use direct-substitution method
d_iterations=0;
while c_error>errormax && d_iterations<=d_maxiterations
c_error=0;
functiontemp=Ftotal(n,m);
RetPres=RetPresDrop(R,Rtotal,RetPres);
for j=1:nos
Ptotal(j) = 0;
Rtotal(j) = 0;
P_sum(j) = 0;
for i=1:noc
for n = 1:now
P(j,i) = P_sum(j) + 1/sqrt(pi)*(w(n)*Permf(n,i)); %Summing up the flows from each fiber type(n)
end
Ptotal(j) = Ptotal(j) + P(j,i);
Rtotal(j) = Rtotal(j) + R(j,i);
end
for i=1:noc
y(j,i) = P(j,i)/Ptotal(j);
end
for i=1:noc
if j==1
x(j,i) = F(i)/Rtotal(j);
%y(j,i) = P(j,i)/Ptotal(j);
J(j,i) = Q(i)*A*(RetPres(j)*x(j,i)-PermPres*y(j,i));
R_new= F(i) - J(j,i);
P_new = P(j+1,i) + J(j,i);
% P_new=(Q(i)*A*RetPres(j)*F(i)/Rtotal(j)+P(j+1,i)*(1+Q(i)*A*RetPres(j)/Rtotal(j)))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*RetPres(j)/Rtotal(j));
% R_new=F(i)-(P_new-P(j+1,i));
elseif j==nos
x(j,i) = R(j,i)/Rtotal(j);
% y(j,i) = Sweep(i)/Ptotal(j);
J(j,i) = Q(i)*A*(RetPres(j)*x(j,i)-PermPres*y(j,i));
R_new= R(j-1,i) - J(j,i);
P_new = Sweep(i) + J(j,i);
% P_new=(Q(i)*A*RetPres(j)*R(j-1,i)/Rtotal(j)+Sweep(i)*(1+Q(i)*A*RetPres(j)/Rtotal(j)))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*RetPres(j)/Rtotal(j));
% R_new=R(j-1,i)-(P_new-Sweep(i));
else
x(j,i) = R(j,i)/Rtotal(j);
%y(j,i) = P(j,i)/Ptotal(j);
J(j,i) = Q(i)*A*(RetPres(j)*x(j,i)-PermPres*y(j,i));
R_new= R(j-1,i) - J(j,i);
P_new = P(j+1,i) + J(j,i);
% P_new=(Q(i)*A*RetPres(j)*R(j-1,i)/Rtotal(j)+P(j+1,i)*(1+Q(i)*A*RetPres(j)/Rtotal(j)))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*RetPres(j)/Rtotal(j));
% R_new=R(j-1,i)-(P_new-P(j+1,i));
end
c_error=c_error+(abs(R_new-R(j,i))+abs(P_new-P(j,i)))/Ftotal(n,m);
P(j,i)=P_new;
R(j,i)=R_new;
end
Ptotal(j)=0;
Rtotal(j)=0;
for i=1:noc
Ptotal(j)=Ptotal(j)+P(j,i);
Rtotal(j)=Rtotal(j)+R(j,i);
end
if m==1
if j==nos
for i=1:noc
Retf(n,i)=R(nos,i);
end
Rtotaltemp=Rtotal(nos);
end
% if j == 1
% for i=1:noc
% Permf(n,i)=P(j,i);
% end
% Ptotaltemp=Ptotal(1);
% end
if m==1
dimrec(n,j)=0;
% dimperm(n,j) = 0;
for i=1:noc
dimrec(n,j)=dimrec(n,j)+R(j,i);
% dimperm(n,j) = dimperm(n,j) + P(j,i);
dimdist(j)=j/nos;
end
dimxfrac(n,j)=R(j,1)/dimrec(n,j);
% dimyfrac(n,j)=P(j,1)/dimperm(n,j);
dimrec(n,j)=dimrec(n,j)/Ftotal(n,1);
% theta(n,j) = dimperm(n,j)/Ftotal(n,1);
end
end
end
d_iterations=d_iterations+1;
end
PresDrop(m)=(FeedPres-RetPres(nos))*14.7/76;
Ftotal(n,2)=1.1*Ftotal(n,1);
end
if n==3
xfrac_nv(pcounter)=Retf(3,1)/Rtotaltemp;
recovery_nv(pcounter)=Rtotaltemp/Ftotal(3,1);
end
alpha=0.1*Ftotal(n,1)/(PresDrop(2)-PresDrop(1));
error_f=error_f+abs(alpha*(PresDropSet-PresDrop(1))/Ftotal(n,1));
Ftemp(n)=Ftotal(n,1);
Ftotal(n,1)=Ftotal(n,1)+alpha*(PresDropSet-PresDrop(1));
end
f_iterations=f_iterations+1;
end
RetTotal=0;
% FeedTotal=0;
%PermTotal= 0;
for i=1:noc
Ret(i)=0;
%Feed(i)=0;
Perm(i)=0;
%Permf_nos(i)=0;
for n=1:now
Ret(i)=Ret(i)+1/sqrt(pi)*(w(n)*Retf(n,i));
%Perm(i)=Perm(i)+1/sqrt(pi)*(w(n)*Permf(n,i)); %% total permeate at stage 1
%Feed(i)=Feed(i)+1/sqrt(pi)*(w(n)*Ftemp(n));
end
RetTotal=RetTotal+Ret(i);
%FeedTotal=FeedTotal+Feed(i);
% PermTotal=PermTotal + Perm(i);
end
FeedTotal=0;
for n=1:now
FeedTotal=FeedTotal+1/sqrt(pi)*(w(n)*Ftemp(n));
end
% MBCheck = - PermTotal + FeedTotal - RetTotal;
xfrac_v(pcounter)=Ret(1)/RetTotal;
recovery_v(pcounter)=RetTotal/FeedTotal;
pcounter=pcounter+1;
PresDropSet=PresDropSet*1.1;
end
function RetPres = RetPresDrop(R,Rtotal,RetPres)
global nos noc nof CompMW CompSigma FeedTemp FeedPres ...
length ir functiontemp
for j=1:nos
RetMu=0;
for i=1:noc
CompMu=((FeedTemp*CompMW(i))^0.5)/(CompSigma(i)^2);
RetMu=RetMu+CompMu*R(j,i)/Rtotal(j);
end
RetMu=RetMu*2.6693e-6;
if j==1
deltaP=(8*RetMu*functiontemp*(76/FeedPres)*(FeedTemp/273.15)* ...
(length/nos))/(pi*(ir^4)*nof);
deltaP=deltaP*7.5e-4; %Pa to cmHg
RetPres(j)=FeedPres-deltaP;
else
deltaP=(8*RetMu*Rtotal(j-1)*(76/RetPres(j-1))*(FeedTemp/273.15)* ...
(length/nos))/(pi*(ir^4)*nof);
deltaP=deltaP*7.5e-4; %Pa to cmHg
RetPres(j)=RetPres(j-1)-deltaP;
end
end
end

  0 Comments

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 28 Sep 2020 at 22:17
Around line 62 or so you have
for n=1:now
Within that for loop, around line 140 or so you have
if m ==1
for j = 1:nos
for n = 1:now
for i=1:noc
Permf(n,i)=P(j,i);
end
end
end
end
Notice that you have the for n inside the for n . MATLAB does define what happens when you do this: the change to n will be in effect until the next outer for n starts, at which point n will be set to the "next" value of the outer loop as-if the change to n had not been done. But meanwhile until the bottom of the loop, n has changed. In particular at your line just inside the inner while loop you are having problems with,
functiontemp=Ftotal(n,m);
the n that is going to be used there is the n that was left-over from the inner for n loop, not the n associated with the outer for n loop.
Notice as well that the assignment to Permf is inside a triple for loop but is assigning to a variable using only two indices. The value that Permf(n,i) is going to be left with is going to be the last value assigned to it. You could achieve the same effect with
Permf(1:now, 1:noc) = P(nos, 1:noc)
and since you are not initializing Permf (growing it inside the loop) and it is always the same size, this in turn could be expressed as
Permf = repmat(P(nos, :), now, 1);
with no loop.
Are you sure that is what you want Permf to be assigned?
Any time you have a for loop that is not doing reduction (e.g., totaling the P values) and has more nesting levels than the number of indices of the output variable, then you should be optimizing to do only what would be left in the final result after the overly-nested assignments were complete.

  0 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!