# NIntegrate with singularity does not work

Consider the following integral $$\int_{-\infty}^\infty \int_{-\infty}^\infty \frac{e^{-x^2-y^2}}{|x^2+y^2-1|^{0.6}} dxdy.$$ By changing to polar coordinate analytically and putting the result into Mathematica, the result is 6.589.

However, I tried to numerically integrate this directly:

int = Exp[-x^2 - y^2]/Abs[x^2 + y^2 - 1]^0.6
NIntegrate[int, {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]

Then, it throws the following error:

Next, I tried to handle the singularity as follows:

NIntegrate[int, {x, -Infinity, Infinity}, {y, -Infinity, Infinity},
Exclusions -> x^2 + y^2 == 1]

However, again I obtain an error with inaccurate result 2.56.

Curiously, changing the exponent 0.6 to 0.5 gives an accurate result for both with and without Exclusions option.

Why this happen, and why putting Exclusion option is ineffective? How can I evaluate this integral accurately?

The main contribution to the integral comes from x^2+y^2==1 as you can see:

int = Exp[-x^2 - y^2]/Abs[x^2 + y^2 - 1]^0.6
Plot3D[int, {x, -2, 2}, {y, -2, 2}, PlotPoints -> 50]

Therefore you need to sample mostly around this region. This can be achieved by the options Method-> "LocalAdaptive":

NIntegrate[int, {x, -Infinity, Infinity}, {y, -Infinity, Infinity},
(*6.56931*)
• Thanks for your answer! Then can I understand that if an integral is sharply peaked (either by singularity or not), then LocalAdaptive method is better? I am asking this question since the document reference.wolfram.com/language/tutorial/… says that generally GlobalAdaptive is better than LocalAdaptive, and furthermore if an integrand has singularity then Exclusions options is good. Oct 10, 2021 at 11:05
• Yes. You can read about the strategies used in: the help under: tutorial/NIntegrateIntegrationStrategies Oct 10, 2021 at 11:11
• Well.. it is that document that I found out that "Globaladaptive is generally better than localadaptive", and also Exclusions option is helpful for dealing with singularities. Oct 10, 2021 at 11:24
• More info than the document gives, I can only guess. "Global" seems to imply hat no region is specially treated, whereas "Local" seem to imply that some regions are treated differently from others. But it would be nice if you could ask Wolfram at [email protected] about it and report it here. Oct 10, 2021 at 12:30
• @eigenvalue “Generally better” means “usually better”, not “always better” in this context. Often local adaptive is much worse, but here it seems to be much better. I’ve never really clearly understood how to predict when it will be better. Oct 10, 2021 at 12:39

This is a complicated improper integral: the integration domain is the whole plane $$\mathbb{R}^2$$ and there is a singular line $$x^2+y^2=1$$. I obtained different results by

int = Exp[-x^2 - y^2]/RealAbs[x^2 + y^2 - 1]^(3/5)]
NIntegrate[int, {x,-Infinity,Infinity}, {y,-Infinity,Infinity},Exclusions -> x^2 + y^2 == 1]

2.56359

and by switching to the polar coordinates

2*Pi*NIntegrate[r*Exp[-r^2]/RealAbs[r^2 - 1]^(3/5), {r, 0, Infinity}]

6.58887

Its exact value is unclear to me in view of

2*Pi*Integrate[r*Exp[-r^2]/RealAbs[r^2 - 1]^(3/5), {r, 0, Infinity}]

(\[Pi] (14 Gamma[2/5] - 25 (-1)^(3/5) Gamma[12/5] + 14 (-1)^(3/5) Gamma[2/5, -1]))/(14 E)

N[%]

6.58887 - 5.86567*10^-16 I

and

2*Pi*Integrate[ r*Exp[-r^2]/RealAbs[r^2 - 1]^(3/5), {r, 0, 1, Infinity}]

(\[Pi] Gamma[2/5])/E

N[%]

2.56359

I don't know good numerical methods for multiple improper integrals.

• It is very mysterious. By the way, I obtained the 6.58887 value by first changing to your form 2*Pi*Integrate[r*Exp[-r^2]/RealAbs[r^2 - 1]^(3/5), {r, 0, Infinity}], and then further analytically performing the substitution $r^2=u$. In this case, regardless of setting the integral range as {u,0,Infinity} or {u,0,1,Infinity}, the value is 6.58887. Oct 10, 2021 at 11:38
• @eigenvalu: Which version do you use? In 12.3.1 on Windows 10 N[Pi*Integrate[Exp[-x]/RealAbs[x - 1]^(3/5), {x, 0, 1, Infinity}]] still results in 2.56359 and N[Pi*Integrate[Exp[-x]/RealAbs[x - 1]^(3/5), {x, 0, Infinity}]] still produces 6.58887 - 5.13246*10^-16 I. Oct 10, 2021 at 14:12
• I use version 12.1 on Window 10. Oct 11, 2021 at 1:15
• @eigenvalue: It's strange: the documentation to Integrate says the command was updated in 12.0 latest time. Oct 11, 2021 at 5:33

First try analytical integration. (For symmetry 4 times the first quadrant.) You get some complicated unsolved integrals of x, you can solve with NIntegrate.

integrand = Exp[-x^2 - y^2]/Abs[x^2 + y^2 - 1]^(6/10);

aint = 4*Integrate[integrand, {x, 0, Infinity}, {y, 0, Infinity}]

(*   4*(Integrate[-((5*Sqrt[Pi]*Gamma[2/5]*
Hypergeometric1F1[1/2, 9/10, -1 + x^2])/
(E^x^2*((1 - x^2)^(1/10)*Gamma[-(1/10)]))),
{x, 0, 1}] +
Integrate[(5*((Gamma[2/5]*Gamma[11/10]*Hypergeometric1F1[
1/2, 9/10, -1 + x^2])/(Sqrt[Pi]*
(1 - x^2)^(1/10)) - Gamma[9/10]*
Hypergeometric1F1[3/5, 11/10, -1 + x^2]))/
E^x^2, {x, 0, 1}] +
Integrate[((1/2)*Sqrt[Pi]*HypergeometricU[3/5, 11/10,
-1 + x^2])/E^x^2, {x, 1, Infinity}])   *)

aint /. Integrate -> NIntegrate

(*   6.58887   *)