Well, I've thought of a strategy, which is the following. First, I create a function to take care of calculating the partial sums over a list of intervals, and checking which are good enough to keep. 416
Second, for those intervals that were not good, a helper function splits them in two. Finally, at each step the good contributions are added to the final result. Here's an implementation in R,
gk_segments <- function(f, segments=list(c(2, 50), c(50, 100)), dimf = 2,
eps = 1e-3, ...){
Nquad <- length(abscissa)
Nsegments <- length(segments)
shifts <- sapply(segments, sum) / 2
scalings <- sapply(segments, diff) / 2
weightsG <- tcrossprod(weights15, scalings)
weightsK <- tcrossprod(weights7, scalings)
x.scaled <- tcrossprod(abscissa, scalings) +
matrix(shifts, ncol=length(shifts), nrow=Nquad, byrow=TRUE)
fvals <- f(c(x.scaled), ...)
## function evaluations are reshaped into a 3D array
## rows correspond to dimf values of integrand
## columns are the Nquad evalution points in the sub-intervals
## slices correspond to the Nsegments sub-intervals
dim(fvals) = c(dimf, Nquad, Nsegments)
ikron <- rep(c(FALSE, TRUE), length.out=Nquad)
gauss <- <PERSON> <- matrix(0, nrow=dimf, ncol=Nsegments)
for (ii in seq.int(Nsegments)){ # integrals for each sub-interval
gauss[,ii] <- fvals[,,ii] %*% weightsG[,ii]
kronrod[,ii] <- fvals[,ikron,ii] %*% weightsK[,ii]
}
errors <- abs(gauss - kronrod) / abs(gauss)
test <- which(apply(errors, 2, max) > eps)
if(!length(test))
return( list(converged = rowSums(gauss),
todo = list()))
list(converged = rowSums(gauss[ ,-test]),
todo = segments[test])
}
# f <- function(x) rbind(exp(x/(x+2))*sin(x), cos(x))
# gk_segments(f, segments=split_interval(c(2, 110), 10))
# split a given interval [x0, x1] into a list of N subintervals
split_interval <- function(x, N=2){
breaks <- seq(x[1], x[2], length=N+1)
lapply(seq.int(N), function(ii) c(breaks[ii], breaks[ii+1]))
}
gk_adaptive <- function(fun, a, b, startN=10, dimf = 2,
eps = 1e-3, maxiter = 100, ...){
## first step
segments <- split_interval(c(a, b), startN)
tmp <- gk_segments(f=fun, segments=segments, dimf=dimf, eps=eps, ...)
result <- tmp$converged # partial sums that have already converged
segments <- unlist(lapply(tmp$todo, split_interval, N=2), recursive=FALSE)
iter <- 1
## we stop if there's no subinterval left to split
## or we've reached the maximum number of iterations
finished <- (!length(tmp$todo)) || (iter > maxiter)
while(!finished){
message("iteration #", iter)
## work on the new subsegments
tmp <- gk_segments(f=fun, segments=segments, dimf=dimf, eps=eps, ...)
result <- result + tmp$converged ## add the good chunks
segments <- unlist(lapply(tmp$todo, split_interval, N=2), recursive=FALSE)
iter <- iter + 1
stopping <- iter > maxiter
if(stopping)
message("Maximum number of iteration reached, result is likely wrong")
finished <- (! 416