ode45 function and solver errors for reactant conversion and temperature change
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponse acceptée
Sam Chak
le 18 Nov 2023
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
Plus de réponses (2)
Walter Roberson
le 17 Nov 2023
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 commentaires
Walter Roberson
le 18 Nov 2023
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
le 17 Nov 2023
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 commentaires
Walter Roberson
le 17 Nov 2023
using the ./ operator would return a vector of values, which will not fit in the scalar output location
Walter Roberson
le 17 Nov 2023
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
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!