Kinetic parameter estimation and initial time setting

10 ビュー (過去 30 日間)
Emmanuel Nzediegwu
Emmanuel Nzediegwu 2021 年 5 月 15 日
編集済み: Star Strider 2021 年 5 月 17 日
Hi: how can I set the initial time on this code to be at 0 min at an initial concentration of 100 mg/ml. This is a case study of a degradation reaction for which at 0 min, 100 mg/ml of the reactant was all present, but at 20 min of the reaction, the reactant started degrading into other products. As such, the X-axis should start at 0 min and not 20 min. Please I need assistance on this. The code is as below.
function HtwoT_How
% 2016 12 03
% NOTES:
%
% 1. The ‘theta’ (parameter) argument has to be first in your
% ‘kinetics’ funciton,
% 2. You need to return ALL the values from ‘DifEq’ since you are fitting
% all the values
function C=kinetics(theta,t)
c0=[100;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(4,1);
dcdt(1)=-theta(1).*c(1);
dcdt(2)= theta(1).*c(1)-theta(2).*c(2)-theta(3).*c(2);
dcdt(3)= theta(2).*c(2)-theta(4).*c(3);
dcdt(4)= theta(4).*c(3);
dC=dcdt;
end
C=Cv;
end
t=[25
35
45
55
65];
c=[20.76 5.93 2.77 69.54
16.37 5.72 0.6 77.31
19.88 3.9 0.19 76.03
8.01 2 0.18 89.81
17.73 0.15 0.16 81.96];
theta0=[1;1;1;1];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time (min)')
ylabel('Concentration (mg/ml)')
legend(hlp, 'C1', 'C2', 'C3', 'C4', 'Location','N')
end

採用された回答

Emmanuel Nzediegwu
Emmanuel Nzediegwu 2021 年 5 月 15 日
編集済み: Star Strider 2021 年 5 月 17 日
_

その他の回答 (1 件)

Alex Sha
Alex Sha 2021 年 5 月 15 日
Hi, see the results below:
Root of Mean Square Error (RMSE): 4.00513564326846
Sum of Squared Residual: 320.822230419589
Correlation Coef. (R): 0.788386977278981
R-Square: 0.621554025943088
Parameter Best Estimate
-------------------- -------------
theta1 0.0529114059211692
theta2 0.337616166932669
theta3 0.0180685434553373
theta4 52.0236365257367
  1 件のコメント
Emmanuel Nzediegwu
Emmanuel Nzediegwu 2021 年 5 月 15 日
Hi Alex, Thank you, please can you share the code you used as I have other data to attempt accordingly

サインインしてコメントする。

Community Treasure Hunt

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

Start Hunting!

Translated by