ode45 function and solver errors for reactant conversion and temperature change
显示 更早的评论
I made a MATLAB code based on an example that is related to the problem I am doing, but I had to tweak it a little due to my problem being slightly different. I am still new to Matlab, so I apologize in advanced if I missed something obvious. I got the errors
Error in ch11solver (line 13)
y_p1 = F(3)/FT;
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);
Error in ch11pt1 (line 11)
[V,x] = ode45(@ch11solver,V,xO); %integration
In my solver file:
clc;
clear all;
CaO = 0.1;
VO = 2;
FaO = CaO*VO;
TO = 300;
xO = [FaO TO];
V = [0, 100]; %reactor vol
[V,x] = ode45(@ch11func,V,xO); %integration
figure(1)
plot(V,x(:,1:end-1),'.'); %individual flow rates kg/s
xlabel('reactor volume (m^3)');
ylabel('flow rate (kmol/s)');
figure(2)
plot(V,x(:,end),','); %reactor temp in K
xlabel('reactor volume (m^3)');
ylabel('temperature (K)');
In my function file:
function dF = ch11func(V,F)
dF = zeros(4,1);
T = F(end); %temperature
T1 = 300;
k1 = 0.01;
E = 10000;
R = 1.987;
FT = sum(F(1:end-1)); %calculate total flow rate
y_r1 = F(1)/FT;
y_r2 = F(2)/FT;
y_p1 = F(3)/FT;
CpO = [15 15 30];
HO = [-20 -15 -41];
DeltaH_RX = [-1 -1 1]*HO; %heat of reaction
k = k1*exp(E/R*((1/T1)-(1/T)));
r1 = k*(y_r1*y_r2);
dF(1) = -r1;
dF(2) = -r1;
dF(3) = r1;
dF(4) = (-r1*DeltaH_RX)/((F(1:end-1))*CpO); %DT calc
采纳的回答
Hi @Maureen
This should get the code running. Please edit the initial values as needed. If you are not very familiar with some special functions or the matrix rules in MATLAB, feel free to use ordinary math notations to describe the differential equations.
CaO = 0.1;
VO = 2;
FaO = CaO*VO;
TO = 300;
xO = [FaO; 1; 0.5; TO]; % initial values
V = [0, 100]; % reactor vol
[V, x] = ode45(@ch11func, V, xO); % integration
figure(1)
plot(V, x(:,1:end-1), '.'), grid on % individual flow rates kg/s
xlabel('reactor volume (m^3)');
ylabel('flow rate (kmol/s)');
figure(2)
plot(V, x(:,end), '.'), grid on % reactor temp in K
xlabel('reactor volume (m^3)');
ylabel('temperature (K)');
function dF = ch11func(V, F)
dF = zeros(4,1);
T = F(4); % temperature
T1 = 300;
k1 = 0.01;
E = 10000;
R = 1.987;
FT = F(1) + F(2) + F(3); % sum(F(1:end-1)); %calculate total flow rate
y_r1 = F(1)/FT;
y_r2 = F(2)/FT;
y_p1 = F(3)/FT;
CpO = [ 15 15 30]';
HO = [-20 -15 -41]';
DeltaH_RX = [ -1 -1 1]*HO; % heat of reaction
k = k1*exp(E/R*(1/T1 - 1/T));
r1 = k*(y_r1*y_r2);
dF(1) = -r1;
dF(2) = -r1;
dF(3) = r1;
dF(4) = - (r1*DeltaH_RX)/([F(1) F(2) F(3)]*CpO); % DT calc
end
更多回答（2 个）
Walter Roberson
2023-11-17
xO = [FaO TO];
those are both scalars so x0 is a vector of length 2. You have two state variables.
dF = zeros(4,1);
but you are returning four state variables
y_p1 = F(3)/FT;
and using 3 state variables on input
3 个评论
Walter Roberson
2023-11-17
HO = [-20 -15 -41];
DeltaH_RX = [-1 -1 1]*HO; %heat of reaction
You have a 1x3 vector matrix multiplied by a 1x3 vector. That is an error for algebraic matrix multiplication. What size of output were you expecting?
Walter Roberson
2023-11-18
You must have the same number of state inputs and outputs. (You are not required to use all of the inputs in your calculations)
Pratyush
2023-11-17
Hi Maureen,
I understand that you are getting error in your MATLAB script, this could be because of last line in your function file.
You are trying to calculate dF(4) using F(1:end-1) (which is a vector) divided by CpO (which is also a vector). MATLAB does not support element-wise division of vectors using the "/" operator. You can use the "./" operator to perform element-wise division. Modify last line as follows:
dF(4) = (-r1*DeltaH_RX)./(F(1:end-1).*CpO); % DT calc
I hope this should resolve the error.
3 个评论
Maureen
2023-11-17
Walter Roberson
2023-11-17
using the ./ operator would return a vector of values, which will not fit in the scalar output location
Walter Roberson
2023-11-17
Depending how some of the other code gets repaired, the line might be intended as a 3x1 vector / a 3x1 vector. If so then the / operation should return a scalar result of least-square fitting
类别
在 帮助中心 和 File Exchange 中查找有关 Chemistry 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
