0
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ The integral with respect to r is 1/4 (\[Pi] - 2 ArcTan[1 - 10 x^2]). Integrating it with respect to x does not converge. Does that seem right? $\endgroup$ – Michael E2 Nov 3 '14 at 0:31
  • $\begingroup$ I'll post an answer in a while which should hopefully clear things up, but I think you might need to recheck your integrand for correctness; it does not decay as $1/x^2$ like you say it does. $\endgroup$ – DumpsterDoofus Nov 3 '14 at 0:44
  • $\begingroup$ Minor note: you write an $\omega$ in the denominator, but it later turns into an $\epsilon=0.1$. Is that correct? $\endgroup$ – DumpsterDoofus Nov 3 '14 at 0:53
  • $\begingroup$ I have set, $\omega = 0.1$. You are right, I did not realize that the integrand was not decaying $1\x^2$ while I was typing this, sorry. Let me carefully check my calculations, this is supposed to be a simple problem. I think the problem is with my calculations, maybe I mistyped something... $\endgroup$ – gurluk Nov 3 '14 at 1:18
4
$\begingroup$

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}]

enter image description here

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:

enter image description here

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}];
$\endgroup$
  • $\begingroup$ You are right that this integral is divergent, I did not realize that while I was typing this, sorry. I am going to check my calculations... Thank you for your attention, and shame on me! $\endgroup$ – gurluk Nov 3 '14 at 1:23
  • $\begingroup$ @gurluk: It's fine, don't worry. Maybe there's a sign problem, and maybe it's supposed to be $\frac{\epsilon r}{\epsilon^2 + (\omega+z^2 + r^2)^2}$? That would decay like $1/r^2$ and $1/z^2$. $\endgroup$ – DumpsterDoofus Nov 3 '14 at 1:23
  • $\begingroup$ Yes, it was a sign typo :) Both in code and here... I have corrected and added my comments in parantheses. $\endgroup$ – gurluk Nov 3 '14 at 1:44
  • $\begingroup$ But, for my real calculation, the function inside the parentheses behaves like x^2, and my real calculation is still giving me some huge, non-sense results. $\endgroup$ – gurluk Nov 3 '14 at 1:47
  • 1
    $\begingroup$ @gurluk: On my 7-year-old laptop (LOL), when I execute NIntegrate[(1/10)/((1/10 - x^2 - r^2)^2 + 1/100) r, {x, -\[Infinity], \[Infinity]}, {r, 0, \[Infinity]}] // AbsoluteTiming, it gives {0.236107, 1.0915}, which means it took 1/4 of a second to evaluate, and it exactly matches the symbolic result. So the problems you are encountering when numerically integrating your complicated expression do not exist in the toy problem you gave us, so we can't help you unless you describe further what you're trying to integrate. $\endgroup$ – DumpsterDoofus Nov 3 '14 at 2:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.