# Module 6.4: Runge-Kutta 4
# Run the sections in order while coding along with the Quarto page.
# Each change is a full-step trial change, not a separate population update.

# One large step for dP/dt = 0.1 P, starting at P(0) = 100.
old_population <- 100
growth_fraction <- 0.1       # per year
delta_t <- 8                 # years

change_1 <- growth_fraction * old_population * delta_t
change_2 <- growth_fraction * (old_population + change_1 / 2) * delta_t
change_3 <- growth_fraction * (old_population + change_2 / 2) * delta_t
change_4 <- growth_fraction * (old_population + change_3) * delta_t

data.frame(stage = c("start", "first midpoint", "second midpoint", "end"),
           trial_population = c(old_population,
                                old_population + change_1 / 2,
                                old_population + change_2 / 2,
                                old_population + change_3),
           full_step_change = c(change_1, change_2, change_3, change_4))

weighted_change <- (change_1 + 2 * change_2 +
                    2 * change_3 + change_4) / 6
rk4_estimate <- old_population + weighted_change
rk4_estimate

# Compare one-step estimates with the exact growth curve.
smooth_time <- seq(0, 8, by = 0.1)
exact_curve <- 100 * exp(growth_fraction * smooth_time)

plot(smooth_time, exact_curve, type = "l", lwd = 2, col = "steelblue",
     xlab = "Years", ylab = "Population", ylim = c(90, 235))
points(c(8, 8, 8), c(180, 212, rk4_estimate),
       col = c("firebrick", "darkgreen", "purple4"), pch = c(19, 17, 15))
legend("topleft", c("Exact curve", "Euler", "RK2", "RK4"),
       col = c("steelblue", "firebrick", "darkgreen", "purple4"),
       lty = c(1, NA, NA, NA), pch = c(NA, 19, 17, 15), bty = "n")

# Repeat RK4 from year 0 through year 8, using two-year steps.
# Try delta_t = 8, 1, and 0.5 after the first run.
initial_population <- 100
end_time <- 8
delta_t <- 2
time <- seq(0, end_time, by = delta_t)
population_rk4 <- numeric(length(time))
population_rk4[1] <- initial_population

for (i in 2:length(time)) {
  old_population <- population_rk4[i - 1]
  change_1 <- growth_fraction * old_population * delta_t
  change_2 <- growth_fraction * (old_population + change_1 / 2) * delta_t
  change_3 <- growth_fraction * (old_population + change_2 / 2) * delta_t
  change_4 <- growth_fraction * (old_population + change_3) * delta_t

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

data.frame(time, population_rk4,
           exact = initial_population * exp(growth_fraction * time))

# The Euler and RK2 numbers come from earlier modules at the same step size.
exact_at_8 <- initial_population * exp(growth_fraction * end_time)
data.frame(method = c("Euler", "RK2", "RK4", "Exact"),
           estimate = c(207.36, 221.5335,
                        population_rk4[length(time)], exact_at_8),
           absolute_error = c(abs(207.36 - exact_at_8),
                              abs(221.5335 - exact_at_8),
                              abs(population_rk4[length(time)] - exact_at_8), 0))