# Why do I receive incorrect result from ode45?

17 views (last 30 days)
Tess Nguyen on 24 Sep 2021 at 0:56
Commented: Tess Nguyen on 24 Sep 2021 at 5:03
Hello, I'm trying to solve this system of differential equations using ode45 function.
The code I designed is written below.
concentration.m file
%This models concentration of CH4 over FBR length
function diffeqs=concentration(~,var)
cH2 = var(1);
cCH4 = var(2); %define variable
cC2H4 = var(3); %Imporant!!
cC2H6 = var(4);
cC6H6 = var(5);
cC7H8 = var(6);
cC10H8 = var(7);
a1 = var (8);
a2 = var (9);
%Rate constants
k1 = 0.01831;
k2 = 1231;
k3 = 11340000;
k4 = 1.697;
k5 = 706300;
k6 = 2917;
%Equilibrium constants
K1 = (3.07941*10^-5);
K2 = 5.094435141; %chane k2 and k6
K3 = 721589.6019;
K4 = 142356.6252;
K5 = 0.007763912;
K6 = 9405.352961;
K1H2 = 0.07581; %Site 1
K1CH4 = 0.09664;
K1C2H4 = 0.07944;
K1C2H6 = 0.08095;
K2H2 = 0.09396; %Site 2
K2CH4 = 0.1189;
K2C2H4 = 0.2269;
K2C2H6 = 0.05123;
K2C6H6 = 0.04881;
K2C7H8 = 0.05219;
K2C10H8 = 0.1175;
%Conversion between concentration, partial pressure and molar flowrate
R = 8.314;
T = 973.15; %in K
PH2 = (cH2*R*T*10^-5); %bar pressure can be drop
PCH4 = (cCH4*R*T*10^-5);
PC2H4 = (cC2H4*R*T*10^-5);
PC2H6 = (cC2H6*R*T*10^-5);
PC6H6 = (cC6H6*R*T*10^-5);
PC7H8 = (cC7H8*R*T*10^-5);
PC10H8 = (cC10H8*R*T*10^-5);
%DEN
DEN1 = (1 + K1H2*PH2 + K1CH4*PCH4 + K1C2H4*PC2H4 + K1C2H6*PC2H6);
DEN2 = (1 + K2H2*PH2 + K2CH4*PCH4 + K2C2H4*PC2H4 + K2C2H6*PC2H6 + K2C6H6*PC6H6 + K2C7H8*PC7H8 + K2C10H8*PC10H8);
%Deactivation
diffeqs(8,1) = -1.5;
diffeqs(9,1) = -0.0091;
%Reaction rate calculation
r1 = (k1*(PCH4^2 - PC2H4*PH2^2/K1)*a1/DEN1);
r2 = (k2*(PC2H4^4 - PC6H6*PC2H6*PH2^2/K2)*a2/DEN2);
r3 = (k3*(PC2H4^3 - PC6H6*PH2^3/K3)*a2/DEN2);
r4 = (k4*(PC6H6*PCH4 - PC7H8*PH2/K4)*a2/DEN2);
r5 = (k5*(PC6H6*PC2H4^2 -PC10H8*PH2^3/K5)*a2/DEN2);
r6 = (k6*(PC2H4*PH2 - PC2H6/K6)*a1/DEN1);
%Constants
us = 0.3068; % Feed velocity meter/h
pB = 500; %kg/m3
%Differential eqs
diffeqs(1,1) = pB*(2*r1 +2*r2 +3*r3 +r4 +3*r5 -r6)/us; %H2
diffeqs(2,1) = pB*(-2*r1 -r4)/us; %CH4
diffeqs(3,1) = pB*(r1 -4*r2 -3*r3 -2*r5 -r6)/us; %C2H4
diffeqs(4,1) = pB*(r2 +r6)/us; %C2H6
diffeqs(5,1) = pB*(r2 + r3 - r4 -r5)/us; %C6H6
diffeqs(6,1) = pB*r4/us; %C7H8
diffeqs(7,1) = pB*r5/us; %C8H10
end
solving.m file
%Script to solve DEs
range =[0 0.39]; %fix to real length
ICs = [0, 0.9, 0, 0, 0, 0, 0, 1, 1]; %mol/m3
[z,conc]=ode45(@concentration,range,ICs);
figure
plot(z,conc(:,1)
%xaxis('z, m')
%ylabel('conc, mol/m3')
From the rate constant values, k3 is much greater than k6 so it expects to get cC6H6 larger than cC2H6 but I got the otherwise result.
Please help me check if the incorrect result coming from the code itself or it might be the logic of differential eqn system.
Many thanks!

Sulaymon Eshkabilov on 24 Sep 2021 at 3:46
Edited: Sulaymon Eshkabilov on 24 Sep 2021 at 3:46
You can try to tighten the relative and absolute error tolerances, e.g.:
OPTS = odeset('reltol', 1e-5, 'abstol', 1e-7);
[z,conc]=ode45(@concentration, Range,ICs, OPTS);
Another point (a good pratice) is not use MATLAB's built in fcn names. In your code, you have range() that should be renamed, e.g.: Range
One typo - missing ) in plot(z,conc(:,1) that was missed due to copy and paste.
Tess Nguyen on 24 Sep 2021 at 5:03
Thanks a lot for your suggestion.
I've tried to change the error tolerances but the result is the same.
I guess the issue coming from other spot.

R2021a

### Community Treasure Hunt

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

Start Hunting!