polyCub.SV  R Documentation 
Product Gauss cubature over polygons as proposed by Sommariva and Vianello (2007).
polyCub.SV(polyregion, f, ..., nGQ = 20, alpha = NULL, rotation = FALSE,
engine = "C", plot = FALSE)
polyregion 
a polygonal domain.
The following classes are supported:

f 
a twodimensional realvalued function to be integrated over

... 
further arguments for 
nGQ 
degree of the onedimensional GaussLegendre quadrature rule
(default: 20) as implemented in function 
alpha 
baseline of the (rotated) polygon at 
rotation 
logical (default: 
engine 
character string specifying the implementation to use.
Up to polyCub version 0.43, the twodimensional nodes and weights
were computed by R functions and these are still available by setting

plot 
logical indicating if an illustrative plot of the numerical integration should be produced. 
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
.
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.
Sommariva, A. and Vianello, M. (2007): Product Gauss cubature over polygons based on Green's integration formula. BIT Numerical Mathematics, 47 (2), 441453. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/s1054300701312")}
Other polyCubmethods:
polyCub.exact.Gauss()
,
polyCub.iso()
,
polyCub.midpoint()
,
polyCub()
## a function to integrate (here: isotropic zeromean 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) # nonvectorized 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))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.