Error in my ode45 equation.
6 views (last 30 days)
Show older comments
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 !
0 Comments
Answers (3)
Walter Roberson
on 29 Dec 2020
clc
Vspan = [0 1];
yo = [4 0 0 1 500.15];
[V,y] = ode45(@fun1, Vspan, yo);
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
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
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.
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);
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
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!