A version of Simpson's rule for multiple integrals

Let M(f) denote the midpoint rule and T(f) the trapezoidal rule for estimating . Then Simpson's rule =[lambda]M(f)+(1-[lambda])T(f), where . We generalize Simpson's rule to multiple integrals as follows. Let Dn be some polygonal region in let P0,...,Pm denote the vertices of Dn, and let Pm+1 equal the center of mass of Dn. Define the linear functionals M(f)=Vol(Dn)f(Pn+1), which generalizes the midpoint rule, and T(f)=Vol(Dn)([1/(m+1)][summation operator]j=0mf(Pj)), which generalizes the trapezoidal rule. Finally, our generalization of Simpson's rule is given by the cubature rule (CR) L[lambda]=[lambda]M(f)+(1-[lambda])T(f), for fixed [lambda], 0[less-than-or-equals, slant][lambda][less-than-or-equals, slant]1. We choose [lambda], depending on Dn, so that L[lambda] is exact for polynomials of as large a degree as possible. In particular, we derive CRs for the n simplex and unit n cube. We also use points Qj[set membership, variant][not partial differential](Dn), other than the vertices Pj, to generate T(f). This sometimes leads to better CRs for certain regions -- in particular, for quadrilaterals in the plane. We use Grobner bases to solve the system of equations which yield the coordinates of the Qj's.