Problem with "roots" command

4 vues (au cours des 30 derniers jours)
Sepanta Gharib
Sepanta Gharib le 14 Juin 2018
Commenté : OCDER le 14 Juin 2018
Hello
In the following code, the "roots" command should give us 3*6 matrix with 18 values, because A and B have 6 values and the equation is 3rd order. But a 3*1 matrix with 3 values have been obtained. It seems that only one step has passed. What is the problem?
clc
clear
close all
k1 = 0.841;
k2 = 0.537;
k3 = 0.942;
k4 = 0.019;
k5 = 0.720;
T = 333.15:10:383.15; % K
P = 2e7:0.2e7:3e7; % Pa
Tc = 304.12; % K
Pc = 73.74e5; % Pa
Vc = 9.407e-5; % m^3/mol
Zc = 0.274;
w = 0.225;
R = 8.314; % m^3*Pa/mol*K
Mw = 44.01; % g/mol
Tr = T/Tc;
Pr = P/Pc;
m1 = 0.2513+(0.4178*w)-(0.0207*w^2);
m2 = -0.1382-(0.5923*w)-(1.0493*w^2);
m3 = 0.1931+(1.3955*w)+(0.1961*w^2);
alpha = 1+(m1*(1-Tr))+(m2*(1-Tr).^2)+(m3*(1-Tr).^3);
n = 0.2155-(0.9069*w)-(0.206*w^2);
beta = 1+(n*(1-Tr));
a = 0.47448*(((R*Tc)^2)/Pc)*alpha;
b = 0.06824*((R*Tc)/Pc)*beta;
zarib = (k1+k2*log(Tr)+k3*log(Pr))./(1+k4*log(Tr).^2+k5*log(Pr).^2);
A = zarib.*(a.*P)./((R.*T).^2);
B = (b.*P)./(R.*T);
for u = 1:6
Root = roots([1 ...
(-1)...
(A(u)-B(u).^2-2*B(u))...
(-A(u).*B(u)-B(u).^2)]);
zmax = 0;
zmin = 1;
for r = 1:3
if imag(Root(r)) == 0
if Root(r) > 0 && Root(r) > zmax && Root(r) < 1
zmax(u) = Root(r);
end
if Root(r) > 0 && Root(r) < zmin
zmin(u) = Root(r);
end
end
end
end
vm = (zmax.*R.*T)./(P); % m^3/mol
ro = (0.001)^2*(Mw./vm); % g/cm^3

Réponse acceptée

OCDER
OCDER le 14 Juin 2018
Root = roots(...) %Not saving output correctly here
This will override the value of Root every iteration, resulting in a 3x1 matrix. To save the results of each iteration, you need to change where the data is saved to in Root.
Root = zeros(3, 6);
for u = 1:6
Root(:, u) = roots(...);
end
  1 commentaire
OCDER
OCDER le 14 Juin 2018
Replace " Root(r) " with " Root(r, u) "

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Functions dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by