How to create a scalar set of values

3 次查看(过去 30 天)
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

回答(1 个)

fred  ssemwogerere
fred ssemwogerere 2020-2-11
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"

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by