Why do I receive incorrect result from ode45?

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;

%Adsorption Equilibrium constants

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!

Answer 1

[Sulaymon Eshkabilov]

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]

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.