Multiple integration over convex polyhedra

This package allows to evaluate multiple integrals like: \[ \int_{-5}^4\int_{-5}^{3-x}\int_{-10}^{6-x-y} f(x,y,z) \,\mathrm{d}z\,\mathrm{d}y\,\mathrm{d}x \]

Using base R only, a possibility is to nest the integrate function to evaluate such an integral:

f <- function(x, y, z) x*(x+1) - y*z^2
integrate(Vectorize(function(x) { 
  integrate(Vectorize(function(y) { 
    integrate(function(z) { 
      f(x,y,z) 
    }, -10, 6 - x - y)$value
   }), -5, 3 - x)$value 
}), -5, 4) 
#> 57892.27 with absolute error < 6.4e-10

This approach works well in general. But it has one default: the estimate of the absolute error it returns is not reliable, because the estimates of the absolute errors of the inner integrals are not taken into account.

Here is how to proceed with the polyhedralCubature package. The domain of integration is defined by the set of inequalities: \[ \left\{\begin{matrix} -5 & \leq & x & \leq & 4 \\ -5 & \leq & y & \leq & 3-x \\ -10 & \leq & z & \leq & 6-x-y \end{matrix} \right. \] which is equivalent to \[ \left\{\begin{matrix} -x & \leq & 5 \\ x & \leq & 4 \\ -y & \leq & 5 \\ x+y & \leq & 3 \\ -z & \leq & 10 \\ x+y+z & \leq & 6 \end{matrix} \right.. \] This set of inequalities defines a convex polyhedron. In order to use polyhedralCubature, you have to construct the matrix A defining the linear combinations of the variables, and the vector b giving the upper bounds of these linear combinations:

A <- rbind(
  c(-1, 0, 0), # -x
  c( 1, 0, 0), # x
  c( 0,-1, 0), # -y
  c( 1, 1, 0), # x+y
  c( 0, 0,-1), # -z
  c( 1, 1, 1)  # x+y+z
)
b <- c(5, 4, 5, 3, 10, 6)

library(polyhedralCubature)
f <- function(x, y, z) {
  x*(x+1) - y*z^2
}
integrateOverPolyhedron(f, A, b)
#> $integral
#> [1] 57892.28
#> 
#> $estAbsError
#> [1] 8.485622e-08
#> 
#> $functionEvaluations
#> [1] 330
#> 
#> $returnCode
#> [1] 0
#> 
#> $message
#> [1] "OK"

Alternatively, you can use the getAb function to get A and b with the help of the user-friendly syntax of the ompr package:

library(ompr)
model <- MIPModel() %>%
  add_variable(x) %>% add_variable(y) %>% add_variable(z) %>%
  add_constraint(-5 <= x) %>% add_constraint(x <= 4) %>%
  add_constraint(-5 <= y) %>% add_constraint(y <= 3 - x) %>%
  add_constraint(-10 <= z) %>% add_constraint(z <= 6 - x - y)
getAb(model)
#> $A
#>      [,1] [,2] [,3]
#> [1,]   -1    0    0
#> [2,]    1    0    0
#> [3,]    0   -1    0
#> [4,]    1    1    0
#> [5,]    0    0   -1
#> [6,]    1    1    1
#> 
#> $b
#> [1]  5  4  5  3 10  6

Observe that the function \(f\) is a polynomial here. Then one can get the exact value of the integral by feeding integrateOverPolyhedron with a spray polynomial instead of a function:

library(spray)
x <- lone(1, 3); y <- lone(2, 3); z <- lone(3, 3)
p <- f(x, y, z)
integrateOverPolyhedron(p, A, b)
#> $integral
#> [1] 57892.28
#> 
#> $functionEvaluations
#> [1] 24

The first step in integrateOverPolyhedron consists in finding the vertices of the polyhedron from the set of linear inequalities. It is possible, and better, to use exact arithmetic for this step, by defining the matrix A and the vector b in character mode, each entry representing an integer or a fraction like "1/3". Here we can simply do:

storage.mode(A) <- "character"
storage.mode(b) <- "character"

as.character(1/3)
#> [1] "0.333333333333333"

Instead, you must directly define A and b in character mode and enter "1/3" for this entry.

Finally, when \(f\) is a polynomial, you can get the exact value of the integral given as a fraction, by using a qspray polynomial instead of a spray polynomial:

library(qspray)
x <- qlone(1); y <- qlone(2); z <- qlone(3)
q <- f(x, y, z)
integrateOverPolyhedron(q, A, b)
#> [1] "2315691/40"