Exact integral of a polynomial over a simplex

General Usage
,
1

Hello,

I wrote a function which computes the exact value of a multivariate polynomial over a simplex. I’m open to any comment on my code.

I also implemented this method in Python. I’m sure it is correct, I multi-checked.

Integration on simplices is important because any convex polyhedron can be partitionned into simplices, thanks to the Delaunay tessellation. Hence integration on simplices allows integration on any convex polyhedron.

By the way, does there exist a Julia package for integration of an arbitrary function over a simplex?

using TypedPolynomials
using LinearAlgebra

function integratePolynomialOnSimplex(P, S)
    gens = variables(P)
    n = length(gens)
    v = S[end]    
    B = Array{Float64}(undef, n, 0)
    for i in 1:n
        B = hcat(B, S[i] - v)
    end
    Q = P(gens => v + B * vec(gens))
    s = 0.0
    for t in terms(Q)
        coef = TypedPolynomials.coefficient(t)
        powers = TypedPolynomials.exponents(t)
        j = sum(powers)
        if j == 0
            s = s + coef
            continue
        end
        coef = coef * prod(factorial.(powers))
        s = s + coef / prod((n+1):(n+j))
    end
    return abs(LinearAlgebra.det(B)) / factorial(n) * s
end


#### --- EXAMPLE --- ####

# define the polynomial to be integrated
@polyvar x y z
P = x^4 + y + 2*x*y^2 - 3*z

#= be careful
if your polynomial does not involve one of the variables, 
e.g. P(x,y,z) = x⁴ + 2xy², it must be defined as a 
polynomial in x, y, and z; you can do:
@polyvar x y z
P = x^4 + 2*x*y^2 + 0.0*z
then check that variables(P) returns (x, y, z)
 =#

# simplex vertices
v1 = [1.0, 1.0, 1.0] 
v2 = [2.0, 2.0, 3.0] 
v3 = [3.0, 4.0, 5.0] 
v4 = [3.0, 2.0, 1.0]
# simplex
S = [v1, v2, v3, v4]

# integral
integratePolynomialOnSimplex(P, S)
2

Here is an example of the calculation of an integral over a polyhedron. This is the Julia equivalent of this R example.

using Polyhedra
using MiniQhull # to use delaunay()
using TypedPolynomials

# define the polyhedron Ax <= b and get its vertices
# see https://laustep.github.io/stlahblog/posts/integralPolyhedron.html
A = [
    -1  0  0;
     1  0  0;
     0 -1  0;
     1  1  0;
     0  0 -1;
     1  1  1.0
]
b = [5; 4; 5; 3; 10; 6.0]
H = hrep(A, b)
P = polyhedron(H)
pts = points(P)
vertices = collect(pts)

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

# polynomial to be integrated
@polyvar x y z
f = x + y*z
# integral is sum of integrals over the tetrahedra
integral = 0
for i in 1:ntetrahedra
    I = integratePolynomialOnSimplex(f, tetrahedra[i])
    global integral = integral + I
end
integral
3

By the way, does there exist a Julia package for integration of an arbitrary function over a simplex?

I found GrundmannMoeller.jl but it seems to me this package has a problem.

4

Well, you could always do it via a coordinate transformation to a hypercube, at which point you could call HCubature.jl or similar (which has the advantage of being an adaptive scheme rather than fixed-order, though there are published works on adaptive quadrature for simplices).

See also Integral of a function inside non-rectangular domain - #9 by stevengj

5

I wrote some code for the integration of an arbitrary function over a simplex. I will package it one day. This is a port of the R package SimplicialCubature.

6

Blog post.

2 Likes