# Integration over a polyhedron with Julia

The Julia function `integrateOnSimplex`

provided in this gist allows to compute the integral of a function over a simplex or more generally a union of simplices (a simplex is a triangle in dimension 2, a tetrahedron in dimension 3). I will package it soon, under the name **SimplicialCubature**.

The Delaunay tessellation of a convex polyhedron in \(\mathbb{R}^n\) provides a partition of this convex polyhedron into \(n\)-dimensional simplices. Therefore, once you have a way to integrate a function over a simplex, you are able to integrate over a convex polyhedron.

Suppose for example you have to evaluate the integral \[ \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 \] for a certain function \(f\).

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. We can get the vertices of this polyhedron with the `Polyhedra`

package:

```
# define the polyhedron Ax <= b and get its vertices
using Polyhedra
= [
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.0 # x + y + z
]= [5; 4; 5; 3; 10; 6.0]
b = hrep(A, b)
H = polyhedron(H)
PH = collect(points(PH))
vertices ## 8-element Vector{Vector{Float64}}
## [4.0, -1.0, 3.0]
## [-5.0, 8.0, 3.0]
## [4.0, -1.0, -10.0]
## [-5.0, 8.0, -10.0]
## [4.0, -5.0, 7.0]
## [-5.0, -5.0, 16.0]
## [4.0, -5.0, -10.0]
## [-5.0, -5.0, -10.0]
```

Then, to evaluate the integral, we proceed as follows:

split the polyhedron into simplices (tetrahedra) with the Delaunay algorithm, implemented in the

`MiniQhull`

package;evaluate the integral over the union of these simplices with the

`integrateOnSimplex`

function.

Let’s go. We split the polyhedron into tetrahedra:

```
# decompose the polyhedron into simplices (tetrahedra)
using MiniQhull # to use delaunay()
= delaunay(hcat(vertices...))
indices , ntetrahedra = size(indices)
_= Vector{Vector{Vector{Float64}}}(undef, ntetrahedra)
tetrahedra for j in 1:ntetrahedra
= vec(indices[:, j])
ids = vertices[ids]
tetrahedra[j] end
```

Now we are ready to use the `integrateOnSimplex`

function. We take as example the function \(f(x, y, z) = x + yz\).

```
# function to be integrated
function f(x)
return x[1] + x[2]*x[3]
end
# integral
= integrateOnSimplex(f, tetrahedra)
I
I.integral## -5358.300000000004
```

Note that for this example we can do better. The function \(f\) is polynomial. So we can get the exact value of its integral over a simplex; see the previous post.