Application Center - Maplesoft

App Preview:

Biological populations and ecological models

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

chaos.mws

Biological Populations and Ecological Models

by Prof. Matt Miller, Department of Mathematics, University of South Carolina
email: miller@math.sc.edu

This Maple worksheet is designed to illustrate two papers.

May, RM 1974. Biological populations with nonoverlapping generations: stable points, stable cycles and chaos. Science 186: 645-647.

May, RM & GF Oster. 1976. Bifurcations and dynamic complexity in simple ecological models. American Naturalist 110: 573-599.

May examined the behavior of a difference equation version of the logistic equation, and found that it exhibited very complex dynamics, depending upon the value of the growth rate `r'. He found monotonic damping to an equilibrium, damped oscillations, stable oscillations, and chaos. The equation he used was

N[t+1] = N[t]*exp(r*(1-N[t]))

This equation can never have negative population sizes, and has a carrying capacity of 1, like the dimensionless logistic that we developed earlier in the course (look what happens when N(t) = 1).

> restart; with(plots):

Warning, the name changecoords has been redefined

> timeplot:=NULL:

> for r from 1 to 3.0 by 0.2 do

> n[0]:=0.01;

> for t from 0 to 30 do

> n[t+1]:=n[t]*exp(r*(1-n[t]));

> od;

> timeplot:=timeplot,plot([ [[x,n[x]]$x=0..30], [x,r,x= -3..-1]], color=[red, blue], thickness=2, xtickmarks=4, ytickmarks=4, title=cat(`N vs t for r=`, convert(r, string)));

> od:

The following graph is an animation of the dynamics of the population at different values of `r'. The short blue line on the left side of the graph represents `r', and the red trajectory is the population size plotted as a function of time. Note that the population shows logistic behavior at low values of r, shows damped oscillations, stable 2-cycles, stable 4-cycles, 8-cycles and then chaos. The behavior is predicted by the plots of N[t+8] , N[t+4] , and N[t+2] vs N[t] above. Compare this graph to the figures in May 1974. Then change the initial condition from n[0]:= 0.01 to n[0]=0.81; what changes do you obseve?

To set the animation in play, click on the graph, and then use the play and stop buttons on the toolbar. The ->| button does single frame advance, the << and >> decrease and increase the speed of playback. The -> and <- buttons change the direction of playback. The short purple line at the left side of the graph is the value of r (read on vertical axis scale). Continuous loop play is also possible.

> display([timeplot],insequence=true);

[Maple Plot]

>

The bifurcation plot shows the stable points as a function of r. When r is less than 2, there is a single stable point. When r reaches 2, the single stable point bifurcates into 2 stable points with a stable period-2 cycle. When r reaches 2.5 the two point cycle bifurcates into a stable period-4 cycle with 4 stable points. This is Figure 8a of May and Oster, who provide an explanation for the changes in behavior of the system as r changes.

> bifurcplot:=NULL:

> for r from 1 to 3.0 by 0.01 do;

> n[0]:=0.99;

> for t from 0 to 70 do;

> n[t+1]:=n[t]*exp(r*(1-n[t]));

> od;

> bifurcplot:=bifurcplot, plot( [ [r,n[x]]$x =40..70 ], color=black, style=point, symbol=POINT):

> od:

> display([bifurcplot]);

[Maple Plot]

>

We will consider equation 1 of May and Oster, which is the same equation used by May 1974.

N[t+1] = F( N[t] )

where

F(N) = N[t]*exp(r*(1-N[t]))

This is analogous to the logistic equation, but is better behaved, because it cannot predict negative population sizes. By deciding on the value of the growth rate `r', we can plot a graph of N[t+1] versus N[t] , and see the shape of the function F(N). The red line is the function F(N), and the black line is a plot of N[t+1] = N[t] . This is figure 2a of May and Oster, using r = 0.5.

> plot([n*exp(0.5*(1-n)),n],n=0..2,color=[red,black]);

[Maple Plot]

>

The intersections of F(N) with the line N(t+1)=N(t) are the equilibria or fixed points, N*. The stability of the fixed points is easily determined because in this situation, the eigenvalue of F at a fixed point N* is the slope of F at the intersection with the 45 degree line. If the eigenvalue is less than 1 (slope <1), the fixed point is stable (attractor), if the eigenvalue is greater than 1 (slope > 1), the fixed point is unstable (repulsor).

The analogous procedure can be used to determine the relation between N(t+2) and N(t), or between N(t+4) and N(t). In these cases, if N(t+4)=N(t), there is a stable period-4 oscillation, and if N(t+2) = N(t), there is a stable period-2 oscillation. The stability of these fixed points is determined in the same way as described above. The population at time t+4 is found by recursively applying the function F . (Maple uses %1 to abbreviate a repeated chunk of stuff.)

> F:= n -> n*exp(0.5*(1-n)); t:= 't';

F := proc (n) options operator, arrow; n*exp(.5-.5*...

t := 't'

> N[t+1]:= F(N[t]);

N[t+1] := N[t]*exp(.5-.5*N[t])

> N[t+2]:= F(N[t+1]);

N[t+2] := N[t]*exp(.5-.5*N[t])*exp(.5-.5*N[t]*exp(....

> N[t+3]:= F(N[t+2]);

N[t+3] := N[t]*exp(.5-.5*N[t])*exp(.5-.5*N[t]*exp(....

> N[t+4]:= F(N[t+3]);

N[t+4] := N[t]*exp(.5-.5*N[t])*exp(.5-.5*N[t]*exp(....
N[t+4] := N[t]*exp(.5-.5*N[t])*exp(.5-.5*N[t]*exp(....

Next we illustrate the May & Oster cobweb plot method to find bifurcations. Note that as the growth rate r increases, the slope of function F(N) changes as does function F(F(N)). There are conditions where a single fixed point becomes 3 fixed points: an unstable point between two stable points with oscillation period 2. When this occurs, there is a `bifurcation' in the population process, and the population switches from a single stable population to a stable oscillation between the two new fixed points. This process repeats as each of the stable points bifurcates into 2 stable points with oscillation period 4, and repeats again, forming now 8 stable points with oscillation period 8.

> p2_list:=NULL: p_list:=NULL: p4_list:=NULL: p8_list:=NULL:

> for r from 1.0 to 3.1 by 0.15 do

> F:= n -> n*exp(r*(1-n)); # modified logistic model

Look at time 2 in the future to find 2 point cycles

> n1:= F(n): n2:= F(n1): n3:= F(n2):

Look at time 4 in the future to find 4 point cycles

> n4:= F(n3): n5:= F(n4): n6:= F(n5): n7:= F(n6):

Look at time 8 in the future to find 8 point cycles.

> n8:= F(n7):

> p2_list:= p2_list, plot([[n,n2,n=0..2], [n,n,n=0..2], [n,r,n=-0.3..-0.15]], color=[red,black,plum], ytickmarks=4);

> p4_list:=p4_list, plot([[n,n4,n=0..2], [n,n,n=0..2], [n,n2,n=0..2],[n,r,n=-0.3..-0.15]], color=[blue,black,red,plum], ytickmarks=4);

> p8_list:=p8_list,plot([[n,n8,n=0..3],[n,n4,n=0..3],[n,n,n=0..3],[n,r,n=-0.3..-0.15]], color=[sienna,blue,black,plum], ytickmarks=4);

> od:

The following animations display the fixed points of the logistic process and their stabilities. The two-point cycle graph is figure 6 in May and Oster. This graph plots N[t+2] vs N[t] in red. The function is determined by N[t+2] = F( N[t+1] ) = F(F( N[t] )). Note that the single stable point bifurcates into a stable 2-point cycle when r > 2, as noted in the legend of Fig 6 in May and Oster. The parameter r is indicated by a short plum line to the left of the vertical axis.

> display([p2_list],insequence=true,title=`2 point cycles`);

[Maple Plot]

>

This graph plots N[t+4] vs N[t] in blue and replots N[t+2] vs N[t] in red. Notice how the 2-cycle stable points bifurcate into 4-cycle stable points when r = 2.6.

> display([p4_list],insequence=true,title=`4 point cycles`);

[Maple Plot]

>

This graph plots N[t+8] vs N[t] in sienna and replots N[t+4] vs N[t] in blue. Note the number of crossing points of the mapping function and the 45 degree line. Each crossing point represents a fixed point that is either stable or unstable, depending upon the slope of the mapping function.

> display([p8_list],insequence=true,title=`8 point cycles`);

[Maple Plot]

>

Finally we use the F(N) vs N and N vs N plots to produce a sort of phase plane equivalent for the plot of N vs t. On the vetical axis we have N(t+1) and on the horizontal N(t). The first point plotted is (N[0], F(N[0])). The next value N[1] = F(N[0]) can be found graphically by extending the horizontal line at level F(N[0]) across to the 45 degree line. Using this position, which is N[1], on the horizontal axis, we move up (or down) to the curve to find N[2] = F(N[1]). From here we again move across to the line, then up or down to the curve, and so on. Try the following simple example first.

> F:= n -> n*exp(r*(1-n));

F := proc (n) options operator, arrow; n*exp(r*(1-n...

> n[0]:= 0.1; r:= 1.8; iterations:= 12;

n[0] := .1

r := 1.8

iterations := 12

> plotlist:= NULL: x:= n[0]:

> for j from 1 to iterations do

> y:= F(x): plotlist:= plotlist, [x, y], [y, y]:

> x:= y;

> od:

> P:= plot([plotlist], color=black, thickness=2): C:= plot( {F(n), n}, n = 0 .. 2.0):

> display( {P, C} );

[Maple Plot]

>

Change n[0], and see what happens. Record the results of your experiment! Next change the value of r and see what happens for the same sequence of values of n[0] that you used previously.

>