Reaction kinetics parameter determination
24 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
Star Strider
2021 年 3 月 3 日
編集済み: Star Strider
2021 年 3 月 3 日
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 件のコメント
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!