SphericalCubature-package {SphericalCubature}R Documentation

Numerical integration over spheres and balls in n-dimensions; multivariate polar/spherical coordinates

Description

Provides functions to integrate a function f(x)=f(x[1],...,x[n]) over the unit sphere and unit ball in n-dimensional Euclidean space:

\int_S f(s) ds \quad \quad \quad \quad \mathrm{and} \quad \quad \quad \int_B f(x) dx,

where the first integral is over the unit sphere S, an (n-1) dimensional surface, and the second integral is over the unit ball B, an n dimensional solid. Vector valued integrands are allowed in some functions; see functions with an argument named fDim.

Integration over general sphers and balls can be done by adjusting the integrand function: if S2 is sphere with center x0 and radius R and B2 is a ball with same center and radius, a simple change of variable shows that

\int_{S2} f(t) dt = \int_S R^{n-1} f(x0+R s) ds \quad \quad \quad \mathrm{and} \quad \quad \int_{B2} f(y) dy = \int_B R^n f(x0+R x) dx.

See the example in adaptIntegrateBallTri.

The package also includes functions to convert to/from polar coordinates in higher dimensions.

There are three cubature methods:

  1. exact methods for polynomials in any dimension (fast)

  2. a method due to Stroud for smooth integrands on the sphere (in dimensions n=3,4,...,16) (slower)

  3. adaptive methods for integrands with different behavior in different regions (slowest)

Methods 2 and 3 are approximations: like any numerical quadrature algorithm, they may give inaccurate results if the integrand changes abruptly on a small region. This happens even in one dimension, and is more difficult to find and deal with in higher dimensions. (One attempt to handle this difficulty is the 'split' versions of the adaptive methods, functions adaptIntegrateSpherePolarSplit() and adaptIntegrateBallPolarSplit(), where one can split the region of integration based on knowledge of the integrand. This can also be done in functions adaptIntegrateSphereTri() and adaptIntegrateBallTri() by specifying an appropriate partition in the simplices.)

It is expected that these methods will yield several significant digits, but not many digits of precision. This seems to be the state of the art in multivariate integration.

Version 1.1 of this package introduces new methods to integrate over spheres. Earlier versions used only polar coordinate representations of the sphere. Now one can use both polar representaions and triangulations of the sphere. The latter has advantages in some cases: it avoids the problems with polar coordinates giving regions that are sometimes rectangles and sometimes triangles (which occurs at the poles), triangles can be approximately equal area in any dimension, etc. While adding these new routines, names became confusing. Apologies to anyone who has trouble because of this, but it seems better in the long run to explicitly name functions based on their approach. Hence adaptIntegrateSphere() has been renamed adaptIntegrateSpherePolar() to indicate that it uses polar coordinates, while the functions adaptIntegrateSphereTri() and adaptIntegrateBallTri() uses spherical triangles.

An explicit goal was to get beyond the cases where n=2, so some efficiency has been sacrificed. In all the methods, the higher the dimension n, the longer the compute time. For methods 2 and 3, compute times get noticeable when n > 5. One application that motivated this package required the ability to work reliably with integrands that have spikes. That requires some sort of adaptive technique, with the possibility of telling the integration algorithm where the spikes are.

This package is an attempt to provide methods for integrating over spheres and balls in multiple dimensions, not a final answer. One possible improvement is speed: coding routines in C would give a significant increase in speed. Another possible extension is to include other multivariate integration methods, e.g. the package R2cuba. This may provide a way to approximate higher dimensional integrals in some cases, if the integrand is well behaved.

Please let me know if you find any mistakes. Constructive comments for improvements are welcome. Fixing bugs or implementing suggestions will be dependent on my workload.

Version history:

Author(s)

John P. Nolan

Maintainer: John P. Nolan <jpnolan@american.edu>

This research was supported by an agreement with Cornell University, Operations Research & Information Engineering, under contract W911NF-12-1-0385 from the Army Research Development and Engineering Command.

See Also

integrateSpherePolynomial, integrateBallPolynomial, integrateSphereStroud11, sphereArea, ballVolume, polar2rect, rect2polar, adaptIntegrateSpherePolar, adaptIntegrateSpherePolarSplit, adaptIntegrateSphereTri, adaptIntegrateSphereTri3d, adaptIntegrateBallPolar, adaptIntegrateBallPolarSplit, adaptIntegrateBallTri, adaptIntegrateBallRadial

Examples

#  integral should just be the area of sphere in n dimensions
f1 <- function( x ) { return(1.0) }
n <- 3
adaptIntegrateBallTri( f1, n )   # exact answer = volume of ball = (4/3)*pi = 4.18879

# other methods
integrateSphereStroud11( f1, n )
adaptIntegrateSpherePolar( f1, n )$integral 
# exact value for a polynomial
p <- list(coef=1.0,k=matrix( rep(0L,n), nrow=1,ncol=n))
integrateSpherePolynomial( p )

# test of exact polynomial integration
f2 <- function( x ) { return(x[1]^2) }
sphereArea(n)/n # exact answer
integrateSphereStroud11( f2, n )
p <- list(coef=1.0,k=matrix( c(2L,rep(0L,n-1)), nrow=1) )
integrateSpherePolynomial( p )
adaptIntegrateSpherePolar( f2, n )$integral

# integration over a ball
adaptIntegrateSphereTri( f1, n ) # exact answer = surface area of sphere = 4*pi = 12.56637

# vector valued integrands
f3 <- function( x ) { c( x[1]^2, x[2]^3 ) }
adaptIntegrateBallTri( f3, n, fDim=2 )
adaptIntegrateSphereTri( f3, n, fDim=2 )

# radial function
g <- function( x ) { sum(x^2) }
adaptIntegrateBallRadial( g, n=3 )

# for more examples enter:   demo(SphericalCubature)
[Package SphericalCubature version 1.5 Index]