How to create a scalar set of values

4 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"

Tags

Community Treasure Hunt

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

Start Hunting!