subtrap <-
function(f, a, b, n=1, ...)
{
h <- (b-a)/n
if(n==1) return( h*mean( f(c(a,b), ...)))
else return(h*sum( f(seq(a+h,b,by=2*h),...) ))
}
#' @export
trap <-
function(f, a, b, tol=1e-8, max.step=1000, ...)
{
i.old <- subtrap(f,a,b,1,...); n <- 2
for(i in 2:max.step) {
s <- subtrap(f, a, b, n, ...)
i.new <- i.old/2 + s
if(abs(i.new-i.old) < tol) break
i.old <- i.new
n <- n*2
}
i.new
}
# simp
#'
#' Numerical integration
#'
#' Perform numerical integration by Simpson's rule or the trapezoidal
#' rule.
#'
#' @aliases trap
#'
#' @param f The integrand; must be a vectorized function.
#'
#' @param a Lower limit of integration.
#'
#' @param b Upper limit of integration.
#'
#' @param tol Tolerance for choosing the number of grid points.
#'
#' @param max.step Log base 2 of the total number of grid points.
#'
#' @param ... Other arguments passed to the integrand, `f`.
#'
#' @details
#' Iterately doubles the number of grid points for the numerical
#' integral, stopping when the integral decreases by less than
#' `tol`.
#'
#' @export
#' @return
#' The integral of `f` from `a` to `b`.
#'
#' @examples
#' f <- function(x) x*x*(1-x)*sin(x*x)
#' I1 <- trap(f,0,2)
#' I2 <- simp(f,0,2)
#'
#' @seealso
#' [stats::integrate()]
#'
#' @keywords
#' math
simp <-
function(f, a, b, tol=1e-8, max.step=1000, ...)
{
i.old <- subtrap(f, a, b, 1, ...)*2/3
n <- 2; old.s <- 0
for(i in 2:max.step) {
s <- subtrap(f, a, b, n, ...)
i.new <- (i.old/2 + (4*s-old.s)/3)
if(abs(i.new-i.old) < tol) break
i.old <- i.new; old.s <- s
n <- n*2
}
i.new
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.