Coefficient estimation for a system of coupled ODEs

5 views (last 30 days)
pablo cardona
pablo cardona on 30 Aug 2021
Commented: Alex Sha on 9 Sep 2021
Hello to you all, I'm a complete beginner to Matlab and I'm trying to solve a system of 4 ODEs that have unknown parameters. My intention is to estimate the parameters to obtain a solution then use experimental data to optimize the solution of the ODE system.
I modiffied the code provided by user Star Strider in the question: "Parameter Estimation for a System of Differential Equations" to use the ODEs and data I posses since the two problems are very similar. However the resulting fit the program found indicated to me that I am doing something wrong.
I will attach the .m with my code, inside is also contained the experimental data.
For reference the ODE system I'm trying to solve is the one displayed on the attached image.

Accepted Answer

Star Strider
Star Strider on 30 Aug 2021
I do not see that you are doing anything wrong, and the code runs without error.
The changes I made were to estimate the initial conditions as well as the other parameters, and use random values for ‘theta0’. I expanded the format for the ‘Constants:’ output in order to see the full precision, and changed the plot loop so that the fitted lines are the same colours ad the data points. (I made this change in other posts on this topic, and so am updating them here as well.) The fit appears to be reasonably good with these changes.
The notice is simply that lsqcurvefit exceeded the function evaluation limit before it was convinced that the model converged. You can increase the function evaluation limit (see the documentation section on options), however that may not be necessary. (I increased it to be 15000 and the iteration limit to 1500, and the fit appears to be even better, and the estimation converges.) Use the estimated parameters here for ‘theta0’ if you want to run it with changed options, since this appears to be reasonably good approximation.
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. Constants: Theta( 1) = 1.225621526046 Theta( 2) = 8.818523145614 Theta( 3) = 0.832213488120 Theta( 4) = 8864.918733082799 Theta( 5) = 1.234266936058 Theta( 6) = 0.810227130705 Theta( 7) = 0.157698146317 Theta( 8) = 0.023932541443 Theta( 9) = 0.089951988393 Theta(10) = 36.855992659497 Theta(11) = 0.565552316098 Theta(12) = 0.000134265635
function modelacion
function C=kinetics(theta,t)
% c0=[1;0;0;0];
c0 = theta(9:12);
function dC=DifEq(t,c)
dcdt(1)= (theta(1).*(1-(c(1)./theta(2)))-theta(3)).*c(1);
dcdt(2)= -((1./theta(4)).*dcdt(1))-(theta(5).*c(1));
dcdt(3)= theta(6).*dcdt(1)+theta(7).*c(1);
dcdt(4)= theta(8).*dcdt(1);
c=[0.075 36.77419355 0.317906336 0
0.075 36.77419355 0.545179063 0.000434783
0.175 36.4516129 0.71184573 0.001304348
0.325 36.12903226 0.909090909 0.003478261
0.425 36.12903226 1.106060606 0.006086957
0.525 35.48387097 1.242424242 0.010869565
0.675 34.19354839 1.439393939 0.013913043
0.925 33.5483871 1.636363636 0.020869565
1.15 32.25806452 1.985123967 0.026521739
1.475 30.32258065 2.454820937 0.032608696
1.8 28.38709677 3.015702479 0.040869565
2.075 26.12903226 3.54600551 0.047826087];
% theta0=[1;1;1;1;1;1;1;1];
theta0 = rand(12,1);
opts = optimoptions('lsqcurvefit', 'MaxFunctionEvaluations', 15000, 'MaxIterations',1500);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)),Inf(size(theta0)), opts);
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%2d) = %18.12f\n', k1, theta(k1))
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
hold off
legend(hlp, compose('C_{%d}(t)',1:size(c,2)), 'Location','best')
Make appropriate changes to get the result you want.
Alex Sha
Alex Sha on 9 Sep 2021
Hi, Pablo, according to the picture of ODE equations you posted, as above, in your code:
dcdt(2)= -((1./theta(4)).*dcdt(1))-(theta(5).*c(1));
be should:
dcdt(2)= -((1./theta(4)).*dcdt(3))-(theta(5).*c(1));

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!