Search code examples
pythonintegral

Integral on a polyhedral non-rectangular domain with Python

You are a Python user and you want to evaluate this triple integral:

enter image description here

for a certain function f, say f(x,y,z) = x + y*z.

In R, a possibility is to nest one-dimensional integrals:

f <- function(x,y,z) x + y*z

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) 
## -5358.3 with absolute error < 9.5e-11

The point I don't like with this method is that it does not take into account the error estimates on the first two integrals.

I'm interested in the Python way to realize this method. But I'm also interested in a Python way to get a result with a more reliable error estimate.

Solution

  • To give a solution to this problem, we will use four libraries:

    import numpy as np
import pypoman
from scipy.spatial import Delaunay
from pysimplicialcubature.simplicialcubature import integrateOnSimplex

    The domain of integration is defined by the set of inequalities:

    enter image description here

    which is equivalent to:

    enter image description here

    This set of inequalities defines a convex polyhedron. We can get the vertices of this polyhedron with the pypoman package (or we could use pycddlib but this would be more lengthy):

    A = np.array([
  [-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 = np.array([5, 4, 5, 3, 10, 6])

vertices = pypoman.compute_polytope_vertices(A, b)

    Now that we get the vertices of this convex polytope, we decompose it into tetrahedra with the help of the Delaunay algorithm:

    dlnay = Delaunay(vertices)
tetrahedra = np.asarray(vertices)[dlnay.simplices]

    A tetrahedron is a three-dimensional simplex. So it remains to apply the pysimplicialcubature library to evaluate the integral:

    # function to integrate
f = lambda x : x[0] + x[1] * x[2]
# integral of f on the tetrahedra
I_f = integrateOnSimplex(f, tetrahedra)
print(I_f)

    We get:

    {
 'integral': -5358.300000000005, 
 'estAbsError': 1.1879700000000003e-08, 
 'functionEvaluations': 330.0, 
 'returnCode': 0, 
 'message': 'OK'
}

    Actually our function f is a multivariate polynomial. In this case, with pysimplicialcubature, it is possible to get the exact value of the integral (function integratePolynomialOnSimplex).

    Edit

    This method is now implemented in the PyPolyhedralCubature package.