Multidimensional NIntegrate problem of the function decaying as 1/x^2

Question

The function I am trying to integrate is more complicated but I can simply write the function as (I had made a typo error, sorry. The '+' sign in front of the r should be '-'): $f(\omega ) = \int \limits_0^{\infty} r dr \int \limits_{-\infty}^{\infty} dz \frac{\epsilon}{\epsilon^2 + (\omega-z^2 - r^2)^2}$

So function decays as $1/x^2$. Let me choose $\epsilon = 0.1$, so function is close to zero everywhere. Mathematica cannot integrate efficiently even this simple function (it does, for the simple case after the correction, sorry everyone :( ). For my real case, it gave very huge numbers, which does not make sense. I have suspected the loss of precision, but still gave me huge numbers. (But this typo correction does not effect my real calculation, it is still giving huge result) I have tried:

• Set very high precision, no luck
• Tried increasing MaxRecursion and MinRecursion (as function is close to 0 mostly), no change
• Tried Trapezoidal method

• Symbolic Processing is off, of course

  NIntegrate[(1/10)/((1/10 - x^2 + r^2)^2 + 1/100) r, {x, -\[Infinity], \[Infinity]}, {r, 0, \[Infinity]}]
  

Any idea how to make this integration more efficient? Or a help understanding why this huge number problem appears?

Answer 1

Here is a plot of your function $\frac{\epsilon r}{\epsilon^2 + (\omega-z^2 + r^2)^2}$ (code for ComplexPlotR2 at end of answer):

ComplexPlotR2[
 CCompileR2[(1/10)/((1/10 - x^2 + y^2)^2 + 1/100) y], {-10, 10, 
  0.02}, {0, 10, 0.02}]

As you can see, it is nonzero on a pair of lines that extend to infinity, so it is not unexpected that the integral might be divergent.

To prove this, note that solving $$\epsilon-z^2+r^2=\pm\epsilon/2$$ for its positive branch gives $$z=\sqrt{r^2+\epsilon/2},\sqrt{r^2+3\epsilon/2}$$

giving a rough half-width in the $z$-axis of $$w(r)=\sqrt{r^2+3\epsilon/2}-\sqrt{r^2+\epsilon/2}$$

which means that when integrating along $z$, we get a contribution of roughly

$$w(r)\frac{\epsilon r}{\epsilon^2+\epsilon^2/4}\rightarrow\infty$$

as $r\rightarrow\infty$. So the integral is so strongly divergent that even its partial integrals along $z$ diverge!

As a result, you should probably check your math to make sure you didn't goof up somewhere.

Edit: After correction, you can do the integral symbolically:

Integrate[(ϵ)/((ϵ - x^2 - 
        y^2)^2 + ϵ^2) y, {x, -∞, ∞}, {y, 
  0, ∞}, Assumptions -> ϵ > 0]

which gives

1/2 I π (Sqrt[(-1 - I) ϵ] - Sqrt[(-1 + I) ϵ])

which when plotted as a function of $\epsilon$ looks like this:

Additional code:

hue = Compile[{{z, _Complex}}, {(1.0 Arg[-z] + π)/(2 π), 
    Exp[1 - Max[Abs[z], 1]], Min[Abs[z], 1]}, 
   CompilationTarget -> "C", RuntimeAttributes -> {Listable}];
ComplexPlotR2[f_, {x0_, x1_, δx_}, {y0_, y1_, δy_}] := 
  Image[hue[
     f[#[[All, All, 1]], #[[All, All, 2]]] &@
      Outer[List, Range[x0, x1, δx], 
       Range[y0, y1, δy]]]\[Transpose], ColorSpace -> Hue, 
   Magnification -> 1];
CCompileR2[expr_] := 
  Compile[{{x, _Real}, {y, _Real}}, Evaluate[expr], 
   CompilationTarget -> "C", RuntimeAttributes -> {Listable}];