Predator-Prey Models

Part 1. Background

The following loads the necessary package.

<< Graphics`PlotField`

Answer the questions to Part 1 in the cells below:

Answer to question 1:

Answer to question 2:

Answer to question 3:

Part 2. The Lotka-Volterra Model

Answer questions 1-3 in the cells below.

Answer to question 1:

Answer to question 2:

Answer to question 3:

Enter sample values of a, b, c, p:

a  =  1;
b  =  0.03;
c  =  0.4;
p  =  0.01;

The next command draws the slope (direction) field.  Note that the first argument is a list {1, dy/dx}.  This calculation will take a while.

fieldPlot  =  
        PlotVectorField[{a*x - b*x*y, -c*y + p*x*y},  {x, 0, 140}, {y, 
            0, 80},  PlotLabel -> "Lotka-Volterra Model",  
          AxesLabel -> {"x", "y"},  
          ScaleFunction -> (0.25 &),  Axes -> True,  
          PlotPoints -> 15];

Answer to question 4:

The following commands plot some trajectories determined by a list of starting points, each given by x and y coordinates at time 0.

system  =  {x'[t] == 
            a*x[t] - b*x[t]*y[t],  y'[t] == -c*y[t] + 
              p*x[t]*y[t]};
initial  =  {{x[0] == 15, y[0] == 15}, {x[0] == 20, 
            y[0] == 20}, {x[0] == 25, 
            y[0] == 25}, {x[0] == 30, 
            y[0] == 30}};
solutions  =  
        Table[NDSolve[
            Join[system, initial[[k]]],  {x,  y},  {t, 0, 10}], {k, 
            Length[initial]}];
solutionsPlot  =  
        Table[ParametricPlot[
            Evaluate[{x[t],  y[t]} /. solutions[[i]]],  {t,  0,  10}, 
            PlotStyle -> RGBColor[0, 0, 1], 
            DisplayFunction -> Identity], {i, 
            Length[initial]}];

Now we plot the solutions with the direction field shown earlier.  This will also take a while.

Show[solutionsPlot,  fieldPlot,  
        DisplayFunction ->  $DisplayFunction, 
        AxesLabel -> {"x", "y"}];

Answer to question 5:

Enter formulas to calculate the coordinates of the equilibrium point.

xs  =  c/p
ys  =  a/b

Answer the rest of the questions in Part 2 in the cells below.

Answer to question 6:

Answer to question 7:

Answer to question 8:

Answer to question 9:

Part 3. Euler's Method for Systems

We repeat one of our earlier calculations of a trajectory in the xy-plane, but using Euler's Method instead of whatever Mathematica did before.  Note that because Mathematica does not have a built-in command for Euler's Method, so we need to write one ourselves.

slope[x_, y_]  =  {a*x - b*x*y, -c*y + p*x*y};
n  =  40;
Delta  =  10/n;
Clear[t];
t[k_] :=  k*Delta;
Clear[x, y];
x[0]  =  20;
x[k_] :=  x[k]  =  
          x[k - 1]  +  
            N[ slope[x[k - 1], y[k - 1]][[1]] ]*
              Delta;
y[0]  =  20;
y[k_] :=  y[k]  =  
          y[k - 1]  +  
            N[ slope[x[k - 1], y[k - 1]][[2]] ]*
              Delta;
Table[{x[k], y[k]}, {k, 0, n}];
ListPlot[%, PlotStyle -> RGBColor[1, 0, 0], 
        Prolog -> {PointSize[0.015]}, PlotJoined -> True, 
        DisplayFunction -> Identity];
Show[%, fieldPlot, DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"x", "y"}, 
        PlotLabel -> "Euler's Method"];

Answer Part 3 questions below.

Answer to question 1:

Answer to question 2:

Part 4. Populations as Functions of Time

Suppose we would like to plot the prey and predator populations.  Let's choose, for example, the first scenario with x[0] = 15 and y[0] = 15, and look at the population over two periods.  First, we solve the system over two periods.

example  =  
        NDSolve[{system[[1]],  system[[2]],  x[0] == 15, 
            y[0] == 15},  {x,  y},  {t, 0, 20}];
curve  =  Flatten[example];

Now, we plot the prey and predator populations together.  Note that we have extracted the individual populations from the parametric solution returned by Mathematica's NDSolve.  The prey population is in red and that of the predator in blue.  Since we will be plotting graphs with similar names, we suppress the spelling warnings.

Off[General::spell1];

preyplot  =  
        Plot[Evaluate[x[t] /. curve[[1]],  {t,  0,  20}],  
          PlotStyle -> RGBColor[1, 0, 0],  
          DisplayFunction -> Identity];
predplot  =  
        Plot[Evaluate[y[t] /. curve[[2]],  {t,  0,  20}],  
          PlotStyle -> RGBColor[0, 0, 1],  
          DisplayFunction -> Identity];
Show[preyplot,  predplot,  
        DisplayFunction -> $DisplayFunction,  
        AxesLabel -> {"t",  "x"},  
        PlotLabel -> "Population for 2 Periods"];

Answer Part 4 questions in the cells below:

Answer to question 1:

Answer to question 2:

Part 5. Summary