You are trying to fit 5 columns of your differential equation to 4 columns of data. That will throw the error you got.
Change the ‘C’ assignment in the ‘kinetics’ objective function to:
and it works.
So ‘klinetics’ is now:
function C=kinetics(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
function dC=DifEq(t,c)
dcdt=zeros(5,1);
dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
dcdt(5)=theta(1).*c(2)-theta(5);
dC=dcdt;
end
C=Cv(:,1:4);
end
If you then want to estimate Compartment 5, run your differential equation function again with all the estimated parameters to include it and plot it.
I created a version of ‘kinetics’ called ‘kinetics5’ to calculate and plot Compartment 5 as well. That entire additional code (with all the previously estimated parameters, so all this goes after the lsqcurvefit call) is:
function C=kinetics5(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
function dC=DifEq(t,c)
dcdt=zeros(5,1);
dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
dcdt(5)=theta(1).*c(2)-theta(5);
dC=dcdt;
end
C=Cv;
end
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics5(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)', 'C_5(t)', 'Location','SW')
That should do what you want.
Best Answer