MATLAB & Van der Waals

[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... [g)]
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. [a)]

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... [g)]
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' (http://www.sosmath.com/algebra/factor/fac11/fac11.html)

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.

[mmwave]

Originally posted by enigma
I'm having some problems here...

% 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');
end


I think if you print p inside your loop (decrease the stop value before you try it) you will find that it is a vector in matlab and solve does not like p(1:700). I'm pretty sure you wanted something more like
for i=[1:700]
v=solve('(p(i)*1.01 ... etc);

with the vector p filled out before the for loop. Sorry I don't have time to check this out before I reply.

mmwave.