Predator-Prey (Foxes and Rabbits)

We'll examine here the system , .

First note that our friend dsolve doesn't seem to help us much.

>	dsolve({D(R)(t)=3*R(t)-1.4*R(t)*F(t),D(F)(t)=-F(t)+0.5*R(t)*F(t)}, [R(t),F(t)]);

We can solve numerically by requesting the appropriate option in dsolve.

>	g:=dsolve({D(R)(t)=3*R(t)-1.4*R(t)*F(t),D(F)(t)=-F(t)+0.5*R(t)*F(t),               R(0)=2,F(0)=2},          [R(t),F(t)],          numeric);

We get a procedure as an output.  The procedure inputs a value of and outputs three equations, the equations for , and .

>

g(0);

>

g(1);

Using the command odeplot (in the package plots) we can graph the solutions.

>	with(plots):

>	plot1:=odeplot(g,[t,R(t)],0..10): plot2:=odeplot(g,[t,F(t)],0..10,color=blue):

Here are the two functions plotted on the same pair of axes.

>	display([plot1,plot2]);

Here we have a 3D plot.

>	odeplot(g,[t,R(t),F(t)],0..10,axes=boxed,color=red);

Here is an attempt to get a phase plane picture.  Note that the inaccuracies in the numerical picture give us a fuzzy picture.  

>	odeplot(g,[R(t),F(t)],0..10);

If we increase numpoints from the default 50 to 500, the picture smooths out.  It certainly looks like the solution here is periodic.

>	odeplot(g,[R(t),F(t)],0..10,numpoints=500);

The period seems to be a little over 3.6.

>	odeplot(g,[R(t),F(t)],0..3.6,numpoints=500);

Here's another cute command that can get you a nice picture of the phase plane.  Here there are three initial conditions specified.

>	phaseportrait( [D(R)(t)=3*R(t)-1.4*R(t)*F(t),D(F)(t)=-F(t)+0.5*R(t)*F(t)], [R(t),F(t)], t=0..6, [[R(0)=2,F(0)=2],[R(0)=2,F(0)=1],[R(0)=1,F(0)=2]], stepsize=.05 );

Modified Predator Prey - Logistic Growth for Rabbits

Here we make the rabbits growth pattern one of logistic growth with a carrying capacity of 4.

>	phaseportrait( [D(R)(t)=3*R(t)*(1-R(t)/4)-1.4*R(t)*F(t),D(F)(t)=-F(t)+0.5*R(t)*F(t)], [R(t),F(t)], t=0..6, [[R(0)=2,F(0)=2],[R(0)=2,F(0)=1],[R(0)=1,F(0)=2]], stepsize=.05 );

Note that the equilibrium point has moved.  It is now at

>	solve({3*R*(1-R/4)-1.4*R*F=0, -F+0.5*R*F=0},{R,F});

>	g2:=dsolve({D(R)(t)=3*R(t)*(1-R(t)/4)-1.4*R(t)*F(t), D(F)(t)=-F(t)+0.5*R(t)*F(t), R(0)=2,F(0)=1}, [R(t),F(t)], numeric);

>

g2(0);

>

g2(1);

>	plot21:=odeplot(g2,[t,R(t)],0..10): plot22:=odeplot(g2,[t,F(t)],0..10,color=blue):

>	display([plot21,plot22]);