Bifurcation GUI

Question

David Arnold 2012-3-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/31285-bifurcation-gui

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