How Do You Implement Root-Finding for Compressibility Factors in MATLAB?

[ShakeSpee]

Homework Statement

Compute and plot the compressibility factor (y) verses pressure (x) for the (1) Van der Waal’s (2) Redlich-Kwong and (3) Peng-Robinson equations of state.

Compressibility Factor, Z = (P*v)/(R*T); where v is the specific volume (V/v).

Data for n-Butane:
T = 500 K; Tc = 425.2 K; Pc = 37.5 atm; R = 0.08206 L.atm/(mol.K);
for P = 1:31 (atm)
Equations of State:

Hint:
1. Make P, Pc, T, Tc, and R global variables
2. Write three functions. One for each of the equations.
3. Write a script file that calls the functions using a root finding method to determine a root
4. Use the root to calculate the compressibility factor
5. Plot the compressibility factor verses pressure.
Note: Show all three graphs in the same plot window. Properly label the axis etc.

Homework Equations

Van der Waal[/B]:
a = 0.42188*(R^2*Tc^2/Pc)
b = 0.125(RTc/Pc)

Pv^3 - (Pb + RT)*v^2 + (a)v + ab = 0.

Redlich-Kwong:
Tr = T/Tc;
alpha = 1/(Tr^0.5);
a = 0.42748*(((R^2)*(Tc^2))/Pc)*alpha;
b = 0.08664*((R*Tc)/Pc);

(P)*(x^3) - (R*T)*(x^2) + (a - P*(b^2) - R*T*b)*x - a*b;

Peng- Robinson:
w = 0.193;
m = 0.37464 + 1.54226*w - 0.26992*w^2;
Tr = T/Tc;
alpha = (1 + m*(1 - (Tr^0.5)))^2;
a = 0.45724*(((R^2)*(Tc^2))/Pc)*alpha;
b = 0.07780*((R*Tc)/Pc);

(P)*(x^3) + (b*P - R*T)*(x^2) + (a -3*P*(b^2) - 2*R*T*b)*(x) + (P*(b^3) + R*T*(b^2) - a*b);

The Attempt at a Solution

[/B]
I finished the first two steps - I created function scripts for all of the equations, but I'm stuck on the third part, which is finding the root of one of the functions. I can use any method to find the root, and for now, I chose the Newton-Raphson Method, so I also created scripts for the derivatives of each function. Whenever I run this script, it goes into an infinite loop, and I can't figure out why. Here is what I have so far:%This script will take a user input (guess) and use the Newton - Raphson
%Method to determine the root of the function. y = vdw(x) and ypr =
%vdwpr(x), the derivative of vdw(x).
T = 500;
Tc = 425.2;
Pc = 37.5;
R = 0.08206;

%Set converge tolerance.
tol = 0.00001;

%Input initial guess.
x = input('Enter initial guess: ');

y = vdw(x);

while abs(y) > tol;
for P = 1:31
y = vdw(x);
ypr = vdwpr(x);
end
x = x - (y/ypr);
end

 

[TJGilb]

It doesn't look like you're setting abs(y) to anything in your while loop.

 

[Mark44]

It doesn't look like you're setting abs(y) to anything in your while loop.

You don't "set abs(y)" to a value -- you set y, and take its absolute value, which the OP is doing in the while loop.

while abs(y) > tol;
    for P = 1:31
       y = vdw(x);  ;; <---  y gets a new value here
      ypr = vdwpr(x);
    end
    x = x - (y/ypr);
end