adaptIntegrateBallTri: Adaptive integration over the unit ball

In SphericalCubature: Numerical Integration over Spheres and Balls in n-Dimensions; Multivariate Polar Coordinates

Description

Adaptively integrate a function over the ball, specified by a set of spherical triangles to define a ball (or a part of a ball). Function adaptIntegrateBallTri() uses spherical triangles and works in n-dimensions; it uses function adaptIntegrateSimplex() in R package SimplicialCubature, which is based on code of Alan Genz. adaptIntegrateBallRadial() integrates radial functions of the form f(x) = g(|x|) over the unit ball.

Usage

adaptIntegrateBallTri( f, n, S=Orthants(n), fDim=1L, maxEvals=20000L, absError=0.0, 
    tol=1.0e-5, integRule=3L, partitionInfo=FALSE, ...  ) 
adaptIntegrateBallRadial( g, n, fDim=1, maxEvals=20000L, absError=0.0, 
    tol=1e-05, integRule=3L, partitionInfo=FALSE, ... )

Arguments

f	integrand function f defined on the sphere in R^n
g	inegrand function g defined on the real line
n	dimension of the space
S	array of spherical triangles, dim(S)=c(n,n,nS). Columns of S should be points on the unit sphere: sum(S[,i,j]^2)=1. Execution will be faster if every simplex S[,,j] is contained within any single orthant. This will happend automatically if function Orthants is used to generate orthants, or if S is a tessellation coming from function UnitSphere in package mvmesh. If one or more simplices interect multiple orthants, the simplices will automatically be subdivided so that each subsimplex is in a single orthant.
fDim	integer dimension of the integrand function
maxEvals	maximum number of evaluations allowed
absError	desired absolute error
tol	desired relative tolerance
integRule	integration rule to use in call to function adsimp
partitionInfo	if TRUE, return the final partition after subdivision
...	optional arguments to function f(x,...) or g(x,...)

Details

adaptIntegrateBallTri() takes as input a function f defined on the unit sphere in n-dimensions and a list of spherical triangles S and attempts to integrate f over (part of) the unit sphere described by S. It uses the R package SimplicialCubature to evaluate the integrals. The spherical triangles in S should individually be contained in an orthant.

If the integrand is nonsmooth, you can specify a set of spherical triangles that focus the cubature routines on that region. See the example below with integrand function f3.

Value

A list containing

status	a string describing result, ideally it should be "success", otherwise it is an error/warning message.
integral	approximation to the value of the integral
I0	vector of approximate integral over each triangle in K
numRef	number of refinements
nk	number of triangles in K
K	array of spherical triangles after subdivision, dim(K)=c(3,3,nk)
est.error	estimated error
subsimplices	if partitionInfo=TRUE, this gives an array of subsimplices, see function adsimp for more details.
subsimplicesIntegral	if partitionInfo=TRUE, this array gives estimated values of each component of the integral on each subsimplex, see function adsimp for more details.
subsimplicesAbsError	if partitionInfo=TRUE, this array gives estimated values of the absolute error of each component of the integral on each subsimplex, see function adsimp for more details.
subsimplicesVolume	if partitionInfo=TRUE, vector of m-dim. volumes of subsimplices; this is not d-dim. volume if m < n.

Examples

# integrate over ball in R^3
n <- 3
f <- function( x ) { x[1]^2  }
adaptIntegrateBallTri( f, n )

# integrate over first orthant only
S <- Orthants( n, positive.only=TRUE ) 
a <- adaptIntegrateSphereTri( f, n, S )
b <- adaptIntegrateSphereTri3d( f, S )
# exact answer, adaptIntegrateSphereTri approximation, adaptIntegrateSphereTri3d approximation
sphereArea(n)/(8*n); a$integral; b$integral

# integrate a vector valued function
f2 <- function( x ) { c(x[1]^2,x[2]^3) }
adaptIntegrateBallTri( f2, n=2, fDim=2 )

# example of specifiying spherical triangles that make the integration easier
f3 <- function( x ) { sqrt( abs( x[1]-3*x[2] ) ) } # has a cusp along the line x[1]=3*x[2]
a <- adaptIntegrateBallTri( f3, n=2, absError=0.0001, partitionInfo=TRUE )
str(a) # note that returnCode = 1, e.g. maxEvals exceeded

# define problem specific spherical triangles, using direction of the cusp
phi <- atan(1/3); c1 <- cos(phi); s1 <- sin(phi)
S <- array( c(1,0,c1,s1,  c1,s1,0,1,    0,1,-1,0,   -1,0,-c1,-s1,  
            -c1,-s1,0,-1, 0,-1,1,0), dim=c(2,2,6) )
b <- adaptIntegrateBallTri( f3, n=2, S, absError=0.0001, partitionInfo=TRUE )
str(b) # here returnCode=0, less than 1/2 the function evaluations and smaller estAbsError

# integrate x[1]^2 over nested balls of radius 2 and 5 (see discussion in ?SphericalCubature)
f4 <- function( x ) { c( 2^2 * (2*x[1])^2, 5^2 * (5*x[1])^2 ) }
bb <- adaptIntegrateBallTri( f4, n=2, fDim=2 )
str(bb)
bb$integral[2]-bb$integral[1] # = integral of x[1]^2 over the annulus 2 <= |x| <= 5

SphericalCubature documentation built on Jan. 13, 2021, 1:04 p.m.