adaptive Gauss-Kronrod quadrature with vector-valued integrand
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. 65