adaptIntegrateBallTri function

Adaptive integration over the unit ball

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

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

Returns

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