adaptIntegrateSphereTri: Adaptive integration over spherical triangles

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

Description

Adaptively integrate a function over a set of spherical triangles. Function adaptIntegrateSphereTri() uses spherical triangles and works in n-dimensions; it uses function adaptIntegrateSimplex() in package SimplicialCubature, which is based on code of Alan Genz. adaptIntegrateSphereTri3d() works only in 3-dimensions and is described in the paper by N. Boal and F-J. Sayas; it is not as sophisticated and is slower than adaptIntegrateSphereTri(), but can be useful and is self contained.

Usage

adaptIntegrateSphereTri( f, n, S, fDim=1L, maxEvals=20000L, absError=0.0, 
    tol=1.0e-5, integRule=3L, partitionInfo=FALSE, ...  ) 
adaptIntegrateSphereTri3d( f, S, maxRef=50, relTol=0.001, maxTri=50000, gamma=1.5 ) 
adaptIntegrateSphereTriI0( f, K )
adaptIntegrateSphereTriI1( f, K )
adaptIntegrateSphereTriSubdivideK( K ) 

Arguments

f	function f defined on the sphere
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.
n	dimension of the ball
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,...)
maxRef	maximum number of refinements allowed
relTol	desired relative tolerance
maxTri	maximum number of triangles allowed while refining
gamma	threshold parameter for when a triangle is subdivided, should be >= 1
K	a single spherical triangle in 3 space, dim(K)=c(3,3)

Details

adaptIntegrateSphereTri3d() takes as input a function f defined on the unit sphere in 3-dimensions and a list of spherical triangles K0 and attempts to integrate f over the part of the unit sphere described by K0. The triangles in K0 should individually contained in an orthant. The algorithm estimates the integral over each triangle, then estimates error on each triangle, and subdivides triangles with large errors. This is repeated until either the desired accuracy is achieved, or there are too many subdivisions.

Functions adaptIntegrateSphereI0(), adaptIntegrateSphereI1(), and adaptIntegrateSphereSubdivideK() are internal functions that compute the integral I0 for a spherical triangle, integral I1 for a spherical triangle, and subdivide spherical triangle K into 4 smaller spherical triangles (using midpoints of each side).

adaptIntegrateSphereTri() 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. The algorithm estimates the integral over each triangle, then estimates error on each triangle, and subdivides triangles with large errors. This is repeated until either the desired accuracy is achieved, or there are too many subdivisions. This function is more general than adaptIntegrateSphereTri3d in two ways: (1) it works in dimension n >= 2, and (2) it allows vector integrands f. It also is generaly faster than adaptIntegrateSphereTri(). Use the function Orthants() to get a rough triangulation of the sphere; see the examples below. For finer resolution triangulation can be obtained from UnitSphere() in R package mvmesh. For description of adaptIntegrateSphereTri3d() see www.unizar.es/galdeano/actas_pau/PDFVIII/pp61-69.pdf.

Thoughtful choice of the spherical trianges can result in faster and more accurate results, see the example below with 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

# test of polynomial integration f(s)=s[1]^2
f <- function( s ) { return( s[1]^2 ) }
n <- 3

# integrate over whole sphere
S <- Orthants( n )
a <- adaptIntegrateSphereTri( f, n, S )
b <- adaptIntegrateSphereTri3d( f, S )
# exact answer, adaptIntegrateSphereTri approximation, adaptIntegrateSphereTri3d approximation
sphereArea(n)/n; a$integral; b$integral

# 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

# 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 <- adaptIntegrateSphereTri( f3, n=2, partitionInfo=TRUE )
str(a)

# 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 <- adaptIntegrateSphereTri( f3, n=2, S, partitionInfo=TRUE )
str(b) # estAbsError here is 1/6 of above with approx. the same number of function evaluations

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