# Handing singularity of numerical integral to high precision

Consider the integral

    int[z] = Integrate[(2 x^2 - 1) (\[Pi] Sqrt[1 - x^2])^-1
(2 y^2 -  1) (\[Pi] Sqrt[1 - y^2])^-1 Log[z Abs[x - y]], {x, -1, 1}, {y, -1, 1}]


This integral can be done analytically, and gives the same answer

int[z] = -1/4

for any z. Let us now try to compute it numerically

  Nint[z] :=   NIntegrate[(2 x^2 - 1) (\[Pi] Sqrt[1 - x^2])^-1 (2 y^2 -
1) (\[Pi] Sqrt[1 - y^2])^-1 Log[z Abs[x - y]], {x, -1, 1}, {y,-1,1}]


for z=1 gives the low precision answer

Nint = -0.250002


while for z>1 it fails to evaluate. Increasing precision does not seem to drastically improve the situation. I suspect the problem is the singularities when

x,y=1,-1


but I have not been able to successfully implement this yet. Any ideas?

• Using "DuffyCoordinates" might be one thing to try to weaken the singularities. – J. M.'s ennui Oct 16 '18 at 12:45

Brute force, but only takes a second:

Block[{z = 1},
NIntegrate[(2 x^2 - 1) (π Sqrt[1 - x^2])^-1 (2 y^2 - 1) *
(π Sqrt[1 - y^2])^-1 Log[z Abs[x - y]],
{x, -1, 1}, {y, -1, (*x,*) 1}, (* specify singularity y == x: unnecesary *)
MaxRecursion -> 20, Method -> "GaussKronrodRule" (* uses open sampling *)
]
]


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.

(*  -0.25  *)


NIntegrate::slwcon is just a warning, not an indication of error. It should be expected in this case.

Response to comment:

Maybe it's worth outlining my thought process in the present case. A general discussion of method selection is beyond the scope of an SE answer. If you get to know the NIntegrate tutorials, as @J.M. suggests, you'll get to know the tools at hand and the examples in which they are effective. Then you're more than halfway to understanding how to approach numerical integration in Mathematica. "There is no royal road to geometry," and however you learn about numerical integration, it takes considerable time.

So first of all, I wasn't sure what would work. I didn't think the inverse square roots and the logarithmic singularity looked too bad, but multivariate integration is still a bit of a mystery to me (in terms of what sort of integrand causes trouble. The first thing that stuck out was the diagonal singularity along y == x. The recursive subdivision divides the integration domain into rectangles which would be unlikely to be successful, since the singularity would cut across the rectangles. So I thought I should explicitly identify the singularity by including x in {y, 0, x, 1} (this turns out to be unnecessary, as I just now discovered). Testing just this alteration leads to overflow, so the change does have an effect.

I know in single-variable numerical integration, such singularities are not that bad and can be handled by a good open-sampling rule such as Gauss-Kronrod. (The function is analytic between the singularities, and previous experience with inverse square root and logarithmic singularities has shown me it might work well enought.) So I decided to try it. I didn't know it would work until it did. Since it did, I stopped thinking about the problem until....

Until the OP commented and I wondered if the automatic symbolic processing might identify the singularity without my help. Removing x from {y, 0, x, 1} shows that manually identifying the singular locus y == x was unnecessary.

• btw, how do you know in each situation whether to use super high precision, multidimensionalrule, GaussKronrodRule, etc? is it just guess and check, or is there some intuition i can learn? – esches Oct 16 '18 at 14:13
• @esches, the advanced documentation for NIntegrate[] is a more-or-less complete guide on how to choose a Method for a particular problem, depending on its special features. – J. M.'s ennui Oct 16 '18 at 14:42