one-dimensional integrals

1
2
3
4
5
6
7
8
9
10

trapezoid <- function(fun, a, b, n=100) {
	# numerical integral of fun from a to b
	# using the trapezoid rule with n subdivisions
	# assume a < b and n is a positive integer
	h <- (b-a)/n
	x <- seq(a, b, by=h)
	y <- fun(x)
	s <- h * (y[1]/2 + sum(y[2:n]) + y[n+1]/2)
	return(s)
}

1
2
3
4
5
6
7
8
9
10
11
12
13

simpson <- function(fun, a, b, n=100) {
	# numerical integral using Simpson's rule
	# assume a < b and n is an even positive integer
	h <- (b-a)/n
	x <- seq(a, b, by=h)
	if (n == 2) {
		s <- fun(x[1]) + 4*fun(x[2]) +fun(x[3])
	} else {
		s <- fun(x[1]) + fun(x[n+1]) + 2*sum(fun(x[seq(2,n,by=2)])) + 4 *sum(fun(x[seq(3,n-1, by=2)]))
	}
	s <- s*h/3
	return(s)
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

simpson_v2 <- function(fun, a, b, n=100) {
	# numerical integral using Simpson's rule
	# assume a < b and n is an even positive integer
	if (a == -Inf & b == Inf) {
		f <- function(t) (fun((1-t)/t) + fun((t-1)/t))/t^2
		s <- simpson_v2(f, 0, 1, n)
	} else if (a == -Inf & b != Inf) {
		f <- function(t) fun(b-(1-t)/t)/t^2
		s <- simpson_v2(f, 0, 1, n)
	} else if (a != -Inf & b == Inf) {
		f <- function(t) fun(a+(1-t)/t)/t^2
		s <- simpson_v2(f, 0, 1, n)
	} else {
		h <- (b-a)/n
		x <- seq(a, b, by=h)
		y <- fun(x)
		y[is.nan(y)]=0
		s <- y[1] + y[n+1] + 2*sum(y[seq(2,n,by=2)]) + 4 *sum(y[seq(3,n-1, by=2)])
		s <- s*h/3
	}
	return(s)
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14

f <- function(x) 1/x
#integrate(f, 0.01, 1) == -log(0.01)
S.trapezoid <- function(n) trapezoid(f, 0.01, 1, n)
S.simpson <- function(n) simpson(f, 0.01, 1, n)
 
n <- seq(10, 1000, by = 10)
St <- sapply(n, S.trapezoid)
Ss <- sapply(n, S.simpson)
 
opar <- par(mfrow = c(2, 2))
plot(n,St + log(0.01), type="l", xlab="n", ylab="error", main="Trapezoidal rule")
plot(n,Ss + log(0.01), type="l", xlab="n", ylab="error",main="Simpson's rule")
plot(log(n), log(St+log(0.01)),type="l", xlab="log(n)", ylab="log(error)",main="Trapezoidal rule")
plot(log(n), log(Ss+log(0.01)),type="l", xlab="log(n)", ylab="log(error)",main="Simpson's rule")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

yIntegrate <- function(fun, a, b, tol=1e-8, method= "simpson", verbose=TRUE) {
	# numerical integral of fun from a to b, assume a < b 
	# with tolerance tol
	n <- 4
	h <- (b-a)/4
	x <- seq(a, b, by=h)
	y <- fun(x)
	yIntegrate_internal <- function(y, h, n, method) {
		if (method == "simpson") {
			s <- y[1] + y[n+1] + 4*sum(y[seq(2,n,by=2)]) + 2 *sum(y[seq(3,n-1, by=2)])
			s <- s*h/3
		} else if (method == "trapezoidal") {
			s <- h * (y[1]/2 + sum(y[2:n]) + y[n+1]/2)
		} else {
		}		
		return(s)
	}
 
	s <- yIntegrate_internal(y, h, n, method)
	s.diff <- tol + 1 # ensures to loop at once.
	while (s.diff > tol ) {
		s.old <- s
		n <- 2*n
		h <- h/2
		y[seq(1, n+1, by=2)] <- y ##reuse old fun values
		y[seq(2,n, by=2)] <- sapply(seq(a+h, b-h, by=2*h), fun)
		s <- yIntegrate_internal(y, h, n, method)
		s.diff <- abs(s-s.old)
	}
	if (verbose) {
		cat("partition size", n, "\n")
	}
	return(s)
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

quadrature <- function(fun, a, b, tol=1e-8) {
	# numerical integration using adaptive quadrature
 
	quadrature_internal <- function(S.old, fun, a, m, b, tol, level) {
		level.max <- 100
		if (level > level.max) {
			cat ("recursion limit reached: singularity likely\n")
			return (NULL)
		}
		S.left <- simpson(fun, a, m, n=2) 
		S.right <- simpson(fun, m, b, n=2)
		S.new <- S.left + S.right
		if (abs(S.new-S.old) > tol) {
			S.left <- quadrature_internal(S.left, fun, a, (a+m)/2, m, tol/2, level+1)
			S.right <- quadrature_internal(S.right, fun, m, (m+b)/2, b, tol/2, level+1)
			S.new <- S.left + S.right
		}
		return(S.new)
	}
 
	level = 1
	S.old <- (b-a) * (fun(a) + fun(b))/2
	S.new <- quadrature_internal(S.old, fun, a, (a+b)/2, b, tol, level+1)
	return(S.new)
}