Error when using ode15

Hello,

I'm a chemical engineering student, so these equations might come difficult to analyze. However, I am receiving following error on my equations

This is the main script

close all
clear
clc

global Fao Tr z a b1 b2 R K30 E2 E3 omega alpha k20 Bo rho_c Phi Ac T0 P0

Fao = 1487.351; % kmol/h


Tr = 273; % K

z = 1;

a = (3 * 1.5) / ((3 * 1.5) + 1);
b1 = (0.5 + (0.5 / 1.5)) / 1;
b2 = (-0.5+(1.5*1.5))/1;

R = 8.314 / 1000; % kJ/mol*K

K30 = 0.0307;
E2 = 72100.82; % kJ/mol
E3 = -80989.98; % kJ/mol

omega = 1.564;
alpha = 0.640;
k20 = 2.12 * (10 ^ 13);

Bo = 0.254626; %atm/m
rho_c = 1600; %kg/m^3
Phi = 0.5;
D = 0.05; %Diameter of tube, m
Ac = pi * (D^2)/4; %m^2

X0 = 0;
T0 = 673; % K
P0 = 150; %atm

x0 = [X0 T0 P0];
Wspan = [0 5];
[W, x] = ode15s(@fun, Wspan, x0);

And this is the function used

function func = fun(x)
%% Identification of global variables
global Fao Tr z a b1 b2 R K30 E2 E3 omega alpha k20 Bo rho_c Phi Ac T0 P0

%% Decleration of the non-linear equations
func = [((k20*(-E2/R*x(2))*(x(3)*(a/3)*(1-b2*z)*((10^(-5.519265*(10^(-5))*x(2)+...
    1.848863*(10^(-7))*(x(2)^2)+(2001.6/x(2))+2.6899)/(x(2)^(5.519265)))^2)-...
    (((x(3)*z)^2)/((x(3)*a*(1-b1*z))^3))))/((1+(K30*exp(-E3/R*x(2)))*(x(3)*z)/...
    ((x(3)*a*(1-b1*z))^omega))^(2*alpha)))/Fao; % dX/dW, x(1) = X
    
    ((2.12*(10^13)*(-E2/R*x(2))*(x(3)*(a/3)*(1-b2*z)*((10^(-5.519265*(10^(-5))*x(2)...
    + 1.848863*(10^(-7))*(x(2)^2)+(2001.6/x(2))+2.6899)/(x(2)^(5.519265)))^2)-...
    (((x(3)*z)^2)/((x(3)*a*(1-b1*z))^3))))/((1+(K30*exp(-E3/R*x(2)))*(x(3)*z)/...
    ((x(3)*a*(1-b1*z))^omega))^(2*alpha)))/(-(-91.63+(-12.96*(x(2)-Tr)+0.00917*((x(2)-Tr)^2)/2)))/...
    (Fao*(6.50+0.00100*x(2))+x(1)*2*(6.70+ 0.00630*x(2))-...
    3*(6.62+0.00081*x(2))-(6.50+0.00100*x(2))); % dT/dW, x(2) = T
    
    -(Bo/(Ac*(1-Phi)*rho_c))*(P0/x(3))*(x(2)/T0)*(1-0.5*x(1)) %dp/dW, x(3) = P
    ];


end

and I receive following error when i start the script:

Error in ode15s (line 150)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Reactor_Project (line 40)
[W, x] = ode15s(@fun, Wspan, x0);

I would be very happy if someone could identify this error. I thank everybody who spends time on this matter.

Answer 1

[David Wilson]

It is good practice to try and avoid globals, but rather pass them as auxilary variables into the function.

Fao = 1487.351; % kmol/h
Tr = 273; % K
z = 1;
a = (3 * 1.5) / ((3 * 1.5) + 1);
b1 = (0.5 + (0.5 / 1.5)) / 1;
b2 = (-0.5+(1.5*1.5))/1;
R = 8.314 / 1000; % kJ/mol*K
K30 = 0.0307;
E2 = 72100.82; % kJ/mol
E3 = -80989.98; % kJ/mol
omega = 1.564;
alpha = 0.640;
k20 = 2.12 * (10 ^ 13);  % should be 2.12e13 I suspect 
Bo = 0.254626; %atm/m
rho_c = 1600; %kg/m^3
Phi = 0.5;
D = 0.05; %Diameter of tube, m
Ac = pi * (D^2)/4; %m^2
X0 = 0;
T0 = 673; % K
P0 = 150; %atm
x0 = [X0 T0 P0];
Wspan = [0 5];
[W, x] = ode15s(@(t,x) fchem(t,x,Fao, Tr, z, a, b1, b2, R, K30, E2, E3, omega, alpha, k20, Bo, rho_c, Phi, Ac, T0, P0) ... 
           , Wspan, x0);
       
plot(W,x)      

Now you function is (which I've renamed to make it more user-friendly)

 function func = fchem(t,x,Fao, Tr, z, a, b1, b2, R, K30, E2, E3, omega, alpha, k20, Bo, rho_c, Phi, Ac, T0, P0)
%% Identification of global variables
%global Fao Tr z a b1 b2 R K30 E2 E3 omega alpha k20 Bo rho_c Phi Ac T0 P0
%% Decleration of the non-linear equations
func = [((k20*(-E2/R*x(2))*(x(3)*(a/3)*(1-b2*z)*((10^(-5.519265*(10^(-5))*x(2)+...
    1.848863*(10^(-7))*(x(2)^2)+(2001.6/x(2))+2.6899)/(x(2)^(5.519265)))^2)-...
    (((x(3)*z)^2)/((x(3)*a*(1-b1*z))^3))))/((1+(K30*exp(-E3/R*x(2)))*(x(3)*z)/...
    ((x(3)*a*(1-b1*z))^omega))^(2*alpha)))/Fao; % dX/dW, x(1) = X
    
    ((2.12*(10^13)*(-E2/R*x(2))*(x(3)*(a/3)*(1-b2*z)*((10^(-5.519265*(10^(-5))*x(2)...
    + 1.848863*(10^(-7))*(x(2)^2)+(2001.6/x(2))+2.6899)/(x(2)^(5.519265)))^2)-...
    (((x(3)*z)^2)/((x(3)*a*(1-b1*z))^3))))/((1+(K30*exp(-E3/R*x(2)))*(x(3)*z)/...
    ((x(3)*a*(1-b1*z))^omega))^(2*alpha)))/(-(-91.63+(-12.96*(x(2)-Tr)+0.00917*((x(2)-Tr)^2)/2)))/...
    (Fao*(6.50+0.00100*x(2))+x(1)*2*(6.70+ 0.00630*x(2))-...
    3*(6.62+0.00081*x(2))-(6.50+0.00100*x(2))); % dT/dW, x(2) = T
    
    -(Bo/(Ac*(1-Phi)*rho_c))*(P0/x(3))*(x(2)/T0)*(1-0.5*x(1)) %dp/dW, x(3) = P
    ];
end

When I plot the solution, it looks pretty boring, but that's a different problem.

[Can Cakiroglu]

Thank you very much for your hard work. The pretty boring part is indeed so, but hopefuly will be fixed since this is our final project calculations. :)