Needs["Graphics`PlotField`"]

Question 12

(i) (a)

Solve[{-x + ((x y)/100) == 0, 2  y - ((2  x y)/25) == 0}, {x, 
        y}]

(ii) (a)

fx1 = D[-x + ((x y)/100), x] /. {x -> 0, y -> 0}

fy1 = D[-x + ((x y)/100), y] /. {x -> 0, y -> 0}

gx1 = D[2  y - ((2  x y)/25), x] /. {x -> 0, y -> 0}

gy1 = D[2  y - ((2  x y)/25), y] /. {x -> 0, y -> 0}

matrixA1=(({{-1, 0}, {0, 2}}))

fx1 = D[-x + ((x y)/100), x] /. {x -> 25, y -> 100}

fy1 = D[-x + ((x y)/100), y] /. {x -> 25, y -> 100}

gx1 = D[2  y - ((2  x y)/25), x] /. {x -> 25, y -> 100}

gy1 = D[2  y - ((2  x y)/25), y] /. {x -> 25, y -> 100}

matrixA2=(({{0, (1/4)}, {-8, 0}}))

(iii) (a)

Eigenvalues[matrixA1]

Eigenvectors[matrixA1]

Eigenvalues[matrixA2]

Eigenvectors[matrixA2]

(iv) (a)

Note
The following procedure for plotting multiple trajectories along with the direction field is adapted from Commbes, K.R. et al (2nd ed 1998) Differential Equations with Mathematica. John Wiley & Sons, Inc., (Revised for Mathematica 3.0), Chapter 12.

nsol1[tval_, a_, b_, 
        t1_] := ({x[t], y[t]} /. 
            First[NDSolve[{x'[
                      t] == -x[t] + ((x[t] y[t])/100), y'[t] == 
                    2  y[t] - ((2  x[t] y[t])/25), x[0] == a, 
                  y[0] == b}, {x[t], y[t]}, {t, 0, t1}]]) /. 
        t -> tval
nphase1[t1_] := 
      ParametricPlot[
        Evaluate[
          Flatten[Table[
              nsol1[t, a, b, t1], {a, 0, 100, 50}, {b, 0, 300, 50}], 1]], {t, 
          0, t1}, AspectRatio -> 1, 
        AxesLabel -> {"x", "y"}]

lines1 = nphase1[5];

arrows1 = 
        PlotVectorField[{-x + ((x y)/100), 2  y - ((2  x y)/25)}, {x, 
            0, 120}, {y, 0, 500}, Axes -> True, ScaleFunction -> (1 &), 
          AspectRatio -> 0.8];

Show[arrows1, lines1, PlotRange -> {0, 500}, 
        AxesLabel -> {"x", "y"}];

Which shows oscillations around the fixed point (x, y) = (25, 100) .

(i) (b)

Solve[{-(1/2) x + ((x y)/((1/4) + y)) == 0, 
        y - (y^2) - ((x y)/((1/4) + y)) == 0}, {x, y}]

(ii) (b)

fx1 = 
      D[-(x/2) + ((x y)/((1/4) + y)), x]  /. {x -> 0, 
          y -> 0}

fy1 = 
      D[-(x/2) + ((x y)/((1/4) + y)), y]  /. {x -> 0, 
          y -> 0}

gx1 = 
      D[y - (y^2) - ((x y)/((1/4) + y)), x] /. {x -> 0, y -> 0}

gy1 = 
      D[y - (y^2) - ((x y)/((1/4) + y)), y] /. {x -> 0, y -> 0}

matrixB1=(({{-(1/2), 0}, {0, 1}}))

fx1 = 
      D[-(x/2) + ((x y)/((1/4) + y)), x]  /. {x -> 0, 
          y -> 1}

fy1 = 
      D[-(x/2) + ((x y)/((1/4) + y)), y]  /. {x -> 0, 
          y -> 1}

gx1 = 
      D[y - (y^2) - ((x y)/((1/4) + y)), x] /. {x -> 0, y -> 1}

gy1 = 
      D[y - (y^2) - ((x y)/((1/4) + y)), y] /. {x -> 0, y -> 1}

matrixB2=(({{(3/10), 0}, {-(4/5), -1}}))

fx1 = 
      D[-(x/2) + ((x y)/((1/4) + y)), x]  /. {x -> (3/8), 
          y -> (1/4)}

fy1 = 
      D[-(x/2) + ((x y)/((1/4) + y)), y]  /. {x -> (3/8), 
          y -> (1/4)}

gx1 = 
      D[y - (y^2) - ((x y)/((1/4) + y)), x] /. {x -> (3/8), 
          y -> (1/4)}

gy1 = 
      D[y - (y^2) - ((x y)/((1/4) + y)), y] /. {x -> (3/8), 
          y -> (1/4)}

matrixB3=(({{0, (3/8)}, {-(1/(2 )), (1/8)}}))

(iii) (b)

Eigenvalues[matrixB1]

Eigenvectors[matrixB1]

Eigenvalues[matrixB2]

Eigenvectors[matrixB2]

Eigenvalues[matrixB3]

Eigenvectors[matrixB3]

(iv) (b)

nsol2[tval_, a_, b_, 
        t1_] := ({x[t], y[t]} /. 
            First[NDSolve[{x'[
                      t] == -(1/2) 
                        x[t] + ((x[t] y[t])/((1/4) + y[t])), y'[t] == 
                    y[t] - y[t]^2 - ((x[t] y[t])/((1/4) + y[t])), 
                  x[0] == a, y[0] == b}, {x[t], y[t]}, {t, 0, t1}]]) /. 
        t -> tval
nphase2[t1_] := 
      ParametricPlot[
        Evaluate[
          Flatten[Table[
              nsol2[t, a, b, t1], {a, 0, 1,  .25}, {b, 0, 1,  .25}], 1]], {t, 
          0, t1}, AspectRatio -> 1, 
        AxesLabel -> {"x", "y"}]

lines2 = nphase2[5];

arrows2 = 
        PlotVectorField[{-(1/2) x + ((x y)/((1/4) + y)), 
            y - (y^2) - ((x y)/((1/4) + y))}, {x, 0, 1.2}, {y, 0, 1}, 
          Axes -> True, ScaleFunction -> (1 &), 
          AspectRatio -> 0.8];

Show[arrows2, lines2, PlotRange -> {0, 1}, 
        AxesLabel -> {"x", "y"}];

Which shows a limit cycle around the fixed point (x, y) = ((3/8), (1/4)) .

(i) (c)

Solve[{x - (x^2) - x y == 0, y - 2  x y  - 2  (y^2) == 0}, {x, 
        y}]

(ii) (c)

fx1 = D[x - (x^2) - x y, x] /. {x -> 0, y -> 0}

fy1 = D[x - (x^2) - x y, y] /. {x -> 0, y -> 0}

gx1 = D[y - 2  x y - 2  (y^2), x] /. {x -> 0, y -> 0}

gy1 = D[y - 2  x y - 2  (y^2), y] /. {x -> 0, y -> 0}

matrixC1=(({{1, 0}, {0, 1}}))

fx1 = D[x - (x^2) - x y, x]  /. {x -> 0, y -> (1/2)}

fy1 = D[x - (x^2) - x y, y]  /. {x -> 0, y -> (1/2)}

gx1 = D[y - 2  x y - 2  (y^2), x]  /. {x -> 0, y -> (1/2)}

gy1 = D[y - 2  x y - 2  (y^2), y]  /. {x -> 0, y -> (1/2)}

matrixC2=(({{(1/2), 0}, {-1, -1}}))

fx1 = D[x - (x^2) - x y, x]  /. {x -> 1, y -> 0}

fy1 = D[x - (x^2) - x y, y]  /. {x -> 1, y -> 0}

gx1 = D[y - 2  x y - 2  (y^2), x]  /. {x -> 1, y -> 0}

gy1 = D[y - 2  x y - 2  (y^2), y]  /. {x -> 0, y -> 0}

matrixC3=(({{-1, -1}, {0, 1}}))

(iii) (c)

Eigenvalues[matrixC1]

Eigenvectors[matrixC1]

Eigenvalues[matrixC2]

Eigenvectors[matrixC2]

Eigenvalues[matrixC3]

Eigenvectors[matrixC3]

(iv) (c)

nsol3[tval_, a_, b_, 
        t1_] := ({x[t], y[t]} /. 
            First[NDSolve[{x'[t] == 
                    x[t] - x[t]^2 - x[t] y[t], y'[t] == 
                    y[t] - 2  x[t] y[t] - 2  y[t]^2, x[0] == a, 
                  y[0] == b}, {x[t], y[t]}, {t, 0, t1}]]) /. 
        t -> tval
nphase3[t1_] := 
      ParametricPlot[
        Evaluate[
          Flatten[Table[
              nsol3[t, a, b, t1], {a, 0, 1,  .25}, {b, 0, 1.5,  .25}], 
            1]], {t, 0, t1}, AspectRatio -> 1, 
        AxesLabel -> {"x", "y"}]

lines3 = nphase3[5];

arrows3 = 
        PlotVectorField[{x - (x^2) - x y, y - 2  x y - 2  (y^2)}, {x, 0, 
            1}, {y, 0, 1.5}, Axes -> True, ScaleFunction -> (1 &), 
          AspectRatio -> 0.8];

Show[arrows3, lines3, PlotRange -> {0, 1.5}, 
        AxesLabel -> {"x", "y"}];

Which shows competing predators at the fixed point (x, y) = (1, 0) .

Question 13

(i)

Critical point found by solving

Solve[{0 == 1.4 (1 - y) x, 0 == 0.6 (1 - 4  y + x) y}, {x, 
        y}]

(ii)

fx1 = D[1.4 (1 - y) x, x]  /. {x -> 0, y -> 0}

fy1 = D[1.4 (1 - y) x, y]  /. {x -> 0, y -> 0}

gx1 = D[0.6 (1 - 4  y + x) y, x]  /. {x -> 0, y -> 0}

gy1 = D[0.6 (1 - 4  y + x) y, y]  /. {x -> 0, y -> 0}

matrixA1=(({{1.4, 0}, {0, 0.6}}))

fx1 = D[1.4 (1 - y) x, x]  /. {x -> 0, y -> 0.25}

fy1 = D[1.4 (1 - y) x, y]  /. {x -> 0, y -> 0.25}

gx1 = 
      D[0.6 (1 - 4  y + x) y, x]  /. {x -> 0, y -> 0.25}

gy1 = 
      D[0.6 (1 - 4  y + x) y, y]  /. {x -> 0, y -> 0.25}

matrixA2=(({{1.05, 0}, {0.15, -0.6}}))

fx1 = D[1.4 (1 - y) x, x]  /. {x -> 3, y -> 1}

fy1 = D[1.4 (1 - y) x, y]  /. {x -> 3, y -> 1}

gx1 = D[0.6 (1 - 4  y + x) y, x]  /. {x -> 3, y -> 1}

gy1 = D[0.6 (1 - 4  y + x) y, y]  /. {x -> 3, y -> 1}

matrixA3=(({{0, -4.2}, {0.6, -2.4}}))

(iii)

Eigenvalues[matrixA1]

Eigenvectors[matrixA1]

Eigenvalues[matrixA2]

Eigenvectors[matrixA2]

Eigenvalues[matrixA3]

Eigenvectors[matrixA3]

nsol13[tval_, a_, b_, 
        t1_] := ({x[t], y[t]} /. 
            First[NDSolve[{x'[t] == 1.4 (1 - y[t]) x[t], y'[t] == 
                    0.6 (1 - 4  y[t] + x[t]) y[t], x[0] == a, 
                  y[0] == b}, {x[t], y[t]}, {t, 0, t1}]]) /. 
        t -> tval
nphase13[t1_] := 
      ParametricPlot[
        Evaluate[
          Flatten[Table[nsol13[t, a, b, t1], {a, 0, 6,  .5}, {b, 0, 2,  .5}], 
            1]], {t, 0, t1}, AspectRatio -> 1, 
        AxesLabel -> {"x", "y"}]

lines13 = nphase13[5];

arrows13 = 
        PlotVectorField[{1.4 (1 - y) x, 0.6 (1 - 4  y + x) y}, {x, 0, 
            6}, {y, 0, 2}, Axes -> True, ScaleFunction -> (1 &), 
          AspectRatio -> 0.8];

Show[arrows13, lines13, PlotRange -> {0, 2}, 
        AxesLabel -> {"x", "y"}];

Illustrates that the system is dominated by equilibrium point (x, y) = (3, 1) .

(iv)

A spreadsheet investigation of the system around point (x, y) = (3, 1) reveals quite a different picture from the one above. There are two reasons for this. First, in the above analysis we have assumed the system is continuous. Discrete systems do not always behave in the same way that continuous systems behave, especially when they are non-linear. Second, the above stability was constructed from a linear approximation in the neighbourhood of the point  (x, y) = (3, 1). How small such a neighbourhood needs to be is not clear. This example illustrates the danger of using continuous linear approximations for non-linear discrete systems.

Question 14

Subtracting (Subscript[x, t]) from the first equation and (Subscript[y, t]) from the second we obtain
	Δ\(TraditionalForm
      \`x\_t = \(x\_\(t + 1\) - x\_t = 
          0.3  x\_t - 0.3  x\_t\%2 - 0.15  x\_t\ y\_t\)\)
	Δ\(TraditionalForm
      \`y\_t = \(y\_\(t + 1\) - y\_t = 
          0.3  y\_t - 0.3  y\_t\%2 - 0.15  x\_t\ y\_t\)\)

(i)

Solve[{x ((3/10) - (3/10) x - (3/20) y) == 0, 
        y ((3/10) - (3/10) y - (3/20) x) == 0}, {x, y}]

(ii)

fx1 = 
      D[x ((3/10) - (3/10) x - (3/20) y), x] /. {x -> 0, 
          y -> 0}

fy1 = 
      D[x ((3/10) - (3/10) x - (3/20) y), y] /. {x -> 0, 
          y -> 0}

gx1 = 
      D[y ((3/10) - (3/10) y - (3/20) x), x] /. {x -> 0, 
          y -> 0}

gy1 = 
      D[y ((3/10) - (3/10) y - (3/20) x), y] /. {x -> 0, 
          y -> 0}

matrixA1=(({{(3/10), 0}, {0, (3/10)}}))

fx1 = 
      D[x ((3/10) - (3/10) x - (3/20) y), x] /. {x -> 0, 
          y -> 1}

fy1 = 
      D[x ((3/10) - (3/10) x - (3/20) y), y] /. {x -> 0, 
          y -> 1}

gx1 = 
      D[y ((3/10) - (3/10) y - (3/20) x), x] /. {x -> 0, 
          y -> 1}

gy1 = 
      D[y ((3/10) - (3/10) y - (3/20) x), y] /. {x -> 0, 
          y -> 1}

matrixA2=(({{(3/20), 0}, {-(3/20), -(3/10)}}))

fx1 = 
      D[x ((3/10) - (3/10) x - (3/20) y), x] /. {x -> (2/3), 
          y -> (2/3)}

fy1 = 
      D[x ((3/10) - (3/10) x - (3/20) y), y] /. {x -> (2/3), 
          y -> (2/3)}

gx1 = 
      D[y ((3/10) - (3/10) y - (3/20) x), x] /. {x -> (2/3), 
          y -> (2/3)}

gy1 = 
      D[y ((3/10) - (3/10) y - (3/20) x), y] /. {x -> (2/3), 
          y -> (2/3)}

matrixA3=(({{-(1/5), -(1/10)}, {-(1/10), -(1/5)}}))

fx1 = 
      D[x ((3/10) - (3/10) x - (3/20) y), x] /. {x -> 1, 
          y -> 0}

fy1 = 
      D[x ((3/10) - (3/10) x - (3/20) y), y] /. {x -> 1, 
          y -> 0}

gx1 = 
      D[y ((3/10) - (3/10) y - (3/20) x), x] /. {x -> 1, 
          y -> 0}

gy1 = 
      D[y ((3/10) - (3/10) y - (3/20) x), y] /. {x -> 1, 
          y -> 0}

matrixA4=(({{-(3/10), -(3/20)}, {0, (3/20)}}))

(iii)

nsol14[tval_, a_, b_, 
        t1_] := ({x[t], y[t]} /. 
            First[NDSolve[{x'[t] == 
                    x[t] (0.3 - 0.3  x[t] - 0.15  y[t]), y'[t] == 
                    y[t] (0.3 - 0.3  y[t] - 0.15  x[t]), x[0] == a, 
                  y[0] == b}, {x[t], y[t]}, {t, 0, t1}]]) /. 
        t -> tval
nphase14[t1_] := 
      ParametricPlot[
        Evaluate[
          Flatten[Table[
              nsol14[t, a, b, t1], {a, 0, 2,  .25}, {b, 0, 2,  .25}], 
            1]], {t, 0, t1}, AspectRatio -> 1, 
        AxesLabel -> {"x", "y"}]

lines14 = nphase14[5];

arrows14 = 
        PlotVectorField[{x (0.3 - 0.3  x - 0.15  y), 
            y (0.3 - 0.3  y - 0.15  x)}, {x, 0, 2}, {y, 0, 2}, 
          Axes -> True, ScaleFunction -> (1 &), 
          AspectRatio -> 0.8];

Show[arrows14, lines14, PlotRange -> {0, 2}, 
        AxesLabel -> {"x", "y"}];

Question 15

(i)

The system takes the form
	(Subscript[x, 1])(t) = 0.42  (Subscript[x, 3])(t)
(Subscript[x, 2])(t) = 0.6  (Subscript[x, 1])(t)
(Subscript[x, 3])(t) = 
        0.75  (Subscript[x, 2])(t) + 0.95  (Subscript[x, 3])(t)

(iii)

mA = {{0, 0, 0.42}, {0.6, 0, 0}, {0, 0.75, 0.95}}

Eigenvalues[mA]

Note also,

R = (Sqrt[(((-0.0774172))^2) + (0.406292^2)])

So the system is stable.

(iv)

From eigenvalue 1.10483 we obtain the growth rate of the bison population of 10.5%.

(v)

Eigenvectors[mA]

sumeig = 0.348901 + 0.189477 + 0.917805

{0.348901, 0.189477, 0.917805}/sumeig

Hence, calves settle down to 24%, yearlings to 13% and adults to 63%.