Message Boards Message Boards

0
|
5310 Views
|
8 Replies
|
0 Total Likes
View groups...
Share
Share this post:

while loop issues

Greetings, I am new to Mathematica. I am trying to do a while loop to do an iteration. As long as the absolute value of Y1 is higher than 1E-4 then Y becomes equal to Y1, a new F, a new Fpr, and a new Y1 are calculated and the process is repeated. Unfortunately it doesn't work and the loop keeps on running until I abort it. I realize that this is a simple issue but I'd really appreciate any help I can get. For general information,the purpose of this code is to estimate the compressibilty factor for natural gases using the Hall-Yarborough method.enter image description here

POSTED BY: Derek du Plessix
8 Replies
Posted 11 years ago

Hi Derek,

I did some syntax correction and divided your code into 3 sections. The first establishes initial values for the iterations. The second performs the iterations 10 times and accumulates the Y1 results. It shows that the iterations quickly take Y1 to a stable value greater than the termination condition. The third it the while loop, but of course it does not terminate. I don't know enough of the fundamentals to suggest more. perhaps other initial conditions would behave differently, but maybe not.

Kind regards, David

PS. In Mathematica, it is good to use variable names starting with lower case, since reserved words all begin with upper.

In[1]:= (* set up initial values *)
gamma = 21/28.96;
P = 2000;
T = 581.67;
Tpc = 168 + 325*gamma - 12.5*gamma^2;
Ppc = 667 + 15*gamma - 37.5*gamma^2;
Tpr = T/Tpc;
t = 1/Tpr;
Ppr = 2000/Ppc;
X1 = -0.06125*Ppr*t*Exp[-1.2*(1 - t)^2];
X2 = (14.76*t - 9.76*t^2 + 4.58*t^3);
X3 = (90.7*t - 242.2*t^2 + 42.4*t^3);
X4 = (2.18 + 2.82*t);
Y = 0.0125*Ppr*t*Exp[-1.2*(1 - t)^2];
F = X1 + (Y + Y^2 + Y^3 + Y^4)/(1 - Y)^3 - X2*Y^2 + X3*Y^X4;
Fpr = (1 + 4*Y + 4*Y^2 - 4*Y^3 + Y^4)/(1 - Y)^4 - 2*X2*Y + 
   X3*X4*Y^(X4 - 1);
Y1 = Y - F/Fpr;

In[17]:= (* do the first 10 iterations and accumulate Y1 in a list *)

y1Values = {};
Do[
 Y = Y1; 
 F = X1 + (Y + Y^2 + Y^3 + Y^4)/(1 - Y)^3 - X2*Y^2 + X3*Y^X4;
 Fpr = (1 + 4*Y + 4*Y^2 - 4*Y^3 + Y^4)/(1 - Y)^4 - 2*X2*Y + 
   X3*X4*Y^(X4 - 1);
 Y1 = Y - F/Fpr;
 AppendTo[y1Values, Y1],
 {10}
 ]

In[19]:= (* appears to stagnate with 0.147417 a stable value for y1 *)
\
    y1Values

Out[19]= {0.14882, 0.147323, 0.147424, 0.147417, 0.147417, 0.147417, \
0.147417, 0.147417, 0.147417, 0.147417}

In[20]:= (* the While will not terminate *)
While[Abs[Y1] > 10^-4,
  Y = Y1; 
  F = X1 + (Y + Y^2 + Y^3 + Y^4)/(1 - Y)^3 - X2*Y^2 + X3*Y^X4;
   Fpr = (1 + 4*Y + 4*Y^2 - 4*Y^3 + Y^4)/(1 - Y)^4 - 2*X2*Y + 
    X3*X4*Y^(X4 - 1) ;
  Y1 = Y - F/Fpr
  ];

Out[20]= $Aborted
POSTED BY: David Keith

I realized that the condition which had to be met was: Abs[Y -Y1] < 10^-4, in other words F/Fpr < 10^-4 not Y1 < 10^-4. Thank you, your post was really helpful because now I know a procedure to follow when I need to figure out what's going on. Best regards,

Derek du Plessix

POSTED BY: Derek du Plessix
Posted 11 years ago

As Daniel says, it would be a big help if you would post code instead of an image. You can find instructions on posting here: How to post

Also, your algorithm appears to be trying to accomplish a fit to some function. You might tell us more about this at the top level, since Mathematica very often has high level functions which serve such purposes better than procedural programming. In my experience, For, Do, and While are seldom the best construct in Wolfram code.

Kind regards, David

POSTED BY: David Keith

Hi David, Thank you for your response. I'm not trying to fit a funtion. I am just trying to continuously calculate Y1 until I reach the condition abs[Y1] < 10^-4.

POSTED BY: Derek du Plessix
Posted 11 years ago

Inside the while loop, you probably want to replace the double equal signs with single equal signs and put the semicolon and the end of those lines.

POSTED BY: Gustavo Delfino
gamma = 21/28.96;
P = 2000 ;
T = 581.67;
Tpc = 168 + 325* gamma - 12.5*gamma^2;
Ppc = 667 + 15*gamma - 37.5*gamma^2;
Tpr = T/Tpc;
t = 1/Tpr;
Ppr = 2000/Ppc;
X1 = -0.06125*Ppr *t* Exp[-1.2* (1 - t)^2] ;
X2 = (14.76 *t - 9.76*t^2 + 4.58*t^3);
X3 = (90.7* t - 242.2*t^2 + 42.4*t^3);
X4 = (2.18 + 2.82*t);
Y = 0.0125*Ppr*t*Exp[-1.2*(1 - t)^2]
F = X1 + (Y + Y^2 + Y^3 + Y^4)/(1 - Y)^3 - X2*Y^2 + X3*Y^X4
Fpr = (1 + 4*Y + 4*Y^2 - 4*Y^3 + Y^4)/(1 - Y)^4 - 2*X2*Y + 
  X3*X4*Y^(X4 - 1)
Y1 = Y - F/Fpr
While[Abs[Y1] > 10^-4,
 Y == Y1
   F == X1 + (Y + Y^2 + Y^3 + Y^4)/(1 - Y)^3 - X2*Y^2 + X3*Y^X4
    Fpr == (1 + 4*Y + 4*Y^2 - 4*Y^3 + Y^4)/(1 - Y)^4 - 2*X2*Y + 
   X3*X4*Y^(X4 - 1)
    Y1 == Y - F/Fpr]
Y1
POSTED BY: Derek du Plessix

What you show is an image, and that cannot be cut and pasted into Mathematica. SO it's not really possible for anyone to experiment with the code.

POSTED BY: Daniel Lichtblau

Thank you for your time, I've added the code.

POSTED BY: Derek du Plessix
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract