Error using Odearguments - vector must return a column

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')

Answer 1

[Star Strider]

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.

.

[Anh Bui]

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]

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>

<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]

Thank you so much !!!

[Star Strider]

As always, my pleasure!

Answer 2

[Torsten]

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:

eps

ans = 2.2204e-16