% 10.34 - Fall 2006 HW#1
% Problem #2
% Equation Solvers
% Rob Ashcraft - Sept. 7th, 2006
% accompanying files (included in this function): fmin_eqn.m and fzero_eqn.m
function Problem2
% This function serves to set up the problem and constants used in solving
% the RK equation of state; generates a plot of the final solution vs. initial guess
% clear the workspace and command window
clear; clc;
% define certain global variables that are available to other functions
% without explicitly passing them (a global statement must also occur in
% the function where these are needed)
global psi omega sigma eps Tc Pc P R w;
% define the variable EOS parameters and physical constants
psi = 0.42748;
omega = 0.08664;
sigma = 1;
eps = 0;
Tc = 647.1; % K
Pc = 217.7; % atm
P = 1; % atm
R = 0.08206; % atm-L/mole-K
w = 0.152; % acentric factor
T = 300; % temp in K
b = omega*R*Tc/Pc; % b parameter in the EOS (used in the plot)
V_0_vec = linspace(1,50,50); %vector of initial guesses in L/mole
% solver options. These can be very important with complex problems. This
% particular option turns off the output display, so it does not have to
% print it to the command window (particularly annoying with iterating over fsolve)
options = optimset('Display','off');
% solve the equation using the three solvers specified
for i=1:length(V_0_vec)
V_fzero_300K_gas(i) = fzero(@fzero_eqn,V_0_vec(i),options,T);
V_fsolve_300K_gas(i) = fsolve(@fzero_eqn,V_0_vec(i),options,T);
V_fminsearch_300K_gas(i) = fminsearch(@fmin_eqn,V_0_vec(i),options,T);
end
disp(['Vapor Volume Root (fzero, V_guess=50) (L/mole): ',num2str(V_fzero_300K_gas(length(V_0_vec)))]);
disp(['Vapor Volume Root (fsolve, V_guess=50) (L/mole): ',num2str(V_fsolve_300K_gas(length(V_0_vec)))]);
disp(['Vapor Volume Root (fminsearch, V_guess=50) (L/mole): ',num2str(V_fminsearch_300K_gas(length(V_0_vec)))]);
% plot the solver solution as a function of initial guess
figure;
plot(V_0_vec,V_fzero_300K_gas,'.',V_0_vec,V_fsolve_300K_gas,'x',V_0_vec,V_fminsearch_300K_gas,'o',[0,50],[b,b],'--')
xlabel('Volume Initial Guess [L/mole]'); ylabel('Resulting Volume [L/mole]');
legend('fzero Result','fsolve Result','fminsearch Result','EOS b-value (lower physical limit)');
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Include other functions in this file.
% This function defines the equation that should be equal to zero as a
% function of Volume and Temperature; this is called when using the
% "fzero" or "fsolve" solver function
function resid = fzero_eqn(V,T)
%inputs:
% V = specific volume in L/mole
% T = temperature of the system in Kelvin
%outputs:
% resid = residual of the EOS equation, which should -> 0 when solved
% correctly
% this defines the global variables available to this function
global psi omega sigma eps Tc Pc P R w;
Tr = T/Tc; %reduced temperature
% Define the EOS parameters
alpha = Tr^(-1/2);
a = psi*alpha*(R*Tc)^2/Pc;
b = omega*R*Tc/Pc;
% Equation to solve (set to zero). Notice the P had to be moved to the RHS
% so that the equation would equal zero. This is a recurring theme with
% many solvers
resid = -P + R*T/(V-b) - a/((V+eps*b)*(V+sigma*b));
return;
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
% This function is used by fminsearch to find the optimum value of the
% variable that minimizes an objective function, Esqr in this case.
function Esqr = fmin_eqn(V,T)
%inputs:
% V = specific volume in L/mole
% T = temperature of the system in Kelvin
%outputs:
% Esqr = the square of the residual of the EOS equation, which should be
% minimized to 0 when solved correctly
global psi omega sigma eps Tc Pc P R w;
Tr = T/Tc;
alpha = Tr^(-1/2);
a = psi*alpha*(R*Tc)^2/Pc;
b = omega*R*Tc/Pc;
resid = -P + R*T/(V-b) - a/((V+eps*b)*(V+sigma*b));
% in this case, the plain residual cannot be used since a minimization is
% being performed, and not trying to set something to zero. The easy way
% to solve this is to simply square the residual to make it always
% positive. If you had multiple equations, you could use the sum of
% squares of errors as the residual in a similar manner.
Esqr = resid^2;
return;