Error when using ode15
3 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
David Wilson
2019-5-24
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!