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
Best Answer