add a week of exercises

sose2021/r/w06/P06-1.R

@@ -0,0 +1,14 @@
+# Josua Kugler, Christian Merten
+
+rec_series <- function(start, f, n, m = length(start)) {
+    if (length(start) >= n) return(start)
+    rec_series(c(start, f(tail(start, m))), f, n, m)
+}
+
+rec_series(rep(1, 2), sum, 16)
+rec_series(rep(1, 3), sum, 16)
+rec_series(rep(1, 4), sum, 16)
+rec_series(c(1,1), function(x) 2*x[1]+x[2], 16)
+rec_series(c(1,1), function(x) x[2]-x[1], 16)
+rec_series(1:5, mean, 16)
+rec_series(1:5, median, 16)

sose2021/r/w06/P06-2.R

@@ -0,0 +1,34 @@
+# Josua Kugler, Christian Merten
+
+make_store <- function(n) {
+    dl <- vector(mode="list", length=n)
+    i <- 1
+    function(item=NULL, reset=FALSE) {
+        if (reset) { dl <<- vector(mode="list", length=n); i <<- 1 }
+        if (is.null(item)) return(dl)
+        dl[[i]] <<- item
+        i <<- i+1
+        if (i > n) i <<- 1
+        invisible(dl)
+    }
+}
+
+store3 <- make_store(3)
+store2 <- make_store(2)
+store3(1)
+store3("asdf")
+str(store3())
+store3(list())
+store3(TRUE)
+store3(1:10)
+str(store3())
+store3(NA, reset=TRUE)
+str(store3())
+store3(1:10)
+str(store3())
+store2(1)
+store2(2)
+store2(3)
+str(store2())
+str(store3())
+

sose2021/r/w06/P06-3.R

@@ -0,0 +1,104 @@
+# Josua Kugler, Christian Merten
+
+# a) ----------------------------------------------------------------------- 
+midpoint <- function(f, a, b) (b-a) * f((a+b)/2)
+trapezoid <- function(f, a, b) (b-a) * (f(a) + f(b))/2
+
+# b) ----------------------------------------------------------------------
+nc_integrate <- function(f, lower, upper, n, rule) {
+    xs <- seq(lower, upper, length.out=n+1)
+    sum(rule(f, xs[-(n+1)], xs[-1]))
+}
+
+# c) ------------------------------------------------------------------------
+closeable_seq <- function(from, to, len, closed=TRUE) {
+    if (closed) seq(from, to, length.out = len)
+    else seq(from, to, length.out = len+2)[-c(1, len+2)]
+}
+
+newton_cotes <- function(coef, closed=TRUE) {
+    m <- length(coef)
+    w <- coef / sum(coef)
+    function(f, a, b) {
+        t <- mapply(closeable_seq, a, b, len=m, closed=closed)
+        (b-a) * colSums(w*matrix(f(t), ncol=length(a)))
+    }
+}
+
+# d) ------------------------------------------------------------------------
+
+# some objects for param_fun-list
+sin1x <- function(x) {
+  y <- suppressWarnings(sin(1/x))
+  y[is.na(y)] <- 0
+  y
+}
+set.seed(0)
+x <- c(0:1, runif(5))
+y <- runif(7)
+
+# list of functions with integral interval
+param_fun <- list(
+  poly = list(f = function(x) x^4 - x^3 - 3*x^2 + x + 2, lower = 0, upper = 2),
+  sin1x = list(f = sin1x, lower = 0, upper = 1),
+  lin = list(f = approxfun(x, y), lower = 0, upper = 1),
+  spline = list(f = splinefun(x, y), lower = 0, upper = 1)
+)
+
+# true values of integrals
+true <- c(
+  poly = 0.4, 
+  sin1x = 0.504067061906928, 
+  lin = 0.472878602164825, 
+  spline = 0.97236924451286)
+
+# options for creating integral-rules
+param_rule <- list(
+  midpoint = list(coef = 1, closed=FALSE),
+  trapezoid = list(coef = c(1,1), closed=TRUE),
+  simpson = list(coef = c(1,4,1), closed=TRUE),
+  boole = list(coef = c(7,32,12,32,7), closed=TRUE),
+  open5 = list(coef = c(611, -453, 562, 562, -453, 611), closed=FALSE)
+)
+
+# how many function evaluations are allowed
+param_n <- 5*(2:20)
+
+library(tidyverse)
+
+param <- as_tibble(
+  expand.grid(n = param_n,
+              fun_name = names(param_fun),
+              rule_name = names(param_rule),
+              stringsAsFactors=FALSE))
+
+nc_rules <- sapply(param_rule, do.call, what=newton_cotes)
+
+# run nc_integrate with given options
+call_nc <- function(n, fun_name, rule_name) {
+  opts <- param_fun[[fun_name]]
+  opts$rule = nc_rules[[rule_name]]
+  # to make things fair:
+  # a rule which evaluates f at m points may only be called n/m times
+  opts$n = round(n / length(param_rule[[rule_name]]$coef))
+  do.call(nc_integrate, opts)
+}
+
+param %>% rowwise() %>%
+    mutate(value=nc_integrate(param_fun[[fun_name]]$f,
+                              lower = param_fun[[fun_name]]$lower,
+                              upper = param_fun[[fun_name]]$upper,
+                              n = n,
+                              rule = nc_rules[[rule_name]]),
+           true=true[[fun_name]],
+           error=abs(true-value)) -> res
+
+# plot results
+plots <- lapply(names(param_fun), function(nm) 
+  res %>%
+    filter(fun_name == nm) %>%
+    ggplot(aes(x = n, y = error, color = rule_name)) +
+    scale_y_log10() +
+    geom_line() + geom_point() + labs(title = nm)
+)
+gridExtra::grid.arrange(grobs = plots, nrow=2)