polyCub.SV: Product Gauss Cubature over Polygonal Domains

Description

Product Gauss cubature over polygons as proposed by Sommariva and Vianello (2007).

Usage

polyCub.SV(polyregion, f, ..., nGQ = 20, alpha = NULL, rotation = FALSE,
  engine = "C", plot = FALSE)

Value

The approximated value of the integral of f over polyregion.

In the case f = NULL, only the computed nodes and weights are returned in a list of length the number of polygons of polyregion, where each component is a list with nodes (a numeric matrix with two columns), weights (a numeric vector of length nrow(nodes)), the rotation angle, and alpha.

Arguments

polyregion: a polygonal domain. The following classes are supported: "owin" from package spatstat.geom, "gpc.poly" from rgeos (or gpclib), "SpatialPolygons", "Polygons", and "Polygon" from package sp, as well as "(MULTI)POLYGON" from package sf. (For these classes, polyCub knows how to get an xylist.)
f: a two-dimensional real-valued function to be integrated over polyregion (or NULL to only compute nodes and weights). As its first argument it must take a coordinate matrix, i.e., a numeric matrix with two columns, and it must return a numeric vector of length the number of coordinates.
...: further arguments for f.
nGQ: degree of the one-dimensional Gauss-Legendre quadrature rule (default: 20) as implemented in function gauss.quad of package statmod. Nodes and weights up to nGQ=60 are cached in polyCub, for larger degrees statmod is required.
alpha: base-line of the (rotated) polygon at \(x = \alpha\) (see Sommariva and Vianello (2007) for an explication). If NULL (default), the midpoint of the x-range of each polygon is chosen if no rotation is performed, and otherwise the \(x\)-coordinate of the rotated point "P" (see rotation). If f has its maximum value at the origin \((0,0)\), e.g., the bivariate Gaussian density with zero mean, alpha = 0 is a reasonable choice.
rotation: logical (default: FALSE) or a list of points "P" and "Q" describing the preferred direction. If TRUE, the polygon is rotated according to the vertices "P" and "Q", which are farthest apart (see Sommariva and Vianello, 2007). For convex polygons, this rotation guarantees that all nodes fall inside the polygon.
engine: character string specifying the implementation to use. Up to polyCub version 0.4-3, the two-dimensional nodes and weights were computed by R functions and these are still available by setting engine = "R". The new C-implementation is now the default (engine = "C") and requires approximately 30% less computation time.
The special setting engine = "C+reduce" will discard redundant nodes at (0,0) with zero weight resulting from edges on the base-line \(x = \alpha\) or orthogonal to it. This extra cleaning is only worth its cost for computationally intensive functions f over polygons which really have some edges on the baseline or parallel to the x-axis. Note that the old R implementation does not have such unset zero nodes and weights.
plot: logical indicating if an illustrative plot of the numerical integration should be produced.

Author

Sebastian Meyer
These R and C implementations of product Gauss cubature are based on the original MATLAB implementation polygauss by Sommariva and Vianello (2007), which is available under the GNU GPL (>=2) license from https://www.math.unipd.it/~alvise/software.html.

References

Sommariva, A. and Vianello, M. (2007): Product Gauss cubature over polygons based on Green's integration formula. BIT Numerical Mathematics, 47 (2), 441-453. tools:::Rd_expr_doi("10.1007/s10543-007-0131-2")

Examples

Run this code

## a function to integrate (here: isotropic zero-mean Gaussian density)
f <- function (s, sigma = 5)
    exp(-rowSums(s^2)/2/sigma^2) / (2*pi*sigma^2)

## a simple polygon as integration domain
hexagon <- list(
    list(x = c(7.33, 7.33, 3, -1.33, -1.33, 3),
         y = c(-0.5, 4.5, 7, 4.5, -0.5, -3))
)

## image of the function and integration domain
plotpolyf(hexagon, f)

## use a degree of nGQ = 3 and show the corresponding nodes
polyCub.SV(hexagon, f, nGQ = 3, plot = TRUE)

## extract nodes and weights
nw <- polyCub.SV(hexagon, f = NULL, nGQ = 3)[[1]]
nrow(nw$nodes)

## manually apply the cubature rule
sum(nw$weights * f(nw$nodes))

## use an increasing number of nodes
for (nGQ in c(1:5, 10, 20, 60))
    cat(sprintf("nGQ = %2i: %.16f\n", nGQ,
                polyCub.SV(hexagon, f, nGQ = nGQ)))

## polyCub.SV() is the default method used by the polyCub() wrapper
polyCub(hexagon, f, nGQ = 3)  # calls polyCub.SV()


### now using a simple *rectangular* integration domain

rectangle <- list(list(x = c(-1, 7, 7, -1), y = c(-3, -3, 7, 7)))
polyCub.SV(rectangle, f, plot = TRUE)

## effect of rotation given a very low nGQ
opar <- par(mfrow = c(1,3))
polyCub.SV(rectangle, f, nGQ = 4, rotation = FALSE, plot = TRUE)
           title(main = "without rotation (default)")
polyCub.SV(rectangle, f, nGQ = 4, rotation = TRUE,  plot = TRUE)
           title(main = "standard rotation")
polyCub.SV(rectangle, f, nGQ = 4,
           rotation = list(P = c(0,0), Q = c(2,-3)), plot = TRUE)
           title(main = "custom rotation")
par(opar)

## comparison with the "cubature" package
if (requireNamespace("cubature")) {
    fc <- function (s, sigma = 5)  # non-vectorized version of f
        exp(-sum(s^2)/2/sigma^2) / (2*pi*sigma^2)
    cubature::hcubature(fc, lowerLimit = c(-1, -3), upperLimit = c(7, 7))
}

Run the code above in your browser using DataCamp Workspace