7
$\begingroup$

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).

Improve this question
$\endgroup$
  • 1
    $\begingroup$ See scicomp.stackexchange.com/questions/159/… for discussion on when to use libraries. I've seen people mess up simple numerical algorithms like a non-adaptive second-order predictor-corrector integrator that can be written in less than 10 lines of code, so I am biased in favor of using libraries whenever possible. $\endgroup$ – Geoff Oxberry Dec 17 '13 at 19:41
  • $\begingroup$ I agree wholeheartedly, and this is really a last resort for me, when I fail to use existing libraries. I was surprised to get promising results with so little code, however, in a matrix-like syntax that I can understand better, and perhaps trust more, than glue code I would have to write to link two independent c++ libraries. $\endgroup$ – baptiste Dec 17 '13 at 19:51
  • $\begingroup$ linking armadillo and cubature may not be so hard after all, there's hope on that front. Extending (properly) this type of embedded quadrature to N>1 dimensions, however, doesn't seem like an easy job (Genz-Malik papers, etc.), unless one trades-off efficiency for ease of implementation with chained 1D integrations. $\endgroup$ – baptiste Dec 18 '13 at 16:11
3
$\begingroup$

I suggest to read the Quadpack book (Quadpack, Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. ISBN 978-3-540-12553-2) and Pedro Gonnet's PhD thesis "Adaptive Quadrature Re-revisited (available as pdf here). Pedro is a contributor to SciComp SE.

What typically is done is a priority queue data structure:

  1. You apply Gauss-Kronrod on the whole initial interval and estimate the error
  2. You split the interval symmetrically in two and again apply Gauss-Kronrod on the intervals, estimate the errors and add a pair (or tuple) "interval, error" to the priority queue
  3. If the total error is smaller than the requested tolerance (or if you reached the maximal amount of subdivisions), you stop
  4. If not, you pop the top element from the queue (highest error), you split it in two and apply Gauss-Kronrod, add those two new elements to the priority queue and go back to 3.
Improve this answer
$\endgroup$
  • $\begingroup$ Thanks, yes, I've been reading the original literature on adaptive cubature, and it's not a trivial undertaking, especially in more than one dimension. For curiosity's sake, I will improve my 1D code above with some kind of priority queue, but beyond 1D I will have to bite the bullet and somehow link my code to the cubature library. I note that picking just one item at a time in the priority queue doesn't seem optimally efficient when the integrand can be called in a vectorised manner, I'd rather deal with a number of "bad" intervals at a time (deciding how many is tricky). $\endgroup$ – baptiste Dec 19 '13 at 12:11
  • $\begingroup$ @baptiste then you'll need some "measure" on the "badness" of the approximation on the (tensorized) intervals. For cubature and references, you can consult nines.cs.kuleuven.be/research/ecf $\endgroup$ – GertVdE Dec 19 '13 at 15:23
0
$\begingroup$

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. 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 <- kronrod <- 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 <- (!length(tmp$todo)) || stopping
  }

    return(result)
}

f <-  function(x) rbind(sin(x), cos(x))
gk_adaptive(f, 2, 110)
gk_adaptive(f, 2, 110, startN=10, eps=1e-10)
gk_adaptive(f, 2, 110, startN=3, eps=1e-10, maxiter=2)

I don't think it will be easy to generalise to 2 or more dimensions, because choosing the right square / cubes to focus on iteratively seems harder to code. As a less efficient alternative, I can just chain multiple integrations in 1D, keeping all but one of the variables constant at a time.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.