MATLAB code for solving system of ODE with unknown parameters

Hello, I was wondering if someone can help me write a code in MATLAB for solving system of nonlinear ODE with unknown parameters that need to fit the experimental data. I tried doing it by using < this > , but it didn't work. Or maybe I didn't do it properly because I didn't really understand what was what and how the variables were defined. My problem include this system:

• dP/dt = (b - r0(N))P + ri(N)Q ;
• dQ/dt = r0(N)P - (ri(N) + μ)Q ;
• dD/dt = μQ - dD;

where

• P,Q,D are function in t (time) ;
• N = P + Q + D ;
• r0(N) = kN/(aN + 1);
• ri(N) = r/(N + m);

and a, b, d, k, m, μ and r are unknown parameters (that can be estimated within an interval, for every parameter; b, μ, d are positive). Initial conditions are:

• P(0) = 0.1 (or 0.01);
• Q(0) = D(0) = 0;

My problem includes finding approximative solution (or just in a form that can be graphed) and fitting parameters to experimental data that includes values for N(t) in some period of time (and then making a plot of P, Q, D, N functions and data on the same graph), for example:

    t   N(t)
    4,6  0,0016308;
    5,7  0,0032148;
    6,7  0,0056140;
    7,9  0,0118598;
    8,8  0,0201500;
    9,8  0,0275380;
   11,0  0,0345460;
   12,0  0,0660800;
   12,9  0,0789320;
   15,2  0,1550000;

Can someone explain me how to work this out in MATLAB? Please start from scratch, write me everything I need to do (like how to define P, Q, D and parameters, what goes in what function, how to use toolboxes and such). I am sorry if this is too much from me to ask, but I struggle with this for weeks now.

Thank you very much.

Answer 1

[Alan Weiss]

I am not going to give you a step-by-step tutorial on how to formulate ODE for solution in MATLAB. I will tell you that you have an optimization problem.

Set the vector

 x = [a, b, d, k, m, μ, r]

I mean, a = x(1), b = x(2), etc. Then you can take a guess x(0) as initial parameters, and try to solve the optimization problem of minimizing the sum of squares of differences of your ODE solution from the data. The optimizer will adjust x to try to achieve the minimal difference.

For examples, see Curve Fitting via Optimization, lsqnonlin with a Simulink model, and nonlinear curve fitting with lsqcurvefit.

For solving an ODE given an initial point, see Ordinary Differential Equations (a very long section, but one that has a bunch of examples).

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Ana]

Thank you, but I don't see how I could use curve fitting without explicit solution to the ODE. That is my main problem here. Okay, I know ode45 returns iterations and plots them, but I don't know how to use fminsearch if my unknown parameters are in function that contains my ODE.

This is the ODE:

function dydt = f(t,y)
a = 1.012; b = 0.35; d = 0.02085; k = 1.25; m = 3.42; mi = 0.058; r = 1.2025;
dydt = [ (b - k*(y(1) + y(2) + y(3))/(a*(y(1) + y(2) + y(3)) + 1))*y(1) + r*y(2)/(y(1) + y(2) + y(3) + m)
          k*(y(1) + y(2) + y(3))*y(1)/(a*(y(1) + y(2) + y(3)) + 1) - (r/(y(1) + y(2) + y(3) + m) + mi)*y(2)
          mi*y(2) - d*y(3)];

and I use this to solve it:

[t, y] = ode45(@f,tspan,y0);

where

y0 = [1.0e-003;0;0];
tspan = 4:1.2:100;

Parameters in function dydt are my estimated values I got from manipulating with the solution plot in Wolfram Mathematica.

[Torsten]

In the sequel, exp_t and exp_y are your measurement data.

a = 1.012;

b = 0.35;

d = 0.02085;

k = 1.25;

m = 3.42;

mi = 0.058;

r = 1.2025;

p0(1)=a;

p0(2)=b;

p0(3)=d;

p0(4)=k;

p0(5)=m;

p0(6)=mi;

p0(7)=r;

p_estimate = lsqnonlin(@(p)odefit(exp_t,exp_y,p),p0);

function err = odefit(exp_t,exp_y,p)

a=p(1);

b=p(2);

d=p(3);

k=p(4);

m=p(5);

mi=p(6);

r=p(7);

y0 = [1.0e-003;0;0];

[t,y] = ode45(@(t,y)f(t,y,a,b,d,k,m,mi,r),exp_t,y0)

err = y-exp_y; % compute error between experimental y and fitted y

end

function dydt = f(t,y,a,b,d,k,m,mi,r)

dydt = [ (b - k*(y(1) + y(2) + y(3))/(a*(y(1) + y(2) + y(3)) + 1))*y(1) + r*y(2)/(y(1) + y(2) + y(3) + m) k*(y(1) + y(2) + y(3))*y(1)/(a*(y(1) + y(2) + y(3)) + 1) - (r/(y(1) + y(2) + y(3) + m) + mi)*y(2) mi*y(2) - d*y(3)]; end

[Torsten]

err = y-exp_y

has to be modified since in your model, there are three instead of one y(i)'s.

I leave it to you to make the necessary changes.

Best wishes

Torsten.