Predator-Prey Models

Part 1. Background

The following line loads the necessary package.

<< Graphics`PlotField`

Answers to the questions in Part 1:

Part 2. The Lotka-Volterra Model

Answers to questions 1-3:

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

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

The next group of commands will draw the direction field.  We name the Lotka-Volterra equations LV1 and LV2.

LV1  =  x'[t] == a*x[t] - b*x[t]*y[t];
LV2  =  y'[t] == -c*y[t] + p*x[t]*y[t];
system  =  {LV1, LV2};

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

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.

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, 12}], {k, 
            Length[initial]}];
solutionsPlot  =  
        Table[ParametricPlot[
            Evaluate[{x[t],  y[t]} /. solutions[[i]]],  {t,  0,  12}, 
            PlotStyle -> RGBColor[0, 0, 1], 
            DisplayFunction -> Identity], {i, 
            Length[initial]}];

Show[fieldPlot, solutionsPlot, 
        DisplayFunction -> $DisplayFunction];

Answer to question 5:

Enter formulas to calculate the coordinates of the equilibrium point.

xs  =  ?;
ys  =  ?;

Answers to the rest of the questions in Part 2:

Part 3. Graphical Representations

The following commands create a plot of the prey population.

example1  =  
        NDSolve[{system[[1]],  system[[2]],  x[0] == 15, 
            y[0] == 15},  {x,  y},  {t, 0, 20}];
curve  =  Flatten[example1];
preyplot  =  
        Plot[Evaluate[x[t] /. curve[[1]],  {t,  0,  20}],  
          PlotStyle -> RGBColor[1, 0, 0], 
          AxesLabel -> {"x", "y"}];

Now we plot the predator population the same way and overlay the two graphs.  Since we will be plotting graphs with similar names, we suppress the spelling warnings.

Off[General::spell1];

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,y"},  
        PlotLabel -> "Population for 2 Periods"];

Answers to questions 3 and 4:

The following commands will construct a 3-dimensional plot of the curve {x[t], y[t], t}.  The options command BoxRatio gives the relative lengths of the sides of the frame.

example2  =  
        NDSolve[{system[[1]],  system[[2]],  x[0] == 15, 
            y[0] == 15},  {x,  y},  {t, 0, 30}];
curve  =  Append[Flatten[example2], t -> t];
spaceplot  =  
        ParametricPlot3D[Evaluate[{x[t], y[t], t} /. curve], {t, 0, 20}, 
          BoxRatios -> {1, 1, 1}, 
          AxesLabel -> {"x", "y", "t"}];

We can change the viewpoint by specifying it as an option using the command Viewpoint->{xscale, yscale, zscale}.  Note that these values do not refer to a Cartesian point.  Think of the center of the box as {0,0,0} with the longest side scaled to one unit.  The ordered triple for ViewPoint is the location of the viewpoint in this scaled coordinate system.  Try experimenting with different viewpoints by changing the value of view.
	There is quick way to find the viewpoint you want.  Go to the toolbar, and click on Input->3D ViewPoint Selector.  Select Cartesian Form, and experiment with the viewpoint until you are satisfied.  Then, click Paste, and the expression Viewpoint->{xvalue, yvalue, zvalue} will automatically be pasted into your open notebook at the location of your cursor.

view  =  {1, 0, 0};   
ParametricPlot3D[Evaluate[{x[t], y[t], t} /. curve], {t, 0, 20}, 
        BoxRatios -> {1, 1, 1}, ViewPoint -> view, 
        AxesLabel -> {"x", "y", "t"}];

To see their projections onto the coordinate planes, we can set the viewpoint very far away, e.g. Viewpoint->{100, 0, 0} will give us the projection onto the y-t plane.

ParametricPlot3D[Evaluate[{x[t], y[t], t} /. curve], {t, 0, 20}, 
        BoxRatios -> {1, 1, 1}, ViewPoint -> {100, 0, 0}, 
        AxesLabel -> {"x", "y", "t"}];

Alternatively, we can use ParametricPlot.

ParametricPlot[
        Evaluate[{x[t], y[t]} /. {curve[[1]], curve[[2]]}], {t, 0, 
          20}, PlotLabel -> "Lotka-Volterra Model",  
        AxesLabel -> {"x", "y"}];
ParametricPlot[
        Evaluate[{x[t], t} /. {curve[[1]], curve[[3]]}], {t, 0, 20}, 
        PlotLabel -> "Lotka-Volterra Model",  
        AxesLabel -> {"x", "t"}];
ParametricPlot[
        Evaluate[{y[t], t} /. {curve[[2]], curve[[3]]}], {t, 0, 20}, 
        PlotLabel -> "Lotka-Volterra Model",  
        AxesLabel -> {"y", "t"}];

Part 4. Varying the Parameters

Use the following commands to answer steps 1-4.  Make sure that you re-compile all Mathematica input cells so that you will have an updated plot to examine.

a  =  ?;
LV1  =  x'[t] == a*x[t] - b*x[t]*y[t];
LV2  =  y'[t] == -c*y[t] + p*x[t]*y[t];
system  =  {LV1, LV2};

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

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, 60}], {k, 
            Length[initial]}];
solutionsPlot  =  
        Table[ParametricPlot[
            Evaluate[{x[t],  y[t]} /. solutions[[i]]],  {t, 0, 60}, 
            PlotStyle -> RGBColor[0, 0, 1], 
            DisplayFunction -> Identity], {i, 
            Length[initial]}];

Show[fieldPlot, solutionsPlot, 
        DisplayFunction -> $DisplayFunction];

Part 5. The Effect of Hunting Predators

The following commands reset the parameters and define the modified differential equations with constant hunting of predators.  Note that the names of the equations have been changed.  Use your choice of graphical representation to determine the effect of hunting.  Be aware that if Mathematica encounters a singularity in NDSolve, the resulting interpolating function will only be valid over the interval [initial, singularity] as indicated in the output (below suppressed) for NDSolve.

a  =  1;
b  =   .03;
c  =   .4;
p  =   .01;
DE1  =  x'[t] == a*x[t] - b*x[t]*y[t];
DE2  =  y'[t] == -c*y[t] + p*x[t]*y[t] - 5;
system  =  {DE1, DE2};

fieldPlot  =  
        PlotVectorField[{a*x - b*x*y, -c*y + p*x*y},  {x, 0, 250}, {y, 
            0, 125},  PlotLabel -> "Hunting Predators",  
          AxesLabel -> {"x", "y"},  
          ScaleFunction -> (0.25 &),  Axes -> True,  
          PlotPoints -> 15, 
          DisplayFunction -> Identity];

solution1  =  
        NDSolve[{system[[1]],  system[[2]],  x[0] == 25, 
            y[0] == 50},  {x,  y},  {t, 0, 60}];
solution2  =  
        NDSolve[{system[[1]],  system[[2]],  x[0] == 60, 
            y[0] == 30},  {x,  y},  {t, 0, 60}];

Here are the commands to plot the solutions over the interval from t = 0 to 15.  Try experimenting with other values for t1max and t2max.  What is the range of valid values for t1max and t2max?

t1max  =  15;
t2max  =  15;
plot1  =  
        ParametricPlot[
          Evaluate[{x[t],  y[t]} /. solution1],  {t, 0, t1max}, 
          PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
plot2  =  
        ParametricPlot[
          Evaluate[{x[t],  y[t]} /. solution2],  {t , 0, t2max}, 
          PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
Show[fieldPlot, plot1, plot2, 
        DisplayFunction -> $DisplayFunction, 
        PlotRange -> {{0, 250}, {0, 125}}];

Part 6. Summary

Answers to the questions in Part 6: