Nothing
# program spuRs/resources/scripts/simpson.r
simpson <- function(ftn, a, b, tol = 1e-8, verbose = FALSE) {
# numerical integral of ftn from a to b
# using Simpson's rule with tolerance tol
#
# ftn is a function of a single variable and a < b
# if verbose is TRUE then n is printed to the screen
# initialise
n <- 4
h <- (b - a)/4
fx <- sapply(seq(a, b, by = h), ftn)
S <- sum(fx*c(1, 4, 2, 4, 1))*h/3
S.diff <- tol + 1 # ensures we loop at least once
# increase n until S changes by less than tol
while (S.diff > tol) {
# cat('n =', n, 'S =', S, '\n') # diagnostic
S.old <- S
n <- 2*n
h <- h/2
fx[seq(1, n+1, by = 2)] <- fx # reuse old ftn values
fx[seq(2, n, by = 2)] <- sapply(seq(a+h, b-h, by = 2*h), ftn)
S <- h/3*(fx[1] + fx[n+1] + 4*sum(fx[seq(2, n, by = 2)]) +
2*sum(fx[seq(3, n-1, by = 2)]))
S.diff <- abs(S - S.old)
}
if (verbose) cat('partition size', n, '\n')
return(S)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.