MATLAB: Error: “Input arguments must be convertible to floating-point numbers.”

MATLABmatlab functionsymbolic

Hello
I get an error from line "max" function in the following line:
Z = max(real(Root));
"Root" is a symbolic variable. The full code is as follows; How can I fix this?
clc
clear
close all
eta1 = -0.585;
eta2 = -2.052;
eta3 = -0.367;
syms y1
y2 = 1-y1; % 1: Nitrogen 2: Methane
T = 328.15; % K


P = [81.1 91.1 101 108.3 120.1 131.5 141.8 158.5 169.4 172.6 187.1 207.6 220.6 234.7 248.3 284]; % bar


R = 82.06; % bar.cm^3/mol.K
Tc1 = 748.4; % K
Tc2 = 190.8; % K
Pc1 = 40.5; % bar
Pc2 = 45.79; % bar
omega1 = 0.302;
omega2 = 0.011;
Tr1 = T/Tc1;
Pr1 = P/Pc1;
Mw = 128.1705; % g/mol
ro = 1.14; % g/cm^3
V = Mw/ro; % cm^3/mol
Psat = 1.004*exp(8.583-(3733.9/(T-0)));
if omega1 <= 0.491 && omega2 <= 0.491
F1 = 0.37464+(1.54226*omega1)-(0.26992*omega1^2);
F2 = 0.37464+(1.54226*omega2)-(0.26992*omega2^2);
else
F1 = 0.379642+(1.48503*omega1)-(0.164423*omega1^2)+(0.016666*omega1^3);
F2 = 0.379642+(1.48503*omega2)-(0.164423*omega2^2)+(0.016666*omega2^3);
end
a1 = 0.45724*(((R^2)*(Tc1^2))/Pc1)*(1+F1*(1-((T/Tc1)^0.5)))^2;
a2 = 0.45724*(((R^2)*(Tc2^2))/Pc2)*(1+F2*(1-((T/Tc2)^0.5)))^2;
b1 = 0.0778*R*Tc1/Pc1;
b2 = 0.0778*R*Tc2/Pc2;
k12 = eta1+eta2*log(Tr1)+eta3*log(Pr1);
k21 = k12;
a11 = a1;
a12 = sqrt(a1*a2)*(1-k12);
a21 = sqrt(a2*a1)*(1-k21);
a22 = a2;
a = (y1*y1*a11)+(y1*y2*a12)+(y2*y1*a21)+(y2*y2*a22);
b = y1*b1+y2*b2;
A = (a.*P)./((R.*T).^2);
B = (b.*P)./(R.*T);
Root = zeros(3,16,'sym');
for u = 1:16
Root(:,u) = roots([1 ...
(-(1-B(u)))...
(A(u)-(3*B(u)^2)-2*B(u))...
(-A(u)*B(u)+B(u)^2+B(u)^3)]);
end
Z = max(real(Root));
phi = exp(((b1./b).*(Z-1))-log(Z-B)-((A./(2.*sqrt(2).*B)).*((2.*(y1.*a11+y2.*a21)./a)-(b1./b)).*log((Z+(1+sqrt(2)).*B)./(Z+(1-sqrt(2)).*B))));
g = y1-((Psat./(P.*phi)).*exp(((P-Psat).*V)./(R.*T)));
for u = 1:16
y(u) = vpasolve(g(u));
end
double(y);
abs(y)
yreal = [0.001313 0.001672 0.00292 0.005464 0.01229 0.02114 0.02544 0.03053 0.03387 0.03473 0.03928 0.04224 0.04366 0.04586 0.04969 0.05382];
erro = 0;
n = size(y);
m = size(yreal);
for ii = 1:n
for jj = 1:m
erro = erro+(abs((y(ii)-yreal(jj))/yreal(jj)));
end
end
error = (erro/16)*100

Best Answer

You cannot do that. Root contains references to the unresolved variable y1, so which element is the max can vary according to the y1 value. You would need to convert Z into a large piecewise() containing all 48 possibilities .
Or switch to numeric.
bar = matlabFunction(real(Root));
Z = @(Y) arrayfun(@(y1) max(reshape(bar(y1),1,[])), Y);
and change phi and g to cell array of function handles in y1 (one cell for each corresponding B and P and a21 value -- that is, do the indexing by u ahead of time. Each of those can then be processed numerically with fsolve() or equivalent.
But do keep in mind that there might not be a zero, or at least not one it can find, so if you use fsolve() itself put a try/catch around it.