# Computing quadrature points on triangle

According to http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF page 7, it is possible to compute coordinate of quadrature points by using the script

n = 2;
m = 2;
p1 = Tuples[{1, Subscript[x, i], Subscript[y, i]}, {n}];
p2 = DeleteDuplicates[p1, Times @@ #1 == Times @@ #2 &];
p3 = Times @@ # & /@ p2;
p4 = Integrate[#, {Subscript[x, i], 0, 1}, {Subscript[y, i], 0,
1 - Subscript[x, i]}] & /@ p3;
sys = And @@
Table[p4[[j]] == Sum[Subscript[w, i] p3[[j]], {i, 1, m}], {j, 1,
6}]
Solve[sys, {Subscript[w, 1], Subscript[w, 2], Subscript[x, 1],
Subscript[x, 2], Subscript[y, 1], Subscript[y, 2]}, Reals]


However, Mathematica gives no solution. There is something wrong with the script or am I missing the real meaning? I know that there should be a solution over two points even this is not acceptable.

Thanks for suggestions

• I think you are missing a factor for the RHS sums (1/2). Does your script work for 1D? Feb 4, 2016 at 14:25

The reference provided in the question recommends m = 3, so that is what I consider. (Edit: m = 2 has no real solution, as can be seen from NSolve[sys, var].) The equations with minor simplification then become

n = 2;
m = 3;
p1 = Tuples[{1, Subscript[x, i], Subscript[y, i]}, {n}];
p3 = DeleteDuplicates[Times @@@ p1];
p4 = Integrate[#, {Subscript[x, i], 0, 1}, {Subscript[y, i], 0,
1 - Subscript[x, i]}] & /@ p3;
sys = Table[p4[[j]] == Sum[Subscript[w, i] p3[[j]], {i, 1, m}]/2, {j, 1, 6}]
var = Flatten@Table[{{Subscript[w, i], Subscript[x, i], Subscript[y, i]}}, {i, m}]
Solve[sys, var]


Also, the right side of the equations in sys is divided by two, a minor correction. As expected, the problem is underdetermined and give a multiplicity of solutions. To obtain the specific one in the reference, replace the definition of sys by

sys = Table[p4[[j]] == Sum[Subscript[w, i] p3[[j]], {i, 1, m}]/2, {j, 1, 6}] /.
{Subscript[w, 1] -> 1/3, Subscript[w, 3] -> 1/3,
Subscript[x, 1] -> 1/6, Subscript[y, 1] -> 1/6}


which yields the desired

{{Subscript[w, 2] -> 1/3, Subscript[x, 2] -> 1/6, Subscript[y, 2] -> 2/3,
Subscript[x, 3] -> 2/3, Subscript[y, 3] -> 1/6},
{Subscript[w, 2] -> 1/3, Subscript[x, 2] -> 2/3, Subscript[y, 2] -> 1/6,
Subscript[x, 3] -> 1/6, Subscript[y, 3] -> 2/3}}

• Then I misunderstood the reference. "In theory,Ng= 2 may work since in this case we are going to have 6 unknowns... but the resulting quadrature is NOT symmetric". I explored this possibility thinking real solutions exist... but apparently the reference is not "well" written. Feb 4, 2016 at 15:56
• @Fabio A better answer (than the one I deleted) to your comment is that the lack of symmetry made me concerned that the roots would not be real, so I tried first Solve without Reals and, when that returned with no answer, NSolve without Reals, which gave complex solutions. Then, I went on to m = 3. I agree that the sentence you quoted is not very informative. Feb 4, 2016 at 16:35