## Improve accuracy of NIntegrate with GlobalAdaptive over ImplicitRegion?

Let’s say that I want to integrate some arbitrarily “nice” function (uniformly $C^{\infty}$ -smooth, for example) over an ImplicitRegion in more than three dimensions. For example, let’s consider the integral $\int_{D\times D} 1\,\mathrm{d}\mu$ , where $D$ is the unit disk in two dimensions. We can approximate this in Mathematica as follows:

disk2 = ImplicitRegion[x1^2 + y1^2 <= 1 && x2^2 + y2^2 <= 1, {x1, y1, x2, y2}]
val[n_] := NIntegrate[1, Element[{x1, y1, x2, y2}, disk2], WorkingPrecision -> n]


The exact value of this is $\pi^2$ , but the numerical value returned by Mathematica 11.1 doesn’t converge to $\pi^2$ for increasing $n$ . Namely, val[50]-Pi^2 and val[100]-Pi^2 both give me the same error, which is roughly $1.38\times 10^{-18}$ . With slightly more complicated functions (that are still $C^1$ , at least), the error we get can be much worse; consider the following integral:

NIntegrate[Sqrt[(x1 - x2)^2 + (y1 - y2)^2], Element[{x1, y1, x2, y2}, disk2], WorkingPrecision -> n]


We should expect this integral to have a value of $\frac{128\pi}{45}$ , but the output differs from that value by a little more than $4\times 10^{-5}$ , regardless of the setting for WorkingPrecision. Additionally, it doesn’t seem that DiscretizeRegion can really be applied to regions that can’t be expressed as the product of $n$ -dimensional regions for $n\leq 3$ , which in the case of a more general region rules out the strategy here. Are there any general strategies I can use to get with arbitrary precision and accuracy the values of numerical integrals of $C^1$ functions on ImplicitRegions in more than three dimensions? Thanks for your time.

Recent Questions – Mathematica Stack Exchange