I am plotting the logistic growth model using ode45,But I am confused because I am getting oscillation while I should get a constant line so do you think there is another routine could I use it or I need to change something to get the right plot??

10 次查看（过去 30 天）
Avan Al-Saffar
Avan Al-Saffar 2014-10-22
评论： Matt Tearle 2014-10-23
The model is :
dxdt=N0*x*(1-x/k)
What I am doing is :
tspan=0:0.001:100;
x0=0.1;
[t,x]=ode45('funname',tspan,x0)
figure
plot(t,x)
  3 个评论
显示 1更早的评论
Avan Al-Saffar
Avan Al-Saffar 2014-10-23
function Runlogisticconstant1
p0=0.1;
tspan=0:0.001:100;
K=10;
N0=1;
[t,p]=ode45(@logisticconstant,tspan,p0,[],N0,K);
plot(t,p)
1;
% function dpdt = logisticconstant(t,p,N0,K)
% dpdt= N0*p*(1-p./K);
% end

请先登录，再进行评论。

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

回答（1 个）

Matt Tearle
Matt Tearle 2014-10-22
The oscillation you're seeing is a standard numerical artifact that comes from using an explicit RK method. (Look at the scale of the oscillation -- it's small.) The solution to the ODE settles very quickly to equilibrium, which causes stability issues for a numerical solver (it's very easy to overshoot the solution, then have to come back, overshoot again, etc., which is the behavior you're seeing).
Check the doc for ode45 to see the other solvers that are available. A stiff solver such as ode23s will give you better behavior (as far as the oscillation is concerned, at least).
  1 个评论
Matt Tearle
Matt Tearle 2014-10-23
With the values of K and N0 you gave in your comment, I got pretty nice results with ode23s and ode15s right off the bat. You can do even better by just requiring tighter tolerances:
p0=0.1;
tspan=0:0.001:100;
K=10;
N0=1;
[t,p]=ode23s(@(t,p) N0*p*(1-p./K),tspan,p0,odeset('AbsTol',1e-8,'RelTol',1e-10'));
plot(t,p)
(Takes about 1.5 seconds to run on my machine.)
But tighter tolerances also help with ode45. Using the above settings, the oscillation is still there, but it's barely noticeable.

请先登录，再进行评论。

类别

在 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