Generating the Figures in Chapter 11 of Otto and Day

Clear["Global`*"]

Figure 11.2

r=2.1;K=100;n[1]=((0.1 K)/((-1+r)/r));
For[i=1,i≤50,i++,n[i+1]=n[i]+r n[i] (1-((n[i])/K))]
((-1+r)/r)
ListPlot[{Table[{i,n[i]},{i,1,50,2}],Table[{i,n[i]},{i,2,50,2}],Table[{i,n[i]},{i,1,50,1}]},PlotStyle->{GrayLevel[0],Dashing[{Dot,Dash}]},PlotMarkers->{{★,15},{■,15},""},Joined->{False,False,True}]
Clear[r]

Export["Fig11-2.eps",%];

Figure 11.3

gz=1;gp=1;gzp=.5;c=.25;b=2;a=5;tmax=80;zinit=.002;pinit=.009;
result=NDSolve[{z'[t]==gz (-c z[t]-c ((z[t])^3)+a p[t])-gzp (b p[t]+b ((p[t])^3)),p'[t]==gzp (-c z[t]-c ((z[t])^3)+a p[t])-gp (b p[t]+b ((p[t])^3)),z[0]==zinit,p[0]==pinit},{z,p},{t,0,tmax},MaxSteps->100000]
Plot[Evaluate[{z[t],p[t]}/.result],{t,0,tmax},PlotRange->{-1,1},PlotStyle->{Dashing[{0}],Dashing[{.025}]}]

Export["fig11-2a.eps",%];

f[z_,p_]:=gz (-c z-c (z^3)+a p)-gzp (b p+b (p^3))
g[z_,p_]:=gzp (-c z-c (z^3)+a p)-gp (b p+b (p^3))

Show[Plot[0,{x,0,0.00001},PlotRange->{{-1.5,1.5},{-0.5,0.5}},AspectRatio->1],
ParametricPlot[Evaluate[{z[t],p[t]}/.result],{t,0,tmax},AspectRatio->1,PlotRange->{{-1.5,1.5},{-.5,.5}}],ContourPlot[{f[z,p]==0,g[z,p]==0},{z,-1.5,1.5},{p,-0.5,0.5},PlotRange->{{-1.5,1.5},{-.5,.5}},AxesOrigin->{0,0},Frame->False,ContourStyle->{Dashing[{0.05}]}]
]

Export["fig11-3b.eps",%];

gz=1;gp=1;gzp=.5;c=.25;b=2;a=5;tmax=80;zinit=3;pinit=2;
result=NDSolve[{z'[t]==gz (-c z[t]-c ((z[t])^3)+a p[t])-gzp (b p[t]+b ((p[t])^3)),p'[t]==gzp (-c z[t]-c ((z[t])^3)+a p[t])-gp (b p[t]+b ((p[t])^3)),z[0]==zinit,p[0]==pinit},{z,p},{t,0,tmax},MaxSteps->100000]

Show[Plot[0,{x,0,0.00001},PlotRange->{{-1.5,1.5},{-0.5,0.5}},AspectRatio->1],
ParametricPlot[Evaluate[{z[t],p[t]}/.result],{t,0,tmax},AspectRatio->1,PlotRange->{{-1.5,1.5},{-.5,.5}},PlotPoints->100],
ContourPlot[{f[z,p]==0,g[z,p]==0},{z,-1.5,1.5},{p,-0.5,0.5},PlotRange->{{-1.5,1.5},{-.5,.5}},AxesOrigin->{0,0},Frame->False,ContourStyle->{Dashing[{0.05}]}]
]

Export["fig11-3c.eps",%];

Clear[gz,gp,gzp,c,b,a,tmax,zinit,pinit]

Figure 11.4

rhB[x_]:=x (1-sh)+(1-x)
rhb[x_]:=x+(1-x) (1-sh)
rpA[x_]:=x+(1-x) (1-sp)
rpa[x_]:=x (1-sp)+(1-x)

sh=.65;sp=.5;qinit=.6;pinit=.505;tmax=100;result=NDSolve[{q'[t]==q[t] (1-q[t]) (rhB[p[t]]-rhb[p[t]]),p'[t]==p[t] (1-p[t]) (rpA[q[t]]-rpa[q[t]]),q[0]==qinit,p[0]==pinit},{q,p},{t,0,tmax},MaxSteps->100000000]
Plot[Evaluate[{p[t],q[t]}/.result],{t,0,tmax},PlotRange->{0,1},PlotStyle->{Dashing[{.01}],Dashing[{0}]}]

Export["fig11-4a.eps",%];

ParametricPlot[Evaluate[{q[t],p[t]}/.result],{t,0,tmax},AspectRatio->Automatic,PlotRange->{{0,1},{0,1}}]

Export["fig11-4b.eps",%];

sh=.65;sp=.5;qinit=.9;pinit=.505;tmax=100;result=NDSolve[{q'[t]==q[t] (1-q[t]) (rhB[p[t]]-rhb[p[t]]),p'[t]==p[t] (1-p[t]) (rpA[q[t]]-rpa[q[t]]),q[0]==qinit,p[0]==pinit},{q,p},{t,0,tmax},MaxSteps->100000000]
Plot[Evaluate[{p[t],q[t]}/.result],{t,0,tmax},PlotRange->{0,1},PlotStyle->{Dashing[{.01}],Dashing[{0}]}]

Export["fig11-4c.eps",%];

ParametricPlot[Evaluate[{q[t],p[t]}/.result],{t,0,tmax},AspectRatio->Automatic,PlotRange->{{0,1},{0,1}}]

Export["fig11-4d.eps",%];

Figure 11.5

f[q_,p_]:=a q (1-q) (2 p-1)
g[q_,p_]:=b p (1-p) (2 q-1)

Solve[{f[q,p]==0,g[q,p]==0},{q,p}]

matrix={{D[f[q,p],q],D[f[q,p],p]},{D[g[q,p],q],D[g[q,p],p]}}

Simplify[matrix/.{{q->0,p->0},{q->0,p->1},{q->(1/2),p->(1/2)},{q->1,p->0},{q->1,p->1}}]

Eigenvalues[{{-a,0},{0,-b}}]
Eigenvalues[{{a,0},{0,b}}]
Eigenvalues[{{0,(a/2)},{(b/2),0}}]
Eigenvalues[{{a,0},{0,b}}]
Eigenvalues[{{-a,0},{0,-b}}]

Simplify[D[((q[t])^((-b/a))) (((1-q[t]))^((-b/a))) p[t] (1-p[t]),t]]

or better yet

Simplify[D[(((q[t] (1-q[t])))^b) (((p[t] (1-p[t])))^(-a)),t]]

Simplify[(%)/.{q'[t]->a q[t] (1-q[t]) (2 p[t]-1),p'[t]->b p[t] (1-p[t]) (2 q[t]-1)}]

sh=0.65;sp=0.5;
a=Exp[-sh]-1;b=1-Exp[-sp];
Plot3D[(((q (1-q)))^b) (((p (1-p)))^(-a)),{q,0,1},{p,0,1},PlotPoints->50,PlotRange->{0,1}]
ContourPlot[(((q (1-q)))^b) (((p (1-p)))^(-a)),{q,0,1},{p,0,1},PlotPoints->50,FrameLabel->{"q","p"}]

Figure 11.6

sh=.009;sp=.5;
a=-sh;b=sp;
Plot3D[(((q (1-q)))^b) (((p (1-p)))^(-a)),{q,0,1},{p,0,1},PlotPoints->50,PlotRange->{0,1}]
ContourPlot[(((q (1-q)))^b) (((p (1-p)))^(-a)),{q,0,1},{p,0,1},PlotPoints->50,FrameLabel->{"q","p"}]

Figure 11.7 and 11.8

Clear["Global`*"]

rhB[x_]:=x Exp[-γh-ch]+(1-x) Exp[-ch]
rhb[x_]:=x Exp[-γh]+(1-x) Exp[-γh]
rpA[x_]:=x Exp[-cp]+(1-x) Exp[-cp]
rpa[x_]:=x Exp[-γp]+(1-x)

γh=.3;γp=.5;ch=.2;cp=.4;qinit=.15;pinit=.7;tmax=1000;result=NDSolve[{q'[t]==q[t] (1-q[t]) (rhB[p[t]]-rhb[p[t]]),p'[t]==p[t] (1-p[t]) (rpA[q[t]]-rpa[q[t]]),q[0]==qinit,p[0]==pinit},{q,p},{t,0,tmax},MaxSteps->100000000]
Plot[Evaluate[{p[t],q[t]}/.result],{t,0,tmax},PlotRange->{{0,770},{0,1}},PlotStyle->{Dashing[{.01}],Dashing[{0}]}]

γh=.3;γp=.5;ch=.2;cp=.4;qinit=.65;pinit=.7;tmax=1000;result=NDSolve[{q'[t]==q[t] (1-q[t]) (rhB[p[t]]-rhb[p[t]]),p'[t]==p[t] (1-p[t]) (rpA[q[t]]-rpa[q[t]]),q[0]==qinit,p[0]==pinit},{q,p},{t,0,tmax},MaxSteps->100000000]
Plot[Evaluate[{p[t],q[t]}/.result],{t,0,tmax},PlotRange->{{0,770},{0,1}},PlotStyle->{Dashing[{.01}],Dashing[{0}]}]

γh=.3;γp=.3;ch=.2;cp=.4;qinit=.65;pinit=.7;tmax=1000;result=NDSolve[{q'[t]==q[t] (1-q[t]) (rhB[p[t]]-rhb[p[t]]),p'[t]==p[t] (1-p[t]) (rpA[q[t]]-rpa[q[t]]),q[0]==qinit,p[0]==pinit},{q,p},{t,0,tmax},MaxSteps->100000000]
Plot[Evaluate[{p[t],q[t]}/.result],{t,0,tmax},PlotRange->{{0,100},{0,1}},PlotStyle->{Dashing[{.01}],Dashing[{0}]}]

?Contours

γh=.3;γp=.5;ch=.2;cp=.4;
a=Exp[-γh] (Exp[-ch]-1)
b=Exp[-ch]-Exp[-γh]
c=Exp[-cp]-Exp[-γp]
d=Exp[-cp]-1
Plot3D[(q^(-d)) (((1-q))^c) (p^b) (((1-p))^(-a)),{q,.001,.999},{p,.001,.999},PlotPoints->50]
ContourPlot[(q^(-d)) (((1-q))^c) (p^b) (((1-p))^(-a)),{q,.001,.999},{p,.001,.999},PlotPoints->200,FrameLabel->{"q","p"}]

Export["Fig11-8a.eps",%];

γh=.3;γp=.3;ch=.2;cp=.4;
a=Exp[-γh] (Exp[-ch]-1)
b=Exp[-ch]-Exp[-γh]
c=Exp[-cp]-Exp[-γp]
d=Exp[-cp]-1
Plot3D[(q^(-d)) (((1-q))^c) (p^b) (((1-p))^(-a)),{q,.001,.999},{p,.001,.999},PlotPoints->50]
ContourPlot[(q^(-d)) (((1-q))^c) (p^b) (((1-p))^(-a)),{q,.001,.999},{p,.001,.999},PlotPoints->200,FrameLabel->{"q","p"}]

Export["Fig11-8b.eps",%];

Clear["Global`*"]

f[q_,p_]:=q (1-q) (a p+b (1-p))
g[q_,p_]:=p (1-p) (c q+d (1-q))

Solve[{f[q,p]==0,g[q,p]==0},{q,p}]

matrix={{D[f[q,p],q],D[f[q,p],p]},{D[g[q,p],q],D[g[q,p],p]}}

Simplify[matrix/.{{q->0,p->0},{q->0,p->1},{q->1,p->0},{q->1,p->1},{q->(d/(-c+d)),p->(b/(-a+b))}}]

Eigenvalues[{{b,0},{0,d}}]
Eigenvalues[{{a,0},{0,-d}}]
Eigenvalues[{{-b,0},{0,c}}]
Eigenvalues[{{-a,0},{0,-c}}]
Eigenvalues[{{0,-(((a-b) c d)/(((c-d))^2))},{-((a b (c-d))/(((a-b))^2)),0}}]

Simplify[D[((q[t])^(-d)) (((1-q[t]))^c) ((p[t])^b) (((1-p[t]))^(-a)),t]]

Simplify[(%)/.{q'[t]->q[t] (1-q[t]) (a p[t]+b (1-p[t])),p'[t]->p[t] (1-p[t]) (c q[t]+d (1-q[t]))}]

Here is an example with γh<ch, which like figure 11.8b shows a discontinuity (here at p = 0):

γh=.3;γp=.2;ch=.5;cp=.4;
a=Exp[-γh] (Exp[-ch]-1)
b=Exp[-ch]-Exp[-γh]
c=Exp[-cp]-Exp[-γp]
d=Exp[-cp]-1
Plot3D[(q^(-d)) (((1-q))^c) (p^b) (((1-p))^(-a)),{q,.001,.999},{p,.001,.999},PlotPoints->50]
ContourPlot[(q^(-d)) (((1-q))^c) (p^b) (((1-p))^(-a)),{q,.001,.999},{p,.001,.999},PlotPoints->200,FrameLabel->{"q","p"}]

Figure 11.9

eqn=n Exp[r (1-n/K)]

Simplify[Solve[eqn==n,n]]

eqn/.n->eqn//Simplify

eqnf=(E^(((2 K+(-1-(E^(r-((n r)/K)))) n) r)/K)) ;

Show[
Plot[eqnf/.K->100/.r->1,{n,0,400},PlotRange->{Automatic,{0,3}}],
Plot[eqnf/.K->100/.r->2.4,{n,0,400},PlotRange->{Automatic,{0,3}},PlotStyle->Dashing[{0.02}]]
]