swh:1:snp:81eadaa089e8253d8469bcef66aa332632c6c669
Tip revision: 113add4c7a9bfd37466a9b41cf2dc060b7b87091 authored by HwB on 06 December 2012, 00:00:00 UTC
version 1.3.1
version 1.3.1
Tip revision: 113add4
quadgr.R
##
## q u a d g r . R Gauss-Richardson Quadrature
##
quadgr <- function(f, a, b, tol = .Machine$double.eps^(1/2), ...) {
stopifnot(is.numeric(a), length(a) == 1,
is.numeric(b), length(b) == 1)
fun <- match.fun(f)
f <- function(x) fun(x, ...)
# check order of limits
if (a == b) return(list(value = 0, rel.err = 0))
else if (a > b) {
tmp <- a; a <- b; b <- tmp
rev_p <- TRUE
} else
rev_p <- FALSE
# Check infinite limits
if (is.infinite(a) || is.infinite(b)) {
if (is.finite(a) && is.infinite(b)) {
f1 <- function(t) f(a + t/(1-t)) / (1-t)^2
Q <- quadgr(f1, 0, 1, tol = tol)
} else if (is.infinite(a) && is.finite(b)) {
f2 <- function(t) f(b + t/(1+t)) / (1+t)^2
Q <- quadgr(f2, -1, 0, tol = tol)
} else if (is.infinite(a) && is.infinite(b)) {
f1 <- function(t) f(t/(1-t)) / (1-t)^2
f2 <- function(t) f(t/(1+t)) / (1+t)^2
Q1 <- quadgr(f1, 0, 1, tol = tol/2)
Q2 <- quadgr(f2, -1, 0, tol = tol/2)
Q <- list(value = Q1$value + Q2$value,
rel.err = Q1$rel.err + Q2$rel.err)
}
if (rev_p) Q$value <- -Q$value
return(Q)
} # else ...
# 12-point Gauss-Legendre quadrature
xq <- c(0.12523340851146894, 0.36783149899818018, 0.58731795428661748,
0.76990267419430469, 0.9041172563704748, 0.98156063424671924)
wq <- c(0.24914704581340288, 0.23349253653835478, 0.20316742672306584,
0.16007832854334636, 0.10693932599531818, 0.047175336386511842)
xq <- matrix(c(xq, -xq), ncol = 1)
wq <- c(wq, wq)
nq <- length(xq)
# Initiate vectors
maxit <- 17 # max number of iterations
Q0 <- zeros(maxit,1) # quadrature
Q1 <- zeros(maxit,1) # first Richardson extrapolation
Q2 <- zeros(maxit,1) # second Richardson extrapolation
# One interval
hh <- (b - a)/2 # half interval length
x <- (a + b)/2 + hh*xq # nodes
Q0[1] = hh * wq %*% fun(x) # quadrature
for (k in 3:maxit) {
hh <- hh/2
x <- cbind(x + a, x + b) / 2
Q0[k] <- hh * wq %*% apply(f(x), 1, sum)
# Richardson extrapolation
if (k >= 5) {
Q1[k] <- .rich(Q0,k)
Q2[k] <- .rich(Q1,k)
} else if (k >= 3) {
Q1[k] <- .rich(Q0,k)
}
# Estimate absolute error
if (k >= 6) {
Qv <- c(Q0[k], Q1[k], Q2[k])
Qw <- c(Q0[k-1], Q1[k-1], Q2[k-1])
} else if (k >= 4) {
Qv <- c(Q0[k], Q1[k])
Qw <- c(Q0[k-1], Q1[k-1])
} else {
Qv <- Q0[k]
Qw <- Q0[k-1]
}
err <- min(abs(Qv - Qw))
j <- which.min(abs(Qv - Qw))
Q <- Qv[j]
# Convergence
if (err < tol || !is.finite(Q)) break
}
if (rev_p) Q <- -Q
return(list(value = Q, rel.err = err))
}
.rich <- function(Q, k) {
if (Q[k] != Q[k-1]) {
cc <- (Q[k-1] - Q[k-2]) / (Q[k] - Q[k-1]) - 1
} else {
cc <- 1
}
cc <- max(cc, 0.07)
return(Q[k] + (Q[k] - Q[k-1])/cc)
}