Multiple integration on convex polytopes.
Project description
pypolyhedralcubature
Multiple integration over a convex polytope.
This package allows to evaluate a multiple integral whose integration bounds are some linear combinations of the variables, e.g.
$$\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 & \leq & x & \leq & 4 \\ -5 & \leq & y & \leq & 3-x \\ -10 & \leq & z & \leq & 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 & \leq & 5 \\ x & \leq & 4 \\ -y & \leq & 5 \\ x+y & \leq & 3 \\ -z & \leq & 10 \\ x+y+z & \leq & 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:
import numpy as np
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])
The function getAb
provided by the package allows to get $A$ and $b$ in a
user-friendly way:
from pypolyhedralcubature.polyhedralcubature import getAb
from sympy.abc import x, y, z
# linear inequalities defining the integral bounds
i1 = (x >= -5) & (x <= 4)
i2 = (y >= -5) & (y <= 3 - x)
i3 = (z >= -10) & (z <= 6 - x - y)
# get the matrix-vector representation of these inequalities
A, b = getAb([i1, i2, i3], [x, y, z])
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:
from pypolyhedralcubature.polyhedralcubature import integrateOnPolytope
# function to integrate
f = lambda x, y, z : x*(x+1) - y*z**2
# integral of f over the polytope defined by the linear inequalities
g = lambda v : f(v[0], v[1], v[2])
I_f = integrateOnPolytope(g, A, b)
I_f["integral"]
# 57892.275000000016
In the case when the function to be integrated is, as in our current example,
a polynomial function, it is better to use the integratePolynomialOnPolytope
function provided by the package. This function implements a procedure
calculating the exact value of the integral. Here is how to use it:
from pypolyhedralcubature.polyhedralcubature import integratePolynomialOnPolytope
from sympy import Poly
from sympy.abc import x, y, z
# polynomial to integrate
P = Poly(x*(x+1) - y*z**2, domain = "RR")
# integral of P over the polytope
integratePolynomialOnPolytope(P, A, b)
# 57892.2750000001
Actually the exact value of the integral is $57892.275$, so there is a slight numerical error in the procedure. We can get this exact value by using the field of rational numbers as the domain of the polynomial:
# polynomial to integrate
P = Poly(x*(x+1) - y*z**2, domain = "QQ")
# integral of P over the polytope
integratePolynomialOnPolytope(P, A, b)
# 2315691/40
Acknowledgments
I am grateful to the StackOverflow user @Davide_sd for the help he provided
regarding the getAb
function.
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distribution
Built Distribution
Hashes for pypolyhedralcubature-0.2.0.tar.gz
Algorithm | Hash digest | |
---|---|---|
SHA256 | 7f8672601c0272a38c8a1ddc77cf0c78540f22f6cc86dd71e31cc89be73eda07 |
|
MD5 | 38c35d411d5232637c78cea366a64503 |
|
BLAKE2b-256 | 3a2e7ddaa4b1c3fa2143169aa3b2cb3f647883f8e1b08e848af20e3aceafe94d |
Hashes for pypolyhedralcubature-0.2.0-py3-none-any.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | d3c44a062b355beadd68255bff02e684c9e5d9f8a0a5a6dd1c360d039f1c5490 |
|
MD5 | abba5d72f041c67e9b03ac1028b63977 |
|
BLAKE2b-256 | 3b9eb07c0421d16859d423a78425c109e7462b6289e771aeb6e2082965b8ebf2 |