function model4
% 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=[7.00E-02;0.00E+00;0.00E+00;0.00E+00];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(4,1);
dcdt(1)=(-theta(1)*c(1)*0.00273-theta(2)*c(1)*0.00273)/(1+(theta(6)*c(1)+theta(7)*c(2)+theta(8)*c(3)+theta(9)*c(4)+theta(10)*0.00273))^3;
dcdt(2)= (theta(1)*c(1)*0.00273-theta(3).*c(2)*0.00273-theta(5)*c(2))/(1+(theta(6)*c(1)+theta(7)*c(2)+theta(8)*c(3)+theta(9)*c(4)+theta(10)*0.00273))^3;
dcdt(3)= (theta(2)*c(1)*0.00273-theta(4)*c(2)*0.00273+theta(5)*c(2))/(1+(theta(6)*c(1)+theta(7)*c(2)+theta(8)*c(3)+theta(9)*c(4)+theta(10)*0.00273))^3;
dcdt(4)= (theta(3)*c(2)*0.00273+theta(4)*c(3)*0.00273)/(1+(theta(6)*c(1)+theta(7)*c(2)+theta(8)*c(3)+theta(9)*c(4)+theta(10)*0.00273))^3;
dC=dcdt;
end
C=Cv;
end
t=[20
40
60
90
120
180
240];
c=[5.69E-02 1.36E-03 8.25E-03 3.48E-03
4.28E-02 9.79E-03 9.50E-03 7.94E-03
2.93E-02 2.85E-03 1.53E-02 2.26E-02
1.46E-02 4.69E-03 2.00E-02 3.28E-02
1.19E-03 3.52E-03 1.10E-02 5.43E-02
1.75E-03 0.00E+00 0.00E+00 6.10E-02
1.00E-09 1.00E-09 1.00E-09 7.00E-02];
theta0=[1.86e-1,1.70e-1,1.28,1.51e-1,1.0e-1,15.90,2.60e-1,1.87e-2,3.25,2.33e-2];
[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')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')
end
THE PROBLEM IS ORANGE STAR AND RED LINE ARE NOT FITTING THE TOP SHOULD BE AROUND 0.02
