Error in my ode45 equation.

9 views (last 30 days)
Zara Freeman
Zara Freeman on 29 Dec 2020
Commented: Star Strider on 29 Dec 2020
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 !

Sign in to answer this question.

Answers (3)

Walter Roberson
Walter Roberson on 29 Dec 2020
Ran in:
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.
Mischa Kim
Mischa Kim on 29 Dec 2020
Edited: Mischa Kim on 29 Dec 2020
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.
  3 Comments
Show 1 older comment
Mischa Kim
Mischa Kim on 29 Dec 2020
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
Zara Freeman on 29 Dec 2020
I got it to work! Thank you !

Sign in to comment.

Star Strider
Star Strider on 29 Dec 2020
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
.
  4 Comments
Show 2 older comments
Zara Freeman
Zara Freeman on 29 Dec 2020
Thank you ! I managed to get it to work

Sign in to comment.

Categories

Find more on Startup and Shutdown in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!