adaptIntegrateSpherePolar: Adaptive integration over sphere or ball in n-dimensions

Description

Approximate the integral over the sphere or ball in n-dimensions using polar coordinates. Can also integrate over sectors of the sphere/ball, see details. These functions will be slow, but may be necessary to get accurate answers if the integrand function f(x) is not smooth. If the integrand changes rapidly in certain regions, the basic routines adaptIntegrateSpherePolar() and codeadaptIntegrateBallPolar() will likely miss these abrupt changes and give inaccurate results. For cases where the location of the rapid changes are known, the functions adaptIntegrateSpherePolarSplit() and adaptIntegrateBallPolarSplit() allow you to split the region of integration and capture those changes.

Usage

adaptIntegrateSpherePolar(f, n, lowerLimit = rep(0, n - 1), 
    upperLimit = c(rep(pi, n - 2), 2 * pi), tol = 1e-05, ...) 
adaptIntegrateSpherePolarSplit(f, n, xstar, width = 0, lowerLimit = rep(0, n - 1), 
    upperLimit = c(rep(pi, n - 2), 2 * pi), tol = 1e-05, ...)
adaptIntegrateBallPolar(f, n, lowerLimit = rep(0, n - 1), 
    upperLimit = c(rep(pi, n - 2), 2 * pi), R = c(0, 1), tol = 1e-05, ...)
adaptIntegrateBallPolarSplit(f, n, xstar, width = 0, lowerLimit = rep(0, n - 1), 
    upperLimit = c(rep(pi, n - 2), 2 * pi), R = c(0, 1), tol = 1e-05, ...)

Arguments

Integrand function f(x)=f(x[1],...,x[n]).

dimension of the space. The sphere is an (n-1) dimensional manifold inside n-space, the ball is an n-dimensional solid.

lowerLimit

Polar angular coordinates of lower limit

upperLimit

Polar angular coordinates of upper limit

tol

tolerance, the desired accuracy of the result. The functions try to get abs(exact-integral) < tol

…

optional arguments passed to f. If used, these should be specified with a tag, e. g. param1=7

a numeric vector of length 2, integration is performed over the region with R[1] < radius < R[2].

xstar

(n x m) matrix whose columns give the directions where the integrand changes quickly, where the region will be subdivided to focus on that region. (The length of a column vector is not used, just it's direction.)

width

width of 'splitting regions', a vector of length m. If it is of length 1, then that value is repeated for each j in 1:m. If width[j]=0, the angular region is split just at the points given by the columns xstar[ ,j]. If width[j] > 0, then angular region is split at an angle plus and minus width[j].

Value

For adaptIntegrateSpherePolar() and adaptIntegrateBallPolar(), the function returns a list containing several fields. There is always a field

integral

Giving the approximate value of the integral.

The other fields depend on the dimension: when n=2, the other fields are what is returned by the function integrate() in base ; when n > 2, the other fields are the fields returned by package cubature.

For adaptIntegrateSpherePolarSplit() and adaptIntegrateBallPolarSplit(), a single value is returned. (This is because these functions make multiple calls to the adaptive integration routine and the results of each call are not saved.

Details

Approximate the integral of f(x) over all/part of the sphere or ball in n-space. The approach is simplistic: reparameterize the region in polar coordinates. For the sphere, this makes the region of integration a rectangle in dimension (n-1) in the angle space (here the radius is fixed: R=1). For the ball, the polar representation in terms of angles and radius gives an region of integration that is an n dimensional rectangle. The R package cubature is used to evaluate the transformed integral.

The region of integration can be a subset of the sphere/ball by specifying a patch/sector in polar coordinates. To integrate over a subregion, bounds for the polar integration have be specified. For example, in two dimensions, you can integrate over the top half of the circle by specifying lowerLimit=0.0 and upperLimit=pi to adaptIntegrateSpherePolar(). Likewise for the ball, to integrate over the part of the annulus with inner radius .2 and outer radius .7 that is in the first quadrant, specify lowerLimit=0.0, upperLimit=pi/2, R=c(.2,.7).

Examples

Run this code

# NOT RUN {
f1 <- function( x ) { return(x[1]^2+3*x[2]+exp(x[3])) }
n <- 3
adaptIntegrateSpherePolar( f1, n )
adaptIntegrateSpherePolarSplit( f1, n, xstar=matrix(c(1,1,1),nrow=3) )
adaptIntegrateBallPolar( f1, n )
adaptIntegrateBallPolarSplit( f1, n, xstar=matrix(c(1,1,1),nrow=3) )

# test of adaptive integration with deliberate splitting
# function f3 has a sharp spike in the direction (1,2),
# elsewhere it has value 1
f3 <- function( x ) {
  x0 <- c(1.0,2.0)/sqrt(5.0)
  dist <- sqrt(sum( (x-x0)^2) )
  y <- 10-5000*dist
  y <- 1 + max(y,0)
  return(y)  }

# no splitting: this straightforward attempt at integration misses
# the spike and sees the integrand as =1 everywhere, so returns 2*pi
n <- 2
adaptIntegrateSpherePolar( f3, n ) 

# deliberate splitting at specified points, but still misses spike
# default width=0 splits the region of integration from [0,2*pi] to [0,a] union [a,2*pi],
# where tan(a)=2/1.
xstar <- matrix( c(1.0,2.0,-1.0,1.0), nrow=2 )
adaptIntegrateSpherePolarSplit( f3, n, xstar=xstar )

# deliberate splitting around specified points, 'smart' choice of width gets the spike
# Here the region of integration is split into [0,a-.01] U [a-.01,a+.01] U [a+.01,2*pi]
adaptIntegrateSpherePolarSplit( f3, n, xstar=xstar, width=c(0.01,0.01) )
# }