Exact integral of a polynomial over a tetrahedron or a simplex, with Python.

integral_on_simplex.py

from math import factorial
from sympy import Poly
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    n = len(monom)
    return coef / np.prod(list(range(n+1, n+j+1)))

def integral(P, S):
    gens = P.gens
    n = len(gens)
    dico = {}
    v = S[n,:]
    columns = []
    for i in range(n):
        columns.append(S[i,:] - v)    
    B = np.column_stack(tuple(columns))
    for i in range(n):
        newvar = v[i]
        for j in range(n):
            newvar = newvar + B[i,j]*Poly(gens[j], gens, domain="RR")
        dico[gens[i]] = newvar.as_expr()
    Q = P.subs(dico, simultaneous=True).as_expr().as_poly(gens)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / factorial(n) * s

###############################################################################
# the simplex
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])
S = np.array([v1, v2, v3, v4])

# polynomial to integrate
from sympy.abc import x, y, z
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, S)
print(val)

integral_on_tetrahedron.py

from math import factorial
from sympy import Poly
from sympy.abc import x, y, z
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    return coef / np.prod(list(range(4, 4+j)))

def integral(P, v1, v2, v3, v4):
    X = Poly(x, x, y, z, domain="RR")
    Y = Poly(y, x, y, z, domain="RR")
    Z = Poly(z, x, y, z, domain="RR")
    B = np.column_stack((v1-v4, v2-v4, v3-v4))
    newx = B[0,0]*X + B[0,1]*Y + B[0,2]*Z + v4[0]
    newy = B[1,0]*X + B[1,1]*Y + B[1,2]*Z + v4[1]
    newz = B[2,0]*X + B[2,1]*Y + B[2,2]*Z + v4[2]
    let = {x: newx.as_expr(), y: newy.as_expr(), z: newz.as_expr()}
    Q = P.subs(let, simultaneous=True).as_expr().as_poly(x, y, z)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / 6 * s

###############################################################################
# tetrahedron vertices
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])

# polynomial to integrate    
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, v1, v2, v3, v4)
print(val)