See Parameter Estimation for a System of Differential Equations for a working example of how to do this correctly.
Unfortunately, since you have 7 data rows (governed by the number of elements in the ‘Time’ vector, since that is the independent variable), and are estimating 12 parameters (my version estimates 21 parameters, however the last 9 are the initial conditions of the differential equations, so they don’t count), this will fail.
This runs without error, however because of the mismatch between the data and the parameters being estimated, the results will not be reliable.
Adapting that approach to your code —
function c = kinetics(k,t)
c0 = k(13:21);
[t,Cv] = ode15s(@DifEq,t,c0);
function dcdt=DifEq(t,c)
dcdt = zeros(9,1);
dcdt(1)=-k(1)*c(1)-k(2)*c(1)-k(3)*c(1);
dcdt(2)=k(2)*c(1)-k(6)*c(2)-k(9)*c(2);
dcdt(3)=k(3)*c(1)-k(5)*c(4)-k(7)*c(6);
dcdt(4)=k(1)*c(1)-k(4)*c(4)-k(5)*c(4);
dcdt(5)=k(6)*c(2)+k(4)*c(4)-k(10)*c(5);
dcdt(6)=k(7)*c(3)-k(8)*c(6);
dcdt(7)=k(8)*c(6)+k(9)*c(2)+k(10)*c(5)-k(11)*c(8)-k(12)*c(9);
dcdt(8)=k(11)*c(7);
dcdt(9)=k(12)*c(8);
end
c = Cv;
end
Time=[0 2 4 10 15 25 40];
c= [100 0 0 0 0 0 0
0 4 37 15 5 4 2
0 10 20 8 4 0.2 0.1
0 5 10 16 4 2 1
0 4 5 11 5 4 3
0 5 6 5 1 0.1 0.04
0 6 12 32 64 44 31
0 0.0071 0.06 0.14 0.3000 3.8500 7.4440
0 0 0 0.0340 0.340 0.5000 0.630];
c = c.';
k0 = rand(21,1)*1E-3;
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,Time,c,zeros(size(k0)));
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(k)
fprintf(1, '\t\tk(%2d) = %8.5f\n', k1, k(k1))
end
tv = linspace(min(Time), max(Time));
Cfit = kinetics(k, tv);
figure(1)
hd = plot(Time, 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')
lgd = legend(hlp,compose('C_{%d}(t)',1:9), 'Location','N');
lgd.NumColumns = 3;
This will work best if you put all of it inside another function, and then call the external function. See the link I referred to for details.
However more data (at least 12 rows of the transposed ‘c’ matrix) are absolutely necessary if you want to estimate all 12 parameters accurately.
EDIT — (3 Mar 2021 at 4:02)
Added plot —
.
