# Module 2.2: Unconstrained Growth
#
# This example uses Euler's method to simulate a population whose growth rate
# is proportional to its current size. Time is measured in hours.

# 1. Time, parameter, and initial condition

start_time <- 0             # hours
end_time <- 8               # hours
delta_t <- 1                # hours per simulation step

growth_rate <- 0.10         # per hour
initial_population <- 100   # individuals

# 2. Simulated time points

time <- seq(from = start_time, to = end_time, by = delta_t)

# 3. Stock variable

population <- numeric(length(time))
population[1] <- initial_population

# 4. Euler simulation
# new amount = old amount + rate of change * time passed

for (i in 2:length(time)) {
  population_old <- population[i - 1]
  growth <- growth_rate * population_old
  change <- growth * delta_t

  population[i] <- population_old + change
}

# 5. Organize and verify the results

results <- data.frame(time, population)
results

# The first two populations should be 100 and 110.
results[1:2, ]

# 6. Compare Euler's method with the analytical solution
# P(t) = P_0 * exp(r * t)

results$exact <- initial_population * exp(growth_rate * results$time)
results$error <- results$population - results$exact
results

# 7. Graph both trajectories

plot(
  results$time,
  results$exact,
  type = "l",
  lwd = 3,
  col = "blue",
  xlab = "Time (hours)",
  ylab = "Population (individuals)",
  main = "Unconstrained population growth"
)

lines(
  results$time,
  results$population,
  lwd = 3,
  lty = 2,
  col = "red"
)

legend(
  "topleft",
  legend = c("Analytical solution", "Euler simulation"),
  col = c("blue", "red"),
  lty = c(1, 2),
  lwd = 3,
  bty = "n"
)

# Try changing delta_t to 0.5 or 0.25. A smaller time step should make the
# Euler simulation closer to the analytical solution.