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

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.

Answer 1

[Alex Sha]

Hi, Vira Roy, attached your data please, if possible.

[Vira Roy]

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.

Answer 2

[Alex Sha]

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

Adjusted R-Square: 0.999930030352741

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]

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)));