ode45 function and solver errors for reactant conversion and temperature change
조회 수: 3 (최근 30일)
이전 댓글 표시
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
채택된 답변
Sam Chak
2023년 11월 18일
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월 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
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
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!