Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

Three Singularities

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:enter image description here

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

share|improve this question

1 Answer 1

up vote 2 down vote accepted

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.)

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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