# NIntegrate(with singularities) gives different results, which one should I trust?

I am doing an Numerical integral with sigularities.

tk = Cos[Sqrt/2 kx + 1/2 ky] + Cos[-Sqrt/2 kx + 1/2 ky] + Cos[ky]
+ I (Sin[Sqrt/2 kx + 1/2 ky] + Sin[-Sqrt/2 kx + 1/2 ky] - Sin[ky])


The integral I want do is:

int1 = NIntegrate[1/Norm[tk], {kx, 0, 2 Sqrt Pi/3}, {ky, 0, 4 Pi/3}]


Without excluding singularities, I got:

57.6751


This integral has singularities in 3 Places:

{4 Pi/(3 Sqrt), 0}, {2 Pi/(3 Sqrt), 2 Pi/3}, {4 Pi/(3 Sqrt), 2 Pi/3}} So I will Exclude them using:

int1 = NIntegrate[1/Norm[tk], {kx, 0, 2 Sqrt Pi/3}, {ky, 0, 4 Pi/3},
Exclusions ->{{4 Pi/(3 Sqrt), 0}, {2 Pi/(3 Sqrt), 2 Pi/3},
{4 Pi/(3 Sqrt), 2 Pi/3}}]


This result is quite different from which I didn't exclude any singular points. However, it gives me a result:

 13.6216 + 0. I


Q:Why a complex term occurs? It looks quiet strange to me. However, excluding other points which are not singularities also gives nearly the same answer. For example:

int1 = NIntegrate[1/Norm[tk], {kx, 0, 2 Sqrt Pi/3}, {ky, 0, 4 Pi/3},
Exclusions -> {{1, 1}}]


Gives:

13.6215


Excluding {0,0} gives:

57.6751


This is even stranger for me. What happened? Can you explain to me? This is very untrustworthy for me. How can I know this numerical integral is right or wrong?

Another integral I want to do is:

int1 = NIntegrate[1/Abs[Norm[tk]-1], {kx, 0, 2 Sqrt Pi/3}, {ky, 0, 4 Pi/3}]


This has singular lines: With an without excluding this line, the integral is a finite number. However, I expect this integral to be infinite. What is wrong?

For your first integral, NIntegrate gives warnings or messages about failing to converge when not excluding the singular points.

When excluding the singular points I trust the result 13.6216 + 0. I because some alternative methods agree. Monte Carlo sampling:

In:= NIntegrate[
1/Abs[tk], {kx, 0, 2 Sqrt Pi/3}, {ky, 0, 4 Pi/3},
Method -> "MonteCarlo", PrecisionGoal -> 2]

Out= 13.7307 + 0. I


Outer product of one-dimensional rules together with higher precision:

In:= int1 =
NIntegrate[1/Abs[tk], {kx, 0, 2 Sqrt Pi/3}, {ky, 0, 4 Pi/3},
Exclusions -> {{4 Pi/(3 Sqrt), 0}, {2 Pi/(3 Sqrt),
2 Pi/3}, {4 Pi/(3 Sqrt), 2 Pi/3}},
Method -> "GaussKronrodRule", WorkingPrecision -> 20,
PrecisionGoal -> 12]

Out= 13.621641914220150760


Your second question was about why there is the tiny imaginary part + 0. I. This is because, when NIntegrate compiles your integrand to a CompiledFunction for faster execution, Compile notices that the expression contains at least one thing that could be complex-valued (namely, it contains I). Therefore it compiles to a numerical function returning a machine-precision complex value rather than a machine-precision real value.

Of course, in theory Mathematica could try to prove (easily in this case) that the expression always has imaginary part == 0, but NIntegrate doesn't bother doing this because the existence of an imaginary part does not affect NIntegrate's choice of methods to apply.

You can see this (at the cost of slower execution) by disabling automatic compilation:

In:= NIntegrate[
1/Abs[tk], {kx, 0, 2 Sqrt Pi/3}, {ky, 0, 4 Pi/3},
Exclusions -> {{4 Pi/(3 Sqrt), 0}, {2 Pi/(3 Sqrt),
2 Pi/3}, {4 Pi/(3 Sqrt), 2 Pi/3}}, Compiled -> False]

Out= 13.6216


If it bothers you, you could also pre-compile the numerical function yourself, or just take the Re part afterward.

int1 = NIntegrate[1/Abs[Norm[tk]-1], {kx, 0, 2 Sqrt Pi/3}, {ky, 0, 4 Pi/3}]