Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Infinite Loop Problem

Subject: Infinite Loop Problem

From: Young Hong

Date: 21 May, 2010 01:46:03

Message: 1 of 3

Hi,

I'm trying to plot a 3D graph of Pr,z,Tr. My program goes into infinite loops and doesn't plot a graph either. I am wondering if something gone wrong or made an error some where.

Main function file:
%main program, Cubic Equation of State
clear all; clc
global A B R T P
R= 0.08206; Tc= 512.6; Pc= 80.9; w=0.559;
 
 
%loop for pressure
countA=0;
for j=[1:2:200];
    countA=countA+1;
    countB=0;
   for i = [300:5:800];
    countB = countB+1;
    press(countA,countB)= j;
    P=press(countA, countB);
    Temp(countA,countB) = i;
    T = Temp(countA, countB);
    
    m = 0.37464+ 1.54226*w-0.26992*(w^2);
    alpha = (1+m*(1-(T/Tc)^0.5))^2;
    a = 0.45724*(((R*Tc)^2)/Pc);
    b = 0.07780*R*Tc/Pc;
    A= a*alpha*P/((R^2)*(T^2));
    B= b*P/R*T;
    
    Pr=P./Pc;
    Tr=T./Tc;
    
    zseed = P/(T*R);
  z(countA, countB) = fzero('algebraic',zseed);
  
   end
end

surf(Tr,z,Pr)


%function file
function y=algebraic(z)
global A B ;
y= z^3-(1-B)*z^2+(A-2*B-3*B^2)*z+(-A*B+B^2+B^3);
end

Thank you.

Subject: Infinite Loop Problem

From: Roger Stafford

Date: 21 May, 2010 06:54:04

Message: 2 of 3

"Young Hong" <hong_young11@yahoo.com> wrote in message <ht4okr$gvp$1@fred.mathworks.com>...
> Hi,
>
> I'm trying to plot a 3D graph of Pr,z,Tr. My program goes into infinite loops and doesn't plot a graph either. I am wondering if something gone wrong or made an error some where.
>
> Main function file:
> %main program, Cubic Equation of State
> clear all; clc
> global A B R T P
> R= 0.08206; Tc= 512.6; Pc= 80.9; w=0.559;
>
>
> %loop for pressure
> countA=0;
> for j=[1:2:200];
> countA=countA+1;
> countB=0;
> for i = [300:5:800];
> countB = countB+1;
> press(countA,countB)= j;
> P=press(countA, countB);
> Temp(countA,countB) = i;
> T = Temp(countA, countB);
>
> m = 0.37464+ 1.54226*w-0.26992*(w^2);
> alpha = (1+m*(1-(T/Tc)^0.5))^2;
> a = 0.45724*(((R*Tc)^2)/Pc);
> b = 0.07780*R*Tc/Pc;
> A= a*alpha*P/((R^2)*(T^2));
> B= b*P/R*T;
>
> Pr=P./Pc;
> Tr=T./Tc;
>
> zseed = P/(T*R);
> z(countA, countB) = fzero('algebraic',zseed);
>
> end
> end
>
> surf(Tr,z,Pr)
>
>
> %function file
> function y=algebraic(z)
> global A B ;
> y= z^3-(1-B)*z^2+(A-2*B-3*B^2)*z+(-A*B+B^2+B^3);
> end
>
> Thank you.

  There was no infinite loop on my computer for this program - it took only a couple of minutes. I can tell you why you got no surf plot, though. You didn't save Tr and Pr in arrays to match the mesh array of z. The single values that were present in Tr and Pr at the time you called 'surf' were the last ones calculated in your for-loops. Naturally 'surf' doesn't like that.

  The cubic you define in 'algebraic' has three roots, apparently all real, which depend on A and B, and your call to 'fzero' finds only one of them each time, based on the estimate given by 'zseed'. The surface plot using it looked very smooth and consistent without any discontinuities. Do you have some way of knowing that this will produce the root you are seeking?

Roger Stafford

Subject: Infinite Loop Problem

From: Darren Rowland

Date: 21 May, 2010 07:01:07

Message: 3 of 3

I encountered a couple of errors when I ran the program. Here is a version which runs for me in a finite time (maybe 20 seconds on my machine)

--------------------------------------
function eqstate

R= 0.08206; Tc= 512.6; Pc= 80.9; w=0.559;
 
%loop for pressure
countA=0;
for j = 1:2:200
    countA=countA+1;
    countB=0;
   for i = 300:5:800
    countB = countB+1;
    press(countA,countB)= j;
    P=press(countA, countB);
    Temp(countA,countB) = i;
    T = Temp(countA, countB);
    
    m = 0.37464+ 1.54226*w-0.26992*(w^2);
    alpha = (1+m*(1-(T/Tc)^0.5))^2;
    a = 0.45724*(((R*Tc)^2)/Pc);
    b = 0.07780*R*Tc/Pc;
    A= a*alpha*P/((R^2)*(T^2));
    B= b*P/R*T;
    
    Pr=P./Pc;
    Tr=T./Tc;
    
    zseed = P/(T*R);
    f = @(x) (x^3-(1-B)*x^2+(A-2*B-3*B^2)*x+(-A*B+B^2+B^3));
  z(countA, countB) = fzero(f,zseed);
  
   end
end

surf(Temp,z,press)

--------------------------------------

Notice the use of an anonymous function for input to fzero.
What you should do next is look at the M-lint report (listed under Tools in the menu)
and see that there are some warnings. By preallocating the matrices Temp, z and press you can remove the warnings and the code will run faster as a result.

There are some other efficiency improvements you can make, particularly with respect to how you generate the Temp and press matrices.

Hth
Darren

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us