Use fsolve to solve matrix

2 次查看（过去 30 天）
Rodrigo Blas
Rodrigo Blas 2020-5-8
评论： Walter Roberson 2020-5-8
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);
>>

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2020-5-8
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.
  2 个评论
Rodrigo Blas
Rodrigo Blas 2020-5-8
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
Walter Roberson 2020-5-8
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Argument Definitions 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by