My question is a follow up to this one.
The setup is the same, we are looking for a third order accurate (i.e. one that works for polynomials of degree 2) quadrature method on the same triangle.
However, on that answer the points are prescribed. In my formulation you must find the ideal points to sample as well.
I found this other post that claims that for an order 3 quadrature method you need 7 points, the 3 vertices, the 3 midpoints and the center of mass obeying:
$$ \int_T f \approx \frac{|T|}{60} \left(27f(a^T) + 3 \sum_{j=1}^3 f(a_j) + 8\sum_{j=1}^3 f(a^j) \right).$$
I am going to assume $|T|$ is the area of the triangle. But I cannot find the book and it doesn't seem obvious to me. I especially don;t understand why we would need 7 points instead of 3 or 6.
I actually immediately found the book with a Google search of the authors and title. I’m not including links because I don’t know whether that might create copyright issues for the site. I found both the first and second edition in PDF; the page number $415$ refers to the second edition. That book is mainly concerned with quadrature on triangulations, e.g. for applying the finite element method, and in that context one prefers quadrature points on the boundary of the elements to ensure continuity of the function across elements. Since you don’t seem to have any preference for points on the boundary, that book and the cited formula aren’t all that relevant.
How many points we need to achieve a certain order of quadrature on a triangle can be worked out using representation theory. Since quadrature is linear, we can use whatever triangle is most convenient. I’ll use an equilateral triangle for the symmetry considerations and then an isosceles right triangle for the actual calculation of the quadrature points.
The symmetry group of an equilateral triangle is $D_3 \cong S_3$, the dihedral group of $6$ elements: The identity, three reflections and two rotations through $\pm\frac{2\pi}3$. Its character table is:
$$ \begin{array}{c|cc} &(1)(2)(3)&(12)(3)&(123)\\\hline \chi_1&1&1&1\\ \chi_2&1&-1&1\\ \chi_3&2&0&-1 \end{array} $$
The coordinate functions $x$ and $y$ afford an irreducible two-dimensional representation with character $\chi_3$, as you can check by taking the traces of the transformation matrices: With suitable coordinate axes, the diagonal elements are $(1,1)$ for the identity, $(1,-1)$ for one of the reflections and $\left(\cos\frac{2\pi}3,\cos\frac{2\pi}3\right)=\left(-\frac12,-\frac12\right)$ for the rotations.
The monomials $x^n,x^{n-1}y,\ldots,xy^{n-1},y^n$ form a basis for the homogeneous polynomials of degree $n$ and afford an $(n+1)$-dimensional representation of $D_3$. This is the symmetrized $n$-fold tensor product of the above two-dimensional coordinate function representation; its character is obtained by substituting $\chi_3(g^k)$ for $a_k$ in the cycle index of $S_n$. For instance, for $n=3$, this is
$$ \frac16\left(1\chi_3(g^3)+3\chi_3(g)\chi_3(g^2)+2\chi_3(g^3)\right)\;. $$
Computing the character and projecting it onto the irreducible characters tells us how many copies of each of the irreducible representations are contained in each of these symmetrized tensor products. The results up to $n=7$ are:
\begin{array}{c|ccc|cc} n&(1)(2)(3)&(12)(3)&(123)&\chi_1&\chi_2&\chi_3\\\hline 0&1&1&1&1&0&0\\ 1&2&0&-1&0&0&1\\ 2&3&1&0&1&0&1\\ 3&4&0&1&1&1&1\\ 4&5&1&-1&1&0&2\\ 5&6&0&0&1&1&2\\ 6&7&1&1&2&1&2\\ 7&8&0&-1&1&1&3\\ \end{array}
If we choose quadrature points that respect the symmetries of the triangle, functions that transform under the non-identity representations $\chi_2$ and $\chi_3$ will automatically be integrated to the correct value $0$. This is analogous to symmetric one-dimensional quadrature schemes automatically correctly integrating odd functions to $0$.
To derive the order of quadrature schemes with a given number of points, let’s first look at how the order is determined in the one-dimensional case, e.g. for Gaussian quadrature. For an odd number $m=2k+1$ of points, there are $2k+1$ free parameters (given the symmetry): $k+1$ weights and $k$ positions (since one position must be $0$). For an even number $m=2k$ of points, there are $2k$ free parameters: $k$ weights and $k$ positions. Thus, in both cases, there are $m$ free parameters, so $m$ linearly independent even functions can be correctly integrated. This leads to the order $2m$ of one-dimensional Gaussian quadrature with $m$ points.
In our case, we can also either include the centroid as a quadrature point or not. If we do, there are $3k+1$ points and $2k+1$ free parameters: $k+1$ weights and $k$ positions. If we don’t, there are $3k$ points and $2k$ parameters: $k$ weights and $k$ positions. (This is for groups of $3$ points on the medians. If we use groups of $6$ points in general symmetric position instead, they only have three free parameters (two coordinates and one weight) per $6$ points, whereas two groups of $3$ points on the medians have four free parameters (two coordinates and two weights).) Thus, as in the one-dimensonal case, for any given number of linearly independent $\chi_1$ components (corresponding to even functions in the one-dimensional case), we can derive a quadrature scheme that integrates them exactly; we just need $3k$ or $3k+1$ points where the one-dimensional scheme needs $2k$ or $2k+1$, respectively.
The simplest scheme is of course the one that uses only the centroid (with the triangle’s area as the weight). It integrates one $\chi_1$ component exactly, and the table says that there’s one such component from constant functions, none from linear functions and one from quadratic functions, so this scheme correctly integrates constant and linear functions and is thus of order $2$.
The next stage uses $3$ points and can correctly integrate two $\chi_1$ components. According to the table, that makes it of order $3$. As I said, I’ll perform the calculation for an isosceles right triangle because that’s more convenient. So consider the triangle with points $(0,0)$, $(1,0)$ and $(0,1)$ (and thus centroid $\left(\frac13,\frac13\right)$). Three symmetrically arranged points must be on the medians in the equilateral triangle, and the linear transformation to the isosceles right triangle preserves the medians, so they must be on the medians here, too. If they’re at a fraction $\lambda$ of the distance from the corner to the centroid, in this triangle that puts them at $\left(\frac\lambda3,\frac\lambda3\right)$, $\left(\frac\lambda3,1-\frac{2\lambda}3\right)$ and $\left(1-\frac{2\lambda}3,\frac\lambda3\right)$. Note that the admissible range for $\lambda$ is $\left[0,\frac32\right]$, since the points can be on either side of the centroid.
The weight must be $\frac16$ (one third of the area of the triangle) to get the constant functions right. That leaves only the parameter $\lambda$ to be determined. Since all quadratic monomials have a component along the quadratic $\chi_1$ component, it doesn’t matter which one we use to determine $\lambda$; let’s use $x^2$. Then the equation for $\lambda$ is
$$ \frac16\left(\left(\frac\lambda3\right)^2+\left(\frac\lambda3\right)^2+\left(1-\frac{2\lambda}3\right)^2\right)=\int_0^1\mathrm dy\int_0^{1-y}\mathrm dxx^2\;, $$
or
$$ \frac1{18}\left(2\lambda^2-4\lambda+3\right)=\frac1{12}\;. $$
There are two solutions, $\lambda=\frac32$ and $\lambda=\frac12$. The first one reproduces the solution mentioned in Philip Roe’s answer to the first question you linked to, where the quadrature points are the midpoints of the sides. In the other solution, with $\lambda=\frac12$, the quadrature points are halfway between the corners and the centroid. In the present triangle, they’re at $\left(\frac16,\frac16\right)$$\left(\frac16,\frac23\right)$ and $\left(\frac23,\frac16\right)$. We can check that $xy$ is also correctly integrated:
$$ \frac16\left(\frac16\cdot\frac16+\frac16\cdot\frac23+\frac23\cdot\frac16\right)=\int_0^1\mathrm dy\int_0^{1-y}\mathrm dxxy=\frac1{24}\;. $$
As regards the $7$-point formula you quote, the question you linked to misinterpreted the order of that formula. The book says that its “degree of exactness” is $3$, which means that its order is $4$. But as we saw, if you don’t limit yourself to points on the perimeter, you can achieve order $4$ with fewer points. For order $4$, according to the table, we need to be able to integrate one more $\chi_1$ component exactly, so we need one more point at the centroid. Then we have two weights, $w_\text{c}$ for the centroid and $w_\lambda$ for the other three points, and again the positional parameter $\lambda$. Again, we can determine the parameters using any of the cubic monomials since they all have a component along the cubic $\chi_1$ component; let’s use $x^3$. The resulting equations are:
\begin{eqnarray} w_\text{c}\left(\frac13\right)^0+w_\lambda\left(\left(\frac\lambda3\right)^0+\left(\frac\lambda3\right)^0+\left(1-\frac{2\lambda}3\right)^0\right)&=&\int_0^1\mathrm dy\int_0^{1-y}\mathrm dxx^0\;,\\ w_\text{c}\left(\frac13\right)^2+w_\lambda\left(\left(\frac\lambda3\right)^2+\left(\frac\lambda3\right)^2+\left(1-\frac{2\lambda}3\right)^2\right)&=&\int_0^1\mathrm dy\int_0^{1-y}\mathrm dxx^2\;,\\ w_\text{c}\left(\frac13\right)^3+w_\lambda\left(\left(\frac\lambda3\right)^3+\left(\frac\lambda3\right)^3+\left(1-\frac{2\lambda}3\right)^3\right)&=&\int_0^1\mathrm dy\int_0^{1-y}\mathrm dxx^3\;, \end{eqnarray}
or
\begin{eqnarray} w_\text{c}+3w_\lambda&=&\frac12\;,\\ \frac19\left(w_\text{c}+3w_\lambda\left(2\lambda^2-4\lambda+3\right)\right)&=&\frac1{12}\;,\\ \frac1{27}\left(w_\text{c}+3w_\lambda\left(-2\lambda^3+12\lambda^2-18\lambda+9\right)\right)&=&\frac1{20}\;. \end{eqnarray}
Solving the first equation for $w_\text{c}$ and substituting into the other two yields
\begin{eqnarray} w_\lambda(1-\lambda)^2&=&\frac1{24}\;,\\ w_\lambda(1-\lambda)^2(4-\lambda)&=&\frac{17}{120} \end{eqnarray}
and thus $\lambda=\frac35$.
Substituting this into the equations and solving for $w_\lambda$ and $w_\text{c}$ yields $w_\lambda=\frac{25}{96}$ and $w_\text{c}=-\frac{27}{96}$. Thus, the weight at the centroid is negative, which is a somewhat undesirable property of a quadrature scheme. The quadrature points for the present triangle are the centroid $\left(\frac13,\frac13\right)$ and the three points $\left(\frac15,\frac15\right)$, $\left(\frac15,\frac35\right)$ and $\left(\frac35,\frac15\right)$. Again, you can check that $x^2y$ is also integrated correctly:
$$ -\frac{27}{96}\left(\frac13\right)^3+\frac{25}{96}\left(\left(\frac\lambda3\right)^3+\left(\frac\lambda3\right)^3+\left(1-\frac{2\lambda}3\right)^3\right)=\int_0^1\mathrm dy\int_0^{1-y}\mathrm dxx^2y=\frac1{60}\;. $$
To avoid the negative weight, we’d have to go on to $6$ points symmetrically arranged. Those could yield a scheme of order $5$, or we could use the extra free parameter to try to get a scheme of order $4$ with positive weights. The calculations become more tedious, though, and I doubt the results can be expressed in closed form, so I’ll leave it at that. The weights and points for quadrature schemes of this type for higher orders are tabulated in Appendix II of High degree efficient symmetrical Gaussian quadrature rules for the triangle by D.A. Dunavant. (To translate the weights and the normalized barycentric coordinates specified there to the present isosceles right triangle, take the $\alpha$ and $\beta$ values as the $x$ and $y$ coordinates, respectively, and multiply the weights by the area $\frac12$ of the triangle.