ode45 and additional parameters - Growth Equation

12 Ansichten (letzte 30 Tage)
Victor Yuan
Victor Yuan am 2 Dez. 2016
Bearbeitet: Victor Yuan am 2 Dez. 2016
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

Akzeptierte Antwort

Mischa Kim
Mischa Kim am 2 Dez. 2016
Bearbeitet: Mischa Kim am 2 Dez. 2016
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 Kommentar
Victor Yuan
Victor Yuan am 2 Dez. 2016
Bearbeitet: Victor Yuan am 2 Dez. 2016
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

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (2)

Torsten
Torsten am 2 Dez. 2016
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 am 2 Dez. 2016
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.

Community Treasure Hunt

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

Start Hunting!

Translated by