MATLAB Answers

Coupled rate ODEs with ode45

65 views (last 30 days)
Ryan McLeod
Ryan McLeod on 24 Feb 2021
Commented: Star Strider on 26 Feb 2021
Dear all,
I have been having trouble trying to model the concentrations of 8 species in the attached file showing the reactions. I have the rate constant values, initial conditions and a time frame but for some reason the plot seems to display no changes in concentration. I have attached my function and script. Any pointers or tips would be greatly appreciated.
Thanks in advance.

Accepted Answer

Star Strider
Star Strider on 24 Feb 2021
The concentrations change appropriately, however they don’t change much and the concentrations are vanishingly small. That’s the reason they don’ appear to change.
Increasing the final time (by 15 times) demonstrates saturation kinetics:
[t,state]=ode15s(@reactions,[0 51.5*15],[2.71E-5;7.8E-6;9.2E-6;0;1.0614E-4;3.359E-3;3.9278E-4;4.521E-4]);
reactants = {'TAG','DAG','MAG','GLY','FAEE','ET','FFA','W'};
figure
for k = 1:8
subplot(4,2,k)
plot(t,state(:,k))
xlabel('Time')
ylabel('Concentration')
title(sprintf('$C_{%s}(t)$',reactants{k}), 'Interpreter','latex')
grid
end
producing:
Adding:
pos = get(gcf,'Position')
set(gcf, 'Position',pos+[0 -500 200 500])
after the loop enlarges the figure making the subplot plots easier to see. (I did not use that in the posted figure.)
  7 Comments
Star Strider
Star Strider on 26 Feb 2021
As always, my pleasure!
The problem with the latest approach is at least in part the problem of estimating four parameters with only two data points (‘Inlet’ and ‘Outlet’). If you have the time evolution of the various reactants and products over at least one more time instants than the number of parameters to be estimated (more is better), you can use the approach I used in the Answer I linked to. (I’ve updated that code since, so if you want to use it and if you post the necessary data here, as well as the appropriate differential equations, I’ll post back the updated code with an attempt at estimating the parameters, since lsqcurvefit may not be optimal and a genetic algorithm approach may be necessary.)

Sign in to comment.

More Answers (1)

Alan Stevens
Alan Stevens on 24 Feb 2021
It can all be done in one script as follows. Because of the orders of magntude difference between the various concentrations they just look constant when plotted on a single graph. They do vary as can be seen by running the following:
tspan = [0 51.5];
Y0 = [2.71E-5;7.8E-6;9.2E-6;0;1.0614E-4;3.359E-3;3.9278E-4;4.521E-4];
[t,Y]=ode15s(@reactions,tspan,Y0);
TG = Y(:,1); FAEE = Y(:,5);
DG = Y(:,2); ET = Y(:,6);
MF = Y(:,3); FFA = Y(:,7);
G = Y(:,4); W = Y(:,8);
subplot(4,2,1)
plot (t,TG),grid
ylabel('TG')
subplot(4,2,3)
plot (t,DG),grid
ylabel('DG')
subplot(4,2,5)
plot (t,MF),grid
ylabel('MF')
subplot(4,2,7)
plot (t,G),grid
xlabel('t'),ylabel('G')
subplot(4,2,2)
plot (t,FAEE),grid
ylabel('FAEE')
subplot(4,2,4)
plot (t,ET),grid
ylabel('ET')
subplot(4,2,6)
plot (t,FFA),grid
ylabel('FFA')
subplot(4,2,8)
plot (t,W),grid
xlabel('t'),ylabel('W')
function dydt = reactions(~,y)
% y1=TG y5=FAEE
% y2=DG y6=ET
% y3=MF y7=FFA
% y4=G y8=W
k1=0.6;
keq1=6.504E-6;
k2=0.8403;
keq2=0.02386;
k3=0.83524;
keq3=2.464;
k4=0.381;
keq4=12.182;
r1 = k1*((y(1)*y(6)*y(7))-((1/keq1)*y(2)*y(5)*y(7)));
r2 = k2*((y(2)*y(6)*y(7))-((1/keq2)*y(3)*y(5)*y(7)));
r3 = k3*((y(3)*y(6)*y(7))-((1/keq3)*y(4)*y(5)*y(7)));
r4 = k4*((y(7)*y(6)*y(7))-((1/keq4)*y(8)*y(5)*y(7)));
dy1dt=-r1;
dy2dt=r1-r2;
dy3dt=r2-r3;
dy4dt=r3;
dy5dt=r1+r2+r3+r4;
dy6dt=-(r1+r2+r3+r4);
dy7dt=-r4;
dy8dt=r4;
dydt=[dy1dt;dy2dt;dy3dt;dy4dt;dy5dt;dy6dt;dy7dt;dy8dt];
end
You can save and run this as a single script file.
  2 Comments
Ryan McLeod
Ryan McLeod on 24 Feb 2021
Too add more contect these are the sorts of plots I was expecting to see within the time domain. I assume the concentrations that I have for my reactor are the reason for the varying plots but as can be seen the FAEE concentration is expected to increase with the time.

Sign in to comment.

Categories

Community Treasure Hunt

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

Start Hunting!