Error when using ode15

Question

Can Cakiroglu 2019-5-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/463712-error-when-using-ode15

评论： Can Cakiroglu 2019-5-24

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.

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

请先登录，再进行评论。

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

Answer 1

David Wilson 2019-5-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/463712-error-when-using-ode15#answer_376453

在 MATLAB Online 中打开

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Can Cakiroglu 2019-5-24

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. :)

请先登录，再进行评论。