2
$\begingroup$

All,

I am looking for Gauss quadrature formulas for a particular geometric setting.

That is, I am integrating functions over the standard simplex (triangle in dimension 2 and tetrahedron in dimension 3). I would like the usual requirements to be satisfied: symmetric quadrature points lying inside the simplex with positive weights. Additionally, I would like integrals of products of multi-linear functions (bi-linear in dimension 2 and tri-linear in dimension 3) to be computed exactly.

In other words, I would like the quadrature formula to give exact results for these integrals: \begin{alignat*}{3} \int_0^1 \int_0^{1-x} x^a y^b dx dx &\qquad a,b \in \{0,1,2\} &\qquad \mbox{(dim = 2)} \\ \int_0^1 \int_0^{1-x} \int_0^{1-x-y} x^a y^b z^c dz dy dx &\qquad a,b,c \in \{0,1,2\} &\qquad \mbox{(dim = 3)} \end{alignat*}

Does anyone know if such quadrature formulas exist? Is there a nice reference which shows how one derives quadrature formulas on simplices? (The equations satisfied by the quadrature points and weights are straightforward to write down, but not so straightforward to tackle without a bit of magic that I could perhaps learn from studying some known derivations!)

Of course, there are formulas in the literature for integration of degree-six polynomials on tetrahedra and degree-four polynomials on triangles which I could use. But since I do not need to reproduce all of these polynomials, rather only the polynomials which arise as products of multilinear polynomials, I am wondering if there's a simpler formula with fewer quadrature points that I can apply or derive.

Thanks for your help.

Adrian

$\endgroup$
3
$\begingroup$

It looks to me like you are searching for what are called monomial cubature rules on the simplex in the literature on numerical integration. As an example, a 4th-degree monomial in 2-D $(x,y)$ would be any of: $x^4, x^3y, x^2y^2, xy^3, y^4$. A 4th-degree monomial cubature rule will exactly integrate any of those, or any linear combination of them.

There are monomial cubature rules over a variety of domains, with various conditions on the points (on/off boundary, symmetric, etc.), with or without positive weights. In some domains there are Gaussian rules, though it's not obvious to me only Gaussian rules would meet your needs. (Interestingly, in the course of digging up references for this answer, I learned that Gauss-Lobatto rules, i.e., Gaussian rules including boundary points, do not exist on the triangle; see arXiv:1007.4848.)

(My apologies for not providing links to everything cited here; I don't yet have the reputation on MO to post more than two links.)

If you use MATLAB, search for "n-dimensional simplex quadrature" at Mathworks.com and you'll find a contributed MATLAB script for Gaussian cubature on an $n$-D simplex. Unfortunately the author does not provide a reference for the algorithm.

Three widely-cited compendiums of cubature rules (albeit somewhat dated) are these by Ronald Cools:

"Monomial cubature rules since “Stroud”: a compilation" (1993; DOI:10.1016/0377-0427(93)90027-9)

"Monomial cubature rules since “Stroud”: a compilation — part 2" (1999; DOI:10.1016/S0377-0427(99)00229-0)

An encyclopaedia of cubature formulas (2003)

These include Gaussian rules on various simplexes, such as rules from "High degree efficient symmetrical Gaussian quadrature rules for the triangle" (1985; DOI:10.1002/nme.1620210612).

Cools maintains an online Encyclopaedia of Cubature Formulas that compiles tables with the nodes and weights for many known rules, including links to the literature. It includes work since the two articles cited above, but I don't know how thorough its recent coverage is. To access the rules, you'll need a password from Cools, but he is generous with those. The site currently has dozens of rules for the triangle and tetrahedron, for degrees as high as 30, and with varying numbers of points. It has about two dozen rules for $n$-D simplexes up to degree 5, and a more general $n$-D "embedded" rule. Each rule has a "quality" label that identifies whether weights are positive and whether nodes are interior; the PI label indicates interior nodes with positive weights. There are many PI rules for the triangle and tetrahedron.

Many of the rules are from a famous earlier compendium of Stroud. If you do a Google search for "gaussian monomial rules on simplex" (no quotes), you'll find free software in C, C++, and various Fortran versions, implementing rules from Stroud's collection, including simplex rules.

In Python, the quadpy package (on GitHub) includes a selection of simplex rules, including Gaussian rules. These are in the nsimplex module. The source file has a link to an unpublished technical report by Noel Walkington detailing the derivation of those rules.

With collaborators, Cools has produced a package of cubature algorithms in Fortran 95 for $n$-D simplices and parellellepids; see "Algorithm 824: CUBPACK: a package for automatic cubature; framework description" (2003; 10.1145/838250.838253). A Google search for "CUBPACK" will take you to the package page. It's not open-source, although source code is available for research purposes only.

Finally, if you need high precision, you may be interested in work Cools did with fellow cubature expert Alan Genz on adaptive simplex cubature. Their algorithm cleverly subdivides a simplex into additional simplices to achieve a targeted precision. See "An adaptive numerical cubature algorithm for simplices" (2003; DOI:10.1145/838250.838254). This algorithm was my own introduction to cubature on simplices.

By the way, you may get more answers to this question if you move or duplicate it to Computational Science Stack Exchange.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for your most thorough and detailed answer to my question. The goldmine of information you have pointed me to will surely contain all the answers that I need. I sincerely hope that the above answer has a very positive impact on your stackexchange reputation. Your reputation in my eyes is certainly very high :) $\endgroup$ – user2698883 Oct 25 '17 at 4:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.