# Module 5.2: Errors
# Run sections in order while following the Quarto page.

# Meeting 1: absolute and relative error
reference <- 6.239
approximation <- 6.24
absolute_error <- abs(approximation - reference)
relative_error <- absolute_error / abs(reference)
absolute_error
100 * relative_error  # percent

# Binary floating-point numbers are not exact decimal fractions.
0.1 + 0.2
0.3
0.1 + 0.2 == 0.3
sprintf("%.17f", 0.1 + 0.2)
sprintf("%.17f", 0.3)
(0.1 + 0.2) - 0.3

# Compare with a tolerance chosen for this small example.
tolerance <- 1e-12
abs((0.1 + 0.2) - 0.3) < tolerance

# The order of subtraction can leave a tiny floating-point remainder.
9.1 - 9 - 0.1
9.1 - 0.1 - 9

# Repeated addition can accumulate a small error.
number_of_steps <- 100000
delta_t <- 0.1
time_by_addition <- 0
for (i in 1:number_of_steps) {
  time_by_addition <- time_by_addition + delta_t
}
time_from_step_number <- number_of_steps * delta_t
time_by_addition
time_from_step_number
time_by_addition - time_from_step_number

# Meeting 2 begins with hand-work on the book's normalized notation:
# 0.0004500 = 0.4500 * 10^-3, precision 4, magnitude 10^-3.
# R's ordinary numeric display need not preserve the written trailing zeros.

# Finite range
1e308 * 10      # overflow: Inf
1e-300 / 1e100  # underflow: 0

# Equivalent mathematical expressions need not produce the same result.
a <- 1e16
b <- -1e16
c <- 1
(a + b) + c
a + (b + c)

# Compare Euler's finite-step approximation with an analytical solution.
# P(0) = 100 and the growth rate is 10% of the current population per hour.
delta_t <- 1
times <- seq(0, 8, by = delta_t)
population <- numeric(length(times))
population[1] <- 100
for (i in 2:length(times)) {
  old_population <- population[i - 1]
  growth_rate <- 0.1 * old_population
  population[i] <- old_population + delta_t * growth_rate
}
one_hour_result <- population[length(population)]
exact_result <- 100 * exp(0.1 * 8)
one_hour_result
exact_result
abs(one_hour_result - exact_result)

# Repeat the same model with half-hour steps.
delta_t <- 0.5
times <- seq(0, 8, by = delta_t)
population <- numeric(length(times))
population[1] <- 100
for (i in 2:length(times)) {
  old_population <- population[i - 1]
  growth_rate <- 0.1 * old_population
  population[i] <- old_population + delta_t * growth_rate
}
half_hour_result <- population[length(population)]
half_hour_result
abs(half_hour_result - exact_result)

# If time permits, change delta_t above to 0.1 and then 0.01.