5
$\begingroup$

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:

enter image description here

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. enter image description here

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?

Improve this question
$\endgroup$
0

3 Answers 3

Reset to default
8
$\begingroup$

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]

enter image description here

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}, Method -> "LocalAdaptive"] (*6.56931*)
Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Oct 10, 2021 at 11:05
  • 1
    $\begingroup$ Yes. You can read about the strategies used in: the help under: tutorial/NIntegrateIntegrationStrategies $\endgroup$ Commented Oct 10, 2021 at 11:11
  • $\begingroup$ 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. $\endgroup$ Commented Oct 10, 2021 at 11:24
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Oct 10, 2021 at 12:30
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Oct 10, 2021 at 12:39
4
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$ Commented Oct 10, 2021 at 11:38
  • $\begingroup$ @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. $\endgroup$ Commented Oct 10, 2021 at 14:12
  • $\begingroup$ I use version 12.1 on Window 10. $\endgroup$ Commented Oct 11, 2021 at 1:15
  • $\begingroup$ @eigenvalue: It's strange: the documentation to Integrate says the command was updated in 12.0 latest time. $\endgroup$ Commented Oct 11, 2021 at 5:33
4
$\begingroup$

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 *)
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.