data fitting starting from a coupled system of differential equations

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

Answer 1

[Richard Brown]

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]

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?

thank you in advance!

Answer 2

[owr]

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.

Answer 3

[Richard Brown]

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]

you can use lsqcurvefit for this to estimate the parameters correct?