Solve a system of equations
Show older comments
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
Accepted Answer
Star Strider
on 2 Mar 2018
0 votes
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.
10 Comments
Jacob Barrett-Newton
on 3 Mar 2018
Star Strider
on 3 Mar 2018
‘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
on 3 Mar 2018
Star Strider
on 3 Mar 2018
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
on 3 Mar 2018
Star Strider
on 3 Mar 2018
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
on 3 Mar 2018
Star Strider
on 3 Mar 2018
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
on 4 Mar 2018
Star Strider
on 4 Mar 2018
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.
More Answers (1)
Jacob Barrett-Newton
on 5 Mar 2018
16 Comments
Star Strider
on 5 Mar 2018
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
on 5 Mar 2018
Edited: Jacob Barrett-Newton
on 5 Mar 2018
Star Strider
on 5 Mar 2018
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
on 9 Mar 2018
Star Strider
on 9 Mar 2018
MCMC?
I’m not familiar with this acronym.
Jacob Barrett-Newton
on 9 Mar 2018
Star Strider
on 9 Mar 2018
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
on 11 Mar 2018
Star Strider
on 11 Mar 2018
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
on 11 Mar 2018
Jacob Barrett-Newton
on 11 Mar 2018
Star Strider
on 11 Mar 2018
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
on 11 Mar 2018
Star Strider
on 11 Mar 2018
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
on 23 Mar 2018
Star Strider
on 23 Mar 2018
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.
Categories
Find more on Get Started with Curve Fitting Toolbox in Help Center and File Exchange
on 2 Mar 2018
on 23 Mar 2018
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
