# Computing Quadratures of Polynomial Functions over k-Simplices $$\triangle^k$$ in $$\mathbb{R}^n$$¶

Contents:

Suppose we are given a k-Simplex $$\triangle^k$$ in $$\mathbb{R}^n$$

$\triangle^k=\left\{ \sum_{i=0}^{k}\lambda_{i}\mathbf{x}_{i}\vert\lambda_{i}\ge0,\,0\leq i\le k,\,\sum_{i=0}^{k}\lambda_{i}=1\right\}$

where $$\mathbf{x}_i$$ represents the ith vertex of the simplex.

$\begin{split}\mathbf{x}_{1}=\left(\begin{array}{c} x_{1}\\ y_{1}\\ z_{1}\\ \vdots \end{array}\right)\qquad\mathbf{x}_{2}=\left(\begin{array}{c} x_{2}\\ y_{2}\\ z_{2}\\ \vdots \end{array}\right)\dots\qquad\mathbf{x}_{i}=\left(\begin{array}{c} x_{i}\\ y_{i}\\ z_{i}\\ \vdots \end{array}\right)\dots\qquad\mathbf{x}_{k}=\left(\begin{array}{c} x_{k}\\ y_{k}\\ z_{k}\\ \vdots \end{array}\right)\end{split}$

We are interested in finding an explicit expression for the quadrature of a homogeneous polynomial over such a simplex:

$\int_{\triangle_{k}}x^{\alpha}y^{\beta}\dots$

and the moments, which are the above expression divided by the volume of the simplex:

$\langle x^\alpha y^\beta \dots \rangle = \frac{1}{\text{vol} (\triangle_k)} \int_{\triangle_{k}}x^{\alpha}y^{\beta}\dots$

We have implemented in python a simple algorithm that computes exactly the above expressions:

Compute the binomial coefficients $$\binom nk = \frac{n!}{k!(n-k)!}$$

>>> nCk(0,0)
1
>>> nCk(1,0)
1
>>> nCk(4, 2)
6
>>> [[nCk(n,k) for k in range(n+1)] for n in range(7)] == [
...          [1],
...        [1, 1],
...       [1, 2, 1],
...      [1, 3, 3, 1],
...     [1, 4, 6, 4, 1],
...    [1, 5, 10, 10, 5, 1],
...  [1, 6, 15, 20, 15, 6, 1]                              ]
True


Compute the quadrature of a homogenous polynomial given by term=’xxy’ over the vertices verts=(0,1,2) of a simplex.

>>> quadrature('', (0,1))
([((), 1)], Fraction(1, 1))
([((('x', 0),), 1), ((('x', 1),), 1)], Fraction(1, 2))
([((('x', 0),), 1), ((('x', 1),), 1), ((('x', 2),), 1)], Fraction(1, 3))
...    [((('x', 0), ('x', 0)), 1),
...     ((('x', 0), ('x', 1)), 1),
...     ((('x', 1), ('x', 1)), 1)], fractions.Fraction(1, 3))
True
...    [((('x', 0), ('y', 0)), 2),
...      ((('x', 0), ('y', 1)), 1),
...      ((('x', 0), ('y', 2)), 1),
...      ((('x', 1), ('y', 0)), 1),
...      ((('x', 1), ('y', 1)), 2),
...      ((('x', 1), ('y', 2)), 1),
...      ((('x', 2), ('y', 0)), 1),
...      ((('x', 2), ('y', 1)), 1),
...      ((('x', 2), ('y', 2)), 2)], fractions.Fraction(1, 12))
True


# Description of Algorithm¶

Our algorithm is general, but for simplicity we will just demonstrate it for the specific case of a triangle in two dimensions. Assume you have a triangle in the plane defined by its three vertices at $$(x_0,y_0)$$, $$(x_1,y_1)$$, and $$(x_2,y_2)$$ - it is a 2-simplex so $$k=2$$. We will compute the moment of the q-homoqenous polynomial $$\langle xy \rangle$$ with $$q=2$$.

The inputs to the algorithm are two tuples - the first is $$\text{term}=(x,y)$$ describing the polynomial and the other is $$\text{verts}=(0,1,2)$$ describing the labels for the vertices. With this notation, the orders of the simplex $$k$$ and the degree of the polynomial $$q$$ are given by

$\begin{split}k &= \text{len}(\text{verts})-1 = 2 \\ q &= \text{len}(\text{term}) = 2\\\end{split}$

Now take all the permutations of the $$\text{term}$$ tuple.

$P = \mathcal{P}(\text{term}) = (xy, yx)$

Compute all the combinations (with repetitions allowed) of size $$q=2$$ of the coordinates of the vertices

$V = \mathcal{C}^q(\text{verts}) = (00, 01, 02, 11, 12, 22)$

Take the product of $$V$$ and $$P$$

$P \times V = \left ( (xy,00), (xy,01),(xy,02),(xy,11),(xy,12),(xy,22), (yx,00),(yx,01),(yx,02),(yx,11),(yx,12),(yx,22) \right)$

For each element in the tuple zip together the variable with the simplex vertices (e.g. $$(xy,00) \mapsto x_0y_0, \quad (xy,01) \mapsto x_0y_1$$ and so on):

$\text{zip}(P \times V) = \left ( x_0 y_0, x_0 y_1, x_0 y_2, x_1 y_1, x_1y_2, x_2y_2, y_0x_0, y_0x_1,y_0x_2,y_1x_1,y_1x_2,y_2x_2 \right)$

Treat each element as a multiplication and sum all the multiples together

$\sum(\text{zip}(P \times V)) = \left ( 2 x_0 y_0 + x_0 y_1 + x_0 y_2 + x_1 y_0 + 2 x_1 y_1 + x_1 y_2 + x_2 y_0 + x_2 y_1 + 2 x_2 y_2 \right)$

Divide the result by $$\binom{k+q}{q} q!$$ to obtain the final expression

$\frac{1}{\binom{k+q}{q} q!} \sum(\text{zip}(P \times V)) = \frac{1}{12} \left ( 2 x_0 y_0 + x_0 y_1 + x_0 y_2 + x_1 y_0 + 2 x_1 y_1 + x_1 y_2 + x_2 y_0 + x_2 y_1 + 2 x_2 y_2 \right)$

which is indeed what one would obtain for $$\langle xy \rangle$$ by manual calculation.

# Pseudocode¶

In pseudocode the algorithm can be concisely given as:

$\begin{split}& \text{def quadrature}(\text{term},\,\text{verts}): \\ & \qquad k = \text{len}(\text{verts})-1 \\ & \qquad q = \text{len}(\text{term}) \\ & \qquad P = \mathcal{P}(\text{term}) \\ & \qquad V = \mathcal{C}^q(\text{verts}) \\ & \qquad \text{return} \quad \frac{1}{\binom{k+q}{q} q!} \sum(\text{zip}(P \times V))\end{split}$

# Validation¶

To validate the algorithm we consider a triangle, a line segment, and a tetrahedron, and we compute a number of low order moments by integrating explicitly. The Mathematica file where all the integration is done is provided for download here.

## Triangle $$\triangle^2$$ in $$\mathbb{R}^2$$¶

The moments of all the homogeneous polynomials over a triangle area given by

$\langle x^{n}y^{m}\rangle=\frac{1}{A}\int_{\triangle}x^{n}y^{m}\,\mathbf{d}x\wedge\mathbf{d}y$

where $$\mathbf{d}x\wedge\mathbf{d}y$$ is the area form for $$\mathbb{R}^2$$ and $$A=\int_{\triangle}\mathbf{d}x\wedge\mathbf{d}y$$ is the total area of the triangle.

For example, $$\langle x \rangle$$ and $$\langle y \rangle$$ are the coordinates of the center of mass of the triangle.

Change of variables $$x,y \mapsto s,t$$

$\begin{split}x &= (1-s-t)x_{0}+sx_{1}+tx_{2}=s(x_{1}-x_{0})+t(x_{2}-x_{0}) \\ y &= (1-s-t)y_{0}+sy_{1}+ty_{2}=s(y_{1}-y_{0})+t(y_{2}-y_{0})\end{split}$

so that $$(s,t)=(0,0)$$ corresponds to vertex $$(x_{0},y_{0})$$, $$(s,t)=(1,0)$$ corresponds to vertex $$(x_{1},y_{1})$$, and $$(s,t)=(0,1)$$ corresponds to vertex $$(x_{2},y_{2})$$. The area form transforms as follows:

$\begin{split}dx\wedge dy &= \left((x_{1}-x_{0})(y_{2}-y_{0})-(x_{2}-x_{0})(y_{1}-y_{0})\right)\, ds\wedge dt \\ &= 2A\, ds\wedge dt\end{split}$

Thus we have an explicit expression for the quadrature

$\langle x^{n}y^{m}\rangle=2\int_{0}^{1}ds\int_{0}^{1-s}dt\left(s(x_{1}-x_{0})+t(x_{2}-x_{0})\right)^{n}\left(s(y_{1}-y_{0})+t(y_{2}-y_{0})\right)^{m}dt$

which can be evaluated by hand to obtain the first few moments:

$\begin{split}\frac{1}{6} \left(x_{0}^{2} + x_{0} x_{1} + x_{0} x_{2} + x_{1}^{2} + x_{1} x_{2} + x_{2}^{2}\right) \\ \frac{1}{10} \left(x_{0}^{3} + x_{0}^{2} x_{1} + x_{0}^{2} x_{2} + x_{0} x_{1}^{2} + x_{0} x_{1} x_{2} + x_{0} x_{2}^{2} + x_{1}^{3} + x_{1}^{2} x_{2} + x_{1} x_{2}^{2} + x_{2}^{3}\right) \\ \frac{1}{12} \left(2 x_{0} y_{0} + x_{0} y_{1} + x_{0} y_{2} + x_{1} y_{0} + 2 x_{1} y_{1} + x_{1} y_{2} + x_{2} y_{0} + x_{2} y_{1} + 2 x_{2} y_{2}\right) \\ \frac{1}{21} \left(x_{0}^{5} + x_{0}^{4} x_{1} + x_{0}^{4} x_{2} + x_{0}^{3} x_{1}^{2} + x_{0}^{3} x_{1} x_{2} + x_{0}^{3} x_{2}^{2} + x_{0}^{2} x_{1}^{3} + x_{0}^{2} x_{1}^{2} x_{2} + x_{0}^{2} x_{1} x_{2}^{2} + x_{0}^{2} x_{2}^{3} + x_{0} x_{1}^{4} + x_{0} x_{1}^{3} x_{2} + x_{0} x_{1}^{2} x_{2}^{2} + x_{0} x_{1} x_{2}^{3} + x_{0} x_{2}^{4} + x_{1}^{5} + x_{1}^{4} x_{2} + x_{1}^{3} x_{2}^{2} + x_{1}^{2} x_{2}^{3} + x_{1} x_{2}^{4} + x_{2}^{5}\right) \\\end{split}$
$\begin{split}\langle 1 \rangle &= 1 \\ \langle x \rangle &= \frac{1}{3}( x_0 + x_1 + x_2 ) \\ \langle x^2 \rangle &= \frac{1}{6}( x_0^2 + x_0 x_1 + x_1^2 + x_0 x_2 + x_1 x_2 + x_2^2 ) \\ \langle x^3 \rangle &= \frac{1}{10}( x_0^3 + x_0^2 x_1 + x_0 x_1^2 + x_1^3 + x_0^2 x_2 + x_0 x_1 x_2 + x_1^2 x_2 + x_0 x_2^2 + x_1 x_2^2 + x_2^3 ) \\ \langle x y \rangle &= \frac{1}{12}( 2 x_0 y_0 + x_1 y_0 + x_2 y_0 + x_0 y_1 + 2 x_1 y_1 + x_2 y_1 + x_0 y_2 + x_1 y_2 + 2 x_2 y_2 ) \\\end{split}$

Those match exactly with the results from our algorithm.

>>> def q(t): return tostr(quadrature(t, (0,1,2)))
>>> q('x')
'1/3 (x0 + x1 + x2)'
>>> q('xx')
'1/6 (x0x0 + x0x1 + x0x2 + x1x1 + x1x2 + x2x2)'
>>> q('xxx')
'1/10 (x0x0x0 + x0x0x1 + x0x0x2 + x0x1x1 + x0x1x2 + x0x2x2 + x1x1x1 + x1x1x2 + x1x2x2 + x2x2x2)'
>>> q('xy')
'1/12 (2 x0y0 + x0y1 + x0y2 + x1y0 + 2 x1y1 + x1y2 + x2y0 + x2y1 + 2 x2y2)'


## Line Segment $$\triangle^1$$ in $$\mathbb{R}^2$$¶

The same be done for line segments defined by their endpoints $$(x_0,y_0)$$ and $$(x_1,y_1)$$. Again the first few moments in this case are:

$\begin{split}\langle 1 \rangle &= 1 \\ \langle x \rangle &= \frac{1}{2}(x_0 + x_1) \\ \langle x^2 \rangle &= \frac{1}{3}(x_0^2 + x_0 x1 + x_1^2) \\ \langle x^3 \rangle &= \frac{1}{4}(x_0^3 + x_0^2 x_1 + x_0 x_1^2 + x_1^3) \\ \langle x y \rangle &= \frac{1}{6}(2 x_0 y_0 + x_1 y_0 + x_0 y_1 + 2 x_1 y_1) \\ \langle x^2 y \rangle &= \frac{1}{12}(3 x_0^2 y_0 + 2 x_0 x_1 y_0 + x_1^2 y_0 + x_0^2 y_1 + 2 x_0 x_1 y_1 + 3 x_1^2 y_1) \\\end{split}$

And our algorithm gives

>>> def q(t): return tostr(quadrature(t, (0,1)))
>>> q('x')
'1/2 (x0 + x1)'
>>> q('xx')
'1/3 (x0x0 + x0x1 + x1x1)'
>>> q('xxx')
'1/4 (x0x0x0 + x0x0x1 + x0x1x1 + x1x1x1)'
>>> q('xy')
'1/6 (2 x0y0 + x0y1 + x1y0 + 2 x1y1)'
>>> q('xxy')
'1/12 (3 x0x0y0 + x0x0y1 + 2 x0x1y0 + 2 x0x1y1 + x1x1y0 + 3 x1x1y1)'


## Tetrahedron $$\triangle^3$$ in $$\mathbb{R}^3$$¶

$\begin{split}\langle 1 \rangle &= 1 \\ \langle x \rangle &= \frac{1}{4}(x_0 + x_1 + x_2 + x_3) \\ \langle x^2 \rangle &= \frac{1}{20}(2 x_0^2 + 2 x_0 x_1 + 2 x_1^2 + 2 x_0 x_2 + 2 x_1 x_2 + 2 x_2^2 + 2 x_0 x_3 + 2 x_1 x_3 + 2 x_2 x_3 + 2 x_3^2) \\ \langle x^3 \rangle &= \frac{1}{20}(x_0^3 + x_0^2 x_1 + x_0 x_1^2 + x_1^3 + x_0^2 x_2 + x_0 x_1 x_2 + x_1^2 x_2 + x_0 x_2^2 + x_1 x_2^2 + x_2^3 + x_0^2 x_3 + x_0 x_1 x_3 + x_1^2 x_3 + x_0 x_2 x_3 + \\ & + x_1 x_2 x_3 + x_2^2 x_3 + x_0 x_3^2 + x_1 x_3^2 + x_2 x_3^2 + x_3^3) \\ \langle xy \rangle &= \frac{1}{20}(2 x_0 y_0 + x_1 y_0 + x_2 y_0 + x_3 y_0 + x_0 y_1 + 2 x_1 y_1 + x_2 y_1 + x_3 y_1 + x_0 y_2 + x_1 y_2 + 2 x_2 y_2 + x_3 y_2 + x_0 y_3 + x_1 y_3 + x_2 y_3 + 2 x_3 y_3) \\\end{split}$
>>> def q(t): return tostr(quadrature(t, (0,1,2,3)))
>>> q('x')
'1/4 (x0 + x1 + x2 + x3)'
>>> q('xx')
'1/10 (x0x0 + x0x1 + x0x2 + x0x3 + x1x1 + x1x2 + x1x3 + x2x2 + x2x3 + x3x3)'
>>> q('xxx')
'1/20 (x0x0x0 + x0x0x1 + x0x0x2 + x0x0x3 + x0x1x1 + x0x1x2 + x0x1x3 + x0x2x2 + x0x2x3 + x0x3x3 + x1x1x1 + x1x1x2 + x1x1x3 + x1x2x2 + x1x2x3 + x1x3x3 + x2x2x2 + x2x2x3 + x2x3x3 + x3x3x3)'
>>> q('xy')
'1/20 (2 x0y0 + x0y1 + x0y2 + x0y3 + x1y0 + 2 x1y1 + x1y2 + x1y3 + x2y0 + x2y1 + 2 x2y2 + x2y3 + x3y0 + x3y1 + x3y2 + 2 x3y3)'