# Module 6.2: Euler's Method
# Run the sections in order while coding along with the Quarto page.

# One very large step: P(0) = 100 and dP/dt = 0.1 P.
initial_population <- 100
growth_fraction <- 0.1       # per year
end_time <- 8                # years
delta_t <- 8                 # one step

old_population <- initial_population
growth_rate <- growth_fraction * old_population
euler_estimate <- old_population + growth_rate * delta_t
exact_value <- initial_population * exp(growth_fraction * end_time)

euler_estimate
exact_value
abs(euler_estimate - exact_value)

# Draw the exact curve and the one-step tangent estimate.
smooth_time <- seq(0, end_time, by = 0.1)
exact_curve <- initial_population * exp(growth_fraction * smooth_time)
plot(smooth_time, exact_curve, type = "l", lwd = 2, col = "steelblue",
     xlab = "Years", ylab = "Population", ylim = c(90, 230))
segments(0, initial_population, end_time, euler_estimate,
         col = "firebrick", lwd = 2)
points(c(0, end_time), c(initial_population, euler_estimate),
       col = "firebrick", pch = 19)
legend("topleft", c("Exact curve", "One Euler step"),
       col = c("steelblue", "firebrick"), lty = 1, lwd = 2, bty = "n")

# Recalculate the rate at the start of each two-year step.
# Try 8, 1, and 0.5 in place of 2 after the first run.
delta_t <- 2
time <- seq(0, end_time, by = delta_t)
population <- numeric(length(time))
population[1] <- initial_population

for (i in 2:length(time)) {
  old_population <- population[i - 1]
  growth_rate <- growth_fraction * old_population
  population[i] <- old_population + growth_rate * delta_t
}

data.frame(time, population)
final_estimate <- population[length(population)]
absolute_error <- abs(final_estimate - exact_value)
final_estimate
absolute_error

# Compare the stepped result with the exact values.
exact_at_times <- initial_population * exp(growth_fraction * time)
plot(time, exact_at_times, type = "l", lwd = 2, col = "steelblue",
     xlab = "Years", ylab = "Population", ylim = c(90, 230))
lines(time, population, type = "b", col = "firebrick", pch = 19)
legend("topleft", c("Exact values", "Two-year Euler steps"),
       col = c("steelblue", "firebrick"), lty = 1,
       pch = c(NA, 19), lwd = 2, bty = "n")

# Apply the same method to a different rate: dP/dt = 10 + P/5.
delta_t <- 0.1
time <- seq(0, 0.2, by = delta_t)
population <- numeric(length(time))
population[1] <- 500

for (i in 2:length(time)) {
  old_population <- population[i - 1]
  rate_now <- 10 + old_population / 5
  population[i] <- old_population + rate_now * delta_t
}

data.frame(time, population)