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 write a function called berrycur[kx,ky] which I will give at the end of the question, and want to numerically integrate this function over {kx, -2π, 2π]}, {ky, 0, 4π/Sqrt[3]}. The plot of berrycur is shown as follows:

enter image description here

I have checked that all values of berrycur in this region is positive. But the numerical integration

NIntegrate[ berrycur[kx, ky][[1]], {kx, -2π, 2π]}, {ky, 0, 4π/Sqrt[3]}]

gives result is zero !!!! This is absolutely wrong!!

Actually, the NIntegrate result of {kx, -2π, 0]}, {ky, 0, 4π/Sqrt[3]} is opposite to NIntegrate result of {kx, 0, 2π]}, {ky, 0, 4π/Sqrt[3]}, this is strange!

What is wrong here?


the definition of berrycur is

Clear[h]
h[kx_, ky_] := {{0.1` (-4 Cos[(Sqrt[3] ky)/2] Sin[kx/2] + 2 Sin[kx]), 
   E^((I ky)/Sqrt[3]) + 
    E^(-(1/6) I (3 kx + Sqrt[3] ky)) (1 + E^(I kx))}, {E^(-((I ky)/
     Sqrt[3])) + 
    E^(-(1/6) I (3 kx - Sqrt[3] ky)) (1 + E^(
       I kx)), -0.1` (-4 Cos[(Sqrt[3] ky)/2] Sin[kx/2] + 2 Sin[kx])}}
dim = Length@h[1, 1];
Clear[hpar1, hpar2];
hpar1[kx_, ky_] = D[h[kx, ky], kx];
hpar2[kx_, ky_] = D[h[kx, ky], ky];
Clear[purifyeigs];
purifyeigs[eigs_] := 
  Transpose@Sort@Transpose@{Re[eigs[[1]]], eigs[[2]]};

Clear[berrycur];
berrycur[kxkx_?NumericQ, kyky_?NumericQ] := Module[{eigs},
  eigs = purifyeigs@Eigensystem[h[kxkx, kyky]];
  Table[Im@
    Sum[((Conjugate[eigs[[2, i]]].hpar1[kxkx, kyky].eigs[[2, 
            j]])*(Conjugate[eigs[[2, j]]].hpar2[kxkx, kyky].eigs[[2, 
            i]]) - (Conjugate[eigs[[2, i]]].hpar2[kxkx, kyky].eigs[[2,
             j]])*(Conjugate[eigs[[2, j]]].hpar1[kxkx, kyky].eigs[[2, 
            i]]))/(eigs[[1, i]] - eigs[[1, j]])^2, {j, 
      DeleteCases[Range[1, dim], i]}], {i, 1, dim}]
  ]
share|improve this question

2 Answers 2

up vote 1 down vote accepted

I believe your use of Part to extract a single value from the result of berrycur is causing NIntegrate to symbolically integrate your function. Try evaluating the indefinite integral for a clue as to what is going on, e.g.

Integrate[berrycur[x, y][[1]], x, y]

(x^2 y)/2

If you clean-up your function definition to return a single value:

berrycur2[kxkx_?NumericQ, kyky_?NumericQ] := 
 Module[{eigs, i = 1, j = 2}, 
  eigs = purifyeigs@Eigensystem[h[kxkx, kyky]];
  Im@((Conjugate[eigs[[2, i]]].hpar1[kxkx, 
          kyky].eigs[[2, j]])*(Conjugate[eigs[[2, j]]].hpar2[kxkx, 
          kyky].eigs[[2, i]]) - (Conjugate[eigs[[2, i]]].hpar2[kxkx, 
          kyky].eigs[[2, j]])*(Conjugate[eigs[[2, j]]].hpar1[kxkx, 
          kyky].eigs[[2, i]]))/(eigs[[1, i]] - eigs[[1, j]])^2]

Then you obtain

NIntegrate[berrycur2[x, y], {y, 0, 4 π/Sqrt[3]}, {x, -2 π, 2 π}]

12.5664

Which looks more reasonable.

share|improve this answer
    
Hi, MikeLimaOscar! Great explanation! How did you come up with the idea to test Integrate[berrycur[x, y][[1]], x, y]? This explains the unreasonable result. But I still want to know why Integrate gives (x^2 y)/2??? It seems that (x^2 y)/2 is the integration of x y in terms of x. –  matheorem Jun 26 at 0:56
    
What is more, the integration of berrycur is rather slow, it tooks me 40 seconds to give result 12.5664. Do you know how to speed it up? –  matheorem Jun 26 at 0:57
    
@matheorem Try Integrate[berrycur[x,y][[1]],x,y]//Trace. You'll see that Mathematica ends up with an integrand of x. –  MikeLimaOscar Jun 26 at 8:15
    
...as for speeding it up, all I can suggest is you look at the documentation for NIntegrate and experiment with the various rules. –  MikeLimaOscar Jun 26 at 8:24
    
Thank you very much! Trace is so useful! –  matheorem Jun 26 at 9:04

Delete the Table in berrycur,and then change to berrycur[kxkx_?NumericQ, kyky_?NumericQ, i_]

NIntegrate[berrycur[kx, ky, 1], {kx, -2 \[Pi], 2 \[Pi]}, {ky, 0, 4 \[Pi]/Sqrt[3]}]

12.5664

share|improve this answer

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.