Error: "Input arguments must be convertible to floating-point numbers."

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

Answer 1

[Walter Roberson]

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.