So I'm trying to implement a Gauss-Kronrod adaptive quadrature. That is, I want to calculate
$$\int_a^b f(x) dx = \sum_i f(x_i) w_i$$
where f(x) is evaluated at multiple points at once for efficiency. A key point to bear in mind is that $f$ is also vector-valued (possible risk of confusion here with the double meaning, and use of, vectorised); in practice, I'm integrating at once the real part and imaginary part of a complex integrand, or the integrand could be a 3D vector (electric field), etc.
I've worked out the basic quadrature step in R (ultimately to be ported to C++ using the Armadillo library, any Matlab-like matrix syntax will do for the sake of discussion), where I first divide the integral into N portions, and evaluate the weighted sums corresponding to the Gauss and Kronrod nodes.
## nodes and weights borrowed from Octave source code
abscissa = c(-0.9914553711208126e+00, -0.9491079123427585e+00,
-0.8648644233597691e+00, -0.7415311855993944e+00,
-0.5860872354676911e+00, -0.4058451513773972e+00,
-0.2077849550078985e+00, 0.0000000000000000e+00,
0.2077849550078985e+00, 0.4058451513773972e+00,
0.5860872354676911e+00, 0.7415311855993944e+00,
0.8648644233597691e+00, 0.9491079123427585e+00,
0.9914553711208126e+00)
weights15 = c(0.2293532201052922e-01, 0.6309209262997855e-01,
0.1047900103222502e+00, 0.1406532597155259e+00,
0.1690047266392679e+00, 0.1903505780647854e+00,
0.2044329400752989e+00, 0.2094821410847278e+00,
0.2044329400752989e+00, 0.1903505780647854e+00,
0.1690047266392679e+00, 0.1406532597155259e+00,
0.1047900103222502e+00, 0.6309209262997855e-01,
0.2293532201052922e-01)
weights7 = c(0.1294849661688697e+00, 0.2797053914892767e+00,
0.3818300505051889e+00, 0.4179591836734694e+00,
0.3818300505051889e+00, 0.2797053914892767e+00,
0.1294849661688697e+00)
gauss_kronrod <- function(f, a, b, startN=10, dimf = 2, eps = 1e-3, ...){
Nquad <- length(abscissa)
intervals <- seq(a, b, length=startN)
## change of variable from [a, b] to [-1, 1]
shifts <- (intervals[-length(intervals)] + intervals[-1]) / 2
scalings <- diff(intervals) / 2
## scale the weights accordingly
weightsG <- tcrossprod(weights15, scalings)
weightsK <- tcrossprod(weights7, scalings)
## transformed node positions
x.scaled <- tcrossprod(abscissa, scalings) +
matrix(shifts, ncol=length(shifts), nrow=Nquad, byrow=TRUE)
## evaluate the function at all points
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 (startN - 1) sub-intervals
dim(fvals) = c(dimf, Nquad, startN-1)
## select which subset of fvals are used for the 7-points kronrod sum
ikron <- rep(c(FALSE, TRUE), length.out=Nquad)
gauss <- kronrod <- matrix(0, nrow=dimf, ncol=startN-1)
for (ii in seq.int(startN-1)){ # integrals for each sub-interval
gauss[,ii] <- fvals[,,ii] %*% weightsG[,ii]
kronrod[,ii] <- fvals[,ikron,ii] %*% weightsK[,ii]
}
## relative difference between the two rules
errors <- abs(gauss - kronrod) / abs(gauss)
test <- which(apply(errors, 2, max) > eps)
## net integral, just sum the partial sums...
rowSums(gauss)
}
# a vector-valued integrand...
f <- function(x) rbind(exp(x/(x+2))*sin(x), cos(x))
gauss_kronrod(f, 2, 110, startN=10)
My problem comes when trying to implement an adaptive scheme; basically, at each iteration I would like to compare the values of gauss
vs kronrod
for each of the subintervals, and from there decide where to subdivide further. What would be a good strategy to store these intermediate results during the iteration, and keep track of which subintervals sums go where etc. in the final summation?
More general advice on this endeavour will also be appreciated, I'm not usually keen to implement such generic operations myself, but it seems it might be easier than linking different libraries together (I was contemplating cubature, initially, because my end goal is 2D integration).