Find kinetic constants from differential equations

Question

Marina Batlló RIus 2023-3-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1933060-find-kinetic-constants-from-differential-equations

评论： Star Strider 2023-3-29

Hey! So a escription on my problem:

I have a compartimental model that contains 4 different elements that have a kinetic behaviour between thm. The differential equations of this model can be described as:

dydt(1) = k(1)*y(4)+k(2)*y(2)+k(3)*y(1)-(k(4)*y(1)+k(5)*y(1)+k(6)*y(1)+1/25.028)*y(1);

dydt(2) = k(2)*y(1)-(k(5)*y(2)+1/25.028)*y(2);

dydt(3) = k(3)*y(1)-(k(4)*y(3)+1/25.028)*y(3);

dydt(4) = k(1)*y(1)-(k(4)*y(4)+1/25.028)*y(4);

I have some experimental data that I want to use to predict th different k, but I have been struggeling with solving them with the ode45 function or similars. Could someone help me?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-3-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1933060-find-kinetic-constants-from-differential-equations#answer_1198175

Some examples: Monod kinetics and curve fitting , Parameter Estimation for a System of Differential Equations , and there are several others.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Star Strider 2023-3-22

在 MATLAB Online 中打开

My pleasure!

The ‘B0’ vector are the initial parameter estimates for lsqcurvefit.

The ‘x0’ vector are the initial conditions to differential equations in the ode45 call

The code in the second link are likely more appropriate for your problem. I have also updated it to:

t=[0.1

0.2

0.4

0.6

0.8

1

1.5

2

3

4

5

6];

c=[0.902 0.06997 0.02463 0.00218

0.8072 0.1353 0.0482 0.008192

0.6757 0.2123 0.0864 0.0289

0.5569 0.2789 0.1063 0.06233

0.4297 0.3292 0.1476 0.09756

0.3774 0.3457 0.1485 0.1255

0.2149 0.3486 0.1821 0.2526

0.141 0.3254 0.194 0.3401

0.04921 0.2445 0.1742 0.5277

0.0178 0.1728 0.1732 0.6323

0.006431 0.1091 0.1137 0.7702

0.002595 0.08301 0.08224 0.835];

theta0 = rand(10,1);

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));

Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.

fprintf(1,'\tRate Constants:\n')

Rate Constants:

for k1 = 1:numel(theta)

fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))

end

Theta( 1) = 0.76485 Theta( 2) = 0.23489 Theta( 3) = 0.20888 Theta( 4) = 0.49193 Theta( 5) = 0.62211 Theta( 6) = 0.00000 Theta( 7) = 0.90288 Theta( 8) = 0.07145 Theta( 9) = 0.02841 Theta(10) = 0.00000

tv = linspace(min(t), max(t));

Cfit = kinetics(theta, tv);

figure

hd = plot(t, c, 'p');

for k1 = 1:size(c,2)

CV(k1,:) = hd(k1).Color;

hd(k1).MarkerFaceColor = CV(k1,:);

end

hold on

hlp = plot(tv, Cfit);

for k1 = 1:size(c,2)

hlp(k1).Color = CV(k1,:);

end

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')

function C=kinetics(theta,t)

c0=theta(end-3:end);

[T,Cv]=ode45(@DifEq,t,c0);

function dC=DifEq(t,c)

dcdt=zeros(4,1);

dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);

dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);

dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);

dcdt(4)= theta(5).*c(2)-theta(6).*c(4);

dC=dcdt;

end

C=Cv;

end

NOTE — In this version, the initial conditions for the differential equationos are parameters to be estimated. The plot colours also match (data and estimates).

.

Marina Batlló RIus 2023-3-29