How to fit data on only one ordinary differential equations (out of multiple equations) with multiple unknown parameters
3 views (last 30 days)
Show older comments
Chinwendu Madubueze
on 12 Nov 2019
Commented: Star Strider
on 13 Nov 2019
Hi,
I am trying to estimate 5 unknown parameters that exist in 10 ordinary differential equations (ODEs) by fitting it to experimental data of only one equation of the ODE.
I used the Monod Kinetics and curve fitting code to estimate the parameters but I got negative parameters instead of postive. This is the output.
Maximum number of function evaluations exceeded;
increase options.MaxFunEvals
Rate Constants:
beta(1) = 1.22325
beta(2) = -0.04436
beta(3) = 0.23296
beta(4) = -0.47970
beta(5) = -0.56506
This is my code
function S = Lassa(beta, t)
x0 = [100000;2000;400;400;5;50;200;500;50;20];
[T,Sv] = ode45(@DifEq, t, x0);
function dS = DifEq(t, x)
xdot = zeros(10,1);
piH = 0.000215; muH = 0.0000548;d1=0.01; piR=0.56;eR=0.85;rhoR=0.3;muR=0.6;d2=0.001;rho3=0.2;
rho6=0.0232;rho2=0.05;rho1=0.9;rho7=0.45;rho5=0.85;omega=0.45;rho4=0.06;gamma=0.03;
NR=x(8)+x(9)+x(10); NH=x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7);
xdot(1)=piH - beta(4)*x(9)*x(1)/NR -(beta(1)*x(5)+beta(2)*x(6)+beta(3)*x(7))*x(1)/NH -muH*x(1)+rho3*x(3)+rho6*x(7);
xdot(2)=beta(4)*x(9)*x(1)/NR+(beta(1)*x(5)+beta(2)*x(6)+beta(3)*x(7))*x(1)/NH -((1-omega)*rho2+ omega*rho1+muH)*x(2);
xdot(3)=(1-omega)*rho2*x(2)-(rho3*rho4+muH)*x(3);
xdot(4)=omega*rho1*x(2)-(rho5+muH)*x(4);
xdot(5)=rho5*x(4)-(rho7+d1+muH)*x(5);
xdot(6)=rho7*x(5)+rho4*x(3)-(rho6+d2+muH)*x(6);
xdot(7)=d1*x(5)+d2*x(6) - gamma*x(7);
xdot(8)=piR - beta(5)*x(9)*x(8)/NR - (muR+rhoR)*x(8);
xdot(9)=beta(5)*x(9)*x(8)/NR -(muR+rhoR+eR)*x(9);
xdot(10)=eR*x(9)-(muR+rhoR)*x(10);
dS = xdot;
end
S = Sv(:,5);
end
t=(1:37);
xdata = [1;5;2;5;7;4;1;0;2;9;3;1;0;1;0;0;0;0;1;0;0;0;0;0;0;0;0;0;3;0;0;0;0;0;0;0;1];
beta0=[0.5921;0.0162;0.01792;0.6506;0.7598];
[beta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Lassa,beta0,t,xdata);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(beta)
fprintf(1, '\t\tbeta(%d) = %8.5f\n', k1, beta(k1))
end
tv = linspace(min(t), max(t));
Cfit = Lassa(beta, tv);
0 Comments
Accepted Answer
Star Strider
on 12 Nov 2019
The lsqcurvefit function allows you to set upper and lower bounds on the parameters it estimates with the ‘lb’ and ‘ub’ vectors. Set the lower bound to:
zeros(size(beta0))
and do not specify the upper bounds (set them equal to the empty array []) or set them equal to
Inf(size(beta0))
The original Monod Kinetics code did not need to constrain the parameters, so I did not constrain them.
7 Comments
More Answers (0)
See Also
Categories
Find more on Get Started with Curve Fitting Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!