ode45 and additional parameters - Growth Equation

5 次查看（过去 30 天）
Victor Yuan
Victor Yuan 2016-12-2
编辑： Victor Yuan 2016-12-2
Hi,
I am trying to solve the differential equation
dP/dt = (b*(P(t))) - (k*(P(t)).^2)
which describes the behavior of a population P(t) over time given an initial population. My goal is to predict the future population given previous populations. Since my equation has more then 2 parameters (t, b, k, and the array P which contains all the previous population counts over the years) I can't just use the regular ode45 since it only allows t and y. Having looked at the example located here:https://www.mathworks.com/help/matlab/ref/ode45.html I have produced the code below:
tspan = [30 55]; %I have data from years 1-30, I want to predict from 30 to 55
y0 = [b 4683139]; %I am not sure this is correct here. More on this below
[t,y] = ode45(@(t,y) dPdt(t,y,k,P),tspan,y0);
and my function dPdt is here:
function f = dPdt(t,y,A,B)
f = zeros(2,1);
f(1) = y(2);
first = y*P(t);
second = k*P(t);
f(2) = first - second.^2;
end
However, when I run this code, I the errors:
In an assignment A(:) = B, the number of elements in A and B must be the same.
Error in dPdt (line 6)
f(2) = first - second.^2;
Error in A8>@(t,y)dPdt(t,y,k_est,C)
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 A8 (line 34)
[t,y] = ode45(@(t,y) dPdt(t,y,k_est,P),tspan,y0);
The only two things I think that may be causing this error is me attempting to pass the matrix P as a parameter or that I am understanding y0 incorrectly. I have y0 = [b 4683139] because
  1. 4683139 is the last recorded population (my starting one in the prediction)
  2. I am trying to pass b as a constant parameter in the equation
Additionally, if my approach to this equation is completely off, can someone give me a hint in how to tackle it?
Thanks in advance.
Victor Yuan

请先登录，再回答此问题。

采纳的回答

Mischa Kim
Mischa Kim 2016-12-2
编辑：Mischa Kim 2016-12-2
Victor, check out the following:
function myODE()
tspan = [0 25]; %
y0 = [1 1]; % these are the initial conditions; essentially you need
% P(30) and P'(30) in here
opts = []; % here you could set integration parameters
b = -1; % replace accordingly
k = -1; % replace accordingly
param = [b k];
[T,Y] = ode45(@dPdt,tspan,y0,opts,param);
plot(T,Y)
end
function dP = dPdt(~,P,param)
b = param(1);
k = param(2);
dP = [P(2); b*P(1) - k*P(1).^2];
end
  1 个评论
Victor Yuan
Victor Yuan 2016-12-2
编辑：Victor Yuan 2016-12-2
Thanks so much for the help, but if you wouldn't mind, I have some further questions.
  1. I don't need to have myODE as a function do I? I can just run the code normally correct?
  2. Could you explain how you passed the additional parameters in param? I don't see that syntax on the ode45 page.
Thanks for all of your help.
Victor Yuan

请先登录，再进行评论。

更多回答（2 个）

Torsten
Torsten 2016-12-2
What you do above does not make much sense.
The way you will have to proceed is to fit parameters b and k to the data from 1-30 and then use these parameters to solve your differential equation from 30-55.
Best wishes
Torsten.
balandong
balandong 2016-12-2
Is this for research purpose, I would suggest that you used structure array to carry your constant. Structure array are more easy to deal as well as help in term of readability.
What you can do is define a function, 1. function s = structure_array () s.constant_A = SOMEVALUE; s.constant_B = SOMEVALUE;.
Then, in your main function, call the structure array s = structure_array ();
then pass this parameter into your ODE 2. [T,Y] = ode45(@dPdt,tspan,y0,opts, s)
3. function dP = dPdt(~,P,s) dP = [P(2); s.constant_A*P(1) - s.constant_A*P(1).^2];
I simplified thing, but hope you get the idea. But take home message, structure array help you allot whenever you have lots of constant to deal.

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by