# Exact integral of a polynomial on a simplex

**Saturn Elephant**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

The paper
Simple formula for integration of polynomials on a simplex
by Jean B. Lasserre provides a method to calculate the exact value of
the integral of a multivariate polynomial on a simplex (i.e. a
tetrahedron in dimension three). I implemented it in
**Julia**, **Python**, and **R**.

Integration on simplices is important, because any convex polyhedron can
be decomposed into simplices, thanks to the Delaunay tessellation.
Therefore one can integrate over convex polyhedra once one can integrate
over simplices (I wrote an
example
of doing so with **R**).

## Julia

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

### Julia example

We define the polynomial to be integrated as follows:

using TypedPolynomials @polyvar x y z P = x^4 + y + 2*x*y^2 - 3*z

*Be careful*. If the expression of your polynomial does not
involve one of the variables, e.g. \(P(x, y, z) = x^4 + 2xy^2\), you must define a polynomial involving this variable:

P = x^4 + 2*x*y^2 + 0.0*z

Now we define the simplex as a matrix whose rows correspond to the vertices:

# 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]

And finally we run the function:

integratePolynomialOnSimplex(P, S)

## Python

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 integratePolynomialOnSimplex(P, S): gens = P.gens n = len(gens) S = np.asarray(S) v = S[n,:] columns = [] for i in range(n): columns.append(S[i,:] - v) B = np.column_stack(tuple(columns)) dico = {} 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) print(Q) 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

### Python example

# 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] # polynomial to integrate from sympy import Poly from sympy.abc import x, y, z P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR") # integral integratePolynomialOnSimplex(P, S)

## R

library(spray) integratePolynomialonSimplex <- function(P, S) { n <- ncol(S) v <- S[n+1L, ] B <- t(S[1L:n, ]) - v gens <- lapply(1L:n, function(i) lone(i, n)) newvars <- vector("list", n) for(i in 1L:n) { newvar <- v[i] Bi <- B[i, ] for(j in 1L:n) { newvar <- newvar + Bi[j] * gens[[j]] } newvars[[i]] <- newvar } Q <- 0 exponents <- P[["index"]] coeffs <- P[["value"]] for(i in 1L:nrow(exponents)) { powers <- exponents[i, ] term <- 1 for(j in 1L:n) { term <- term * newvars[[j]]^powers[j] } Q <- Q + coeffs[i] * term } s <- 0 exponents <- Q[["index"]] coeffs <- Q[["value"]] for(i in 1L:nrow(exponents)) { coef <- coeffs[i] powers <- exponents[i, ] d <- sum(powers) if(d == 0L) { s <- s + coef next } coef <- coef * prod(factorial(powers)) s <- s + coef / prod((n+1L):(n+d)) } abs(det(B)) / factorial(n) * s }

### R example

library(spray) # variables x <- lone(1, 3) y <- lone(2, 3) z <- lone(3, 3) # polynomial P <- x^4 + y + 2*x*y^2 - 3*z # simplex (tetrahedron) vertices v1 <- c(1, 1, 1) v2 <- c(2, 2, 3) v3 <- c(3, 4, 5) v4 <- c(3, 2, 1) # simplex S <- rbind(v1, v2, v3, v4) # integral integratePolynomialonSimplex(P, S)

## Note

The functions do not check whether the given matrix
`S`

defines a non-degenerate simplex. This is equivalent to
the invertibility of the matrix `B`

constructed in the
functions.

**leave a comment**for the author, please follow the link and comment on their blog:

**Saturn Elephant**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.