Exploiting symmetry makes a 2D NIntegrate harder?

Question

NIntegrate is having trouble converging, for reasons that I don't understand. Here is a simple example.

I can compute the area of the curvilinear triangle

RegionPlot[ArcTan[x]+ArcTan[y] <= Pi/6, {x,0,1}, {y,0,1}]

easily by saying

NIntegrate[Boole[ArcTan[x]+ArcTan[y] <= Pi/6], {x,0,1}, {y,0,1}]

and the answer comes back 0.150728 with no difficulty.

You would think that I would make things easier for Mathematica if I instead exploit the symmetry around the main diagonal by saying

2 NIntegrate[Boole[ArcTan[x]+ArcTan[y] <= Pi/6 \[And] x <= y], {x,0,1}, {y,0,1}]

But no! In this case, Mathematica first complains that

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

It then gives up on the integration, commenting that

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.07504438521372334 and 0.0000162228158507992 for the integral and error estimates.

and returning the fairly poor approximation 0.150089.

Why is the latter integration having trouble? The integrand is piecewise constant. Raising the working precision doesn't seem to help, nor does raising the parameter MaxErrorIncreases.

I am eager to understand the source of this convergence problem, since I am getting similar convergence errors on intriguing integrals for which no symmetry provides a workaround.

Answer 1

I dont know a fix, but at least its intersting to see whats going on: Start with the problematic case:

pts = Reap[
  2 NIntegrate[
    Boole[ArcTan[x] + ArcTan[y] <= Pi/6 && x <= y], {x, 0, 1}, {y, 0, 
     1}, EvaluationMonitor :> Sow[{x, y}]]]
Show[{ListPlot[pts[[2, 1]]], Plot[x, {x, 0, 1}], 
  ContourPlot[ArcTan[x] + ArcTan[y] == Pi/6 , {x, 0, 1}, {y, 0, 1}, 
   ContourStyle -> Red]}]

You see the sampling scheme treats your Boole as a black box and does a huge number of evaluations hunting for the boundary.

Now look at the simpler expression:

pts = Reap[
     NIntegrate[
     Boole[ArcTan[x] + ArcTan[y] <= Pi/6 ], {x, 0, 1}, {y, 0, 1}, 
     EvaluationMonitor :> Sow[{x, y}]]];
Show[{ListPlot[pts[[2, 1]]], Plot[x, {x, 0, 1}], 
  ContourPlot[ArcTan[x] + ArcTan[y] == Pi/6 , {x, 0, 1}, {y, 0, 1}, 
   ContourStyle -> Red]}]

Evidently NIntegrate analytically recognizes the simple Boole expression and so automatically only samples over the appropriate region. Since the value is 1 everywhere it returns quickly with very few samples.

The same issue arises for the full area case if you set Method -> {Automatic, "SymbolicProcessing" -> False}