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

# ======================================================================
# Build the simulation in R
# ======================================================================

# 1. Define time, the parameter, and the initial condition

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

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

# 2. Create the simulated times

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

# 3. Create and initialize the stock

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

# 4. Advance through time

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 inspect the results

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

# 6. Graph the trajectory

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

# ======================================================================
# Compare with the analytical solution
# ======================================================================

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

plot(
  results$time,
  results$exact,
  type = "l",
  lwd = 3,
  col = "#214f73",
  xlab = "Time (hours)",
  ylab = "Population (individuals)",
  main = "Euler simulation and analytical solution"
)
lines(results$time, results$population, lwd = 3, lty = 2, col = "#b65335")
legend(
  "topleft",
  legend = c("Analytical solution", "Euler simulation"),
  col = c("#214f73", "#b65335"),
  lty = c(1, 2),
  lwd = 3,
  bty = "n"
)

# ======================================================================
# Refine the time step
# ======================================================================

small_delta_t <- 0.25
small_time <- seq(from = start_time, to = end_time, by = small_delta_t)

small_population <- numeric(length(small_time))
small_population[1] <- initial_population

for (i in 2:length(small_time)) {
  population_old <- small_population[i - 1]
  growth <- growth_rate * population_old
  change <- growth * small_delta_t

  small_population[i] <- population_old + change
}

small_results <- data.frame(
  time = small_time,
  population = small_population
)

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

exact_at_end <- initial_population * exp(growth_rate * end_time)

data.frame(
  delta_t = c(delta_t, small_delta_t),
  number_of_steps = c(length(time) - 1, length(small_time) - 1),
  final_population = c(
    tail(results$population, 1),
    tail(small_results$population, 1)
  ),
  absolute_error = c(
    abs(tail(results$population, 1) - exact_at_end),
    abs(tail(small_results$population, 1) - exact_at_end)
  )
)

# ======================================================================
# Unconstrained decay
# ======================================================================

decay_start <- 100
decay_rate <- -0.10
decay_end_time <- 20
decay_delta_t <- 0.5

decay_time <- seq(from = 0, to = decay_end_time, by = decay_delta_t)
decay_amount <- numeric(length(decay_time))
decay_amount[1] <- decay_start

for (i in 2:length(decay_time)) {
  amount_old <- decay_amount[i - 1]
  decay_flow <- decay_rate * amount_old
  change <- decay_flow * decay_delta_t

  decay_amount[i] <- amount_old + change
}

decay_results <- data.frame(
  time = decay_time,
  amount = decay_amount
)

decay_results$exact <- decay_start *
  exp(decay_rate * decay_results$time)

plot(
  decay_results$time,
  decay_results$exact,
  type = "l",
  lwd = 3,
  col = "#214f73",
  xlab = "Time (hours)",
  ylab = "Amount",
  main = "Unconstrained decay"
)
lines(
  decay_results$time,
  decay_results$amount,
  lwd = 3,
  lty = 2,
  col = "#b65335"
)
legend(
  "topright",
  legend = c("Analytical solution", "Euler simulation"),
  col = c("#214f73", "#b65335"),
  lty = c(1, 2),
  lwd = 3,
  bty = "n"
)

# ======================================================================
# Use the model to answer questions
# ======================================================================

# When will the population double?

doubling_positions <- which(small_results$population >= 200)
first_double <- doubling_positions[1]

small_results[c(first_double - 1, first_double), ]

exact_doubling_time <- log(2) / 0.10
exact_doubling_time

# When will a decaying amount fall below 25?

below_25_positions <- which(decay_results$amount <= 25)
first_below_25 <- below_25_positions[1]

decay_results[c(first_below_25 - 1, first_below_25), ]

exact_time_below_25 <- log(25 / 100) / -0.10
exact_time_below_25