inst/resources/scripts/simpson.r In spuRs: Functions and Datasets for "Introduction to Scientific Programming and Simulation Using R"

```# 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)
}
```

Try the spuRs package in your browser

Any scripts or data that you put into this service are public.

spuRs documentation built on May 21, 2018, 9:03 a.m.