Error in my ode45 equation.

I am unsure why my code is wrong (pasted below) as using a very similar example this code worked. When running the code below, these errors are displayed:

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in reactorassignmentpart1plot (line 5)

[V,y] = ode45(@fun1, Vspan, yo);

This is the code used:

clc
Vspan = [0 1];
yo = [4 0 0 1 500.15];
[V,y] = ode45(@fun1, Vspan, yo);
 
subplot(2,1,1)
plot(V,y(:,1),V,y(:,2),V,y(:,3));
legend('Fa','Fb','Fc');
ylabel('Molar flowrate, mol/min');
xlabel('Volume,dm3');
subplot(2,1,2)
plot(V,y(:,4));
legend('Temperature (K)');
ylabel('Temperature (K)');
xlabel('Volume,dm3');
 
% Compose the function
function f = fun1(V,Y)
% Define the differential equations that need to be solved 
% Y is the concentration and V is the PFR volume
Fa = Y(1);
Fb = Y(2);
Fc = Y(3);
T = Y(4);
 
% Define initial conditions 
deltaH1 = -25000; %kJ/molA
deltaH2 = 35000; %kJ/molB
deltaH2T = 35000-(80*(T-298)); %kJ/molB
CTo = 0.3996; % mol/L
To = 500.15; % K
Fio = 1; % mol/min
FTo = 4; % mol/min
Cpa = 20; % J/molK
Cpb = 80; % J/molK
Cpc = 100; % J/molK
Cpi = 20; % J/molK
Ua = 150; %J/dm3.min.K
FT = Fa + Fb + Fc + Fio;
Ca = CTo*((Fa/FT)*(To/T));
Cb = CTo*((Fb/FT)*(To/T));
Cc = CTo*((Fc/FT)*(To/T));
Ci = CTo*((Fio/FT)*(To/T));
K1(T) = 50*exp((8000/8.314)*((1/315)-(1/T))); % dm3/mol.min
Kc(T) = 10*exp((-25/8.314)*((1/315)-(1/T))); % dm3/mol.min
K2(T) = 400*exp((4000/8.314)*((1/310)-(1/T))); % dm3/mol.min-1
Ta = 523.15; % K
ra = (K1*((Ca^2)-((1/Kc)*Cb))) + (K2*Ca*(Cb^2));
rb = 1/2*((K1*((Ca^2)-((1/Kc)*Cb)))+(2*K2*(Cb^2)*Ca));
rc = K2*Ca*(Cb^2);
% Differential equations
dFadV = -ra;
dFbdV = -rb;
dFcdV = rc;
dTdV = (Ua*(Ta-T)-((ra-rc)*(deltaH1))-(2*rc*(deltaH2)))/(Cpa*Fa+Cpb*Fb+Cpc*Fc+Cpi*Fio);
f = [dFadV; dFbdV; dFcdV; dTdV];
end

Thank you !

Answer 1

[Walter Roberson]

clc
Vspan = [0 1];
yo = [4 0 0 1 500.15];
[V,y] = ode45(@fun1, Vspan, yo);

Error using odearguments (line 95)
FUN1 returns a vector of length 4, but the length of initial conditions vector is 5. The vector returned by FUN1 and the initial conditions vector must have the same number of elements.

Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Sure enough, you are passing in a yo of length 5,

 
subplot(2,1,1)
plot(V,y(:,1),V,y(:,2),V,y(:,3));
legend('Fa','Fb','Fc');
ylabel('Molar flowrate, mol/min');
xlabel('Volume,dm3');
subplot(2,1,2)
plot(V,y(:,4));
legend('Temperature (K)');
ylabel('Temperature (K)');
xlabel('Volume,dm3');
 
% Compose the function
function f = fun1(V,Y)
% Define the differential equations that need to be solved 
% Y is the concentration and V is the PFR volume
Fa = Y(1);
Fb = Y(2);
Fc = Y(3);
T = Y(4);
 
% Define initial conditions 
deltaH1 = -25000; %kJ/molA
deltaH2 = 35000; %kJ/molB
deltaH2T = 35000-(80*(T-298)); %kJ/molB
CTo = 0.3996; % mol/L
To = 500.15; % K
Fio = 1; % mol/min
FTo = 4; % mol/min
Cpa = 20; % J/molK
Cpb = 80; % J/molK
Cpc = 100; % J/molK
Cpi = 20; % J/molK
Ua = 150; %J/dm3.min.K
FT = Fa + Fb + Fc + Fio;
Ca = CTo*((Fa/FT)*(To/T));
Cb = CTo*((Fb/FT)*(To/T));
Cc = CTo*((Fc/FT)*(To/T));
Ci = CTo*((Fio/FT)*(To/T));
K1(T) = 50*exp((8000/8.314)*((1/315)-(1/T))); % dm3/mol.min
Kc(T) = 10*exp((-25/8.314)*((1/315)-(1/T))); % dm3/mol.min
K2(T) = 400*exp((4000/8.314)*((1/310)-(1/T))); % dm3/mol.min-1
Ta = 523.15; % K
ra = (K1*((Ca^2)-((1/Kc)*Cb))) + (K2*Ca*(Cb^2));
rb = 1/2*((K1*((Ca^2)-((1/Kc)*Cb)))+(2*K2*(Cb^2)*Ca));
rc = K2*Ca*(Cb^2);
% Differential equations
dFadV = -ra;
dFbdV = -rb;
dFcdV = rc;
dTdV = (Ua*(Ta-T)-((ra-rc)*(deltaH1))-(2*rc*(deltaH2)))/(Cpa*Fa+Cpb*Fb+Cpc*Fc+Cpi*Fio);
f = [dFadV; dFbdV; dFcdV; dTdV];

but your f is only putting together 4 values.

end

Every yo input value is a boundary condition. You have to return the derivative for each of the entries.

You do not seem to use Y(5) in your fun1, so perhaps you should only be passing in 4 initial values instead of 5.

[Zara Freeman]

Thank you

Answer 2

[Mischa Kim]

Hi Zara,

your vector of initial conditions has 5 components

yo = [4 0 0 1 500.15];

however, you only have 4 differential equations:

f = [dFadV; dFbdV; dFcdV; dTdV];

The numbers need to match, so, you are either missing a differential equation or your yo vector is too long.

Also, in your ode file you probably want to change K1(t), Kc(t), and K2(t) to K1, Kc, and K2.

[Zara Freeman]

Hi Thank you for your help! I have made the relevant corrections, however this error is now displayed:

Warning: Failure at t=1.109564e-01. Unable to meet integration tolerances without reducing the step

size below the smallest value allowed (2.220446e-16) at time t.

> In ode45 (line 360)

In reactorassignmentpart1plot (line 5)

%% Chemical Systems Engineering Reactor Assignment 
clc
Vspan = [0:0.1:1];
yo = [4 0 0 1 500.15];
[V,y] = ode45(@fun1, Vspan, yo);
 
subplot(2,1,1)
plot(V,y(:,1),V,y(:,2),V,y(:,3));
legend('Fa','Fb','Fc');
ylabel('Molar flowrate, mol/min');
xlabel('Volume,dm3');
subplot(2,1,2)
plot(V,y(:,4));
legend('Temperature (K)');
ylabel('Temperature (K)');
xlabel('Volume,dm3');
 
% Compose the function
function f = fun1(V,Y)
% Define the differential equations that need to be solved 
% Y is the concentration and V is the PFR volume
Fa = Y(1);
Fb = Y(2);
Fc = Y(3);
Fi = Y(4);
T = Y(5);
 
% Define initial conditions 
deltaH1 = -25000; %kJ/molA
deltaH2 = 35000; %kJ/molB
deltaH2T = 35000-(80*(T-298)); %kJ/molB
CTo = 0.3996; % mol/L
To = 500.15; % K
Fio = 1; % mol/min
FTo = 4; % mol/min
Cpa = 20; % J/molK
Cpb = 80; % J/molK
Cpc = 100; % J/molK
Cpi = 20; % J/molK
Ua = 150; %J/dm3.min.K
FT = Fa + Fb + Fc + Fio;
Ca = CTo*((Fa/FT)*(To/T));
Cb = CTo*((Fb/FT)*(To/T));
Cc = CTo*((Fc/FT)*(To/T));
Ci= CTo*((Fio/FT)*(To/T));
K1 = 50*exp((8000/8.314)*((1/315)-(1/T))); % dm3/mol.min
Kc = 10*exp((-25/8.314)*((1/315)-(1/T))); % dm3/mol.min
K2 = 400*exp((4000/8.314)*((1/310)-(1/T))); % dm3/mol.min
Ta = 523.15; % K
ra = (K1*((Ca^2)-((1/Kc)*Cb))) + (K2*Ca*(Cb^2));
rb = 1/2*((K1*((Ca^2)-((1/Kc)*Cb)))+(2*K2*(Cb^2)*Ca));
rc = K2*Ca*(Cb^2);
% Differential equations
dFadV = -ra;
dFbdV = -rb;
dFcdV = rc;
dFidV = 0;
dTdV = (Ua*(Ta-T)-((ra-rc)*(deltaH1))-(2*rc*(deltaH2T)))/(Cpa*Fa+Cpb*Fb+Cpc*Fc+Cpi*Fio);
f = [dFadV; dFbdV; dFcdV; dTdV; dFidV];
end

[Mischa Kim]

Difficult to say what the issue is without digging further into the differential equation and your code. You can try to play with different solvers. And I also suggest double-checking your code. Sometimes it is one little plus vs minus sign (or vice versa) that makes all the difference. Feel free to share a screen shot of the equations.

[Zara Freeman]

I got it to work! Thank you !

Answer 3

[Star Strider]

The most significant problem is that you need to define the functions as anonymous functions:

K1 = @(T) 50*exp((8000/8.314)*((1/315)-(1/T))); % dm3/mol.min
Kc = @(T) 10*exp((-25/8.314)*((1/315)-(1/T))); % dm3/mol.min
K2 = @(T) 400*exp((4000/8.314)*((1/310)-(1/T))); % dm3/mol.min-1

and then call them with the appropriate arguments:

ra = (K1(T)*((Ca^2)-((1/Kc(T))*Cb))) + (K2(T)*Ca*(Cb^2));
rb = 1/2*((K1(T)*((Ca^2)-((1/Kc(T))*Cb)))+(2*K2(T)*(Cb^2)*Ca));
rc = K2(T)*Ca*(Cb^2);

See the documentation section on Anonymous Functions for details.

The other problem is that you have 4 differential equations:

f = [dFadV; dFbdV; dFcdV; dTdV];

and 5 initial conditions:

yo = [4 0 0 1 500.15];

Correct these and your code runs, however it gives a Warning:

Warning: Failure at t=7.629620e-01.  Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (1.776357e-15) at time t. 
> In ode45 (line 360)

meaning that it has found a singularity and cannot go further. (I may not have edited ‘yo’ correctly to make the intiial conditions equal to the number of differential equations, so check that.)

The corrected ‘fun1’ code is:

function f = fun1(V,Y)
% Define the differential equations that need to be solved 
% Y is the concentration and V is the PFR volume
Fa = Y(1);
Fb = Y(2);
Fc = Y(3);
T = Y(4);
 
% Define initial conditions 
deltaH1 = -25000; %kJ/molA
deltaH2 = 35000; %kJ/molB
deltaH2T = 35000-(80*(T-298)); %kJ/molB
CTo = 0.3996; % mol/L
To = 500.15; % K
Fio = 1; % mol/min
FTo = 4; % mol/min
Cpa = 20; % J/molK
Cpb = 80; % J/molK
Cpc = 100; % J/molK
Cpi = 20; % J/molK
Ua = 150; %J/dm3.min.K
FT = Fa + Fb + Fc + Fio;
Ca = CTo*((Fa/FT)*(To/T));
Cb = CTo*((Fb/FT)*(To/T));
Cc = CTo*((Fc/FT)*(To/T));
Ci = CTo*((Fio/FT)*(To/T));
K1 = @(T) 50*exp((8000/8.314)*((1/315)-(1/T))); % dm3/mol.min
Kc = @(T) 10*exp((-25/8.314)*((1/315)-(1/T))); % dm3/mol.min
K2 = @(T) 400*exp((4000/8.314)*((1/310)-(1/T))); % dm3/mol.min-1
Ta = 523.15; % K
ra = (K1(T)*((Ca^2)-((1/Kc(T))*Cb))) + (K2(T)*Ca*(Cb^2));
rb = 1/2*((K1(T)*((Ca^2)-((1/Kc(T))*Cb)))+(2*K2(T)*(Cb^2)*Ca));
rc = K2(T)*Ca*(Cb^2);
% Differential equations
dFadV = -ra;
dFbdV = -rb;
dFcdV = rc;
dTdV = (Ua*(Ta-T)-((ra-rc)*(deltaH1))-(2*rc*(deltaH2)))/(Cpa*Fa+Cpb*Fb+Cpc*Fc+Cpi*Fio);
f = [dFadV; dFbdV; dFcdV; dTdV];
end

.

[Zara Freeman]

Hi Thank you for your help! I have made the relevant corrections, however this error is now displayed:

Warning: Failure at t=1.109564e-01. Unable to meet integration tolerances without reducing the step

size below the smallest value allowed (2.220446e-16) at time t.

> In ode45 (line 360)

In reactorassignmentpart1plot (line 5)

%% Chemical Systems Engineering Reactor Assignment 
clc
Vspan = [0:0.1:1];
yo = [4 0 0 1 500.15];
[V,y] = ode45(@fun1, Vspan, yo);
 
subplot(2,1,1)
plot(V,y(:,1),V,y(:,2),V,y(:,3));
legend('Fa','Fb','Fc');
ylabel('Molar flowrate, mol/min');
xlabel('Volume,dm3');
subplot(2,1,2)
plot(V,y(:,4));
legend('Temperature (K)');
ylabel('Temperature (K)');
xlabel('Volume,dm3');
 
% Compose the function
function f = fun1(V,Y)
% Define the differential equations that need to be solved 
% Y is the concentration and V is the PFR volume
Fa = Y(1);
Fb = Y(2);
Fc = Y(3);
Fi = Y(4);
T = Y(5);
 
% Define initial conditions 
deltaH1 = -25000; %kJ/molA
deltaH2 = 35000; %kJ/molB
deltaH2T = 35000-(80*(T-298)); %kJ/molB
CTo = 0.3996; % mol/L
To = 500.15; % K
Fio = 1; % mol/min
FTo = 4; % mol/min
Cpa = 20; % J/molK
Cpb = 80; % J/molK
Cpc = 100; % J/molK
Cpi = 20; % J/molK
Ua = 150; %J/dm3.min.K
FT = Fa + Fb + Fc + Fio;
Ca = CTo*((Fa/FT)*(To/T));
Cb = CTo*((Fb/FT)*(To/T));
Cc = CTo*((Fc/FT)*(To/T));
Ci= CTo*((Fio/FT)*(To/T));
K1 = 50*exp((8000/8.314)*((1/315)-(1/T))); % dm3/mol.min
Kc = 10*exp((-25/8.314)*((1/315)-(1/T))); % dm3/mol.min
K2 = 400*exp((4000/8.314)*((1/310)-(1/T))); % dm3/mol.min
Ta = 523.15; % K
ra = (K1*((Ca^2)-((1/Kc)*Cb))) + (K2*Ca*(Cb^2));
rb = 1/2*((K1*((Ca^2)-((1/Kc)*Cb)))+(2*K2*(Cb^2)*Ca));
rc = K2*Ca*(Cb^2);
% Differential equations
dFadV = -ra;
dFbdV = -rb;
dFcdV = rc;
dFidV = 0;
dTdV = (Ua*(Ta-T)-((ra-rc)*(deltaH1))-(2*rc*(deltaH2T)))/(Cpa*Fa+Cpb*Fb+Cpc*Fc+Cpi*Fio);
f = [dFadV; dFbdV; dFcdV; dTdV; dFidV];
end

[Star Strider]

My pleasure!

Changing solvers could work:

[V,y] = ode23(@fun1, Vspan, yo);

.

[Zara Freeman]

Thank you ! I managed to get it to work

[Star Strider]

My pleasure!