# difficulty in numerical estimate of double integral

I entered the following simple integral in mathematica:

NIntegrate[1/(0.01)^2*Exp[-(R - 1)^2/(0.01)^2] (1/Abs[z]^(1/3) - 1/(R^2 +z^2)^(1/6)),
{R,0, 10000}, {z, -100000, 100000}]


and it returns 0, which is really strange. The result should not be zero. Now, instead of $\exp(-(R-1)^2/(0.01)^2)$, I try $\delta(R-1)$, after which the integral is

NIntegrate[DiracDelta[R-1]* (1/Abs[z]^(1/3) - 1/(R^2 +z^2)^(1/6)), {R,0, 10000},
{z, -100000, 100000}]


which should definitely yield a non-trivial answer, but Mathematica still returns 0. Any idea what is going on?

• Are your integration limits actually 100000, or is that supposed to be ∞? – J. M.'s ennui Mar 31 '16 at 3:05
• That is infact supposed to Infinity, but to avoid problems at branches, I just use some large finite number instead. Shouldn't matter any ways as the integrand is rapidly decaying with increasing $z$. – titanium Mar 31 '16 at 3:07
• You may have great loss of precision in your integration. Have you tried replacing the machine precision numbers (e.g. 0.01) with arbitrary precision ones, and increasing the WorkingPrecision setting? – MarcoB Mar 31 '16 at 3:39
• At least the DiracDelta one should be relatively straightforward, right? But Mathematica just gives zero even for that which is what troubles me. – titanium Mar 31 '16 at 3:41
• For the second integral: I don't think DiracDelta is meant to work with NIntegrate. If you Integrate it instead, you get an expression that evaluates to 1.29355 after applying N, and this agrees with the integral over z with 1 in place of R. – march Mar 31 '16 at 3:46

If we take the freedom - as J.M. suggested in a comment - of chosing infinity instead of large numbers as limits of integration and allow for a general parameter "a" in the exponent (and in front) the double integral can be done explicitly:

a Integrate[Exp[-a (R - 1)^2] (1/Abs[z]^(1/3) - 1/(R^2 + z^2)^(1/6)), {R, 0, \[Infinity]}, {z, -\[Infinity], \[Infinity]}, Assumptions -> a > 0]

(* Out[8]= - a ((Sqrt[\[Pi]]
Gamma[-(1/3)] (Gamma[5/6] Hypergeometric1F1[-(1/3), 1/2, -a] +
2 Sqrt[a] Gamma[4/3] Hypergeometric1F1[1/6, 3/2, -a]))/(
2 a^(5/6) Gamma[1/6])) *)


You can then insert a = 10^4 to get

%/. a -> 10^4 // N

(* Out[11]= 229.275 *)


This value should be very close to the original integral with the finite Limits of integration.

In your second integral you have tried to replace the exponential by DiracDelta[]. But the replacement is not consistent, because the factor in front should be chosen accordingly:

2 Sqrt[ a/\[Pi] ]
Integrate[Exp[-a x^2], {x, 0, \[Infinity]}, Assumptions -> a > 0]

(* Out[18]= 1 *)


Furthermore with that intention it is logical to let R->0 and calculate the z-Integral which gives

Integrate[(1/Abs[z]^(1/3) - 1/(R^2 + z^2)^(1/6)) /.
R -> 1, {z, -\[Infinity], \[Infinity]}]

(* Out[21]= -((Sqrt[\[Pi]] Gamma[-(1/3)])/Gamma[1/6]) *)

% // N

(* Out[22]= 1.29355 *)