Error using Odearguments - vector must return a column

Question

Anh Bui 2024-4-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2108256-error-using-odearguments-vector-must-return-a-column

评论： Star Strider 2024-4-17

I kept on trying to reword and rearrange my code however it asked me to do it so that the line return a column vector. I tried reshape(), [],1, (does not work). Also tried .' at the end which also doesnt work. I'm not sure where I can continue with this?

%r_H2= -(-Rds - Rwf - Rh);
%Rd = kds*(CH20*((theta_O2-X)/(1+eps*X)))*(CH20*((1-X)/(1+eps*X)))^2
%Rwf = kwf* (CH20*((theta_O2-(b/a)*X)/(1+wf_eps*X)))^.5*CH20.*((1-X)./(1+wf_eps*X))
%Rh = kh* (CH2O2*(theta_H2-(d_b/d_a)*X)) *CH2O2*(1-X)
wspan = [0 1000];
dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + ...
    kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2.*(theta_H2-(d_b/d_a).*X)) .*CH2O2.*(1-X));
IC = [0 0 0 0 0 0 0 0 0 0 0];
[W,X] = ode45(dxdw,wspan,IC);
plot(W,X)
title ('Catalyst v.s Conversion')
xlabel('Weight of Catalyst (kg)')
ylabel('Conversion of H2')

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Star Strider 2024-4-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2108256-error-using-odearguments-vector-must-return-a-column#answer_1443241

在 MATLAB Online 中打开

You only have one differential equation, so only one initial condition is necessary.

%r_H2= -(-Rds - Rwf - Rh);
%Rd = kds*(CH20*((theta_O2-X)/(1+eps*X)))*(CH20*((1-X)/(1+eps*X)))^2
%Rwf = kwf* (CH20*((theta_O2-(b/a)*X)/(1+wf_eps*X)))^.5*CH20.*((1-X)./(1+wf_eps*X))
%Rh = kh* (CH2O2*(theta_H2-(d_b/d_a)*X)) *CH2O2*(1-X)
wspan = [0 1000];
dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + ...
    kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2.*(theta_H2-(d_b/d_a).*X)) .*CH2O2.*(1-X));
% IC = [0 0 0 0 0 0 0 0 0 0 0];
IC = eps;
[W,X] = ode45(dxdw,wspan,IC);
Unrecognized function or variable 'kds'.

Error in solution>@(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2+kwf.*(CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+kh.*(CH2O2.*(theta_H2-(d_b/d_a).*X)).*CH2O2.*(1-X)) (line 6)
dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + ...

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 104)
    odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
plot(W,X)
title ('Catalyst v.s Conversion')
xlabel('Weight of Catalyst (kg)')
ylabel('Conversion of H2')

All the missing constants prevent me from running this.

.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Anh Bui 2024-4-17

在 MATLAB Online 中打开

my X variable is a vector [0:.1:1];

hence the 0s, and it needed to come back as a vector still.

Reaction Kinetics

%% Temperature range
    T = 313; % Kelvin
    %T = 40+273; %K
%% Pressure 
    % P = [20 30 50]; % Bar
    % P = P.*(10^5); % Pa
    P = 5.1*(10^6); %Pa
%% Constants (A, Ea, and R)
    A_ds = 4.5*(10^7); 
    A_wf = 6.5*(10^7);
    A_d = 1.005*(10^7);
    A_h = 280000;
    Ea_ds = 6.9*(10^3); % J/mol
    Ea_wf = 20*(10^3); % J/mol
    Ea_d = 53*(10^3); % J/mol
    Ea_h = 44.3*(10^3); % J/mol
    R = 8.314; % J/molK or Pa*m3/molK
%% Solve for kinetic constants (ki)
    kds = A_ds.*exp((-Ea_ds)./(R.*T));
    kwf = A_wf.*exp((-Ea_wf)./(R.*T));
    kd = A_d.*exp((-Ea_d)./(R.*T));
    kh = A_h.*exp((-Ea_h)./(R.*T));

Direct Synthesis

%initial molar flow mol/s
F_H2_in  = 15; 
F_O2_in = 87.85;
F_CO2_in = 278.82;
yH2o = .04;
theta_O2= F_O2_in/F_H2_in;
theta_CO2 = F_CO2_in/F_H2_in;
theta_H2O2 = 0;
delta = -1;
eps = yH2o*delta;
%Temperature change
T0T = 1;
%Pressure Drop
PP0 = 1;
%initial mass flow g/s
M_H2 = F_H2_in/1.02;
M_O2 = F_O2_in/16;
M_CO2 = F_CO2_in/44;
%density of gas mixtures g/L
rho_mix = 77.555;
%volumetric flow L/s
v_H2 = M_H2/rho_mix;
v_O2 = M_O2/rho_mix;
v_CO2 = M_CO2/rho_mix;
v0 = v_CO2+v_O2+v_H2;
%initial concentration of H2 mol/L
CH20 = F_H2_in/v0;
%mass balance for direct synthesis
X = [0:0.1:1];
CH2 = CH20.*((1-X)./(1+eps.*X))*T0T*PP0;
CO2 = CH20.*((theta_O2-X)./(1+eps*X))*T0T*PP0;
CH2O2 = CH20.*((theta_H2O2+X)./(1+eps*X))*T0T*PP0;
Rds = kds.*CO2.*CH2.^2;
plot (X,CH2,X,CO2,X,CH2O2)
legend('C_H2','C_O2','C_H2O2',Location='northwest')
title('Direct Synthesis Conversion')
xlabel('Conversion')
ylabel('Concentration (mol/L)')
%rate of direct synthesis

Water Formation

%2H2+O2->2H2O 
%coefficient of the chemical equation
a = 2; %H2 coeff
b = 1; %O2 coeff
c = 2;%water coeff
wf_yH20 = .04;
wf_delta = -a-(b/a)+(c/a);
wf_eps = wf_yH20*delta;
wf_C_H2 = CH20.*((1-X)./(1+wf_eps*X))*T0T*PP0;
wf_C_O2 = CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))*T0T*PP0;
wf_C_H2O = CH20.*((0+(c/a).*X)/(1+wf_eps.*X))*T0T*PP0;
wf_C_CO2 =((CH20*theta_CO2)./(1+eps*X))*T0T*PP0;
Rwf = kwf.*(wf_C_O2.^.5).*wf_C_H2;
plot (X,wf_C_H2,X,wf_C_CO2,X,wf_C_H2O,X,wf_C_CO2)
legend('C_H2','C_O2','C_H2O2','C_CO2',Location='northwest')
title('Water Formation Conversion')
xlabel('Conversion')
ylabel('Concentration (kmol/L)')

Decomposition

%H2O2(l)-> .5 O2 (g) + H2O(l)
%coeff of the chemical equation
d_a = 1;%H2O2 coeff
d_b = .5;%O2 coeff
d_c = 1; %water coeff
d_C_H2O2= CH2O2.*(1-X);
d_C_O2 = CH2O2.*(0+(d_b/d_a).*X);
d_C_H2O = CH2O2.*(0+(d_c/d_a).*X);
Rd= kd.*(d_C_H2O2);
plot (X,d_C_H2O2,X,d_C_O2,X,d_C_H2O)
legend('C_H2O2','C_O2','C_H2O',Location='northwest')
title('Decomposition Conversion')
xlabel('Conversion')
ylabel('Concentration (kmol/L)')

Hydrogenation

%H2O2(l)+H2(g)-> 2H2O(l)
%coeff of the chemical equation
h_a = 1;%coeff of h2o2
h_b = 1; %coeff of h2
h_c = 2; %coeff of water
%molar flow rate kmol/hr
F_H2O2 = 33.95;
F_H2 = 10.72;
h_yH20 = .04;
h_delta = 0;
h_eps = h_yH20*h_delta;
theta_H2  = F_H2/F_H2O2;
h_C_H2O2 = CH2O2.*(1-X);
h_C_H2 = CH2O2.*(theta_H2-(h_b/h_a).*X);
h_C_H2O = CH2O2.*(0+(h_c/h_a).*X);
Rh = kh.*h_C_H2.*h_C_H2O2;
plot (X,h_C_H2O2,'-o',X,h_C_H2,'ob',X,h_C_H2O,'-r')
legend('C_H2O2','C_H2','C_H2O',Location='northwest')
title('Hydrogenation Conversion')
xlabel('Conversion')
ylabel('Concentration (kmol/L)')

Rate vs Conversion

r_H2O2 = Rds-Rd-Rh;
r_H2O= Rwf+Rd+2.*Rh;
r_O2 = -(-Rds-0.5.*Rwf-0.5.*Rh);
r_H2= -(-Rds - Rwf - Rh);
plot(X,r_H2,'-o',X,r_H2O,'ob',X,r_O2,'b',X,r_H2O2,'r')
legend('Rate of H2','Rate of water','Rate of O2','Rate of H2O2')
title('Rate vs Conversion')
xlabel('Conversion')
ylabel('Rate of Components')

Finding weight of Catalyst From H2

%r_H2= -(-Rds - Rwf - Rh);
%Rd = kds*(CH20*((theta_O2-X)/(1+eps*X)))*(CH20*((1-X)/(1+eps*X)))^2
%Rwf = kwf* (CH20*((theta_O2-(b/a)*X)/(1+wf_eps*X)))^.5*CH20.*((1-X)./(1+wf_eps*X))
%Rh = kh* (CH2O2*(theta_H2-(d_b/d_a)*X)) *CH2O2*(1-X)
wspan = [0 1000];
dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2.*(theta_H2-(d_b/d_a).*X)) .*CH2O2.*(1-X));
IC = eps;
[W,X] = ode45(dxdw,wspan,IC);
plot(W,X)
title ('Catalyst v.s Conversion')
xlabel('Weight of Catalyst (kg)')
ylabel('Conversion of H2')

Star Strider 2024-4-17

在 MATLAB Online 中打开

Your ‘dxdw’ function produces an (11x11) matrix of identical (matching) columns in each iteration. That was the reason for the ‘column vector’ error. I wrote a ‘helper’ function, ‘dxdwc’ that returns only the first column of the matrix, and changed the solver to ode15s because the system takes a while to integrate, and like every kinetic equation I’ve ever encountered, is ‘stiff’ because of the wide range of parameter magnitudes, so this is appropriate.

There may still be problems. (I have no idea what the result is supposed to look like, however it makes sense that the more catalyst added the more quickly the system reaches equilibrium.)

Try this —

%% Temperature range

T = 313; % Kelvin

%T = 40+273; %K

%% Pressure

% P = [20 30 50]; % Bar

% P = P.*(10^5); % Pa

P = 5.1*(10^6); %Pa

%% Constants (A, Ea, and R)

A_ds = 4.5*(10^7);

A_wf = 6.5*(10^7);

A_d = 1.005*(10^7);

A_h = 280000;

Ea_ds = 6.9*(10^3); % J/mol

Ea_wf = 20*(10^3); % J/mol

Ea_d = 53*(10^3); % J/mol

Ea_h = 44.3*(10^3); % J/mol

R = 8.314; % J/molK or Pa*m3/molK

%% Solve for kinetic constants (ki)

kds = A_ds.*exp((-Ea_ds)./(R.*T));

kwf = A_wf.*exp((-Ea_wf)./(R.*T));

kd = A_d.*exp((-Ea_d)./(R.*T));

kh = A_h.*exp((-Ea_h)./(R.*T));

%initial molar flow mol/s

F_H2_in = 15;

F_O2_in = 87.85;

F_CO2_in = 278.82;

yH2o = .04;

theta_O2= F_O2_in/F_H2_in;

theta_CO2 = F_CO2_in/F_H2_in;

theta_H2O2 = 0;

delta = -1;

eps = yH2o*delta;

%Temperature change

T0T = 1;

%Pressure Drop

PP0 = 1;

%initial mass flow g/s

M_H2 = F_H2_in/1.02;

M_O2 = F_O2_in/16;

M_CO2 = F_CO2_in/44;

%density of gas mixtures g/L

rho_mix = 77.555;

%volumetric flow L/s

v_H2 = M_H2/rho_mix;

v_O2 = M_O2/rho_mix;

v_CO2 = M_CO2/rho_mix;

v0 = v_CO2+v_O2+v_H2;

%initial concentration of H2 mol/L

CH20 = F_H2_in/v0;

%mass balance for direct synthesis

X = [0:0.1:1];

CH2 = CH20.*((1-X)./(1+eps.*X))*T0T*PP0;

CO2 = CH20.*((theta_O2-X)./(1+eps*X))*T0T*PP0;

CH2O2 = CH20.*((theta_H2O2+X)./(1+eps*X))*T0T*PP0;

Rds = kds.*CO2.*CH2.^2;

plot (X,CH2,X,CO2,X,CH2O2)

legend('C_H2','C_O2','C_H2O2',Location='east')

title('Direct Synthesis Conversion')

xlabel('Conversion')

ylabel('Concentration (mol/L)')

%rate of direct synthesis

%2H2+O2->2H2O

%coefficient of the chemical equation

a = 2; %H2 coeff

b = 1; %O2 coeff

c = 2;%water coeff

wf_yH20 = .04;

wf_delta = -a-(b/a)+(c/a);

wf_eps = wf_yH20*delta;

wf_C_H2 = CH20.*((1-X)./(1+wf_eps*X))*T0T*PP0;

wf_C_O2 = CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))*T0T*PP0;

wf_C_H2O = CH20.*((0+(c/a).*X)/(1+wf_eps.*X))*T0T*PP0;

wf_C_CO2 =((CH20*theta_CO2)./(1+eps*X))*T0T*PP0;

Rwf = kwf.*(wf_C_O2.^.5).*wf_C_H2;

plot (X,wf_C_H2,X,wf_C_CO2,X,wf_C_H2O,X,wf_C_CO2)

legend('C_H2','C_O2','C_H2O2','C_CO2',Location='east')

title('Water Formation Conversion')

xlabel('Conversion')

ylabel('Concentration (kmol/L)')

%H2O2(l)-> .5 O2 (g) + H2O(l)

%coeff of the chemical equation

d_a = 1;%H2O2 coeff

d_b = .5;%O2 coeff

d_c = 1; %water coeff

d_C_H2O2= CH2O2.*(1-X);

d_C_O2 = CH2O2.*(0+(d_b/d_a).*X);

d_C_H2O = CH2O2.*(0+(d_c/d_a).*X);

Rd= kd.*(d_C_H2O2);

plot (X,d_C_H2O2,X,d_C_O2,X,d_C_H2O)

legend('C_H2O2','C_O2','C_H2O',Location='northwest')

title('Decomposition Conversion')

xlabel('Conversion')

ylabel('Concentration (kmol/L)')

%H2O2(l)+H2(g)-> 2H2O(l)

%coeff of the chemical equation

h_a = 1;%coeff of h2o2

h_b = 1; %coeff of h2

h_c = 2; %coeff of water

%molar flow rate kmol/hr

F_H2O2 = 33.95;

F_H2 = 10.72;

h_yH20 = .04;

h_delta = 0;

h_eps = h_yH20*h_delta;

theta_H2 = F_H2/F_H2O2;

h_C_H2O2 = CH2O2.*(1-X);

h_C_H2 = CH2O2.*(theta_H2-(h_b/h_a).*X);

h_C_H2O = CH2O2.*(0+(h_c/h_a).*X);

Rh = kh.*h_C_H2.*h_C_H2O2;

plot (X,h_C_H2O2,'-o',X,h_C_H2,'ob',X,h_C_H2O,'-r')

legend('C_H2O2','C_H2','C_H2O',Location='northwest')

title('Hydrogenation Conversion')

xlabel('Conversion')

ylabel('Concentration (kmol/L)')

r_H2O2 = Rds-Rd-Rh;

r_H2O= Rwf+Rd+2.*Rh;

r_O2 = -(-Rds-0.5.*Rwf-0.5.*Rh);

r_H2= -(-Rds - Rwf - Rh);

plot(X,r_H2,'-o',X,r_H2O,'ob',X,r_O2,'b',X,r_H2O2,'r')

legend('Rate of H2','Rate of water','Rate of O2','Rate of H2O2')

title('Rate vs Conversion')

xlabel('Conversion')

ylabel('Rate of Components')

%r_H2= -(-Rds - Rwf - Rh);

%Rd = kds*(CH20*((theta_O2-X)/(1+eps*X)))*(CH20*((1-X)/(1+eps*X)))^2

%Rwf = kwf* (CH20*((theta_O2-(b/a)*X)/(1+wf_eps*X)))^.5*CH20.*((1-X)./(1+wf_eps*X))

%Rh = kh* (CH2O2*(theta_H2-(d_b/d_a)*X)) *CH2O2*(1-X)

wspan = [0 1000]/250;

% dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2.*(theta_H2-(d_b/d_a).*X)) .*CH2O2.*(1-X));

IC = ones(size(X))*1E-8;

[W,X] = ode15s(@(W,X)dxdwc(W,X, kds,CH20,theta_O2,eps,kwf,b,a,wf_eps,kh,CH2O2,theta_H2,d_b,d_a),wspan,IC);

Out = [W X]

Out = 150x12

0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0008 0.0008 0.0008 0.0008 0.0008 0.0008 0.0008 0.0008 0.0008 0.0008 0.0008 0.0000 0.0016 0.0016 0.0016 0.0016 0.0016 0.0016 0.0016 0.0016 0.0016 0.0016 0.0016 0.0000 0.0024 0.0024 0.0024 0.0024 0.0024 0.0024 0.0024 0.0024 0.0024 0.0024 0.0024 0.0000 0.0071 0.0071 0.0071 0.0071 0.0071 0.0071 0.0071 0.0071 0.0071 0.0071 0.0071 0.0000 0.0118 0.0118 0.0118 0.0118 0.0118 0.0118 0.0118 0.0118 0.0118 0.0118 0.0118 0.0000 0.0164 0.0164 0.0164 0.0164 0.0164 0.0164 0.0164 0.0164 0.0164 0.0164 0.0164 0.0000 0.0210 0.0210 0.0210 0.0210 0.0210 0.0210 0.0210 0.0210 0.0210 0.0210 0.0210 0.0000 0.0450 0.0450 0.0450 0.0450 0.0450 0.0450 0.0450 0.0450 0.0450 0.0450 0.0450 0.0000 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678

<mw-icon class=""></mw-icon>

figure

plot(W,X)

Ax = gca;

Ax.XScale = 'log';

grid

title ('Catalyst v.s Conversion')

xlabel('Weight of Catalyst (kg)')

ylabel('Conversion of H2')

function d = dxdwc(W,X, kds,CH20,theta_O2,eps,kwf,b,a,wf_eps,kh,CH2O2,theta_H2,d_b,d_a)

dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2.*(theta_H2-(d_b/d_a).*X)) .*CH2O2.*(1-X));

% Q = dxdw(W,X)

for k = 1:numel(X)

dv(k,:) = dxdw(W,X(k));

end

d = dv(:,1);

end

The ‘helper’ function ‘dxdwc’ extracts the first column from the matrix created by ‘dxdw’ (all the columns appear to be identical) and returns it to ode15s. It does not otherwise change the code in any way.

There may still be problems with the differential equation. I defer to you to address any problems that remain.

.

Anh Bui 2024-4-17

Thank you so much !!!

Star Strider 2024-4-17

As always, my pleasure!

请先登录，再进行评论。

Answer 2

Torsten 2024-4-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2108256-error-using-odearguments-vector-must-return-a-column#answer_1443301

移动：Torsten 2024-4-17

在 MATLAB Online 中打开

This works - although I'm not sure if it is what you want:

wspan = [0 1000];
CH2O2 = @(X)CH20.*((theta_H2O2+X)./(1+eps*X))*T0T*PP0;
dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2(X).*(theta_H2-(d_b/d_a).*X)) .*CH2O2(X).*(1-X));
IC = eps*ones(11,1);
[W,X] = ode15s(dxdw,wspan,IC);

And better use a variable name different from eps. It usually is a predefined constant is MATLAB: