next up previous
Next: Fun with competition Up: No Title Previous: No Title

Two-species models

In this short tutorial, we will examine some of the features of XPPAUT for two-dimensional systems of interacting species. We will start with the Lotka-Volterra model:

\begin{displaymath}x' = ax-bxy \qquad y'=-cx+dxy

A possible ODE file lv.ode , has the following simple form:
# lotka-volterra eqns
par a=1,b=1,c=1,d=1
init x=1,y=.1
Run XPPAUT by typing xpp lv.ode. Click into the Main Window and tap \fbox{{\bf i}} \fbox{{\bf g}} to integrate the equations (On the menu, \fbox{{\tt Initconds}} \fbox{{\tt Go}}.) Tap \fbox{{\bf x}} and \fbox{{\bf Enter}} to plot the variable x against time. (On the menu it is \fbox{{\tt Xvst}}.) Note the periodic behavior. Add y(t) to the picture by choosing \fbox{{\tt Graphic stuff}} \fbox{{\tt Add curve}} ( \fbox{{\bf g}} \fbox{{\bf a}}) and fill the dialog box as
\fbox{\begin{tabular}{ll\vert ll\vert ll }
\par {\tt X-axis:}& t & {\tt Z-axis:}...
... type: & 0 \\
{\tt Y-axis:}& y & {\tt Color:}& 1 & {} &
You should see the prey x in white and the predator y in red. Note the predator and prey undergo large oscillations and prey leads the predator. In the Main Window click on \fbox{{\tt Erase}} ( \fbox{{\bf e}}) to erase the screen. Now, in the Initial Data Window change the initial condition for y to .8 and click on \fbox{{\tt Go}} in the Initial Data Window . Note there are still oscillations, but they have a much smaller amplitude. This is because the Lotka-Volterra model is integrable and has infinitely many periodic orbits. None of them are limit cycles (which are isolated periodic solutions.)

Now lets look at the phase-plane. In the Main Window , click on \fbox{{\tt Graphic stuff}} \fbox{{\tt Remove all}} ( \fbox{{\bf g}} \fbox{{\bf r}}). Change the view by clicking on \fbox{{\tt View}} \fbox{{\tt 2D}} ( \fbox{{\bf v}} \fbox{{\bf 2}}). Fill in the dialog box as follows:

\fbox{\begin{tabular}{l\vert l}
{\tt X-axis:} X & {\tt Xmax:} 3 \\
...label:} prey \\
{\tt Ymin:} -.1 & {\tt Ylabel:} predator
You should see a small circle representing the periodic orbit. Click on \fbox{{\tt Initconds}} \fbox{{\tt mIce}} ( \fbox{{\bf i}} \fbox{{\bf i}}) and use the mouse to put choose a bunch of initial conditions. Click on \fbox{{\bf Esc}} when you get bored. Notice every positive initial condition goes to a periodic orbit.

Click on \fbox{{\tt Dir.field/flow}} \fbox{{\tt Direct Field}} ( \fbox{{\bf d}} \fbox{{\bf d}}) and accept the default of 10. A 10x10 array of small lines will be drawn showing the direction field. Now Click on \fbox{{\tt Nullcline}} \fbox{{\tt New}} ( \fbox{{\bf n}} \fbox{{\bf n}}) to draw the nullclines. There are 2 fixed points, (0,0) and (1,1). Let's check the stability of the (1,1) fixed point. Click on \fbox{{\tt Sing pts}} \fbox{{\tt Mouse}} ( \fbox{{\bf s}} \fbox{{\bf m}}) and click the mouse near the (1,1) fixed point. Click on No (don't print the eigenvalues) and a new window will appear showing the fixed point and the structure of the eigenvalues. ( c+ is the number of complex eigenvalues with positive real parts, and so on, with im the number of pure imaginary eigenvalues.) Note the equilibrium point has neutral stability. Quit XPPAUT by tapping \fbox{{\bf f}} \fbox{{\bf q}} \fbox{{\bf y}}.

Now let's solve the next simplest predator-prey model in which the prey has a finite carrying capacity:

\begin{displaymath}x' = ax(1-x/k -by) \qquad y'=y(-c+d*x).

The XPPAUT file, pp.ode has the form:
# predator-prey with carrying capacity
x'=a*x*(1-x/k -b*y)
par a=1,b=1,c=1,d=.4,k=2
init x=1.9,y=.2
@ xlo=-.1,xhi=3,ylo=-.1,yhi=3,xp=x,yp=y
This file has an added line that starts with the @ symbol. This tells XPPAUT that some control commands follow. The ones here tell XPPAUT the dimensions of the plotting window and which variables to plot. There are many other control commands that eliminate the necessity of calling up dialog boxes.

Start the program by typing xpp pp.ode. Here are some things to try:

Draw the nullclines (Tap \fbox{{\bf n}} \fbox{{\bf n}}).
Draw the direction fields (Tap \fbox{{\bf d}} \fbox{{\bf d}}.)
Choose several initial data with the mouse (Tap \fbox{{\bf i}} \fbox{{\bf i}}) and type \fbox{{\bf Esc}} when you get tired of this.
What does all this tell you? Note that there are only two visible fixed points. The predators get very little from each prey since d is small (only 0.4). Recall that c/d<k for stable coexistence of the predators. In the Parameter Window change d to 1 and then click on Ok in the Parameter Window . Click on the Main Window and erase the screen. Repeat the above three things. Notice all positive solutions spiral into the coexistent fixed point with non-zero predators. There are no periodic solutions - ever. Note that all of Kolmogorov's assumptions hold except the condition on the slope of the nullclines. The condition for periodic orbits is that the coexistent equilibrium should occur in a region where the prey nullcline has a positive slope.

Here is another example for you to try. This is a predator prey model that has growth of the prey increasing at low populations and then decreasing at high rates. The model is similar to our other models in other respects:

\begin{displaymath}x'=x(1+ax-x^2-by) \qquad y'=y(-c+dx)

We will look at this model as the parameter d increases similarly to what we did above. Here is the XPPAUT file, pplc.ode:
# a simple model that obeys Kolmogorov's assumptions
par a=2,b=1,c=1,d=.35
@ xlo=-.1,ylo=-.1,xhi=3,yhi=3,xp=x,yp=y
@ total=40
As in the previous file, I have set up the file so that it comes up in the phaseplane. Draw the nullclines and the direction fields. Draw some representative trajectories. Here d is so small that the y-nullcline x=c/d does not intersect the x-nullcline in the positive quadrant. Thus, the predators die out and the prey go to their maximum value of around 2.5. Increase d to 0.6. Erase the screen and draw the nullclines and direction fields. Draw some representative trajectories. Note that now, the prey and predators co-exist. Finally, increase d to 1.2. Now erase the screen and look at the nullclines, etc? The nullclines indicate that this system now satisfies Kolmogorov's criterion for limit cycles.

At what value of d do limit cycles appear? The way to figure this out is to see at what value of d the coexistent state loses stability. This will happen when the y-nullcline intersects the x-nullcline to the left of the maximum and the x-nullcline has a positive slope. The x-nullcline is

y = (1+ax-x2)/b

which has a maximum at x=a/2. Thus, if c/d<a/2, then the coexistent state will be unstable.

next up previous
Next: Fun with competition Up: No Title Previous: No Title
G. Bard Ermentrout