# 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")
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).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)