Solving VdW's Eq of State with MATLAB

[chronicals]

Homework Statement

Van der Waals’s equation of state for an imperfect gas is

( P + a / v ^ 2) ( v - b ) = R T

where P is the pressure (atm), v is the molar volume (liters/mole), T is the absolute temperature (K), and a (liter
2.atm/mol2), and b (liter/mol) are constants that depend on the particular gas. Find
molar volume of carbon dioxide at 3 atm and 320 K. For carbon dioxide, you may use a=3.592

and b=0.04267. (Hint: Use ideal gas law to make your initial guess.)

a) Write a generic MATLAB m-file that uses Newton’s (Newton-Raphson) method to solve for x
given the function f(x), initial guess x0, and error tolerance es. Make use of the attached
pseudocode given for the fixed-point iteration method.
b) Write a MATLAB m-file that uses the generic m-file you have developed in part (a) to solve for
the molar volume of any gas, given the P, T, R, a, and b. Use a relative error tolerance of
0.01% for convergence of molar volume.
c) Run the m-file you have developed in part (b) along with the generic m-file of part (a) for the
case of carbon dioxide given above. Make sure you display the result at each iteration.
d) Plot error as a function of the iteration step. Do not forget to label the axes.
e) MATLAB has a built-in function called fzero to find the real root of a single equation. A
simple representation of its syntax is
fzero(function, x0)
where function is the name of the function being evaluated, and x0 is the initial guess.
Solve the above problem using fzero. Comment on the result by comparing with the result
of (c).

Homework Equations

The Attempt at a Solution

How can i solve this question?

 

[Päällikkö]

As required by the rules of the homework questions forum, please show your own attempt at a solution.

 

[Mark44]

The rules are here: https://www.physicsforums.com/showthread.php?t=5374.

Homework Help:
On posting questions: Any and all high school and undergraduate homework assignments or textbook style exercises for which you are seeking assistance are to be posted in the appropriate forum in our Homework & Coursework Questions area--not in blogs, visitor messages, PMs, or the main technical forums. This should be done whether the problem is part of one's assigned coursework or just independent study. The reason for this is that the scientific and mathematical sections of Physics Forums are to be reserved for discussions and not academic assistance. Since graduate level assignments are meant to be more thought provoking (and hence more worthy of discussion), graduate level questions will be allowed in the relevant part of the main section of PF, provided that the graduate student attempts the problem and shows his work. NOTE: You MUST show that you have attempted to answer your question in order to receive help. You MUST make use of the homework template, which automatically appears when a new topic is created in the homework help forums.

We do not support cheating in any form: Do not ask for solution manuals, answers to exams, or instructor's manuals. Every school and instructor has their own policies or honor codes on what constitutes cheating, and it is up to the individual student to adhere to those policies when seeking help here. If you are in doubt as to whether you are permitted to seek help, consider erring on the side of caution and not asking for help.

On helping with questions: Any and all assistance given to homework assignments or textbook style exercises should be given only after the questioner has shown some effort in solving the problem. If no attempt is made then the questioner should be asked to provide one before any assistance is given. Under no circumstances should complete solutions be provided to a questioner, whether or not an attempt has been made.

 

[chronicals]

This is my m-file,can you check it,please? I can't see the plot?Why? I use [root]=Newtoniandream( ' y ', 1e-3)


function [f, derf ]=y(v)
a=3.592; b=0.04267; P=3; T=320; R=0.082;
f= (P+a/v^2)*(v-b)-R*T;
derf=P-a/v^2+2*a*b/v^3;


function [root]=Newtoniandream(v0,tol)
v=v0;
n=0;
err=3;
fprintf('iteration x percent relative error\n')
while err>tol
n=n+1;
[f, derf ]=feval('y',v);
vnew=v-f/derf;
err(n)=(abs((vnew-v)/vnew))*100;
fprintf('%2d %f %f\n',n,vnew,err(n))
v=vnew;
plot(n,err)
hold on
end
root=vnew

 

[Mark44]

This is my m-file,can you check it,please? I can't see the plot?Why? I use [root]=Newtoniandream( ' y ', 1e-3)


function [f, derf ]=y(v)
a=3.592; b=0.04267; P=3; T=320; R=0.082;
f= (P+a/v^2)*(v-b)-R*T;
derf=P-a/v^2+2*a*b/v^3;


function [root]=Newtoniandream(v0,tol)
v=v0;
n=0;
err=3;
fprintf('iteration x percent relative error\n')
while err>tol
n=n+1;
[f, derf ]=feval('y',v);
vnew=v-f/derf;
err(n)=(abs((vnew-v)/vnew))*100;
fprintf('%2d %f %f\n',n,vnew,err(n))
v=vnew;
plot(n,err)
hold on
end
root=vnew

You're not using the plot command correctly. It should be used on two arrays, not on a single pair of coordinates, which means it should probably not be inside your while loop.

Also, I don't think you are using err correctly. In one line, you set err to 3, but inside your while loop, you treat err as an array.

Some documentation that might be helpful is here: http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/getstart.pdf.

Ch. 2 deals with data types, including matrices. See Working with Matrices, p. 2-16, and the following pages.

Ch. 3 deals with plotting. See Using Basic Plotting Functions on page 3-56.

Ch. 4 deals with Flow Control, including functions. See Functions, page 4-22.