polyCub.SV {polyCub}

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)

Arguments

polyregion

a polygonal domain. The following classes are supported: "owin" from package spatstat.geom, "gpc.poly" from 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.

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.

Author(s)

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. doi:10.1007/s10543-007-0131-2

See Also

Other polyCub-methods: polyCub(), polyCub.exact.Gauss(), polyCub.iso(), polyCub.midpoint()

Examples

## 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))
}