1. [More computation problem] In lab on Monday, professor Aaron discussed the inherent inaccuracies and imprecisions when you use computers to solve mathematical problems. For most ecological models we don’t have analytical solutions, so we must rely upon our sometimes-innaccurate and -imprecise computational tools. In professor Aaron’s lecture notes he showed pseudocode on how to build a solver on page 4 from here: link. What I woudl like to do is for you to use an equation for which we know the analytic solution, the logistic equation, and then write your own solver to see how accurate it is. The solution to logistic equation is \[N(t) = \frac{KN(0)e^{rt}}{K + N(0)\left(e^{rt} - 1\right)}.\] Write a solver like professor Aaron and use \(\Delta t = 1, 0.1, 0.01,0.001\). For the sake of standardization across groups, use the following parameter values: \(K = 10, r = 1, N(0) = 0.1\).
# Parms
K <- 10
r = 1
N0 <- 0.1

t_vec <- seq(from = 0, to = 10, length.out = 1000)
Nt <- K*N0*exp(r*t_vec) / (K + N0*(exp(r*t_vec) - 1))
#Nt <- N0*exp(r*t_vec)
# Time series
plot(x = t_vec, y = Nt, type = "l")

# Solver.simulation
simLength <- 10
initPop <- N0
initr <- r
initK <- K

deltat_a <- 1
deltat_b <- 0.1
deltat_c <- 0.01
deltat_d <- 0.001

numiter_a <- simLength/deltat_a
numiter_b <- simLength/deltat_b
numiter_c <- simLength/deltat_c
numiter_d <- simLength/deltat_d

pop_vec_a <- rep(x = NA, times = numiter_a)
pop_vec_b <- rep(x = NA, times = numiter_b)
pop_vec_c <- rep(x = NA, times = numiter_c)
pop_vec_d <- rep(x = NA, times = numiter_d)

pop_vec_a[1] <- initPop
pop_vec_b[1] <- initPop
pop_vec_c[1] <- initPop
pop_vec_d[1] <- initPop

# Simulation a
for (i in 2:numiter_a) {
  popt <- pop_vec_a[i - 1]
  pgr <- initr*popt*(initK - popt)/initK
  deltaPop <- pgr*deltat_a
  poptp1 <- popt + deltaPop
  pop_vec_a[i] <- poptp1
}

# Simulation b
for (i in 2:numiter_b) {
  popt <- pop_vec_b[i - 1]
  pgr <- initr*popt*(initK - popt)/initK
  deltaPop <- pgr*deltat_b
  poptp1 <- popt + deltaPop
  pop_vec_b[i] <- poptp1
}

# Simulation c
for (i in 2:numiter_c) {
  popt <- pop_vec_c[i - 1]
  pgr <- initr*popt*(initK - popt)/initK
  deltaPop <- pgr*deltat_c
  poptp1 <- popt + deltaPop
  pop_vec_c[i] <- poptp1
}

# Simulation d
for (i in 2:numiter_d) {
  popt <- pop_vec_d[i - 1]
  pgr <- initr*popt*(initK - popt)/initK
  deltaPop <- pgr*deltat_d
  poptp1 <- popt + deltaPop
  pop_vec_d[i] <- poptp1
}

par(mfrow = c(2,2), mar = c(4.5, 4.5, 1, 0.5))
plot(x = t_vec, y = Nt, type = "l", las = 1, xlab = "Time", ylab = "Density (N)", main = bquote(Delta~"t " == 1))
  lines(x = seq(from = deltat_a, to = simLength, by = deltat_a), y = pop_vec_a, col = "red")
plot(x = t_vec, y = Nt, type = "l", las = 1, xlab = "Time", ylab = "Density (N)", main = bquote(Delta~"t " == 0.1))
  lines(x = seq(from = deltat_b, to = simLength, by = deltat_b), y = pop_vec_b, col = "red")
plot(x = t_vec, y = Nt, type = "l", las = 1, xlab = "Time", ylab = "Density (N)", main = bquote(Delta~"t " == 0.01))
  lines(x = seq(from = deltat_c, to = simLength, by = deltat_c), y = pop_vec_c, col = "red")
  plot(x = t_vec, y = Nt, type = "l", las = 1, xlab = "Time", ylab = "Density (N)", main = bquote(Delta~"t " == 0.001))
  lines(x = seq(from = deltat_d, to = simLength, by = deltat_d), y = pop_vec_d, col = "red")

  1. [More ecological problem] 1. Please recreate Figure 9.7 from your text. Make the figures look tidy, and include a vector field and nullclines. Make up your own parameter values but also look for some guidance within the book if you don’t know where to start.
    If you want to create a multi-panel plot, then use the function par() and set the argument mfrow = c(n, m) with n being the number of rows and m being the number of columns (e.g., par(mfrow = c(2, 3)) creates a plot with 2 panels for rows and 3 columns).
    There is often a lot of white space as a default around plot(). You can control that with mar = c(a, b, c, d)) as an argument in par(). a—d respectfully adjust the whitespace at the bottom, left, top, and right of each panel.

I write out the nullcline equations here, but using phaseR::nullclines() would work just as well.

n_x <- 1000
N1_seq <- seq(from = 0, to = 10, length.out = n_x)

# Panel A
alpha12 <- alpha21 <- 0.2
K1 <- K2 <- 5
N2_eqn1_A <- -K1/alpha12 + 1/alpha12*N1_seq
N2_eqn2_A <- K2 + alpha21*N1_seq

# Panel B
alpha12 <- alpha21 <- 1
h12 <- h21 <- 0.2
N2_eqn1_B <- h12*(N1_seq - K1)/(alpha12 - N1_seq + K1)
N2_eqn2_B <- K2 + alpha21*N1_seq/(h21 + N1_seq)
N1_seq_pos_B <- which(N2_eqn1_B > 0)

# Panel C
alpha12 <- alpha21 <- 0.75
K1 <- K2 <- 2.5
N2_eqn1_C <- K1/alpha12 + 1/alpha12*N1_seq
N2_eqn2_C <- -K2 + alpha21*N1_seq

# Panel D
K1 <- -3
K2 <- -3
alpha12 <- 9
alpha21 <- 9
h12 <- 1
h21 <- 1
N2_eqn1_D <- h12*(N1_seq - K1)/(alpha12 - N1_seq + K1)
N2_eqn2_D <- K2 + alpha21*N1_seq/(h21 + N1_seq)
N1_seq_pos_D <- which(N2_eqn1_D > 0)

# Panel E
K1 <- -5
K2 <- -5
alpha12 <- 0.2
alpha21 <- 1.5
N2_eqn1_E <- K1/alpha12 + 1/alpha12*N1_seq
N2_eqn2_E <- K2 + alpha21*N1_seq

# Panel F
K1 <- -3
K2 <- 3
alpha12 <- 9
alpha21 <- 2
h12 <- 1
h21 <- 1
N2_eqn1_F <- h12*(N1_seq - K1)/(alpha12 - N1_seq + K1)
N2_eqn2_F <- K2 + alpha21*N1_seq/(h21 + N1_seq)
N1_seq_pos_F <- which(N2_eqn1_F > 0)

# Panel G
alpha12 <- alpha21 <- 2
K1 <- 5
K2 <- -5
N2_eqn1_G <- -K1/alpha12 + 1/alpha12*N1_seq
N2_eqn2_G <- K2 + alpha12*N1_seq

# Panel H
K1 <- -7.25
K2 <- 0.5
alpha12 <- 9
alpha21 <- 7
h12 <- 0.25
h21 <- 2
N2_eqn1_H <- h12*(N1_seq - K1)/(alpha12 - N1_seq + K1)
N2_eqn2_H <- K2 + alpha21*N1_seq/(h21 + N1_seq)
N1_seq_pos_H <- which(N2_eqn1_H > 0)



plot_lims <- c(0, 7)

phase_plot <- function(x_seq, N1, N2, lims) {
  plot(x = x_seq, y = N1, col = "red", type = "l", xlim = lims, ylim = lims, las = 1, xlab = bquote("N"[1]), ylab = bquote("N"[2]))
  lines(x = N1_seq, y = N2, col = "blue")
  segments(x0 = 0, x1 = 0, y0 = -1, y1 = 10, col = "red")
  segments(x0 = -1, x1 = 10, y0 = 0, y1 = 0, col = "blue")
}

  
par(mfrow = c(4, 2), mar = c(4.5, 4.5, 0.5, 0.5))
# Panel A
  phase_plot(x_seq = N1_seq, N1 = N2_eqn1_A, N2 = N2_eqn2_A, lims = plot_lims)
# Panel B
  phase_plot(x_seq = N1_seq[N1_seq_pos_B], N1 = N2_eqn1_B[N1_seq_pos_B], N2 = N2_eqn2_B, lims = plot_lims)
# Panel C
  phase_plot(x_seq = N1_seq, N1 = N2_eqn1_C, N2 = N2_eqn2_C, lims = plot_lims)
# Panel D
  phase_plot(x_seq = N1_seq[N1_seq_pos_D], N1 = N2_eqn1_D[N1_seq_pos_D], N2 = N2_eqn2_D, lims = plot_lims)
# Panel E
  phase_plot(x_seq = N1_seq, N1 = N2_eqn1_E, N2 = N2_eqn2_E, lims = plot_lims)
# Panel F
  phase_plot(x_seq = N1_seq[N1_seq_pos_F], N1 = N2_eqn1_F[N1_seq_pos_F], N2 = N2_eqn2_F, lims = plot_lims)
# Panel G
  phase_plot(x_seq = N1_seq, N1 = N2_eqn1_G, N2 = N2_eqn2_G, lims = plot_lims)
# Panel H
  phase_plot(x_seq = N1_seq[N1_seq_pos_H], N1 = N2_eqn1_H[N1_seq_pos_H], N2 = N2_eqn2_H, lims = plot_lims)