swh:1:snp:81eadaa089e8253d8469bcef66aa332632c6c669
Tip revision: 26e049d70b4a1c237987e260cba68f6a9413736c authored by Hans W. Borchers on 09 April 2019, 04:10:07 UTC
version 2.2.5
version 2.2.5
Tip revision: 26e049d
quadgk.R
##
## g a u s s g k . R Adapitve Gauss-Kronrod
##
quadgk <- function(f, a, b, tol = .Machine$double.eps^0.5, ...) {
stopifnot(is.numeric(a), length(a) == 1,
is.numeric(b), length(b) == 1)
eps <- .Machine$double.eps
fun <- match.fun(f)
f <- function(x) fun(x, ...)
if (a == b) return(0)
else if (a > b) return(-1 * quadgk(f, b, a, tol = tol))
# Nodes and weights for Gauss-Kronrod (7, 15)
n15 <- c(-0.9914553711208126, -0.9491079123427585, -0.8648644233597691,
-0.7415311855993944, -0.5860872354676911, -0.4058451513773972,
-0.2077849550078985, 0.0, 0.2077849550078985,
0.4058451513773972, 0.5860872354676911, 0.7415311855993944,
0.8648644233597691, 0.9491079123427585, 0.9914553711208126)
n7 <- c(-0.9491079123427585, -0.7415311855993944, -0.4058451513773972,
0.0,
0.4058451513773972, 0.7415311855993944, 0.9491079123427585)
w15 <- c(0.02293532201052922, 0.06309209262997855, 0.1047900103222502,
0.1406532597155259, 0.1690047266392679, 0.1903505780647854,
0.2044329400752989, 0.2094821410847278, 0.2044329400752989,
0.1903505780647854, 0.1690047266392679, 0.1406532597155259,
0.1047900103222502, 0.06309209262997855, 0.02293532201052922)
w7 <- c(0.1294849661688697, 0.2797053914892767, 0.3818300505051189,
0.4179591836734694,
0.3818300505051189, 0.2797053914892767, 0.1294849661688697)
.gkadpt <- function(f, a, b, tol = tol) {
# use nodes and weights from the environment
x15 <- 0.5 * ((b - a) * n15 + b + a)
x7 <- 0.5 * ((b - a) * n7 + b + a)
Q7 <- sum(w7 * f(x7)) * (b-a)/2
Q15 <- sum(w15 * f(x15)) * (b-a)/2
if (!is.finite(Q7) || !is.finite(Q15)) {
warning("Infinite or NA function value encountered.")
return(Q15)
} else if (abs(Q15 - Q7) < tol) {
return(Q15)
} else if (abs(b-a) < 16*eps) {
warning("Minimum step size reached; singularity possible.")
return(Q2)
} # else
Q2 <- .gkadpt(f, (a+b)/2, b, tol = tol)
Q1 <- .gkadpt(f, a, (a+b)/2, tol = tol)
return(Q1 + Q2)
}
# start the recursive procedure
.gkadpt(f, a, b, tol = tol)
}