#may.R
#R worksheet to accompany Robert M. May, 1974 , Biological populations with non-overlapping generations:
# stable points, stable cycles and chaos , Science, 186: 645-647.
#This worksheet has three parts.
# In the first two parts we look at the model itself, and especially at the long term behaviour
# of the trajectories and how this may depend on initial conditions, growth parameters, etc.
# In section 1 we consider the dynamics of a single species;
# in section 2 we consider two competing species. (Not done yet)
# The third part introduces the idea of a bifurcation diagram, which is a record of the repeated
# "stable points" that occur for different growth parameter values.
#Section 1. Stable cycles and chaos.
#May examines the behavior of difference equation analogues of the logistic equation.
# He looks at two equations:
# 1. N(t+1) = N(t) * exp[r * (1 - N(t)/K]
# 2. N(t+1) = N(t) * [ 1 + r * (1 - N(t)/K)]
#Both are non-linear equations, and both have a stable equilibrium at N = K, when 2 > r > 0.
# The second, although more clearly "logistic-based" has the unfortunate tendancy to permit overshoot
# in which the population can dip negative under extreme oscillation, even with reasonably small stepsize.
# The reader is invited to experiment with this model as well,
# but in this worksheet all computations in this model have been "commented out".
#We can demonstrate this stability by calculating the trajectory of populations from an initial size of 10 individuals,
# over 30 time steps. We will try several values of r from 0.5 to 2.0.
#We will use equation 1 as May did in his paper.
#To avoid having extraneous "junk" getting printed, we suppress all printing and plotting in the main loop
#(under control of r), give the plots names, and plot them afterwards.
# Two plots are created: plot_list gives N as a function of t; pairs gives N(t+1) against N(t).
# Clustering of points in the latter is an indication of stable cyclic behaviour.
N<-vector(length=31)
N[1]<- 10 ; K<-1000 ; # initial condition and `carrying capacity`
title1<-"N(t+1) vs N(t) when r="
title2<-"N(t) vs t when r= " # title blanks for plots
sample_r<-c (0.4, 0.7, 1.0, 1.3, 1.6, 1.9);
P_frames<- NULL; R_frames<- NULL;
par(mfrow=c(2,3))
for (r in sample_r)
{
for (t in 1:30)
{
#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] ); # logistic model
N[t+1]<- N[t] * exp(r * (1 - N[t] / K)) #ricker model
}
plot(N[1:(length(N)-1)],N[2:length(N)], col="green", main=paste(title1, r),
xlab="N[t]", ylab="N[t+1]")
}
for (r in sample_r)
{
for (t in 1:30)
{
#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] ); # logistic model
N[t+1]<- N[t] * exp(r * (1 - N[t] / K)) #ricker model
}
plot(1:length(N),N, col="blue",type="l",main=paste(title1, r),
xlab="t", ylab="N[t]")
}
#You should now have six graphs of the trajectory of the population reaching a stable equilibrium
# of K, at growth rates r of 0.4, 0.7, 1.0, 1.3, 1.6, and 1.9.
# Let us now try values of r close to 2.0 and see if the behavior indeed changes
# where May says it does.
N[1]<- 10 ; K<-1000 ; # initial condition and `carrying capacity`
title1<-"N(t+1) vs N(t) when r="
title2<-"N(t) vs t when r= " # title blanks for plots
sample_r<-c ( 1.99, 1.999, 2.001, 2.01);
P_frames<- NULL; R_frames<- NULL;
par(mfrow=c(2,2))
for (r in sample_r)
{
for (t in 1:2000)
{
#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] ); # logistic model
N[t+1]<- N[t] * exp(r * (1 - N[t] / K)) #ricker model
}
plot(N[1:(length(N)-1)],N[2:length(N)], col="green", main=paste(title1, r),
xlab="N[t]", ylab="N[t+1]")
}
for (r in sample_r)
{
for (t in 1:2000)
{
#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] ); # logistic model
N[t+1]<- N[t] * exp(r * (1 - N[t] / K)) #ricker model
}
plot(1:length(N),N, col="blue",type="l",main=paste(title2, r),
xlab="t", ylab="N[t]",xlim=c(1900,2000),ylim=c(0,1200))
}
#We have plotted from time step 1900 to time step 2000 on the graphs to see whether
# the populations have reached an equilibrium or whether they have stabilized.
# The model with r=1.999 appears to be a damped oscillation which is still damping,
# whereas the model with r=2.001 appears to have entered into a stable oscillation.
#Let us now try examining the regions where May says there is a 2 point cycle (2.526 > r > 2.0),
# a 4 point cycle (2.656 > r > 2.526), an 8 point cycle (2.685 > r > 2.656),
# a cycle of 16 or 32 points (2.692>r>2.685).
# We will also plot graphs of N(t+1) vs N(t) to see the underlying structure in the data:
N[1]<- 10 ; K<-1000 ; # initial condition and `carrying capacity`
title1<-"N(t+1) vs N(t) when r="
title2<-"N(t) vs t when r= " # title blanks for plots
sample_r<-c ( 2.4,2.6,2.67,2.69);
P_frames<- NULL; R_frames<- NULL;
par(mfrow=c(2,2))
for (r in sample_r)
{
for (t in 1:60)
{
#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] ); # logistic model
N[t+1]<- N[t] * exp(r * (1 - N[t] / K)) #ricker model
}
plot(N[1:(length(N)-1)],N[2:length(N)], col="green", main=paste(title1, r),
xlab="N[t]", ylab="N[t+1]")
}
for (r in sample_r)
{
for (t in 1:60)
{
#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] ); # logistic model
N[t+1]<- N[t] * exp(r * (1 - N[t] / K)) #ricker model
}
plot(1:length(N),N, col="blue",type="l",main=paste(title2, r),
xlab="t", ylab="N[t]",xlim=c(1,60),ylim=c(0,2000))
}
#Now let us examine May's assertion that in the region of chaotic behavior,
# the initial conditions determine the pattern of behavior.
# Let us try four simulations with r = 3.3, each with a different starting population.
# The starting conditions used in the second and fourth graphs correspond to May's conditions
# of No/K = 0.02 (No = 20) and No/K = 1.5 (No = 1500); see Fig 1(d) and Fig 1(e).
# We run the time up to 50 instead of 30 (why stop when you're having fun?).
r<- 3.3; K<- 1000; sample_inits<- c(10, 20, 500, 1500); # parameters and initial conditions
title1<-"N(t+1) vs N(t) when N[0]="
title2<-"N(t) vs t when N[0]= " # title blanks for plots
for (j in sample_inits)
{ plotnum<-plotnum+1
N[1]<-j
for (t in 1:50)
{
#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] ); # logistic model
N[t+1]<- N[t] * exp(r * (1 - N[t] / K)) #ricker model
}
plot(1:length(N),N, col="blue",type="l",main=paste(title1, j),
xlab="t", ylab="N[t]",xlim=c(1,60),ylim=c(0,3000))
}
for (j in sample_inits)
{
N[1]<-j
for (t in 1:50)
{
#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] ); # logistic model
N[t+1]<- N[t] * exp(r * (1 - N[t] / K)) #ricker model
}
plot(N[1:(length(N)-1)],N[2:length(N)], col="green", main=paste(title2, r),
xlab="N[t]", ylab="N[t+1]")
}
# Section 3. Bifurcation analysis.
#We recall that May examines the behavior of difference equation analogues of the
# logistic equation. He looks at two equations:
# 1. N(t+1) = N(t) * exp[r * (1 - N(t)/K]
# 2. N(t+1) = N(t) * [ 1 + r * (1 - N(t)/K)]
#Both are non-linear equations, and both have a stable equilibrium at N = K,
# when 2 > r > 0. In the long term the population N(t) exhibits cyclic behaviour,
# which becomes more pronounced as the value of the growth parameter r increases.
# The extremes of the oscillations, known as stable points, recur periodically.
# Here we examine the distribution of stable points as a function of the growth rate
# parameter r. For the purpose of this exercise we will assume the system becomes
# stable after 20 time steps, so the last 10 time steps will be used to estimate the
# locations of the 'stable points'. This kind of plot is called a bifurcation plot
#
number_of_steps<- 30; record_after_step<- 21;
rvals<-seq(1.0,3.0,0.02)
par(mfrow=c(1,1))
for (r in rvals)
{
for (t in 1:number_of_steps)
{
#N[t+1]<- N[t] + ( r / K ) * N[t] * ( K - N[t] ); # discrete logistic model
N[t+1]<- N[t] * exp(r * (1-N[t] / K)); # May's eq 1
#Save the "stable points" in a data list called plot_list. Only use data from time steps 21 - 30 in each run.
}
if(r==1.0){ plot(rep(r,10),N[record_after_step:number_of_steps],
xlab="r",ylab="N",xlim=c(1,3),ylim=c(0,3000),
main="Bifurcation Map of N vs r")}
if(r>1.0){ points(rep(r,10),N[record_after_step:number_of_steps])}
}