Data fitting to the solution of coupled differential equation with all the different parameters.

15 views (last 30 days)
Vira Roy on 3 Oct 2019
Edited: Vira Roy on 25 Oct 2019
So here is my problem. The general idea is that, I have a set of data for growth (n) and time(t) which i need to fit to the coupled differential equations. Star Strider has an excellent answer here using Monod Kinetics but after modelling my code over his, i am not able to generate correct fitting to my data.
Here is my differential equation and for convenience i have kept the coefficients to be theta itself as in his code.
Now here lies the difference, from the earlier code. I do not have any kind of coupling between the coefficients and all of them appear only once and i have to find their optmum value for the fit and those have to be positive.
Here is how i tried to do. I have defined a function Growth_fn.m and stored the differential equations there and then in code_body.m i modified the code.
Growth_fn.m
function G=Growth_fn(theta,t)
IC=[0.1;0];
[T,Cv]=ode45(@DifEq,t,IC);
function dG=DifEq(t,y)
rkm=zeros(2,1);
rkm(1)= theta(1).*y(1).*(1-(y(1)./theta(2)))+ theta(3).*y(1).*(1-exp(-0.04*t))-theta(4).*y(1).*y(2);
rkm(2)=0.3+((theta(5).*y(2).*y(1))./(theta(6)+y(1)))-theta(7).*y(1).*y(2)- 0.0208.*y(2);
dG=rkm;
end
G=Cv(:,1);
end
code_body.m
clear all; clc
%data calling here
theta0=[0.1;0.1;0;0;0;0;0];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Growth_fn,theta0,t_data,y,zeros(size(theta0)));
fprintf(1,'\t Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t_data), max(t_data));
Cfit = kinetics(theta, tv);
figure(1)
plot(t_data, y, 'p')
hold on
hlp = plot(tv, Cfit);
hold on
grid
xlabel('Time')
ylabel('growth')
The output
Constants:
Theta(1) = 0.08016
Theta(2) = 0.44688
Theta(3) = 0.04996
Theta(4) = 0.02754
Theta(5) = 0.00005
Theta(6) = 0.35045
Theta(7) = 0.34710
The fit i get from the code.
This is what i get when i actually put those theta values into my original differential equation. They do not match at all.
I must tell you that, this code works good enough for other data sets like this one gives a considerably good fit.
t_data=[11,15,18,23,26,31,39,44,54,64,74]';
y=[0.00476,0.0105,0.0207,0.0619,0.337,0.74,1.7,2.45,3.5,4.5,5.09]';
Should i use an entirely different method or a different optimzation technique. Is this code not workable for this particular problem?
Edit: I had added the data and got a reply but didn't solve the problem, and due to certain constraint i have to remove it.
Thank You for any help anyways.

Alex Sha on 6 Oct 2019
Vira Roy on 14 Oct 2019
Alex, I tried a more but couldn't find my mistake. I have attached the data, data.xlsx. Could you see where i am making the mistake.

Alex Sha on 14 Oct 2019
if take initial conditions as: t=30,n=0.1178512,i=0;
Root of Mean Square Error (RMSE): 0.00101378248364856
Sum of Squared Residual: 0.000101747737491112
Correlation Coef. (R): 0.999965728565097
R-Square: 0.999931458304726
Determination Coef. (DC): 0.999931458057279
F-Statistic: 223692.310353222
Parameter Best Estimate
-------------------- -------------
theta1 -0.208539738687139
theta2 0.528354845354833
theta3 0.17419907588849
theta4 -0.0567427601648523
theta5 -17.6206493824257
theta6 0.239479728245077
theta7 -28.4178353966237
Vira Roy on 14 Oct 2019
Thanks for quick reply and may be this is trivial Alex, but i need all the postive values for the coefficients. For that in the above code i used the zeros(size(theta0)). It is a physical constraint in my problem.
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Growth_fn,theta0,t_data,y,zeros(size(theta0)));

R2016a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!