MATLAB Answers

Kinetic parameter estimation and initial time setting

3 views (last 30 days)
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

More Answers (1)

Alex Sha
Alex Sha on 15 May 2021
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 Comment
Emmanuel Nzediegwu
Emmanuel Nzediegwu on 15 May 2021
Hi Alex, Thank you, please can you share the code you used as I have other data to attempt accordingly

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by