How to create a scalar set of values

3 views (last 30 days)
Quincey Jones
Quincey Jones on 10 Feb 2020
Answered: fred ssemwogerere on 11 Feb 2020
function [v,err,count] = Clapeyrons(R,T,a,b,P);
% Calculate molar volume as predicted by the Ideal Gas Law equation
% using the Newton-Raphson algorithm with initial estimate
% determined by ideal gas law.
% Inputs:
R = 0.082057; % Ideal gas constant in L atm / K mol
T = 293; % Temperature in K
Tc = 416.90; % Critical temperature of Cl2 in K
Pc = 78.72918; % Critical Pressure of Cl2 in atm
% a and b are Van der Waals constants for gas considered in
% atm L2/mol2 and L/mol respectively.
a = ((R^2)*(Tc^(5/2)))/(9*Pc*(2^(1/3)-1));
b = (R*Tc*(2^(1/3)-1))/(3*Pc);
P = 1; % Pressure in atm
% Outputs:
% v equals molar volume in L/mol as predicted by
% equation for gas considered at temperature T and pressure P.
% err equals modulus of function evaluated at approximate root.
% count is number of iterations taken by Newton-Raphson algorithm.
% Version 1: created 24/04/19. Author: Paul Curran
% Version 2: created 06/02/20. Author: Paul Curran
if (~isscalar(P)) || (~isreal(P)) || P <= 0
error('Input argument P must be positive real scalar.')
end
Iteration_limit = 20; % maximum number of iterations permitted
Tolerance=10^-7; % maximum acceptable value for modulus of
A = (a*P)/((R^2)*(T^(5/2)));
B = (b*P)/(R*T);
v = R * T / P;
Z = (P*v)/(R*T);
% Employ ideal gas law to obtain initial estimate
for count = 1:Iteration_limit + 1
% Terminate with error message if iteration limit exceeded:
if count == Iteration_limit + 1
error('Iteration limit reached. Iteration did not converge.')
break
end
f = (Z^3) - (Z^2) + (A-B-(B^2))*Z - (A*B);
% Terminate iteration if function is sufficiently small at current
% estimate
if abs(f)<Tolerance
break
end
f_dash = 3*Z^2 - 2*Z + (A - B - (B^2)); % Evaluate derivative of f
Z = Z - (f/f_dash); % Newton-Raphson iteration
v = Z*R*T/P; % molar volume found using Newton-Raphson
end
err=abs(f); % Error is magnitude of f(v) at final root estimate
end
I need to change P=1 so...
P = [1,1.5,2,2.5,3,5,10,15,25,50,100];
R = 0.082057;
T = 293;
Tc = 416.90;
Pc = 78.72918;
V_real =zeros(11,1);
V_ideal =zeros(11,1);
for count = 1:11
V_real(count) = Clapeyrons(P(count));
V_ideal(count) = (R*T)/P(count);
%V_difference = V_ideal - V_waals;
end
plot(P,V_real)
hold on
plot(P,V_ideal,'x')
title('Molar Volume vs Pressure for Cl2')
xlabel('Pressure P (atm)')
ylabel('Molar Volume v (L/mol)')
legend({'Redlich-Kwong Equation','Ideal Gas Law'},'Location','northeast')
...this graph for the Redlich Kwon Equation line will go through the same P values as in the Ideal Gas Law

Answers (1)

fred  ssemwogerere
fred ssemwogerere on 11 Feb 2020
Hello, using arrayfun could be more efficient:
% Taking your vector, say "P" to compute real "v_real" using "arrayfun". "arrayfun" applies the function or formula to each value of "P" to get your outputs
P = [1,1.5,2,2.5,3,5,10,15,25,50,100];
[v_real,err,count]=arrayfun(@(x) Clapeyrons(R,T,a,b,x), P);
% The outputs from the above function will be 3 vectors of "v_real","err" and "count", for each point in "P".
% Computing the ideal "v" values
v_ideal=arrayfun(@(x) (R*T/x),P); % output will be a vector of values for each "P" value
% From here you can proceed with plotting your "v_real" and "v_ideal" vectors with "P"

Categories

Find more on Programming in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!