# Numerical Methods in Practice - COSC/MATH 201
# Compare Euler, RK2, and RK4 on the same logistic-growth model.

start_time <- 0
end_time <- 8
delta_t <- 1                 # Try 0.5 and 2 after the first run.
time <- seq(start_time, end_time, by = delta_t)

initial_population <- 100
growth_fraction <- 0.8       # per year
capacity <- 1000
target <- 500

euler <- numeric(length(time))
rk2 <- numeric(length(time))
rk4 <- numeric(length(time))
euler[1] <- initial_population
rk2[1] <- initial_population
rk4[1] <- initial_population

for (i in 2:length(time)) {
  # Euler: use the current rate for the entire step.
  old_euler <- euler[i - 1]
  euler_rate <- growth_fraction * old_euler *
    (1 - old_euler / capacity)
  euler[i] <- old_euler + euler_rate * delta_t

  # RK2: predict a trial endpoint, then average the two rates.
  old_rk2 <- rk2[i - 1]
  rate_1 <- growth_fraction * old_rk2 *
    (1 - old_rk2 / capacity)
  predicted <- old_rk2 + rate_1 * delta_t
  rate_2 <- growth_fraction * predicted *
    (1 - predicted / capacity)
  rk2[i] <- old_rk2 + (rate_1 + rate_2) / 2 * delta_t

  # RK4: four full-step trial changes from the same old value.
  old_rk4 <- rk4[i - 1]
  change_1 <- growth_fraction * old_rk4 *
    (1 - old_rk4 / capacity) * delta_t

  trial_mid_1 <- old_rk4 + change_1 / 2
  change_2 <- growth_fraction * trial_mid_1 *
    (1 - trial_mid_1 / capacity) * delta_t

  trial_mid_2 <- old_rk4 + change_2 / 2
  change_3 <- growth_fraction * trial_mid_2 *
    (1 - trial_mid_2 / capacity) * delta_t

  trial_end <- old_rk4 + change_3
  change_4 <- growth_fraction * trial_end *
    (1 - trial_end / capacity) * delta_t

  rk4[i] <- old_rk4 +
    (change_1 + 2 * change_2 + 2 * change_3 + change_4) / 6
}

results <- data.frame(time, euler, rk2, rk4)
print(round(results[1:4, ], 3))

# Vector positions are not times: use time[position].
first_euler <- which(euler >= target)[1]
first_rk2 <- which(rk2 >= target)[1]
first_rk4 <- which(rk4 >= target)[1]

threshold_results <- data.frame(
  method = c("Euler", "RK2", "RK4"),
  first_recorded_year = c(time[first_euler],
                          time[first_rk2],
                          time[first_rk4]),
  population_at_that_year = round(c(euler[first_euler],
                                    rk2[first_rk2],
                                    rk4[first_rk4]), 2)
)
print(threshold_results)

plot(time, euler, type = "o", pch = 16, col = "firebrick",
     ylim = range(c(euler, rk2, rk4, capacity)),
     xlab = "Time (years)", ylab = "Population",
     main = "One growth model, three numerical methods")
lines(time, rk2, type = "o", pch = 17, col = "darkgreen")
lines(time, rk4, type = "o", pch = 15, col = "purple4")
abline(h = target, lty = 2, col = "gray40")
abline(h = capacity, lty = 3, col = "gray60")
legend("bottomright", c("Euler", "RK2", "RK4", "Target = 500", "Capacity = 1000"),
       col = c("firebrick", "darkgreen", "purple4", "gray40", "gray60"),
       lty = c(1, 1, 1, 2, 3), pch = c(16, 17, 15, NA, NA), bty = "n")

# Optional: compare the year-8 estimates with the exact logistic solution.
exact_at_8 <- capacity /
  (1 + (capacity / initial_population - 1) *
     exp(-growth_fraction * end_time))

accuracy_results <- data.frame(
  method = c("Euler", "RK2", "RK4"),
  estimate_at_year_8 = c(tail(euler, 1), tail(rk2, 1), tail(rk4, 1)),
  absolute_error = abs(c(tail(euler, 1), tail(rk2, 1),
                         tail(rk4, 1)) - exact_at_8)
)
print(accuracy_results)