Solve a system of equations

So I have written a system of equations and used ode45 to solve it. I was just wondering if there is a more efficient way to do it. I am creating an ODE model and will later use certain methods to find the unknown parameters, but for now I am just guessing random values.

I have written my model as

 function dydt = ODE_Model(t,y,p_1,p_2,p_3,d_1,d_2,d_3,a,b,Tau)
dydt = zeros(3,1); 
dydt(1) = p_1*(a*cos(t*2*pi/365-Tau)+b)-d_1*y(1); 
dydt(2) = p_2*(a*cos(t*2*pi/365-Tau)+b)-d_2*y(2); 
dydt(3) = p_3*y(1)-d_3*y(2)*y(3); 
end

and then solved it using

    tspan = [0 5];
y0 = [4 5 6];
p_1 = 1;
p_2 = 1;
p_3 = 1;
d_1 = 1;
d_2 = 1;
d_3 = 1;
a=1;
b=1;
Tau=1;
SolExp = ode45(@(t,y) ODE_Model(t,y,p_1,p_2,p_3,d_1,d_2,d_3,a,b,Tau), tspan, y0);
TimePoints = 0:0.1:2;
y = deval(SolExp, TimePoints);
clf
plot(TimePoints, y)

I wanted to do P and D as a combined vector, so i would have P(1) .. P(6) and then i could just write something like P = [1 1 1 1 1 1], but reading about ode45 it seems I can only use one variable? I am not sure. Is there a neater way to do what I want to do, or is this as neat it gets?

Thanks

Answer 1

[Star Strider]

These may provide you with an approach appropriate to your problem:

• Monod kinetics and curve fitting (link)
• Solving Coupled Differential Equations (link)
• Parameter Estimation for a System of Differential Equations (link)

All of these involve estimating the parameters of a system of differential equations, using lsqcurvefit to do the regressions.

[Jacob Barrett-Newton]

I have had a look over the links you provided, and I am struggling to see how they solve my problem. Perhaps I didn't explain it correctly, or I just don't understand how it links. p_i and d_i are parameters which for the mean time, I am going to be putting in random values. So as it is, I can just put p_1 = 1, p_2 = 1, etc

But i would like to be able to write p as a vector, pretty much just for aesthetic reasons. E.g, p = [1 1 1 1 1 1], and this would be all my p_i and d_i values.

I think, the links you provided me will help me with my next step where I will be estimating the parameters. Although I could be wrong and your links do actually solve my current problem and I just cant see it.

Thanks in advance

[Star Strider]

‘I am creating an ODE model and will later use certain methods to find the unknown parameter ...’

I am not certain what you are currently doing. My intent is to provide you with help to fit or optimise your ODE, in the event you want to fit your system of differential equations to data, or oprimise the parameters with respect to certain criteria.

If you want to optimize your differential equations instead, see: Optimizing a Simulation or Ordinary Differential Equation (link).

If you want to create ‘p_1’ and the others as a vector, just create them as such: ‘p(1)’, ‘p(2)’, and so for the others, including ‘d_1’ assigned as ‘p(4)’, and pass it and the others as one variable (perhaps naming it ‘pd’) to your ODE function. You can re-assign them inside the function, or just keep track of them. This is the best way to pass such parameters, and the only way to work with them for optimization and parameter estimation.

The ODE solvers, such as ode45, are not the problem. They require that the function arguments be a vector of the independent and dependent variables, in that order. However you can easily pass other parameters to your functions. See: Parameterizing Functions (link). This is also what my earlier links are intended to demonstrate, since I do that extensively in them.

[Jacob Barrett-Newton]

I do want to fit my system of equations to some data, that is my next step. I have some data, and I have plotted a graph and eventually i will want to estimate the parameters to fit that data, which is what I think you are saying those links will help me to do if i have understood it correctly.

I have tried to use p(1) etc, but i get an error. I have written my ode as:

function dydt = prac(t,y,p,a,b,Tau)
dydt = zeros(3,1); 
dydt(1) = p(1)*(a*cos(t*2*pi/365-Tau)+b)-p(4)*y(1); 
dydt(2) = p(2)*(a*cos(t*2*pi/365-Tau)+b)-p(5)*y(2); 
dydt(3) = p(3)*y(1)-p(6)*y(2)*y(3); 
end

And then tried to solve:

tspan = [0:5];
y0 = [0 1 1];
p0 = [1 1 1 1 1 1];
a=3;
b=9;
Tau=4;
SolExp = ode45(@(t,y,p) prac(t,y,p,a,b,Tau), tspan, y0);
TimePoints = 0:0.1:2;
y = deval(SolExp, TimePoints);
clf
plot(TimePoints, y)

But i get a series of errors:

Not enough input arguments.

Error in untitled7>@(t,y,p)prac(t,y,p,a,b,Tau)

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in untitled7 (line 7)
SolExp = ode45(@(t,y,p) prac(t,y,p,a,b,Tau), tspan, y0, p0);

I am sure it is something really simple, but i am quite new to matlab and not great at debugging my errors.

[Star Strider]

The ODE solvers only want the independent and dependent variables to be arguments to your function.

Do this instead:

[T,Y] = ode45(@(t,y) prac(t,y,p,a,b,Tau), tspan, y0);

The anonymous function will pick up the other parameters from your workspace.

If you are going to fit ‘prac’ to data, you need to return the independent variable and the integrated dependent variables vectors, as I did here. Using deval is not recommended in this instance, and will only slow your code.

If you are fitting only one column of ‘Y’ (as I have defined it here) to your data, you need to define it and return it. The ‘Monod Kinetics’ solution gives an example of that.

[Jacob Barrett-Newton]

I still get errors using that. Undefined function or variable p.

[Star Strider]

One option is to change this:

p0 = [1 1 1 1 1 1];

to:

p = [1 1 1 1 1 1];

(I thought I submitted much the same Comment hours ago. It’s not appearing, so this may be a duplicate.)

[Jacob Barrett-Newton]

Nope still doesnt work

[Star Strider]

I cannot reproduce the problem.

This runs for me without error:

      function dydt = prac(t,y,p,a,b,Tau)
          dydt = zeros(3,1);
          dydt(1) = p(1)*(a*cos(t*2*pi/365-Tau)+b)-p(4)*y(1);
          dydt(2) = p(2)*(a*cos(t*2*pi/365-Tau)+b)-p(5)*y(2);
          dydt(3) = p(3)*y(1)-p(6)*y(2)*y(3);
      end

tspan = [0 5];
y0 = [4 5 6];
p = [1 1 1 1 1 1];
a=1;
b=1;
Tau=1;
[T,Y] = ode45(@(t,y) prac(t,y,p,a,b,Tau), tspan, y0);
TimePoints = 0:0.1:2;
% y = deval(SolExp, TimePoints);
plot(T, Y)
grid
legend('y_1(t)', 'y_2(t)', 'y_3(t)')

[Jacob Barrett-Newton]

Yes I just reran my code and it worked. Not sure what the error was, but thank you it looks much more pleasing now.

I am now onto my second problem which is what you initially suggested with those links. I have looked through them, recreated the problems to try and understand them better and I am now looking to apply them to my problem. Apologies for asking so many questions, but I started using matlab a few weeks ago and my supervisor is on strike for weeks so there's not a lot I can do.

My problem is simple, I have a set of data taken over a period of time, with repeats on each taking. I have calculated the means, of all this data, and plotted them. I now have the model that i described above, and wish to estimate the parameters to fit the data that I have. However, from reading through the other links it seems that in those problems some data was already known, such as where the parameter should lie within and what not, but I do not have any of this data. I am simply trying to guess the parameters to fit the graph. Is this similar to those other problems?

[Star Strider]

It seems to be similar. The point is that if you believe that guessing the parameters could eventually fit your data, then a more sophisticated and less stochastic approach will likely work.

I would have to have a representative sample and description of your data, what results of your system of differential equations you are fitting to it (for example, the Igor Moura code fits all the results of the integration to the data, while the Monod Kinetics code fits one result), and a description of what you want to do.

The optimization technique you choose to estimate your parameters can be anything. I have use everything from fminsearch to ga (genetic algorithm), depending on the context.

So it is probably possible to estimate your parameters. I can’t say for certain until I see your data, and know what you want to do.

Answer 2

[Jacob Barrett-Newton]

Although this is causing me a lot of stress, I am a bit hesitant to upload my data and for you to write a code that does it for me like you have done in other questions, as this is for my dissertation and I don't want to inadvertently get done for plagiarism, I'm not really sure how the rules around it work.

But perhaps you can help me with a simialr problem I am having on the same topic. I am trying to make another model that models another aspect of the same data. It uses time dependent variables (I think) and so I have tried to follow the help in ode45 but haven't had much luck.

So far I have for my function:

function dydt = model2(t,y,p,ft,f,gt,g)
f = interp1(ft, f, t);
g = interp1(gt, g, t); 
dydt = p(1)*f-p(2)*g*y(1);

and I then try to solve using:

ft = MeansR{8,:};
f =ft
gt = MeansR{11,:};
g = gt
TSPAN = [1 5];
IC = 1; 
p = [1 1];
[T Y] = ode45(@(t,y) model2(t, y,p, ft, f, gt, g), TSPAN, IC);
plot(T, Y);
grid

But the plot is blank. The output variable Y returns a column of NaN and I cannot work out why. If i replace MeansR with some other numbers, it works. I have saved mat files of the two Means. Any ideas?

[Star Strider]

This is a guess, since I cannot run your code. The NaN values may be coming from the interp1 calls.

Try this:

f = interp1(ft, f, t, 'linear', 'extrap');
g = interp1(gt, g, t, 'linear', 'extrap');

[Jacob Barrett-Newton]

Amazing, it worked. You have been a great help, I'm sure i'll be back with more questions. I Appreciate all your comments.

What does 'linear' and 'extrap' actually do?

[Star Strider]

Thank you!

The interp1 function by default returns NaN for values interpolated beyond the range of the first argument. (The default interpolation method is 'linear'.) In order for interp1 to extrapolate, it is necessary to supply a method (here 'linear' since I assume that is what you want), and tell it specifically to extrapolate, by adding 'extrap' as the third argument.

Note that the other interpolation functions, such as interp2, deal with extrapolation differently, so investigate those options if you use the other functions.

[Jacob Barrett-Newton]

Are any of those methods similar to MCMC methods?

[Star Strider]

MCMC?

I’m not familiar with this acronym.

[Jacob Barrett-Newton]

Markov chain Monte Carlo

[Star Strider]

Not that I’m aware of. I’ve not looked through the code for the interpolation functions. However they appear to use the stated method to extrapolate to new values.

Extrapolation, especially far from the data, is not generally considered appropriate, since the behavior of the data beyond the region-of-fit is by definition unknown.

[Jacob Barrett-Newton]

Well at the moment I am trying to do something similar to your monod kinetics example. I have a model of 3 equations, and I am trying to fit those 3 equations to 3 sets of data. In theory, what I want to do is somehow feed in a matrix of the 3 sets of data and try to fit each equation, if that makes sense?

At the moment, I am feeding in the data sets for one equation and I think matlab is trying to fit all 3 equations to that one set of data. Is it possible to do what I am proposing?

[Star Strider]

It makes sense, and lsqcurvefit (and ga and some other optimisation routines) can do it.

The ‘Parameter Estimation for a System of Differential Equations’ example is likely more applicable to your problem than the ‘Monod kinetics and curve fitting’ example, since the first fits a matrix dependent variable and the second fits a vector dependent variable.

[Jacob Barrett-Newton]

Yeah I tried using that method, and I made a model that works for me and when i used some simple test values it worked. But now I have put in my actual data and nothing happens. I think it is still running, but seeing as how with test values it was almost instantaneous I'm thinking it's stuck in a loop with my values but I am not sure why

[Jacob Barrett-Newton]

Actually I think it might just be taking a long time because of how big the numbers are. Because i divided all the numbers by 100 and it worked fine

[Star Strider]

Scaling the data can be appropriate, although you have to scale the parameter estimates as well afterwards.

The ode45 function can be the reason the convergence is slow, because you could have a ‘stiff’ system. I would do the regression again with the unscaled numbers, and using ode15s instead. The advantage to using unscaled data and ode15s is that the parameter estimates would then be correct.

[Jacob Barrett-Newton]

Right, so ode15s does a bit more than ode45 does, but now when it runs it goes for maybe 10 seconds, goes through a load of values and then comes to a stop with errors:

dcdt =

     NaN
     NaN
     NaN

> In ode15s (line 589)
  In Igor_Moura/kinetics (line 12)
  In lsqcurvefit/objective (line 264)
  In snls (line 337)
  In lsqncommon (line 166)
  In lsqcurvefit (line 257)
  In Igor_Moura (line 54) 
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
Warning: Failure at t=4.837389e+00.  Unable to meet integration tolerances without reducing the step size
below the smallest value allowed (1.421085e-14) at time t. 
> In ode15s (line 668)
  In Igor_Moura/kinetics (line 12)
  In lsqcurvefit/objective (line 264)
  In snls (line 337)
  In lsqncommon (line 166)
  In lsqcurvefit (line 257)
  In Igor_Moura (line 54) 
Matrix dimensions must agree.

Error in lsqcurvefit/objective (line 265)
         F = F - YDATA;

Error in snls (line 337)
            newfvec = feval(funfcn{3},xcurr,varargin{:});

Error in lsqncommon (line 166)
            snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...

Error in lsqcurvefit (line 257)
    lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,caller,...

Error in Igor_Moura (line 54)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

So it must be something to do with my data? Because i can switch the data out for random numbers and it all works fine, there is only an error when I change data

[Star Strider]

It has to do with your data and probably your initial parameter estimates. And, of course, that you are using a gradient-descent solver.

When I have problems such as yours, I go to a derivative-free solver, my personal favourite being the Genetic Algorithm (the ga function). If you have the Global Optimization Toolbox, try something like the approach I used in How to find parameter estimation with fminsearch? And solve it with ode45? (link). It is not difficult to program your own fitness function, that simply takes your current objective function, compares it to your data, and outputs a scalar result that ga can use. In this problem, I began with a very large initial population, since a large initial population will quite likely will converge on an optimal parameter set. (With simple problems the ga defaults work. I rarely use ga with simple problems.)

Note that ga is robust but not fast, so it is likely best to get your optimization working and then go do something else for the few minutes it will take ga to converge. I also always define an 'OutputFcn' (see options in Input Arguments) to avoid anxiety about not knowing what ga is doing. Note that if you have defined parameter limits or other constraints with lsqcurvefit, you can use them in ga.

[Jacob Barrett-Newton]

Hi, i'm wondering if you can help me with this issue. The code I have attached works perfectly, bar the fact that one line doesn't quite fit to the data, but other than that no problem. However, for whatever reason I am very limited in what I can change values to. For example, If i change y = [1;0;0] to y = [1;1;0], I get the error "Matrix dimensions must agree" and I cannot seem to figure out why this error is occurring.

[Star Strider]

When I ran your code with ‘y0(2)~=0’, the warning is:

Warning: Failure at t=2.354295e+02.  Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (4.547474e-13) at time t.

So it’s encountering a singularity there, and fails to complete the ‘dY’ vector. That’s the reason it throws the

Matrix dimensions must agree.

error in lsqcurvefit. I don’t see any obvious reason for that to occur. I defer to you to troubleshoot it, and be certain your kinetics equations in ‘DifEq’ are correct.