ODE23s Solver "Not enough input arguments"

hello everybody,

i adpted this code to my specific problem, but its not working and im running out of ideas how to fix it (unfortunately im rather new to matlab).

I always get the error: Not enough input arguments for my ODE23s solver.

I would be very thankful, if somebody of you could take a look at it and give me some hints.

The problem should be in the integration section of the ODE23s function.

You can more or less skip the first part of the code, because its only defining some constants.

Thank you in advance

global sweep Pp C1 C2 Pt

% constants
r0=2e-5;        % Radius of the tube without membrane, [m]
Ri=0.0225;      % Inner radius of the reactor, [m]
PhiB=0.6;       % Bed porosity
rhocat=2893.4;  % catalyst density [kg m^-3]
Pp=1;           % Permeate pressure [bar]
Pi=3.14159;
R=8.3145;       % [J mol^-1 K^-1]


% kinetic constants
k0=3.3e5;           % Frequency factor [mol(NH3) g(cat)^-1 s^-1]
Ek=107;             % Activation energy [kJ mol^-1]
K1=1.2;             % Ammonia adsorption constant [m^3 mol^-1]
K4=4.5e-2;          % Hydrogen desorption constant [mol m^-3]

% stochiometric coefficients
nu=[-1, 1/2, 3/2];     % NH3, N2, H2

% reactor & membrane
L=0.28;         % Reactor Length, [m]
delta=6e-6;     % membrane thickness [m]
rit=r0+delta;   % inner tube radius [m]

% inlet
Ptot=[2;5;10;20;30];
kmax=length(Ptot);
%Pt=30;         % Total pressure [bar]
Fa0=1;          % NH3 initial feed, [kmol h^-1]
sweep=10;       % sweep factor F0(sweep)/F0(NH3)
T0=350+273.15;  % Inlet temperature [K]
theta=0.001;    % F0(H2)/F0(NH3) [just mathematically needed]

nmax=401;
X1=zeros(nmax,kmax);X2=zeros(nmax,kmax);
X3=zeros(nmax,kmax);

for k=1:kmax
    Pt=Ptot(k);                     % If Pt is the parameter/Matrix
    %Pt=Pt0+(k-1)*1;                % If Pt is an axis
    %L=Lsize(k);                    % If L is the parameter/Matrix
    
    for n=1:nmax
        %delta=delta0+(n-1)*1e-6;   % If delta is the parameter/Matrix
        %rit=r0+delta;              % inner tube radius [m]
        %sweep=s0+(n-1);            % If sweep is the parameter/Matrix
        %T=T0;                       % If T is a constant
        T=T0+(n-1)*1;
        
        % kinetic values
        kr=k0*exp(-Ek/(R*10e-3*T));
        
% other parameters
Ac=Pi*(Ri^2-rit^2);                     % area, [m2]
Qpd=6.3227*(1e-3)*exp(-15630/(R*T));    % H2 permeability in Pd, [m3.um.h.atm0.5/m2] PRÜFEN!
rhoB=rhocat*(1-PhiB);                   % Bed density [kg/m3]
C1=rhoB*Ac*L/Fa0;                       % Constant for diff. eq. 1 and 2, PBR and MR
C2=Qpd*2*Pi*rit*L/(Fa0*delta*22.4);     % Constant for diff. eq. 3, MR

% INTEGRATION
u0=[0 0];                               % initial values
tspan=[0 1];                            % precision
[t,S1]=ode23s(f2,tspan,u0);             % PBR reactor
X1(n,k)=real(S1(end,1));                % Ammonia finale conversion
u1(:,k)=S1(:,1);                        % Ammonia conversion along ksi

u0=[0 0 0];
[t,S2]=ode23s(f3,tspan,u0);             % MR reactor
X2(n,k)=real(S2(end,1));                % Ammonia finale conversion
X3(n,k)=real(S2(end,3));                % hydrogen 
% finale conversion
u1(:,k)=S2(:,1);                        % Ammonia conversion along ksi
    
    
    end
end

Error using /
Matrix dimensions must agree.

Error in solution>f2 (line 95)
S=Pt/(1+theta);

T0=T0-273.15;T=T-273.15;    % Temperature expressed in Celsius
Tspan=(T0:1:T);
Pspan=(Pt0:1:Pt);
%Dspan=(delta0*1e6:1:delta*1e6);
%Sspan=(s0:1:sweep);
D(:,:)=X2(:,:)-X1(:,:);     % Delta = XNH3(MR) - XNH3(PBR)

%%

% Packed bed reactor function, Isothermal
% the variable Y=u(3) is not considered
function diff=f2(t,u)
global C1 Pt theta

% Expression of the partial pressures, functions of ammonia conversion, theta and Pt.
S=Pt/(1+theta);
P(1)=(1-u(1))*S;
P(2)=(1/2)*u(1)*S;
P(3)=(theta+(3/2)*u(1))*S;

% Expression of the reaction rates, functions of Kinetic coefficients, partial pressures and conversions.
r=(k*K1*(P(1)/(R*T)))/((1+K1*(P(1)/(R*T))+sqrt(P(3)/(R*T*K4))^2));

% Differential System
diff1=C1*r;
diff=(diff1)';

end

%%

% Membrane reactor function, Isothermal
function diff=f3(t,u)
global C1 C2 Pt Pp sweep theta

% Expression of the partial pressures, functions of ammonia conversion, hydrogen permeability, theta and Pt.
S=Pt/(1+theta);
P(1)=(1-u(1))*Pt;
P(2)=(1/2)*u(1)*Pt;
P(3)=(theta+(3/2)*u(1)-u(3))*Pt;

% Expression of the reaction rates, functions of Kinetic coefficients, partial pressures and conversions.
r=(k*K1*(P(1)/(R*T)))/((1+K1*(P(1)/(R*T))+sqrt(P(3)/(R*T*K4))^2));

% Pressure of hydrogen in the permeate side
Pmem=Pp*u(3)/(u(3)+sweep);

% Differential system
diff1=C1*r;
diff2=C2*(P(3)^0.5-Pmem^0.5);
diff=[diff1 diff2]';

end

Answer 1

[William Rose]

@ephoi,

Your code:

global sweep Pp C1 C2 Pt theta
% constants
r0=2e-5;        % Radius of the tube without membrane, [m]
Ri=0.0225;      % Inner radius of the reactor, [m]
PhiB=0.6;       % Bed porosity
rhocat=2893.4;  % catalyst density [kg m^-3]
Pp=1;           % Permeate pressure [bar]
Pi=3.14159;
R=8.3145;       % [J mol^-1 K^-1]
% kinetic constants
k0=3.3e5;           % Frequency factor [mol(NH3) g(cat)^-1 s^-1]
Ek=107;             % Activation energy [kJ mol^-1]
K1=1.2;             % Ammonia adsorption constant [m^3 mol^-1]
K4=4.5e-2;          % Hydrogen desorption constant [mol m^-3]
% stochiometric coefficients
nu=[-1, 1/2, 3/2];     % NH3, N2, H2
% reactor & membrane
L=0.28;         % Reactor Length, [m]
delta=6e-6;     % membrane thickness [m]
rit=r0+delta;   % inner tube radius [m]
% inlet
Ptot=[2;5;10;20;30];
kmax=length(Ptot);
%Pt=30;         % Total pressure [bar]
Fa0=1;          % NH3 initial feed, [kmol h^-1]
sweep=10;       % sweep factor F0(sweep)/F0(NH3)
T0=350+273.15;  % Inlet temperature [K]
theta=0.001;    % F0(H2)/F0(NH3) [just mathematically needed]
nmax=401;
X1=zeros(nmax,kmax);X2=zeros(nmax,kmax);
X3=zeros(nmax,kmax);
for k=1:kmax
    Pt=Ptot(k);                     % If Pt is the parameter/Matrix
    %Pt=Pt0+(k-1)*1;                % If Pt is an axis
    %L=Lsize(k);                    % If L is the parameter/Matrix
    
    for n=1:nmax
        %delta=delta0+(n-1)*1e-6;   % If delta is the parameter/Matrix
        %rit=r0+delta;              % inner tube radius [m]
        %sweep=s0+(n-1);            % If sweep is the parameter/Matrix
        %T=T0;                       % If T is a constant
        T=T0+(n-1)*1;
        
        % kinetic values
        kr=k0*exp(-Ek/(R*10e-3*T));
        
        % other parameters
        Ac=Pi*(Ri^2-rit^2);                     % area, [m2]
        Qpd=6.3227*(1e-3)*exp(-15630/(R*T));    % H2 permeability in Pd, [m3.um.h.atm0.5/m2] PRÜFEN!
        rhoB=rhocat*(1-PhiB);                   % Bed density [kg/m3]
        C1=rhoB*Ac*L/Fa0;                       % Constant for diff. eq. 1 and 2, PBR and MR
        C2=Qpd*2*Pi*rit*L/(Fa0*delta*22.4);     % Constant for diff. eq. 3, MR
        % INTEGRATION
        u0=[0 0];                               % initial values
        tspan=[0 1];                            % precision
        [t,S1]=ode23s(f2,tspan,u0);             % PBR reactor
        X1(n,k)=real(S1(end,1));                % Ammonia finale conversion
        u1(:,k)=S1(:,1);                        % Ammonia conversion along ksi
        u0=[0 0 0];
        [t,S2]=ode23s(f3,tspan,u0);             % MR reactor
        X2(n,k)=real(S2(end,1));                % Ammonia finale conversion
        X3(n,k)=real(S2(end,3));                % hydrogen 
        % finale conversion
        u1(:,k)=S2(:,1);                        % Ammonia conversion along ksi    
    end
end

Not enough input arguments.

Error in solution>f2 (line 75)
P(1)=(1-u(1))*S;

% Packed bed reactor function, Isothermal
% the variable Y=u(3) is not considered
function diff=f2(t,u)
global C1 Pt theta
% Expression of the partial pressures, functions of ammonia conversion, theta and Pt.
S=Pt/(1+theta);
P(1)=(1-u(1))*S;
P(2)=(1/2)*u(1)*S;
P(3)=(theta+(3/2)*u(1))*S;
% Expression of the reaction rates, functions of Kinetic coefficients, partial pressures and conversions.
r=(k*K1*(P(1)/(R*T)))/((1+K1*(P(1)/(R*T))+sqrt(P(3)/(R*T*K4))^2));
% Differential System
diff1=C1*r;
diff=(diff1)';
end

The error displayed is not the error you report. The error displayed happens because you did not declare theta global on line 1. I fixed the error shown in your post by adding theta to the list of global variables on line 1.

Then the "not enough input arguments" error occurs. What is the intended size of u? You define u0 to have size 1x2. However, diff, the output from f2, appears to be a scalar. Also, you define P(2) in function f2, but you never use P(2). Why not? There may be issues in function f2 that need attention.