How should I approach solving differential equation as function of two variables r, k?

5 views (last 30 days)
I have a dataset containing two variables, time (hours) and colony size of yeast growth (10^6 cells/mL). I also have the differential equation dy/dt = r*y(t)*(1-y(t)/k) and have this so far:
function dydt = deq(t, y)
dydt = r * y * (1 - (y/k));
end
And my task is to solve the differential equation as a function of r, k and also compute the least squares error between the solution to the DEQ and dataset. For solving the DEQ, I have this so far:
t0 = 0; %initial time
tf = 40; %final time
y0 = 1; %first colony size at t = 0
h = 0.5; %step size?
k = 1; %initial guesses for constants?
r = 0.5;
tempFile = 'yeastGrowth_01.xlsx'; %dataset given
myData = xlsread(tempFile);
tVec = (t0 : h : tf);
[r, k] = ode45(deq, tVec, y0);
%plot(r, k, 'ko')
  1 Comment
Torsten
Torsten on 6 May 2022
Why don't you solve your ODE analytically ? If you can't do it using paper and pencil, use "dsolve".

Sign in to comment.

Answers (2)

Alan Stevens
Alan Stevens on 6 May 2022
You need to call the ode with something like
[t, y] = ode45(@(t,y) deq(t,y,r,k),tVec,y0,r,k);
where you pass values of r and k to deq and get back the corresponding values of t and y.
Use fminsearch to find the values of r and k that minimize the resulting least squares error - see
doc fminsearch
  1 Comment
Evan Olivas
Evan Olivas on 6 May 2022
So I'm looking at the fminsearch help and how do I determine what to set my x0 as? Is it based off the dataset I have?

Sign in to comment.


John D'Errico
John D'Errico on 6 May 2022
While you COULD use an ODE solver, this is something you would ONLY do if the ODE had no analytical solution. And this ODE has a simple solution. That is:
syms r k y(t)
ysol = dsolve(diff(y) == r * y * (1 - (y/k)), y(0) == 1)
ysol = 
ysol = simplify(ysol)
ysol = 
And now you have data that can be used to estimate r and k. As long as you have multiple data points, that will provide an estimate.
Now you could use a tool like the curve fitting toolbox, or you could use a tool from the optimization toolbox, so lsqcurvefit or lsqnonlin, or you could use nlinfit form the stats TB, or you could use a simple optimzer like fminsearch to minimize the sum of squares of errors of your data.
Only if you cannot solve the ODE as I did here would you use a solver and then apply it on top of an ODE solver like ODE45. Why is that something to be avoided if you can avoid it? There are multiple issues I can think of.
  1. The ODE solver itself is not exact. It has tolerances. So the numbers that come out ot ODE45 may only be accurate to within some tolerance. And that means the optimizer does not see the truly correct values, and that means the objective function it is trying to minimize is not exact. It need not even be a differentiable function And that can cause the optimizer to fail.
  2. An ODE solver can sometimes have convergence problems for some sets of parameters. And since an optimizer will send in all sorts of parameters into the ODE solver, you might see all sorts of garbage coming out. ODE45 is a classic example of a solver that will simply fail to converge if the problem is stiff. And some problems can be stiff for some sets of parameters, but not others.
  3. It is inefficient to solve an ODE many times, compared to solving it once, and then using a curve fitting tool on the analytical solution.
Since this is surely homework, and you may have been instructed to use ODE45 to solve the ODE, and then to fit the data to the solution as you want to do, you should do as you were told, IF that is the case. However, you should understand this is not remotely the solution method you SHOULD take for this problem.
  1 Comment
Evan Olivas
Evan Olivas on 6 May 2022
Thank you for the detailed explanation! And yes, I am instructed to use ode45 to complete this task.
If you wouldn't mind with some input on this as well: I'm trying find the least squares error between my solution to the DEQ and my data by using fminsearch but am getting an error upon doing so. This is what I have:
x0 = [t0; tf];
rkBest = fminsearch(deq(t, y, r, k), x0);
There is an error with using multiplication in the dydt = r * y * ( 1 - (y/k)) line saying it is incorrect dimensions for matrix multiplication and to use ' .* ' instead but then I get another error saying that I no longer am passing a function for fminsearch to calculate

Sign in to comment.

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!