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.
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.
gamma = 21/28.96;
P = 2000 ;
T = 581.67;
Tpc = 168 + 325* gamma - 12.5gamma^2;
Ppc = 667 + 15gamma - 37.5gamma^2;
Tpr = T/Tpc;
t = 1/Tpr;
Ppr = 2000/Ppc;
X1 = -0.06125Ppr t Exp[-1.2* (1 - t)^2] ;
X2 = (14.76 t - 9.76t^2 + 4.58t^3);
X3 = (90.7 t - 242.2t^2 + 42.4t^3);
X4 = (2.18 + 2.82t);
Y = 0.0125PprtExp[-1.2*(1 - t)^2]
F = X1 + (Y + Y^2 + Y^3 + Y^4)/(1 - Y)^3 - X2Y^2 + X3Y^X4
Fpr = (1 + 4Y + 4Y^2 - 4Y^3 + Y^4)/(1 - Y)^4 - 2X2Y +
X3X4Y^(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 - X2Y^2 + X3Y^X4
Fpr == (1 + 4Y + 4Y^2 - 4Y^3 + Y^4)/(1 - Y)^4 - 2X2Y +
X3X4Y^(X4 - 1)
Y1 == Y - F/Fpr]
Y1
Thank you for your time, I’ve added the code.
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.
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
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.
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
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