Let me explain the problem. I am trying to integrate a one dimensional integral: $$\int {g\left( {{k_x}},parameter1,parameter2,...\right)d{\mkern 1mu} {k_x}} $$
for the sake of clarity, I will give function g in the last.
I choose a specific set of parameter, and do the following integration:
In:=NIntegrate[g[kx, 9, -2, 7, 1, 1, 0.05], {kx, -\[Pi]/3,\[Pi]/3}]//AbsoluteTiming
Out:={422.396160, 0.163126 + 0.103155 I}
This takes me 422s on my computer!!!
I use the following code to show the sampling points used during the integration.
sampp = Reap[
NIntegrate[#, {kx, -\[Pi]/(3), \[Pi]/(3)},
EvaluationMonitor :> Sow[{kx, #}]]] &[
g[kx, 9, -2, 7, 1, 1, 0.05]] ;
{Length[#], ListPlot[#, AxesOrigin -> {0, 0}, Filling -> Axis]} &[
relist@sampp[[2, 2, 1]]]
it shows that it used 733 sample points. But should that take 400 seconds?
Actually, the evaluation of the integrand at each sample point is acually fast. I use the most simple interpolation method to do the same integration as below, I use 1000 sample points that is evenly distributed and it only takes 1.3 second. Besides the integration result is quite accurate.
In:=(sampp = Table[{kx,
g[kx, 9, -2, 7, 1, 1, 0.05]}, {kx, -\[Pi]/3, \[Pi]/3,
2 \[Pi]/3/1000}];
Integrate[
Interpolation[sampp, InterpolationOrder -> 2][
x], {x, -\[Pi]/3, \[Pi]/3}]) // AbsoluteTiming
Out:={1.308075, 0.163126 + 0.103155 I}
So what did Mathematica do in the extremely long 400 second? Of course, it needs to do many error estimation, but that should not take too much time.
I scan the documents and related questions in stackexchange. Somebody suggest to set "SymbolicProcessing"->0. OK, I tried
In:=NIntegrate[g[kx, 9, -2, 7, 1, 1, 0.05], {kx, -\[Pi]/3, \[Pi]/3},
Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule",
"SymbolicProcessing" -> 0}]//AbsoluteTiming
Out:={401.791981, 0.163126 + 0.103155 I}
unfortunately, the result is the same.
Even more peculiar, when I use "Trapezoidal" strategy and limit the MaxRecursion
In:=NIntegrate[g[kx, 9, -2, 7, 1, 1, 0.05], {kx, -\[Pi]/3, \[Pi]/3},
Method -> "Trapezoidal", MinRecursion -> 1, MaxRecursion -> 4]
Out:={413.894674, 0.163126 + 0.103155 I}
the sample points shows as follows
it only used 257 sample point, but the time cost is the same??!!It is unbelievable! Anyway, I think the "Trapezoidal" should be the same as evenly sample interpolation method.
I really can't understand what did Mathematica do in the NIntegrate of function g. Can somebody explain it?
Of course, it seems that I can use the simple interpolation integration method manually, because it is fast, and the result seems good in this case. But it is not robust, I insisted using the built in "Adaptive" strategy, because I have to Integrate some other functions which are sharp at specific points.
In the last
the form of function g
g[kx_, xx_, e_, width_, ii_, label_, \[Eta]_] :=
Inverse[(e - I \[Eta]) IdentityMatrix[2*width] -
armchairibbonmat[kx, width]][[2*(ii - 1) + label,
2*(ii - 1) + label]] E^(I kx xx)
(*armchairibbonmat is a function used in function g*)
armchairibbonmat[kx_, n_] :=
(
h = Table[0, {i, 1, 2 (n + 2)}, {j, 1, 2 (n + 2)}];
a[m_] := 2 m - 1;
b[m_] := 2 m;
t1 = 1;
aa = 1;
Do[
h[[b[i], a[i]]] = t1 E^(I kx aa);
h[[a[i], b[i]]] = t1 E^(-I kx aa);
h[[b[i + 1], a[i]]] = t1 E^(-I kx aa/2);
h[[b[i - 1], a[i]]] = t1 E^(-I kx aa/2);
h[[a[i - 1], b[i]]] = t1 E^(I kx aa/2);
h[[a[i + 1], b[i]]] = t1 E^(I kx aa/2);
, {i, 2, n + 1}];
h = ArrayPad[h, -2]
)
