MATLAB Pressure-Density Relationship with Van der Waals Equation in MATLAB

[enigma]

I'm having some problems here...

Part of an assignment is to find how density of a gas changes as a function of pressure (1 - 700 atmospheres) using Van der Waals. Unless I've really missed the boat, there is no closed form solution for specific volume in that equation,

(p+a/v^2)(v-b)=RT

so a numerical method will need to be used.
I tried using a <for> loop to use Matlab's solve function, but I can't get the blasted thing to accept the variable <p> to operate in the loop.

Anybody have any ideas? I REALLY don't want to have to calculate the value for density 700 times for each of 2 gases... no... really, I don't.

% Nitrogen
% density variable = rho_N
% a=179.65
% b=.001398
% m=28
% R=8314
% T=298.15
for p=[1:700] % atmospheres
    FLAG=false; % Catch multiple real answers
    % (p*Pascals/atm + a/v^2)*(v-b) = R/m*t 
    v=solve('(p*1.01*10^5+179.65/v^2)*(v-.001398)=8314/28*298.15');
    for i=[1:length(v)] % find real component
        if isreal(v(i))
            if FLAG==true; %Catch multiple solutions
                error('multiple real solutions for N2 VdW')
            end
            rho_N(p)=1/v(i);
            FLAG==true;
        else % No reals error
            p
            error('no real solutions for N2 VdW')
        end
    end
end

 

[Bystander]

Van der Waals did NOT have access to computers and numerical methods; that's hint #1.

 

[enigma]

Wow... that was helpful, Bystander...
Van der Waals also wasn't asked to churn out 1400 data points by hand for a rocket propulsion systems homework assignment...

Anyway, I couldn't find a way to get the solve function to accept a parameter, so I dumped it and replaced it with a Newton-Raphson iteration.

Graphs came out fine.

EDIT: If anyone knows how to get the Matlab symbolic solver to accept a parameter, I'm still interrested.

 

[Bystander]

Originally posted by enigma
Wow... that was helpful, Bystander...
Van der Waals also wasn't asked to churn out 1400 data points by hand for a rocket propulsion systems homework assignment...

Needed something less subtle? Yes, you missed the boat --- VDW is a cubic EOS --- something engineers love for just the property you've missed. It does require that the user pick one of three roots, but even engineers are capable of that much.

1400 pts.? VdW did more than that a day, by hand --- 'course, it wasn't "rocket science" in his day. Quit yer belly-aching.

 

[enigma]

Is this what you're referring to: 'Cubic Formula'

You've got to be joking. That's 25 steps. This sort of problem kinda does need to be done on a computer... hence me asking a question about Matlab code. Do you honestly think that coding a 25 step algoritm is an effective way to solve something when (assuming I found an answer to my question) it could be solved by typing solve('function'); ?!? At _very_ least 25 lines vs. 5 lines for 'solve' or 12 lines for Newton-Raphson including variable definition and error catches. And you're insulting engineers' problem solving skills? Please...

I don't feel I am bellyaching. I am irritated that I asked for help with a question about computer code and got a deliberately vague, tongue-in-cheek insult about my math skills.

 

[UntchbleConce]

I know this post is really old. But i was having a similar problem. I first went to google, and when that didn't work. I actually had to think! lol. Either way I figured it out. You just need to create a user defined function with the solve function in it. From there you create a for loop with inputs for the function and tada! Since I can never make sense when people explain their code, see below to see it actually done. (i kinda just copied and pasted out of my code without my comments (cause it just looks messy in here with it so i just sort of added another set of comments to try to keep this post short-er) Hope this helps!

% My user defined function: See link below for more info on this, if needed
% http://www.mathworks.com/help/techdoc/ref/function.html
function a2 = ind_factor(lambda)
syms a l
S=solve('(((1-4*a)^2*(1-a))/(1-3*a))=l^2','a');
l=lambda;
% It's kinda long to explain why i did the next step, just know I needed the real part of the
% second row every time.
a2=real(eval(S(2,);

% How I used it in my code
lambda=0.5:0.1:10;
for i=1:length(lambda)
a2(i)=ind_factor(lambda(i));
end