# ---------------------------
# Part 1: Numerical methods
# ---------------------------
# (1) Bisection method - generic implementation
bisection <- function(f, a, b, tol = 1e-6, max.iter = 1000) {
if (f(a) * f(b) > 0) stop("f(a) and f(b) must have opposite signs")
fa <- f(a); fb <- f(b)
iter <- 0
while ((b - a) / 2 > tol && iter < max.iter) {
c <- (a + b) / 2
fc <- f(c)
iter <- iter + 1
if (fc == 0) {
a <- b <- c
break
}
if (fa * fc < 0) {
b <- c
fb <- fc
} else {
a <- c
fa <- fc
}
}
root <- (a + b) / 2
list(root = root, froot = f(root), iterations = iter, interval = c(a, b))
}
# functions for (1)
f1 <- function(x) x^3 - 4*x - 9 # 1a
f2 <- function(x) x^2 - 5 # 1b
# Find sign-changing intervals first (simple scan) and apply bisection
# For f1: scan -10..10
scan_interval <- function(f, from = -50, to = 50, step = 0.5) {
xs <- seq(from, to, by = step)
segs <- list()
for (i in seq_along(xs)[-1]) {
if (f(xs[i-1]) * f(xs[i]) <= 0) {
segs[[length(segs) + 1]] <- c(xs[i-1], xs[i])
}
}
segs
}
cat("Bisection method results:\n")
# f1
segs1 <- scan_interval(f1, -10, 10, 0.5)
cat(" f1 sign-change intervals found:\n"); print(segs1)
res1 <- bisection(f1, segs1[[1]][1], segs1[[1]][2], tol = 1e-8)
print(res1)
# f2
segs2 <- scan_interval(f2, -10, 10, 0.5)
cat(" f2 sign-change intervals found:\n"); print(segs2)
# choose positive root
res2 <- bisection(f2, segs2[[length(segs2)]][1], segs2[[length(segs2)]][2], tol = 1e-8)
print(res2)
# (2) Newton-Raphson for f(x) = x^5 + 3x^2 + x - 10
f_nr <- function(x) x^5 + 3*x^2 + x - 10
fprime_nr <- function(x) 5*x^4 + 6*x + 1 # derivative f'(x)
newton_raphson <- function(f, df, x0, tol = 1e-6, max.iter = 100) {
x <- x0
iter <- 0
while (iter < max.iter) {
fx <- f(x); dfx <- df(x)
if (abs(dfx) < .Machine$double.eps) stop("Derivative near zero.")
xnew <- x - fx/dfx
iter <- iter + 1
if (abs(xnew - x) < tol) break
x <- xnew
}
list(root = xnew, froot = f(xnew), iterations = iter)
}
cat("\nNewton-Raphson (start x0 = 1.25):\n")
nr_res <- newton_raphson(f_nr, fprime_nr, x0 = 1.25, tol = 1e-8, max.iter = 100)
nr_res$root_3dp <- round(nr_res$root, 3)
print(nr_res)
# ---------------------------
# Part 2: Cooling dataset & Newton's Law of Cooling
# ---------------------------
# Provided data
time <- seq(0, 20, by = 2)
temp <- c(85.0, 76.4, 69.1, 63.2, 58.4, 54.3, 50.9, 47.9, 45.4, 43.2, 41.5)
data <- data.frame(time = time, temp = temp)
T_ambient <- 22
T0 <- temp[1] # 85
# (a) Linearised fit: ln(T - T_ambient) = ln(A) - k t
# where A = T0 - T_ambient
data$TminusA <- data$temp - T_ambient
if (any(data$TminusA <= 0)) stop("T - ambient is <= 0 for some points; cannot take log")
data$lnTminusA <- log(data$TminusA)
lm_fit <- lm(lnTminusA ~ time, data = data)
lm_coef <- coef(lm_fit)
intercept_lm <- lm_coef[1] # ln(A_hat)
slope_lm <- lm_coef[2] # -k_hat
A_hat_from_lm <- exp(intercept_lm)
k_hat_from_lm <- -slope_lm
cat("\nLinear regression on ln(T - T_ambient) ~ time results:\n")
cat(sprintf(" intercept (ln A) = %g, so A_hat = %g\n", intercept_lm, A_hat_from_lm))
cat(sprintf(" slope = %g => k_hat = %g per minute\n", slope_lm, k_hat_from_lm))
cat(sprintf(" Implied model: T(t) = %g + %g * exp(-%g * t)\n", T_ambient, A_hat_from_lm, k_hat_from_lm))
# (b) Non-linear least squares (nls) to refine A and k
# Model: temp ~ T_ambient + A * exp(-k * time)
# Starting values from linear fit
start_list <- list(A = A_hat_from_lm, k = k_hat_from_lm)
nls_fit <- nls(temp ~ T_ambient + A * exp(-k * time),
data = data,
start = start_list,
control = nls.control(maxiter = 200, tol = 1e-8, minFactor = 1/1024))
nls_coef <- coef(nls_fit)
A_nls <- unname(nls_coef["A"])
k_nls <- unname(nls_coef["k"])
cat("\nNLS fit results:\n")
cat(sprintf(" A (amplitude) = %g\n k = %g per minute\n", A_nls, k_nls))
cat(sprintf(" NLS model: T(t) = %g + %g * exp(-%g * t)\n", T_ambient, A_nls, k_nls))
# (c) Predictions, % error, RMSE, R^2, residuals
data$pred_lm <- T_ambient + A_hat_from_lm * exp(-k_hat_from_lm * data$time)
data$pred_nls <- T_ambient + A_nls * exp(-k_nls * data$time)
# choose nls predictions for final table (more accurate)
data$pred <- data$pred_nls
data$residual <- data$temp - data$pred
data$percent_error <- abs(data$residual) / data$temp * 100
# RMSE
rmse <- sqrt(mean(data$residual^2))
# R-squared (use 1 - SSE/SST on original temp vs model predictions)
SSE <- sum(data$residual^2)
SST <- sum((data$temp - mean(data$temp))^2)
R2 <- 1 - SSE/SST
cat("\nModel performance (NLS model):\n")
cat(sprintf(" RMSE = %g °C\n", rmse))
cat(sprintf(" R-squared = %g\n", R2))
# Print table comparing actual vs predicted and % error
cat("\nTime | Actual | Predicted | Residual | %Error\n")
print(round(data[, c("time", "temp", "pred", "residual", "percent_error")], 4))
# (d) Residual pattern check (print residuals & simple summary)
cat("\nResidual summary:\n"); print(summary(data$residual))
cat("\nResiduals by time (first few):\n"); print(round(data[, c("time", "residual")], 4))
# (e) Final model & parameters (rounded)
cat("\nFinal reported model (rounded):\n")
cat(sprintf(" T(t) = %g + %g * exp(-%g * t)\n", T_ambient, round(A_nls, 4), round(k_nls, 5)))
cat(sprintf(" k (decay constant) = %g per minute\n", round(k_nls, 5)))
cat(sprintf(" A = %g (should be approx T0 - T_ambient = %g)\n", round(A_nls,4), T0 - T_ambient))
# Also show linear-fit based k for reference
cat(sprintf("\nLinearized fit gave k = %g per minute (for comparison)\n", k_hat_from_lm))
# Save a neat table for potential copy-paste in a report
final_table <- data.frame(
time = data$time,
actual_T = data$temp,
predicted_T = round(data$pred, 4),
residual = round(data$residual, 4),
percent_error = round(data$percent_error, 3)
)
cat("\nFinal table (actual vs predicted):\n")
print(final_table)
# End of script