Strange numerical integration result

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:

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

-

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.

-
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 '14 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 '14 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 '14 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 '14 at 8:24
Thank you very much! Trace is so useful! – matheorem Jun 26 '14 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

-