data fitting starting from a coupled system of differential equations

5 views (last 30 days)
simona on 4 Apr 2012
Commented: Arbol on 5 Jun 2017
I have a coupled system of differential equations, which I defined in the following way:
dx= zeros(2,1);
dx(1)=K1*K2*x(1)*x(2)./(K3+x(2));
dx(2)= -K1*x(1)*x(2)./(K3+x(2));
I have some experimental data for x(1) at a given time grid. (I have about 30 values).
I would like to find the values of the parameters K1,K2,K3 such that the solution of the ODE fits these experimental data.
How can i program this in matlab?? I'd like to use a least square minimisation to get my parameters... how can I use the function ''lsqcurvefit'' in the optimisation toolbox starting from a coupled system of differential equations??
thank you very much in advance.
simona

Richard Brown on 5 Apr 2012
You're not going to be able to use a curve-fitting tool to find your parameters, you're solving an optimisation problem. Essentially you're interested in minimising an objective (e.g. sum of squared error) given input K1, K2, and K3. Here's approximately what I'd do as a first cut:
1. Write a function that for specific K1, K2, K3, solves the ode on your time grid (specify the grid points as the TSPAN option to ode45/23/23s/etc to get the output interpolated to those values) and evaluates your objective (sum of squared deviations at the time points)
2. Find an approximate K1, K2, K3 by trial and error that give you a somewhat correct answer (by eye). You probably want to be graphing your solutions vs your data here.
3. Use something like fminsearch to optimise your function from 1., starting at your approximate solution from 2.
simona on 14 Apr 2012
Thank you very much Mr. Brown! can I ask you something more? I created an M-file function such as:
function [ dx ] = mymodel( t,x )
global K1;
global K2;
global K3;
dx= zeros(2,1);
dx(1)=K1*K2*x(1)*x(2)./(K3+x(2));
dx(2)= -K1*x(1)*x(2)./(K3+x(2));
dx=[dx(1);dx(2)];
end
then i'll have to use this function in the input arguments of ode45:
[t x]=ode45(@mymodel,[timegriddata],[initial conditions])
can you explain me something more about the first point you wrote?
I have to write an M-file function that for specific K1,K2,K3 BOTH solves the ode AND evaluates the sum of squared deviations?
How can i programme this?

owr on 4 Apr 2012
I would start by figuring out how to simulate the system for some fixed values of K1, K2, etc. and sample this numerical solution at exactly the points in time where you have your experimental data sampled. Check out ode45 or something similar.

Richard Brown on 15 Apr 2012
Hi Simona. I'm adding this as a new answer so that I can use code markup. You need a function that takes in K1, K2, K3 as inputs, and returns all the deviations (if you're using lsqnonlin), or the sum of the squared deviations (if you're using something like fminsearch). It may also be useful to have a flag that you can set if you want graphical output.
You can do it fairly neatly using nested functions. I'll give you the overall structure, I'm sure you can fill in the details. T and X are your data values and times.
function r = foo(K1, K2, K3, T, X, is_plotting)
x0 = blah;
[~, Y] = ode45(@myode, T, x0)
if is_plotting
plot(T, Y, 'b', T, X, 'rx');
end
% The format of this line depends on what solver you use
r = sum(X - Y).^2;
function dx = myode(t, x)
% Because myode is nested, it can "see" K1, K2, and K3
dx = [K1*K2*x(1)*x(2)./(K3+x(2));
-K1*x(1)*x(2)./(K3+x(2))];
end
end
end
Arbol on 5 Jun 2017
you can use lsqcurvefit for this to estimate the parameters correct?