Standard Predator/Prey Model

The standard predator/prey model assumes the prey would grow exponentially ( x' = ax) without the interaction effect with predators ( -bxy). The interaction is detrimental to the prey and so the interaction coefficient -b is negative. The predators would "decay" exponentially (y' = -ny) if they had no food source, namely the prey. The interaction nxy between predator and prey is advantageous to the predators and so the interaction coefficient n is positive. These assumptions give rise to the standard predator/prey model.

> restart;with(student):

> `Standard predator/prey model`;`x'` = a*x - b*x*y;`y'` = -n*y + m*x*y; `a, b, n, m all positive`;

It seems impossible to solve the model explicitly for x(t) and y(t), but a calculus trick allows an implicit solution.

> `dy/dx`=`y'`/`x'`;`dy/dx` = (a*x-b*x*y)/(-n*y+m*x*y);

By separation of variables we get

> ((a-b*y)/y)*`dy` = ((m*x-n)/x)*`dx`;

> Int((a-b*y)/y,y)=Int((m*x-n)/x,x);

> int((a-b*y)/y,y)=int((m*x-n)/x,x)+ln(c);

> simplify(exp(lhs(%))=exp(rhs(%)));

> simplify(lhs(%)/rhs(%)*c=c);

> restart:with(plots):with(student):

Thus when the standard predator-prey model is integrated, the resulting orbits are level curves (contours) of the function J of two variables x and y below:

> J:=(x,y)->y^a*x^n*exp(-b*y)*exp(-m*x);

For the case a = 5, b = 3, m=2, n=6 we have:

> a:=5;b:=3;m:=2;n:=6; J(x,y);

The level curves have traditionally been computed by a method called "The Volterra Mapping Technique". Using plot3d however we can easily draw whatever level curves are required.

> plot3d(J(x,y),x=0..10,y=0..10,orientation=[-153,49],axes=BOXED,projection=1,shading=ZHUE,style=PATCHCONTOUR,contours=[.0001635,.001,.005,.05,.1,.154],title=`3D plot of J(x,y)`);

Using the mouse or by setting the orientation at [-90,0] will turn this curve over so that we are looking down the J(x,y) (or z) axis and the plot looks like a contour plot. My guess is that the orbits in the magenta part of the graph are very unstable. Small random events when both predator and prey are near extinction would almost surely perturb the point (x,y) into another orbit.

> plot3d(J(x,y),x=0..10,y=0..10,orientation=[-90,0],axes=BOXED,projection=1,shading=ZHUE,style=PATCHCONTOUR,contours=[.0001635,.001,.005,.05,.1,.154],title=`Looking down the z axis at J(x,y)`);

Maple offers a contourplot command which allows for a more direct way to do this.

> `Initial conditions x(0) = 10 and y(0) = 2 yield a c value of:`;evalf(J(10,2),4);

Below is a plot of this contour.

> contourplot(J(x,y),x=0..20, y=0..10, contours=[.0001635],numpoints=10000,title=`Blue: t=0 x(0)=10.0 y(0)=2.0 Magenta: t=1.3 Coral: del t = 0.1`,color=green);

> pic1:=%:

This gives the shape of this orbit, but it has no time labels, because the calculus procedure used to solve the pred-prey model elliminates time at the first step. The following data points including time were computed using FORTRAN and numerical methods.

>

> L:=[[5.29,5.58],[1.47,5.56],[.61,3.67],[.42,2.22],[.41,1.32],[.50,.79],[.68,.49],[.99,.32],[1.51,.22],[2.35,.18],[3.68,.18],[5.70,.24]];

> plot(L,style=point,symbol=circle,color=coral);

> pic2:=%:

> plot([[10.,2.]],color=blue,style=point,symbol=circle);

> pic3:=%:

> plot([[8.44,.55]],style=point, symbol=circle, color=magenta);

> pic4:=%:

Pics 2, 3 and 4 are needed to make the final graph, but need not be displayed. They are displayed anyway for instructional purposes. The following command puts pics 1, 2, 3 and 4 together into a graph of the orbit with time labels. The orbit starts at the blue point (10,2) at time 0 and proceeds anticlockwise around to the magenta point each coral point representing 1/10th time unit. The magenta point is time 1.3. The position at time 1.4 is almost exactly back at the starting point making the orbital period 1.4 time units.

> display(pic1,pic2,pic3,pic4);