Use fsolve to solve matrix

Would like to use fsolve to solve xeq.

When i multiply with vector 'T' and 'Y' I only get one value.

How would I make it return a vector. THe issue is at the botto of the code

function ydot=Zazarian_f_5(~,x)

X=x(1); y=x(2); T=x(3);

% given & calculation

T0=325;
Ta=300;
Fa0=5;
Fb0=10;
Fi0=10;
Ca0=0.2;
Cp_a=20;
Cp_b=20;
Cp_c=20;
Cp_i=18;
alpha=0.00015;
k1=0.0002;
Ua=320;
pb=1400;
e=25000;
mc=18;
thb=2;
thi=2;
delta_Hrx=-20000;
Cp_cool=18;
r=1.987;
kc1=1000;
sumcp=Cp_a+Cp_b*thb+Cp_i*thi;
%calculation

k2=k1*exp((e/r)*((1/300)-(1/T)));

Kc2=kc1*exp((-delta_Hrx/r)*((1/305)-(1/T)));

CA=Ca0*(1-X)*y*(T0/T);

CB=Ca0*(2-X)*y*(T0/T);

CC=2*Ca0*X*y*(T0/T);

Ci=2*Ca0*y*(T0/T);

ra=-k2*((CA*CB)-((CC^2)/Kc2));

ydot(1)=(-ra)/(Fa0);

ydot(2)=(-alpha/(2*y))*(T/T0);

ydot(3)=(-delta_Hrx*ra)/(Fa0*sumcp);

ydot = ydot';

end 




clear all
clc
y0=[0 1 325]; 

W=[0 1000];
r=1.987;
kc1=1000;
Tkc1 = 305;
T0=325;
delta_Hrx=-20000;

[t,z]=ode45(@Zazarian_f_5, W, y0);
T = z(:,3);
y=z(:,2);
X = z(:,1);

figure
plot(t, T);

To=325;
plot(t, z(:,3));
xlabel('catalyst weight in kg');
ylabel('Temperature');
legend('T');
title('Temperature');

cao=0.2;

Kc=kc1*exp((-delta_Hrx/r)*((1/Tkc1)-(1/T)));

Z=To./T.*y;

cae=@(xeq)cao.*(1-xeq).*Z;
cbe=@(xeq)cao.*(thb-2*xeq).*Z;
cce=@(xeq)cao.*2.*xeq.*Z;
F=@(xeq) cce.^2./(cae.*cbe)-kc1*exp((dhrx/R)*((1/305)-(1./T)));
xeq=fsolve(F,.7);

figure
plot(t, z(:,1),t,xeq);
ylim([0,1.2])
xlabel('catalyst weight in kg');
ylabel('conversion X');
legend('X','Xe');
title('converison X and Xe')

Undefined operator '.^' for input arguments of type
'function_handle'.

Error in
zzz>@(xeq)cce.^2./(cae.*cbe)-kc1*exp((dhrx/R)*((1/305)-(1./T)))

Error in fsolve (line 242)
            fuser = feval(funfcn{3},x,varargin{:});

Error in zzz (line 44)
xeq=fsolve(F,.7);
 
>> 

Answer 1

[Walter Roberson]

F = @(xeq) cce(xeq).^2 ./ (cae(xeq)  .* cbe(xeq)) - kc1*exp((dhrx/R)*((1/305)-(1./T)));

It looks like the part starting from kc1 is constant, so for efficiency calculate the term first and assign to a variable and use the variable.

[Rodrigo Blas]

I get the same undefined operator error

clear all
clc
y0=[0 1 325]; % why a initial pressure drop of 1
W=[0 1000];
r=1.987;
kc1=1000;
Tkc1 = 305;
T0=325;
delta_Hrx=-20000;
% X vs W
[t,z]=ode45(@Zazarian_f_5, W, y0);
T = z(:,3);
y=z(:,2);
X = z(:,1);
%T = T0 + (X*delta_Hrx)/(60);
% Temperature vs W
figure
plot(t, T);
%using other guys
To=325;
plot(t, z(:,3));
xlabel('catalyst weight in kg');
ylabel('Temperature');
legend('T');
title('Temperature');
% step 2 Kc -- equilibrium constant equation at a t
cao=0.2;
%Kc2=kc1*exp((-delta_Hrx/r)*((1/305)-(1/T)));
Kc=kc1*exp((-delta_Hrx/r)*((1/Tkc1)-(1/T)));
% step 6 xe equation
%base = (3*Kc)/4;
%base2 = (Kc/4)-1;
%top =((base)-sqrt(base.^2 - 2.*Kc.*(base2)));
%bot = (2.*(base2));
%Xe = top/bot;
%Xe = ((base)-sqrt(base.^2 - 2.*Kc.*(base2)))/(2.*(base2));
%Xe = (((3*Kc)/4)-sqrt(((3*Kc)/4).^2 - (2*Kc)*((Kc/4)-1)))/(2*((Kc/4)-1));
Z=To./T.*y;
dhrx=-20000;
R=1.987;
cae=@(xeq)cao.*(1-xeq).*Z;
cbe=@(xeq)cao.*(thb-2*xeq).*Z;
cce=@(xeq)cao.*2.*xeq.*Z;
KC=kc1*exp((dhrx/R)*((1/305)-(1./T)));
F=@(xeq) cce.^2./(cae.*cbe)-KC;
xeq=fsolve(F,.7);

figure
plot(t, z(:,1),t,xeq);
ylim([0,1.2])
xlabel('catalyst weight in kg');
ylabel('conversion X');
legend('X','Xe');
title('converison X and Xe')

function ydot=Zazarian_f_5(~,x)

X=x(1); y=x(2); T=x(3);

% given & calculation

T0=325;
Ta=300;
Fa0=5;
Fb0=10;
Fi0=10;
Ca0=0.2;
Cp_a=20;
Cp_b=20;
Cp_c=20;
Cp_i=18;
alpha=0.00015;
k1=0.0002;
Ua=320;
pb=1400;
e=25000;
mc=18;
thb=2;
thi=2;
dhrx=-20000;
Cp_cool=18;
r=1.987;
kc1=1000;
sumcp=Cp_a+Cp_b*thb+Cp_i*thi;
%calculation

k2=k1*exp((e/r)*((1/300)-(1/T)));

Kc2=kc1*exp((-dhrx/r)*((1/305)-(1/T)));

CA=Ca0*(1-X)*y*(T0/T);

CB=Ca0*(2-X)*y*(T0/T);

CC=2*Ca0*X*y*(T0/T);

Ci=2*Ca0*y*(T0/T);

ra=-k2*((CA*CB)-((CC^2)/Kc2));

ydot(1)=(-ra)/(Fa0);

ydot(2)=(-alpha/(2*y))*(T/T0);

ydot(3)=(-dhrx*ra)/(Fa0*sumcp);

ydot = ydot';
end

I get the same undefined operator error.

[Walter Roberson]

You had several problems, with the first one being that you missed the part where I passed xeq into each of the function handles.

You also had thb undefined in a function handle.

Your fsolve was never going to work. Your first part of your equation is scalar, and you subtract from that the vector KC (it is a vector because it involves 1/T and T is a vector). It is not possible to find a single xeq that solves all of those equations simultaneously -- and if you did manage to do so then it would give you a single scalar result where your plot() call shows that you expect xeq to be the same length as t is. So we deduce that what you want is to solve for each t value individually.

y0=[0 1 325]; % why a initial pressure drop of 1
W=[0 1000];
r=1.987;
kc1=1000;
Tkc1 = 305;
T0=325;
thb = 2;
delta_Hrx=-20000;
% X vs W
[t,z]=ode45(@(t,x) Zazarian_f_5(t,x,thb), W, y0);
T = z(:,3);
y=z(:,2);
X = z(:,1);
%T = T0 + (X*delta_Hrx)/(60);
% Temperature vs W
figure
plot(t, T);
%using other guys
To=325;
plot(t, z(:,3));
xlabel('catalyst weight in kg');
ylabel('Temperature');
legend('T');
title('Temperature');
% step 2 Kc -- equilibrium constant equation at a t
cao=0.2;
%Kc2=kc1*exp((-delta_Hrx/r)*((1/305)-(1/T)));
Kc=kc1*exp((-delta_Hrx/r)*((1/Tkc1)-(1/T)));
% step 6 xe equation
%base = (3*Kc)/4;
%base2 = (Kc/4)-1;
%top =((base)-sqrt(base.^2 - 2.*Kc.*(base2)));
%bot = (2.*(base2));
%Xe = top/bot;
%Xe = ((base)-sqrt(base.^2 - 2.*Kc.*(base2)))/(2.*(base2));
%Xe = (((3*Kc)/4)-sqrt(((3*Kc)/4).^2 - (2*Kc)*((Kc/4)-1)))/(2*((Kc/4)-1));
Z=To./T.*y;
dhrx=-20000;
R=1.987;

cae=@(xeq)cao.*(1-xeq).*Z;
cbe=@(xeq)cao.*(thb-2*xeq).*Z;
cce=@(xeq)cao.*2.*xeq.*Z;
KC=kc1*exp((dhrx/R)*((1/305)-(1./T)));

options = optimset('fsolve');
options = optimset(options, 'Algorithm', 'trust-region-reflective', 'display', 'none');

xeq = zeros(size(T));
x0 = 0.7;
for K = 1 : numel(KC)
    F = @(xeq) cce(xeq).^2./(cae(xeq).*cbe(xeq))-KC(K);
    x0 = fsolve(F, x0, options);
    xeq(K) = x0;
end

figure
plot(t, z(:,1), t, xeq);
ylim([0,1.2])
xlabel('catalyst weight in kg');
ylabel('conversion X');
legend('X','Xe');
title('converison X and Xe')

function ydot=Zazarian_f_5(~, x, thb)

X=x(1); y=x(2); T=x(3);

% given & calculation

T0=325;
Ta=300;
Fa0=5;
Fb0=10;
Fi0=10;
Ca0=0.2;
Cp_a=20;
Cp_b=20;
Cp_c=20;
Cp_i=18;
alpha=0.00015;
k1=0.0002;
Ua=320;
pb=1400;
e=25000;
mc=18;
thi=2;
dhrx=-20000;
Cp_cool=18;
r=1.987;
kc1=1000;
sumcp=Cp_a+Cp_b*thb+Cp_i*thi;
%calculation

k2=k1*exp((e/r)*((1/300)-(1/T)));

Kc2=kc1*exp((-dhrx/r)*((1/305)-(1/T)));

CA=Ca0*(1-X)*y*(T0/T);

CB=Ca0*(2-X)*y*(T0/T);

CC=2*Ca0*X*y*(T0/T);

Ci=2*Ca0*y*(T0/T);

ra=-k2*((CA*CB)-((CC^2)/Kc2));

ydot(1)=(-ra)/(Fa0);

ydot(2)=(-alpha/(2*y))*(T/T0);

ydot(3)=(-dhrx*ra)/(Fa0*sumcp);

ydot = ydot';
end