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

Answer 1

[Sam Chak]

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

[Maureen]

Thank you so much! this helped a lot

Answer 2

[Walter Roberson]

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

[Walter Roberson]

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?

[Maureen]

Sorry I meant to reply earlier. So i have two state variables in line 13 in my solver file: xO = [FO TO]; (please correct me if I am using the terms wrong), and I am trying to have three outputs for the first one, hence in my function file i have in lines 22-24, these define FO:

dF(1) = -r1;
dF(2) = -r1;
dF(3) = r1;

and the second state variable is meant to have only one output, hence in my function file I have in line 25:

dF(4) = (-r1*DeltaH_RX)/((F(1:end-1))'*CpO); %DT calc

so technically i have x0 has a vector length of 4?

[Walter Roberson]

You must have the same number of state inputs and outputs. (You are not required to use all of the inputs in your calculations)

Answer 3

[Pratyush]

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.

[Maureen]

This did not work unfortunately

[Walter Roberson]

using the ./ operator would return a vector of values, which will not fit in the scalar output location

[Walter Roberson]

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