Functions

## generate a sequence of N(from N0 to Nt)
d.pop.growth <- function(N0, r, K, t){
    N <- replicate(t+1, N0)
    for (i in 1:t) {
        N[i+1] <- N[i] + r*N[i]*(1-N[i]/K)}
    N
}

## color generator
col.generate <- function(n) {
    rainbow(n, s = 1, v = 1, start = 0, end = max(1, n - 1)/n, alpha = 1)
}

## plot the first x-y diagram
plot_N <- function(N, col, main=NULL){
    plot(N, xlab = "Generation", main=main,ylab = "Pop. size", col=col, type="l", ylim =c(100,800))
}

## plot subsequent diagram
plot_more <- function(N, col){
    lines(N, col=col)
}

# Q3: time step 1 return map
N.by.N <- function(N){
    N <- rep(N, each = 2)
    N.plus <- lead(N)      # lead() is a function in dplyr
    mat <- cbind(N, N.plus)
}

plot_N2 <- function(mat, col, main=NULL){
    plot(mat, main=main, xlab = "Pop. size at t", ylab = "Pop. size at t+1", xlim=c(50, 700),ylim = c(100,700), pch=20,col=col)
    lines(mat, col=col) 
}

plot_more2 <- function(mat, col){
    points(mat, col=col,pch=20)
    lines(mat, col=col)
}

Q1

K=500, N0=100, r=1.8, 2.3, 2.45, 2.56, 2.8

col <- col.generate(5)
par(mfrow=c(2,2))
plot_N(d.pop.growth(100, 1.8, 500, 50), col=col[1])
plot_more(d.pop.growth(100, 2.3, 500, 50),col=col[2])
legend(40,830,legend=c("r = 1.8", "r = 2.3"),col=col[1:2], cex=.45, lty=1)
plot_N(d.pop.growth(100, 2.45, 500, 50),col=col[3])
legend(40,830,legend="r = 2.45",col=col[3], cex=.45, lty=1)
plot_N(d.pop.growth(100, 2.56, 500, 50),col=col[4])
legend(40,830,legend="r = 2.56",col=col[4], cex=.45, lty=1)
plot_N(d.pop.growth(100, 2.8, 500, 50),col=col[5])
legend(40,830,legend="r = 2.8",col=col[5], cex=.45, lty=1)

Q2

K=500, N0=101, r=1.8, 2.3, 2.45, 2.56, 2.8

Blue: N=101

Red: N=100

col <- c("blue", "red")
plot_N(d.pop.growth(101, 2.8, 500, 50),col=col[1], main = "r = 2.8")
plot_more(d.pop.growth(100, 2.56, 500, 50),col=col[2])

par(mfrow=c(2,2))
plot_N(d.pop.growth(101, 1.8, 500, 50), col=col[1],main = "r = 1.8")
plot_more(d.pop.growth(100, 1.8, 500, 50),col=col[2])

plot_N(d.pop.growth(101, 2.3, 500, 50),col=col[1], main = "r = 2.3")
plot_more(d.pop.growth(100, 2.3, 500, 50),col=col[2])
      
plot_N(d.pop.growth(101, 2.45, 500, 50),col=col[1], main = "r = 2.45")
plot_more(d.pop.growth(100, 2.45, 500, 50),col=col[2])

plot_N(d.pop.growth(101, 2.56, 500, 50),col=col[1], main = "r = 2.56")
plot_more(d.pop.growth(100, 2.56, 500, 50),col=col[2])

Q3: Step 1 Return Map

K=500, N0=100, r=1.8, 2.3, 2.45, 2.56, 2.8

col <- replicate(5, col.generate(5)[5])
plot_N2(N.by.N(d.pop.growth(100, 2.8, 500, 200)),col=col[5], main="r=2.8")
eq <- function(N){N + 2.8*N*(1-N/500)}
par(new=TRUE)
plot(eq, xlim=c(50, 700),ylim = c(100,700), col="green", xlab="", ylab="")
par(new=TRUE)
plot(function(x){x}, xlim=c(50, 700),ylim = c(100,700), xlab="", ylab="")

plot_N2(N.by.N(d.pop.growth(100, 2.56, 500, 200)),col=col[4], main="r=2.56")
eq <- function(N){N + 2.56*N*(1-N/500)}
par(new=TRUE)
plot(eq, xlim=c(50, 700),ylim = c(100,700), col="green", xlab="", ylab="")
par(new=TRUE)
plot(function(x){x}, xlim=c(50, 700),ylim = c(100,700), xlab="", ylab="")

plot_N2(N.by.N(d.pop.growth(100, 2.45, 500, 200)),col=col[3], main="r=2.45")
eq <- function(N){N + 2.45*N*(1-N/500)}
par(new=TRUE)
plot(eq, xlim=c(50, 700),ylim = c(100,700), col="green", xlab="", ylab="")
par(new=TRUE)
plot(function(x){x}, xlim=c(50, 700),ylim = c(100,700), xlab="", ylab="")

par(mfrow=c(1,2))
plot_N2(N.by.N(d.pop.growth(100, 1.8, 500, 50)), col=col[1], main="r=1.8")
eq <- function(N){N + 1.8*N*(1-N/500)}
par(new=TRUE)
plot(eq, xlim=c(50, 700),ylim = c(100,700), col="green", xlab="", ylab="")
par(new=TRUE)
plot(function(x){x}, xlim=c(50, 700),ylim = c(100,700), xlab="", ylab="")


plot_N2(N.by.N(d.pop.growth(100, 2.3, 500, 200)),col=col[2], main="r=2.3")
eq <- function(N){N + 2.3*N*(1-N/500)}
par(new=TRUE)
plot(eq, xlim=c(50, 700),ylim = c(100,700), col="green", xlab="", ylab="")
par(new=TRUE)
plot(function(x){x}, xlim=c(50, 700),ylim = c(100,700), xlab="", ylab="")

Q4 Bifurcation diagram

rmax <- 3
plot(-1, -1, xlim = c(1.8, rmax), ylim = c(0, 700), xlab = "r", ylab = "N")
r <- seq(1.9, rmax, by = 0.01)

n <- 200
for (z in 1:length(r)) {
    N <- vector()
    N[1] <- 100
    for (i in 1:n) {
        N[i+1] <- N[i] + r[z]*N[i]*(1-N[i]/500)
    }
    points(rep(r[z], length(N[100:n])), N[100:n], cex = 0.1, pch = 19)
}