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

I am doing an Numerical integral with sigularities.

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


The integral I want do is:

int1 = NIntegrate[1/Norm[tk], {kx, 0, 2 Sqrt[3] 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[3]), 0}, {2 Pi/(3 Sqrt[3]), 2 Pi/3}, {4 Pi/(3 Sqrt[3]), 2 Pi/3}}


So I will Exclude them using:

int1 = NIntegrate[1/Norm[tk], {kx, 0, 2 Sqrt[3] Pi/3}, {ky, 0, 4 Pi/3},
Exclusions ->{{4 Pi/(3 Sqrt[3]), 0}, {2 Pi/(3 Sqrt[3]), 2 Pi/3},
{4 Pi/(3 Sqrt[3]), 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[3] 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[3] 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?

-

## 1 Answer

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[76]:= NIntegrate[
1/Abs[tk], {kx, 0, 2 Sqrt[3] Pi/3}, {ky, 0, 4 Pi/3},
Method -> "MonteCarlo", PrecisionGoal -> 2]

Out[76]= 13.7307 + 0. I


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

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

Out[79]= 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[85]:= NIntegrate[
1/Abs[tk], {kx, 0, 2 Sqrt[3] Pi/3}, {ky, 0, 4 Pi/3},
Exclusions -> {{4 Pi/(3 Sqrt[3]), 0}, {2 Pi/(3 Sqrt[3]),
2 Pi/3}, {4 Pi/(3 Sqrt[3]), 2 Pi/3}}, Compiled -> False]

Out[85]= 13.6216


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

(Your third question about

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


I can't offer anything yet.)

-
Thanks very much about your detailed answer, it really helps me a lot. Especially the idea that check the result using different integral schemes. However I got no warning or message about failing to coverage or something. It simply gives the wrong answer. I am using Mathematica 9 on windows7 –  luming Mar 19 at 3:52
Yeah, I found the error message when execute on Linux MatheKernel. It indeed gives the warnings: NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 18 recursive bisections in ky near {kx, ky} = {1.20921, 2.09441}. NIntegrate obtained 13.6215 and 0.0000840493 for the integral and error estimates. –  luming Mar 19 at 4:47
And in MathKenal, I can directly see that my integral can't coverage: NIntegrate::maxp: The integral failed to converge after 50100 integrand evaluations. NIntegrate obtained 3370.06 + 0. I and 3323.2 for the integral and error estimates. –  luming Mar 19 at 4:59