Bifurcation GUI

4 vues (au cours des 30 derniers jours)
David Arnold
David Arnold le 6 Mar 2012
All, I've written a little gui to show the bifurcation that occurs as H increases in the differential equation involving logistic growth and harvesting.
dP/dt = rP(1 - P/K) - H
I'm no pro, so it's a hack. But it works. There are some annoying warnings that float by in the command window, but it works.
I'd love to see some pros give this a try and submit their attempts. It would be a great learning opportunity for me.
Thanks.
David Arnold College of the Redwoods Eureka, CA 95501 http://msemac.redwoods.edu/~darnold/index.php
function bif1
f=figure(...
'Menubar','none',...
'Toolbar','none',...
'Name','Logistic with Harvesting: dP/dt = rP(1 - P/K) - H',...
'NumberTitle','off',...
'Position',[400,300,800,400]);
a1=axes(...
'Parent',f,...
'Units','normalized',...
'Position',[0.05,0.2,0.4,0.7]);
a2=axes(...
'Parent',f,...
'Units','normalized',...
'Position',[0.55,0.2,0.4,0.7]);
s=uicontrol(...
'Style','slider',...
'Min',0,...
'Max',4,...
'Value',0,...
'SliderStep',[.05,.1],...
'Units','normalized',...
'Position',[0.1,0.01,.8,.1],...
'Callback',@sliderCallback);
r=1;
K=10;
H=0;
P=linspace(-5,15,200);
dPdt=r*P.*(1-P/K)-H;
axes(a1)
h=plot(P,dPdt,'b')
line([-1,11],[0,0],'Color','black')
line([0,0],[-2.5,2.5],'Color','black')
m1=line(0,0,'Marker','o',...
'MarkerSize',8,...
'MarkerEdgeColor','red',...
'MarkerFaceColor','white')
m2=line(10,0,'Marker','o',...
'MarkerSize',8,...
'MarkerEdgeColor','red',...
'MarkerFaceColor','red')
axis([-1,11,-2.5,2.5])
xlabel('P')
ylabel('dP/dt')
T1=title(['H = ',num2str(H)])
axes(a2)
line([0,0],[-5,15],...
'Color','black')
tspan=[0,5];
set(a2,'NextPlot','add')
for init=-5:15
[t,y]=ode45(@harvest,tspan,init);
plot(t,y)
end
tspan=[0,-5];
for init=-5:15
[t,y]=ode45(@harvest,tspan,init);
plot(t,y)
end
eq1=line([-5,5],[0,0],...
'LineStyle','--',...
'Color','red');
eq2=line([-5,5],[10,10],...
'LineStyle','-',...
'Color','red');
axis([-5,5,-5,15])
xlabel('t')
ylabel('P')
T2=title(['H = ',num2str(H)])
function sliderCallback(hObject,eventdata)
H=get(s,'Value');
set(T1,'String',['H = ', num2str(H)]);
set(T2,'String',['H = ', num2str(H)]);
ydata=r*P.*(1-P/K)-H;
set(h,'YData',ydata);
tspan=[0,5];
axes(a2)
cla
line([0,0],[-5,15],...
'Color','black')
for init=-5:15
[t,y]=ode45(@harvest,tspan,init);
plot(t,y)
end
tspan=[0,-5];
for init=-5:15
[t,y]=ode45(@harvest,tspan,init);
plot(t,y)
end
D=r^2*K^2-4*r*H*K;
if D>0
m1x=(r*K-sqrt(D))/(2*r);
m2x=(r*K+sqrt(D))/(2*r);
set(m1,'XData',m1x);
set(m2,'XData',m2x);
set(m1,'Visible','on')
set(m2,'Visible','on')
line([-5,5],[(r*K-sqrt(D))/(2*r),(r*K-sqrt(D))/(2*r)],...
'LineStyle','--','Color','red');
line([-5,5],[(r*K+sqrt(D))/(2*r),(r*K+sqrt(D))/(2*r)],...
'LineStyle','-','Color','red');
else
set(m1,'Visible','off');
set(m2,'Visible','off');
end
axis([-5,5,-5,15])
end
function Pprime=harvest(t,P)
Pprime=r*P*(1-P/K)-H;
end
end
  1 commentaire
Walter Roberson
Walter Roberson le 6 Mar 2012
What are the warnings that show up? (Some of us don't have access to MATLAB from home to run experiments with.)

Connectez-vous pour commenter.

Réponses (0)

Catégories

En savoir plus sur MATLAB dans Help Center et File Exchange

Tags

Produits

Community Treasure Hunt

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

Start Hunting!

Translated by