That does not appear to me to be the same data.
I get decent results except for 
and 
, leading me to believe that there is something wrong with their differential equations. Please check them.
and , leading me to believe that there is something wrong with their differential equations. Please check them. % clc;
% clear all;
t=[ 0
10
20
30
40
50
60];
c=[ 100 100 0 0 0 0 0
94.77333333 93.73333333 0.400145927 0.27701074 0.107176434 0.158130584 0.000505722
84.48 88.87333333 0.604820791 0.42039699 0.148464431 0.359954181 0.001077632
73.91333333 83.23333333 0.589530533 0.461026019 0.17071451 0.508774341 0.001205872
68.80666667 76.43333333 0.529979 0.48559649 0.180350987 0.635601375 0.001720902
63.56666667 73.92 0.502282172 0.52315421 0.198918834 0.700114548 0.002064255
59.76666667 70.59333333 0.478273783 0.574401193 0.207380132 0.845269187 0.002895749
];
theta0 = [rand(11,1)*1E-3; c(1,:).'];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));
fprintf(1,'\tRate Constants:\n')
for k1 = 1:numel(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure(1)
hd = plot(t, 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')
legend(hlp, 'c(1)', 'c(2)', 'c(3)', 'c(4)', 'c(5)','c(6)','c(7)', 'Location','best');
figure
for k = 1:size(c,2)
subplot(4,2,k)
plot(t,c(:,k),'bp', 'MarkerFaceColor','b')
hold on
plot(tv, Cfit(:,k))
hold off
grid
title("C_"+k)
end
function C=kinetics(theta,t)
% c0=[100; 100; 0; 0; 0; 0; 0];
c0 = theta(end-6:end);
[T,Cv]=ode45(@DifEq,t,c0);
function dcdt=DifEq(t,c)
dcdt=zeros(7,1);
dcdt(1)=-theta(1)*c(1)-theta(2)*c(1)-theta(3)*c(1)-theta(9)*c(1)*c(2);
dcdt(2)=-theta(4)*c(2)*c(3)-theta(5)*c(2)*c(4)-theta(6)*c(2)*c(3)-theta(7)*c(2)*c(4)-theta(8)*c(2)-theta(9)*c(2)*c(1);
dcdt(3)= theta(2)*c(1)-theta(4)*c(2)*c(3)-theta(6)*c(2)*c(3);
dcdt(4)= theta(3)*c(1)-theta(5)*c(2)*c(4)-theta(7)*c(2)*c(4);
dcdt(5)= theta(4)*c(2)*c(3)-theta(10)*c(5)+theta(9)*c(2)*c(1);
dcdt(6)= theta(5)*c(2)*c(4)-theta(11)*c(6);
dcdt(7)= theta(6)*c(2)*c(3)+theta(7)*c(2)*c(4)+theta(8)*c(2);
end
C=Cv;
end
I changed the code slightly to add the initial conditions as parameters to be estimated, since in my experience, that generally produces better results, even if the initial conditions themselves do not change much from their original values. (I added the subplot array because otherwise the last 5 compartments are squashed together at the lower end of the first plot.)
.
