I'm comparing different quadrature formulas over triangular regions in 2D on MATLAB. My professor gave me the article from Xiao-Gimbutas (2010): A numerical algorithm for the construction of efficient quadrature rules in two and higher dimensions and in the Conclusion paragraph they give nodes (x,y) and weights w for some formulas in 2D, with different degrees of precision, over the triangle T with vertices $(-1,-1/\sqrt3), \, (0,2/\sqrt3), \, (1,-1/\sqrt3)$. The first formula they give is, as they say: "A reflectively symmetric quadrature of degree 10 with a total number of 24 nodes on T. Only the generating 14 nodes on either the generator triangle $T_2$
or the median are listed."
$\begin{bmatrix}
x & y & w \\
0.1551275866992061 & 0.3926587334507982 & 0.05567338455800167 \\
0.9252314250365582 & -0.5241973822167591 & 0.00896962036448401 \\
0.1384253518864596 & 0.8231152771873217 & 0.02240118142259882 \\
0.3609812845058779 & 0.0267765158868186 & 0.06089917117558221 \\
0.5812005281770837 & 0.0484050181187623 & 0.03108755935951332 \\
0.3458850891448130 & 0.4561437111102405 & 0.03146650526166210 \\
0.4510519110813880 & -0.3059662347573659 & 0.06246416621868210 \\
0.6893778940829133 & -0.5226600964300047 & 0.02582709676771077 \\
0.7634762995170715 & -0.3046970738749675 & 0.03178319217145733 \\
0.2563805151012611 & -0.5238003914950184 & 0.03640547863919981 \\
0.0000000000000000 & 0.7059408419246426 & 0.03498602106696763 \\
0.0000000000000000 & 1.0463216388186550 & 0.01202119786551837 \\
0.0000000000000000 & 0.0240306520799389 & 0.07686209708916766 \\
0.0000000000000000 & -0.3100722765497620 & 0.07278083120266070 \\
\end{bmatrix}$
With the generator triangle $T_2$ they mean, as they show in the article, the right half of the triangle T, wich would give back the other half after a riflection on the y axis.
As I understood from the following theorem (the demonstration is on the article), the weight of the node obtained by symmetry is the same as the weight of the node that mapped into it (correct me if I'm wrong):
Theorem 3.2: Suppose that $X_B=\{x_1, … , x_n\}$ and $W_B=\{w_1, … , w_n\}$ is an n-point quadrature on B, an S-generator of $\Omega$, such that
$$\int_B \phi_j = \sum_{x_i\in X_B}w_i\phi_j(x_i) $$
is accurate for all $\{\phi_1, … , \phi_m\}$ defined on $\Omega$. Then, the S-symmetric p-point quadrature satisfies the equation
$$\int_{\Omega} \phi_j = \sum_{x_i\in X}w_i\phi_j(x_i) $$
for all $\{\phi_1, … , \phi_m\}$ if the nodes $X=\{x_1, … , x_p\}$ are given by the formula
$$X=\bigcup_{x_i\in X_B}O_{x_i}^S $$
and all weights belonging to nodes in the same orbit are equal.
S-symmetric means symmetric in regard of the symmetry S of $\Omega$; S-generator of $\Omega$ means that B can generate $\Omega$ via action of S; $O_{x_i}^S$ means the orbit of $x_i$ under the action of S. For further explanation on this definitions see the following:
"As is well known, the symmetry group G of $\Omega \subset \mathbb{R}^d$ is the set of all orthogonal linear transformations of $\mathbb{R}^d$ that map $\Omega$ onto itself. A subset $S\subseteq G$ is called a symmetry subgroup of $\Omega$ if S is itself a group. We say that two points $x_1, \, x_2$ in $\Omega$ are S-symmetrically related if there exists an element $s \in S$ such that $x_2 = sx_1$. We call the set of all points that are S-symmetrically related with x the orbit of x, and denote it by $O^S_x$.
If for each point x in $X\subset \Omega$, all S-symmetrically related points are also in X, then we say that X is S-symmetric (or invariant with respect to S). It can be easily seen that any S-symmetric set X is a union of distinct orbits. Suppose that the set B contains one and only one point from each orbit contained in X. Then we say that B is an S-generator of X. Thus, all S-symmetric sets X can be expressed as $X=\bigcup_{x\in B}O_{x}^S $, where $O_{x_j}^S\bigcap O_{x_i}^S=\emptyset$, for all $x_i\in B, x_j \in B$ and $x_i \not= x_j$. Obviously, many possible B's may exist for an arbitrary X. For example, X may be partitioned geometrically based on the location of its elements within."
So in order to obtain the full set of 24 nodes and weights I simply copy-pasted the list of the first ten rows and changed the sign on the x column:
$\begin{bmatrix}
x & y & w \\
0.1551275866992061 & 0.3926587334507982 & 0.05567338455800167 \\
0.9252314250365582 & -0.5241973822167591 & 0.00896962036448401 \\
0.1384253518864596 & 0.8231152771873217 & 0.02240118142259882 \\
0.3609812845058779 & 0.0267765158868186 & 0.06089917117558221 \\
0.5812005281770837 & 0.0484050181187623 & 0.03108755935951332 \\
0.3458850891448130 & 0.4561437111102405 & 0.03146650526166210 \\
0.4510519110813880 & -0.3059662347573659 & 0.06246416621868210 \\
0.6893778940829133 & -0.5226600964300047 & 0.02582709676771077 \\
0.7634762995170715 & -0.3046970738749675 & 0.03178319217145733 \\
0.2563805151012611 & -0.5238003914950184 & 0.03640547863919981 \\
0.0000000000000000 & 0.7059408419246426 & 0.03498602106696763 \\
0.0000000000000000 & 1.0463216388186550 & 0.01202119786551837 \\
0.0000000000000000 & 0.0240306520799389 & 0.07686209708916766 \\
0.0000000000000000 & -0.3100722765497620 & 0.07278083120266070 \\
-0.1551275866992061 & 0.3926587334507982 & 0.05567338455800167 \\
-0.9252314250365582 & -0.5241973822167591 & 0.00896962036448401 \\
-0.1384253518864596 & 0.8231152771873217 & 0.02240118142259882 \\
-0.3609812845058779 & 0.0267765158868186 & 0.06089917117558221 \\
-0.5812005281770837 & 0.0484050181187623 & 0.03108755935951332 \\
-0.3458850891448130 & 0.4561437111102405 & 0.03146650526166210 \\
-0.4510519110813880 & -0.3059662347573659 & 0.06246416621868210 \\
-0.6893778940829133 & -0.5226600964300047 & 0.02582709676771077 \\
-0.7634762995170715 & -0.3046970738749675 & 0.03178319217145733 \\
-0.2563805151012611 & -0.5238003914950184 & 0.03640547863919981 \\
\end{bmatrix}$
Then I copied this entire set of numbers into a matrix A in MATLAB and wrote this simple program to test it:
syms x y
f=@(x,y) x*y;
I=0;
for i=1:24
I=I+A(i,3)*f(A(i,1),A(i,2));
end
disp(I)
The result, wich is $0$ on T, should be exact because the order of precision is said to be $10$ so it should always integrate exactly all polynomials on two variables up to the degree $10$ (and $xy$ is a polynomial of degree $2$) but it gives me back $-8.673617379884035e-19$. I tried with other functions and it always gives me wrong results. Then I thought that it should also integrate exactly the function $1$, but with the function $1$ the quadrature formula gives only the sum of the weights: the integral of $1$ on T is the area of T (wich is $\sqrt3$) but the sum of the weights gives me back $0.930604859102099$.
Sorry for the very long question but I wanted to write down all the materials that I used to get this problem. Hope someone can explain where I made a mistake (I feel it must be a dumb one since I don't have any clue where it could be).
Best Answer
I'm summarizing my comments (which seem to have resolved the problem) as an answer so that the question can be marked as answered.
The $0$ value was rounded to $2^{-60}$, which is about the expected rounding error in this situation, so there was no problem in that regard.
The non-zero results were all off by a factor of $\sqrt[4]{12}$, which seems to be a normalization factor required for applying these weights in the present setup.