#test the accuracy of Simpson's rule

#clear the workspace
rm(list = ls())

source("D:/FaraDars.org/T1365/FVR97041/source/codes/ch2/pr05.r")

ftn <- function(x) return(1/x)
S <- function(n) simpson_n(ftn, 0.01, 1, n)

n.vec <- seq(10, 1000, by = 10)
S.vec <- sapply(n.vec, S)

opar <- par(mfrow = c(1,2), pty = "s", mar = c(4,4,2,1), las = 1)
plot(n.vec, S.vec + log(0.01), type = "l", xlab = "n", ylab = "error")
plot(log(n.vec), log(S.vec + log(0.01)), type = "l", xlab = "log(n)", ylab = "log(error)")
par(opar)
     