MATLAB Answers

Reaction kinetics parameter determination

34 views (last 30 days)
iPanda on 3 Mar 2021
Commented: Star Strider on 3 Mar 2021
Hi, I am new to MATLAB and I have a series of data points to fit and then calculate the rate constants.
here are the equations:
dcdt(1)=-k(1)*c(1)-k(2)*c(1)-k(3)*c(1)
dcdt(2)=k(2)*c(1)-k(6)*c(2)-k(9)*c(2)
dcdt(3)=k(3)*c(1)-k(5)*c(4)-k(7)*c(6)
dcdt(4)=k(1)*c(1)-k(4)*c(4)-k(5)*c(4)
dcdt(5)=k(6)*c(2)+k(4)*c(4)-k(10)*c(5)
dcdt(6)=k(7)*c(3)-k(8)*c(6)
dcdt(7)=k(8)*c(6)+k(9)*c(2)+k(10)*c(5)-k(11)*c(8)-k(12)*c(9)
dcdt(8)=k(11)*c(7)
dcdt(9)=k(12)*c(8)
and the data:
Time=[0 2 4 10 15 25 40]
c= [100 0 0 0 0 0 0
0 4 37 15 5 4 2
0 10 20 8 4 0.2 0.1
0 5 10 16 4 2 1
0 4 5 11 5 4 3
0 5 6 5 1 0.1 0.04
0 6 12 32 64 44 31
0 0.0071 0.06 0.14 0.3000 3.8500 7.4440
0 0 0 0.0340 0.340 0.5000 0.630]
I would really appreciate the help! Thanks in advance
0 CommentsShowHide -1 older comments

Sign in to comment.

Accepted Answer

Star Strider on 3 Mar 2021
Edited: Star Strider on 3 Mar 2021
See Parameter Estimation for a System of Differential Equations for a working example of how to do this correctly.
Unfortunately, since you have 7 data rows (governed by the number of elements in the ‘Time’ vector, since that is the independent variable), and are estimating 12 parameters (my version estimates 21 parameters, however the last 9 are the initial conditions of the differential equations, so they don’t count), this will fail.
This runs without error, however because of the mismatch between the data and the parameters being estimated, the results will not be reliable.
Adapting that approach to your code —
function c = kinetics(k,t)
c0 = k(13:21);
[t,Cv] = ode15s(@DifEq,t,c0);
function dcdt=DifEq(t,c)
dcdt = zeros(9,1);
dcdt(1)=-k(1)*c(1)-k(2)*c(1)-k(3)*c(1);
dcdt(2)=k(2)*c(1)-k(6)*c(2)-k(9)*c(2);
dcdt(3)=k(3)*c(1)-k(5)*c(4)-k(7)*c(6);
dcdt(4)=k(1)*c(1)-k(4)*c(4)-k(5)*c(4);
dcdt(5)=k(6)*c(2)+k(4)*c(4)-k(10)*c(5);
dcdt(6)=k(7)*c(3)-k(8)*c(6);
dcdt(7)=k(8)*c(6)+k(9)*c(2)+k(10)*c(5)-k(11)*c(8)-k(12)*c(9);
dcdt(8)=k(11)*c(7);
dcdt(9)=k(12)*c(8);
end
c = Cv;
end
Time=[0 2 4 10 15 25 40];
c= [100 0 0 0 0 0 0
0 4 37 15 5 4 2
0 10 20 8 4 0.2 0.1
0 5 10 16 4 2 1
0 4 5 11 5 4 3
0 5 6 5 1 0.1 0.04
0 6 12 32 64 44 31
0 0.0071 0.06 0.14 0.3000 3.8500 7.4440
0 0 0 0.0340 0.340 0.5000 0.630];
c = c.';
k0 = rand(21,1)*1E-3;
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,Time,c,zeros(size(k0)));
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(k)
fprintf(1, '\t\tk(%2d) = %8.5f\n', k1, k(k1))
end
tv = linspace(min(Time), max(Time));
Cfit = kinetics(k, tv);
figure(1)
hd = plot(Time, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
lgd = legend(hlp,compose('C_{%d}(t)',1:9), 'Location','N');
lgd.NumColumns = 3;
This will work best if you put all of it inside another function, and then call the external function. See the link I referred to for details.
However more data (at least 12 rows of the transposed ‘c’ matrix) are absolutely necessary if you want to estimate all 12 parameters accurately.
EDIT — (3 Mar 2021 at 4:02)
Added plot —
.
4 CommentsShowHide 3 older comments
Star Strider on 3 Mar 2021
As always, my pleasure!

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!