MATLAB Answers

How to estimate coefficients for a system of ODE's?

10 views (last 30 days)
Hi guys, I need some help with how to process the code in matlab. I have set up a system of ODE's, got the initial values and have values for A at specific times and now I need to estimate the values of k1, k2 and k3. I tried several things, but I still don't have a clue on how to process this into matlab, this is what i got:
dAdt = -k1*A*B-k2*A*C-k3*A*D
dBdt = -k1*A*B
dCdt = k1*A*B-k2*A*C
dDdt = k1*A*B+k2*A*C-k3*A*D
dEdt = k3*A*D
dFdt = k2*A*C+k3*A*D
Inital values
A = 20.09; B = 6.96; C = 0; D = 0; E = 0; F = 0
Specific values for A at different times:
numdata = xlsread('dataset.xlsx')
t = numdata(:,1);
A = numdata(:,2);
eg. numdata =
4.5000 15.4000
8.6700 14.2200
12.6700 13.3500
Is there someone that can help me?
  1 Comment
Nisrina Pargustan
Nisrina Pargustan on 21 Nov 2020
hai danny, i got the exact same problem as you. so i tried to run your code but i got the error message. can you please send me your matlab code for this problem?

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 24 Jan 2019
If you are doing parameter estimation of a system of ordinary differential equations, these will likely provide some guidance:
This is usually a fairly straightforward problem.
Star Strider
Star Strider on 25 Jan 2019
Here it is (a version of previous code, this time using ga), attached.
Since all the parameters are positive (by definition of the model), I constrained them to be greater than or equal to zero. Change that if your parameters are allowed to be negative.
Note that you will have to re-state your objective function in terms of the fitness function (in this code, ‘ftns’):
ftns = @(B) norm(A_data - Kinetics(B, time));
You mentioned that your parameters are on the order of 100, so change the options structure to:
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+2,PopSz,Parms), 'MaxGenerations',2E3, 'PlotFcn','gaplotbestf');
That should make it converge faster, however it will likely still take a few minutes. Also, ‘theta’ in my code is ‘B’ in yours.

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!

Translated by