Coefficient estimation for a system of coupled ODEs
11 views (last 30 days)
Show older comments
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.
我将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
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.
modelacion
functionmodelacion
functionC=kinetics(theta,t)
% c0=[1;0;0;0];
c0 = theta(9:12);
[T,Cv]=ode45(@DifEq,t,c0);
%
functiondC=DifEq(t,c)
dcdt=zeros(4,1);
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);
dC=dcdt;
end
C=Cv;
end
t=[0
0.928
1.984
2.992
3.952
5.008
5.968
6.976
7.98400000000
8.992
9.952
11.008];
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);
fprintf(1,'\tConstants:\n')
fork1 = 1:length(theta)
fprintf(1,'\t\tTheta(%2d) = %18.12f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure(1)
hd = plot(t, c,'p');
fork1 = 1:尺寸(c, 2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
holdon
hlp = plot(tv, Cfit);
fork1 = 1:尺寸(c, 2)
hlp(k1).Color = CV(k1,:);
end
holdoff
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_{%d}(t)',1:size(c,2)),'Location','best')
end
Make appropriate changes to get the result you want.
.
3 Comments
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));
?
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!