Generating the Figures in Chapter 11 of Otto and Day

Inputs (Enter first)

Off[General::"spell"]

Off[General::"spell1"]

<< Graphics`Graphics`

<< Graphics`Graphics3D`

<< Graphics`MultipleListPlot`

<< Graphics`Arrow`

<< Graphics`PlotField`

<< Graphics`ImplicitPlot`

<< Graphics`MultipleListPlot`

Figure 11.2

r = 2.1; K = 100; 
    n[1] =  .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)
MultipleListPlot[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}]}, 
      SymbolShape -> {PlotSymbol[Star], PlotSymbol[Box]}, 
      SymbolStyle -> {GrayLevel[0], GrayLevel[ .5], GrayLevel[ .5]}, 
      PlotJoined -> {False, False, True}]
Clear[r]

Display["graph1.eps", %, "EPS"];

Figure 11.3

<< Graphics`ImplicitPlot`

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}]}]
ParametricPlot[Evaluate[{z[t], p[t]} /. result], {t, 0, tmax}, 
      AspectRatio -> 1, 
      PlotRange -> {{-1.5, 
            1.5}, {- .5,  .5}}]
Clear[gz, gp, gzp, c, b, a, tmax, zinit, pinit]

Display["f2a.eps", %88, "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))

gz = 1; gp = 1; gzp =  .5; c =  .25; b = 2; a = 5; tmax = 80; 
    zinit =  .002; pinit =  .009;
ImplicitPlot[{f[z, p] == 0, g[z, p] == 0}, {z, -1.5, 
        1.5}, AspectRatio -> 1, 
      PlotRange -> {{-1.5, 1.5}, {- .5,  .5}}, 
      PlotStyle -> {Dashing[{ .025}]}]
Clear[gz, gp, gzp, c, b, a, tmax, zinit, pinit]

Show[%58, %44, 
      PlotRange -> {{-1.5, 1.5}, {- .5,  .5}}]

Display["f2b.eps", %92, "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]
ParametricPlot[Evaluate[{z[t], p[t]} /. result], {t, 0, tmax}, 
      AspectRatio -> 1, 
      PlotRange -> {{-1.5, 1.5}, {- .5,  .5}}, 
      PlotPoints -> 100]
Clear[gz, gp, gzp, c, b, a, tmax, zinit, pinit]

Show[%76, %44, 
      PlotRange -> {{-1.5, 1.5}, {- .5,  .5}}]

Display["f2c.eps", %78, "EPS"];

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}]}]
ParametricPlot[Evaluate[{q[t], p[t]} /. result], {t, 0, tmax}, 
      AspectRatio -> Automatic, 
      PlotRange -> {{0, 1}, {0, 1}}]

Display["1.eps", %13, "EPS"]

Display["2.eps", %14, "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}]}]
ParametricPlot[Evaluate[{q[t], p[t]} /. result], {t, 0, tmax}, 
      AspectRatio -> Automatic, 
      PlotRange -> {{0, 1}, {0, 1}}]

Display["3.eps", %16, "EPS"]
Display["4.eps", %17, "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)}]

Simulation

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}]}]
ParametricPlot[Evaluate[{q[t], p[t]} /. result], {t, 0, tmax}, 
      AspectRatio -> Automatic, 
      PlotRange -> {{0, 1}, {0, 1}}]

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}]}]
ParametricPlot[Evaluate[{q[t], p[t]} /. result], {t, 0, tmax}, 
      AspectRatio -> Automatic, 
      PlotRange -> {{0, 1}, {0, 1}}]

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]

Figure 11.7 and 11.8

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 =  .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, 1}, 
      PlotStyle -> {Dashing[{ .01}], Dashing[{0}]}]
ParametricPlot[Evaluate[{q[t], p[t]} /. result], {t, 0, tmax}, 
      AspectRatio -> Automatic, 
      PlotRange -> {{0, 1}, {0, 1}}]

?Contours

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]

Export["Fig.eps", %];

γ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, 1}, 
      PlotStyle -> {Dashing[{ .01}], 
          Dashing[{0}]}]
ParametricPlot[Evaluate[{q[t], p[t]} /. result], {t, 0, tmax}, 
      AspectRatio -> Automatic, 
      PlotRange -> {{0, 1}, {0, 1}}]

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]

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]))}]

γ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]

Figure 11.9

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

SP[Solve[eqn == n, n]]

eqn /. n -> eqn // Simplify

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

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