PolyhedralCubature.jl documentation

This package allows to evaluate multiple integrals like

\[\int_{-5}^4\int_{-5}^{3-x}\int_{-10}^{6-x-y} f(x, y, z)\,\text{d}z\,\text{d}y\,\text{d}x.\]

In other words, the domain of integration is given by a set of linear inequalities:

\[\left\{\begin{matrix} -5 & \leqslant & x & \leqslant & 4 \\ -5 & \leqslant & y & \leqslant & 3-x \\ -10 & \leqslant & z & \leqslant & 6-x-y \end{matrix}\right..\]

These linear inequalities define a convex polytope (in dimension 3, a polyhedron). In order to use the package, one has to get the matrix-vector representation of these inequalities, of the form

\[A {(x,y,z)}' \leqslant b.\]

The matrix $A$ and the vector $b$ appear when one rewrites the linear inequalities above as:

\[\left\{\begin{matrix} -x & \leqslant & 5 \\ x & \leqslant & 4 \\ -y & \leqslant & 5 \\ x+y & \leqslant & 3 \\ -z & \leqslant & 10 \\ x+y+z & \leqslant & 6 \end{matrix}\right..\]

The matrix $A$ is given by the coefficients of $x$, $y$, $z$ at the left-hand sides, and the vector $b$ is made of the upper bounds at the right-hand sides:

A = [
  -1  0  0;   # -x
   1  0  0;   # x
   0 -1  0;   # -y
   1  1  0;   # x + y
   0  0 -1;   # -z
   1  1  1    # x + y + z
]
b = [5; 4; 5; 3; 10; 6]

Integrating an arbitrary function

Now assume for example that $f(x,y,z) = x(x+1) - yz^2$. Once we have $A$ and $b$, here is how to evaluate the integral of $f$ over the convex polytope:

using PolyhedralCubature
function f(x, y, z)
  return x*(x+1) - y*z^2
end
function g(v)
  return f(v[1], v[2], v[3])
end
I_f = integrateOnPolytope(g, A, b)
I_f.integral
# 57892.27499999998

As we shall see, the exact value of this integral is $57892.275$.

Integrating a multivariate polynomial

The integratePolynomialOnPolytope function provided by the package implements a procedure to get the exact value of the integral in the case when the integrand is polynomial, as in our current example.

using TypedPolynomials
@polyvar x y z
poly = x + y*z
integratePolynomialOnPolytope(poly, A, b)
# 2315691//40 (= 57892.275)

Method

The method runs as follows.

• Firstly, the Polyhedra package is used to get the vertices of the convex

polytope from the matrix $A$ and the vector $b$.

• Then, with the help of the MiniQhull package, the Delaunay tessellation

is applied to the vertices and one gets a partition of the convex polytope into simplices.

• Finally, the SimplicialCubature package is used to evaluate the integrals

on the simplicies.

When the element type of A and b is Int64, then the vertex coordinates obtained at the first step will have the Rational{BigInt} type. If possible, always use the element type Int64 or Rational{Int64} or Rational{BigInt} for A and b. This allows to get an exact rational result from integratePolynomialOnPolytope when the polynomial has integer or rational coefficients.

___

Member functions

integrateOnPolytope(f, A, b; dim, maxEvals, absError, tol, rule, info, fkwargs...)

integratePolynomialOnPolytope(poly, A, b)