Exercise 1: Plotting Verhulst dynamics

f <- function(t, p) {
  with(as.list(p), {
    x <- K*x0/(x0+(K-x0)*exp(-r*t))
    out <- data.frame(t = t, x = x, x0 = x0, r = r, K = K)  
    return(out)
  })
}
library(ggplot2)
theme_set(theme_bw())
times <- seq(0, 1, len = 100)
x0Sequence <- seq(0, 1, by = 0.1)

xlist <- lapply(x0Sequence, function(x0){
  p <- c(x0 = x0, K = 0.5, r = 4)
  out <- f(times, p)
  return(out)
})

out <- do.call(rbind, xlist)

P1 <- ggplot(out, aes(x = t, y = x, group = x0, color = x0)) + geom_line()
print(P1)

In addition, vary \(K\) within [0.2,1] and plot the curves in panels separated by their \(K\) values.

times <- seq(0, 1, len = 100)
x0Sequence <- seq(0, 1, by = 0.1)
KSequence <- seq(0.2, 1, len = 5)

x0KComb <- expand.grid(x0 = x0Sequence, K = KSequence)
NComb <- nrow(x0KComb)

xlist <- lapply(1:NComb, function(n){
  p <- c(x0 = x0KComb$x0[n], K = x0KComb$K[n], r = 4)
  out <- f(times, p)
  return(out)
})
out <- do.call(rbind, xlist)

P2 <- ggplot(out, aes(x = t, y = x, group = x0, color = x0)) + facet_wrap(~K) + geom_line()
print(P2)

times <- seq(0, 1, len = 100)
x0Sequence <- seq(0, 1, by = 0.1)
KSequence <- seq(0.2, 1, len = 5)
rSequence <- seq(1, 10, by = 2)
parComb <- expand.grid(x0 = x0Sequence, K = KSequence, r = rSequence)
NComb <- nrow(parComb)

xlist <- lapply(1:NComb, function(n){
  p <- c(x0 = parComb$x0[n], K = parComb$K[n], r = parComb$r[n])
  out <- f(times, p)
  return(out)
})
out <- do.call(rbind, xlist)

P2 <- ggplot(out, aes(x = t, y = x, group = x0, color = x0)) + facet_grid(r~K) + geom_line()
print(P2)

VH <- function(x, p) {
  with(as.list(p), {
    dx <- r*x*(1-x/K)  
    out <- data.frame(x = x, dx = dx)
    return(out)
  })
}
x <- seq(0, 1, len = 100)
p <- c(K = 0.75, r = 4, x0 = 1)
out <- VH(x, p)
P2 <- ggplot(out, aes(x = x, y = dx)) + geom_line()
print(P2)

x2 <- seq(0.01, 1, len = 10) 
out2 <- VH(x2, p)
P2 <- P2 + geom_segment(data = out2, aes(x = x, y = 0, xend = x + 0.1*dx, yend = 0),  arrow = arrow(angle=20, length = unit(0.12, "inches")), color="red")
print(P2)

Exercise 2: Data simulation and error bars

sigma <- 0.05
timesData <- seq(0, 1, len = 10)
p <- c(x0 = 0.1, K = 0.5, r = 4)
outData <- f(timesData, p)
outData$sigma <- sigma 
noise <- rnorm(length(timesData), 0, sigma)
outData$x <- outData$x + noise
times <- seq(0, 1, len = 100)
outSmooth <- f(times, p)
outSmooth$sigma <- sigma
P3 <- ggplot(outSmooth, aes(x = t, y = x, ymin = x - sigma, ymax = x + sigma)) + 
  geom_line() + geom_ribbon(alpha = 0.2) 
P3 <- P3 + geom_point(data = outData) + geom_errorbar(data = outData, width = 0)
print(P3)