1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Root-Finding method in MATLAB

  1. Dec 3, 2016 #1
    1. The problem statement, all variables and given/known data
    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:

    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.

    2. Relevant equations

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

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

    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);

    3. The attempt at a solution

    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);
    x = x - (y/ypr);
  2. jcsd
  3. Dec 4, 2016 #2
    It doesn't look like you're setting abs(y) to anything in your while loop.
  4. Dec 4, 2016 #3


    Staff: Mentor

    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.
    Code (Matlab M):
    while abs(y) > tol;
        for P = 1:31
           y = vdw(x);  ;; <---  y gets a new value here
          ypr = vdwpr(x);
        x = x - (y/ypr);
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted