adaptIntegrateSphereTri {SphericalCubature}

Adaptive integration over spherical triangles

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