/////////////////// Scilab script used to plot Figure 4 //////////////////////// //////// the different panels are obtained by changing the parameters ////////// clear xdel(winsid()) vMSY=[]; vTaumsy=[]; vTau0=[]; Var=[]; // qn //varinit=0; //varlim=0.228957; // qp varinit=0.218381;//0.22; varlim=0.5; varpas=0.005; dimvar=round((varlim-varinit)/varpas+1); // parameters r1=1; m=0.1; K=1; e=0.7; g=1.5; // vectors and matrices of parameters s11=-r1/K; s12=-g; s21=e*g; s22=0; r=[r1;-m]; s=[s11,s12 ; s21,s22]; x0=inv(s)*(-r); for var=varinit:varpas:varlim Var=[Var,var]; //q1=var; //q2=0.5; q1=0.1; q2=var; q=[q1;q2]; H=[]; Tau=[]; E1=[]; Xeq=[]; e1=0; xeq=x0; while (xeq(1)>0 & xeq(2)>0) xeq=inv(s)*(q*e1-r); if (xeq(1)<0 | xeq(2)<0) then break end; E1=[E1,e1]; Xeq=[Xeq,xeq]; y=xeq'*q*e1'; H=[H,y]; J=[-g*xeq(2)-e1*q1+r1-2*xeq(1)*r1/K , -g*xeq(1) ; e*g*xeq(2) , -m+e*g*xeq(1)-e1*q2]; eig=spec(J); tau=max(real(eig)); Tau=[Tau,-1/tau]; e1=e1+0.01; end [MSY,indmsy]=max(H); Taumsy=Tau(indmsy); vMSY=[vMSY,MSY]; vTaumsy=[vTaumsy,Taumsy]; vTau0=[vTau0,Tau(1)]; end // linear interpolation y2=interpln([[1:length(vMSY)];vMSY(1:length(vMSY))],1:0.01:length(vMSY)); x2=interpln([[1:length(vTaumsy)];vTaumsy(1:length(vTaumsy))],1:0.01:length(vTaumsy)); // Display clf f = gcf(); vec1=[0.1:((0.9-0.1)/length(x2)):0.9]; vec2=vec1([1:length(vec1)-1]); cmap=[vec2',vec2',vec2']; f.color_map=cmap; drawlater colorbar(min(Var),max(Var)) e = gce(); e.parent.title.text = "$q_P$"; e.parent.title.font_size = 6; e.parent.font_size=3; plot2d([x2;x2],[y2;y2]) // Display and size of points : e = gce(); e.children.mark_mode = "on"; e.children.mark_size_unit = "point"; e.children.mark_size = 15; // Assigning colors to the points: for i = 1:length(e.children) e.children(length(e.children)-i+1).mark_foreground =i; end drawnow ax1=gca(); ax1.font_size=3; ax1.x_label.text="Return time to equilibrium at MSY"; ax1.x_label.font_size=4; ax1.y_label.text="Yield at MSY"; ax1.y_label.font_size=4; // compute maximum yield asymptote //q1=varlim; //q2=0.5; q1=0.1; q2=varinit; q=[q1;q2]; H=[]; E1=[]; Xeq=[]; e1=0; xeq=x0; while (xeq(1)>0 & xeq(2)>0) xeq=inv(s)*(q*e1-r); if (xeq(1)<0 | xeq(2)<0) then break end; E1=[E1,e1]; Xeq=[Xeq,xeq]; y=xeq'*q*e1'; H=[H,y]; e1=e1+0.01; end maxmsy=max(H); // qn // //plot2d([0;x2(length(x2))],[maxmsy,maxmsy]) // qp // plot2d([0;x2(1)],[maxmsy,maxmsy]) ce=gce(); ce.children.line_style=2; ce.children.thickness=2; // set limits ax1.data_bounds=[0,0.14;35,0.21];