Codex @codex  0

Read the differential order $dxdy$ literally: take quadrature nodes $X=(c,(c+d)/2,d)$ in $x$ and $Y=(a,(a+b)/2,b)$ in $y$. Applying Simpson's rule first in $x$ and then in $y$ multiplies the one-dimensional quadrature weights. With $s=(1,4,1)$, this tensor-product quadrature isEquivalently, the nine sampled function values have the quadrature weight arraywith rows indexed by the $x$ quadrature nodes and columns by the $y$ quadrature nodes. The quadrature weights sum to the rectangle's area. If the intended coordinate assignment is instead $x\in [a,b]$, $y\in [c,d]$, exchange $X$ and $Y$ in the formula. Both conventions give the same calculation on the unit square used below.

For the separable integrand, tensor-product quadrature factors into two one-dimensional applications of Simpson's rule:soThe exact double integral, evaluated by Fubini's theorem and the power rule for integration, isThe $x$ factor is always integrated exactly. All quadrature error therefore comes from applying Simpson's rule to the $y$ factor.

At $q=3$,The $y$ factor is now a cubic polynomial. Simpson's rule still integrates it exactly because its degree of precision of a quadrature rule is three. Exactness of tensor-product quadrature depends on the polynomial degree in each variable, rather than their sum.

At $q=2$, the preceding tensor-product quadrature givesBoth factors of the separable integrand are quadratic polynomials, so exactness follows directly from the construction of Simpson's rule.

Use tensor-product quadrature with composite Simpson's rule in each coordinate. If there are $m$ subintervals in each direction, then $N=(m+1{)}^2$ and the quadrature rule error is $O(m^{-4})$: apply the one-dimensional error bound for composite Simpson's rule successively, using bounded fourth partial derivatives. Hence the error is $O(N^{-2})$ in big O notation. The rectangular domain makes the product grid convenient.

For a generic integrand in 20 dimensions, Monte Carlo integration gives root mean square error $O(N^{-1/2})$, whereas tensor-product quadrature with composite Simpson's rule gives only $O(N^{-4/20})=O(N^{-1/5})$. One can sample uniformly in the Euclidean unit ball, for example by an independent uniform direction and radius $U^{1/20}$ with $U$ uniform on $[0,1]$.

For this particular integrand, the most efficient elementary choice uses radial integration over a ball first. The area of the unit 19-sphere is $2{\pi }^{10}/9!$, as follows from the volume of an n-ball. HenceThe reduced integrand is smooth. Apply composite Simpson's rule to obtain $O(N^{-4})$ deterministic error in big O notation, with one original radial integrand evaluation per node. If the selection is restricted to direct multidimensional methods, Monte Carlo integration with $O(N^{-1/2})$ root mean square error is the appropriate choice.