A version of Simpson's rule for multiple integrals (Q5946722)

The Simpson rule (SR) for the integral \(\int_a^b f(x) dx\) can be written as \(L_\lambda(f)=\lambda M(f)+(1-\lambda)T(f)\) where \(\lambda=2/3\), \(M(f)\) is the midpoint rule and \(T(f)\) is the trapezoidal rule. It is the integral of the interpolating polynomial of degree 2. For integrals in \({\mathbb R}^n\), over a polygonal region \(D_n\) with vertices \(P_0,\ldots,P_m\), and \(P_{m+1}\) the center of mass, these rules are generalized as \(M(f)=\text{vol}(D_n)f(P_{m+1})\), \(T(f)=\text{vol}(D_n)/(m+1)\sum_{j=0}^m f(P_j)\), and the Simpson rule can still be defined as \(L_\lambda(f)\) for some \(0\leq\lambda\leq 1\). A particular property of the SR in 1D is that it is exact for polynomials of degree at most 2. In the \(n\)-dimensional case this is not obvious in general and it will depend on \(D_n\) whether this can be reached or not. The choice of \(\lambda\) is discussed which gives a maximal degree of exactness, i.e., the degree of the polynomials for which the SR is exact is maximized by choosing \(\lambda\). Also cubature rules are constructed that use points that are in the faces, but not at the vertices. Gröbner bases are used to find these points. If possible, a connection with polynomial interpolation is given. In this way, formulas are given for the \(n\)-simplex, the \(n\)-cube, and plane regions like the unit disk and a trapezoid. Some of the formulas are known, others are new.