MATLAB Answers

Bifurcation GUI

3 views (last 30 days)
David Arnold
David Arnold on 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 Comment
Walter Roberson
Walter Roberson on 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.)

Sign in to comment.

Answers (0)

Tags

Products

Community Treasure Hunt

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

Start Hunting!

Translated by