Speed up NIntegrate, which takes eigenvalues from another function

Question

I have simplified my problem here, in my actual problem the matrix is much bigger which makes it impossible to find eigenvalues analytically. So, I chose standard BCS problem (2x2 matrix) to demonstrate my problem with Mathematica. But, with even this simple case, it takes forever to integrate. My code is the following:

hamiltonian[k_, mu_, delta_] := {{k^2 - mu, delta}, {delta, -k^2 + mu}}
eigens[k_, mu_, delta_] := Eigenvalues[hamiltonian[k, mu, delta]]
fermitotal[beta_, k_, mu_, delta_] := Block[{ee = eigens[k, mu, delta]},   1/(1 + Exp[beta ee[[1]]]) + 1/(1 + Exp[beta ee[[2]]])]
nTotal[beta_, mu_, delta_] := NIntegrate[fermitotal[beta, k, mu, delta] k^2, {k, 0, 20}, Method -> {"LocalAdaptive", "SymbolicProcessing" -> False}, AccuracyGoal -> 4, PrecisionGoal -> 4, MinRecursion -> 10, MaxRecursion -> 300, WorkingPrecision -> 13]
nTotal[50,1,1/10]//AbsoluteTiming

Any idea how to speed up the calculation?

I have tried, using Compile for the hamiltonian matrix. Specifying the precision for the numbers,

nTotal[N[50,30],N[1,30],N[1/10,30]]//AbsoluteTiming

Did not help. I am completely stuck here and I don't know what causes this slowness. Any help would be very much appreciated.

Answer 1

The given integral can be integrated exactly and quickly.

Integrate[fermitotal[beta, k, mu, delta] k^2, {k, 0, 20}] // AbsoluteTiming
N@Last[%]
(*
  {0.006876, 8000/3}
  2666.67
*)

But since, presumably, the OP's actual use-case is not, I'll comment on the set up and the relation of OP's choices to speed.

Remarks on the OP's option settings

Since it's a toy example, it's hard to speculate further about the problem. The "LocalAdaptive" strategy is usually slower than "GlobalAdaptive" and sometimes a lot slower. The integrand

fermitotal[beta, k, mu, delta] k^2

does not look problematic. "GlobalAdaptive" should be sufficiently robust and considerably quicker.

Using arbitrary precision, with WorkingPrecision -> 13, will be slower than machine precision, even when the precision is less. Arbitrary precision will usually give more accurate results, but I would use it only when the machine precision result is inaccurate.

The setting MinRecursion -> 10 forces the initial partition of the interval to be subdivided over and over 10 times (so at least 2^10 times the number of initial subintervals), before NIntegrate begins to analyze the error on the subintervals. This saves times only when the integrand would require most intervals to be subdivided 10 times. The OP's example integrand does not require it, and probably many more subdivisions are done than necessary.

MaxRecursion -> 300 seems absurdly large for someone not wanting to wait a long time. That could mean over 10^90 function evaluations (it probably wouldn't). I suppose it might be reasonable in cases where the integrand had a nasty singularity that affected a very small region. One can often handle such singularities in a more efficient way.

"SymbolicProcessing" -> False sometimes saves time and sometimes doesn't. It depends on what you gain by the preprocessing.

Setting the AccuracyGoal and PrecisionGoal lower, as in the OP's code, generally saves time while lowering the quality of the result.

Summary

The given example is quite simple to deal with. It may be that it does not capture enough of the problem to make any potential answer here helpful in the OP's actual use-case. On the other hand maybe it is a case of one or two misunderstood NIntegrate options.

Answer 2

How about this:

fermitotal[beta_?NumericQ, k_?NumericQ, mu_?NumericQ, delta_?NumericQ] :=
With[{ee = Eigenvalues[{{k^2 - mu, delta}, {delta, -k^2 + mu}}]},
    1/(1 + Exp[beta*ee[[1]]]) + 1/(1 + Exp[beta*ee[[2]]])
];


nTotal[beta_?NumericQ, mu_?NumericQ, delta_?NumericQ] := NIntegrate[fermitotal[beta, k, mu, delta] *k^2, {k, 0, 20}]

nTotal[50, 1, 1/10] // AbsoluteTiming

{0.435001, 2666.67}

I hope this is fast enough?

Answer 3

I'm afraid that this is a bit puzzling indeed. If I use your example and run the NIntegrate without all the fancy options I get {0.006944, 2666.67}. So I assume it has something to do with your options.

Leaving out the method, but leaving everything else in place, slows down things dramatically, but it is still doable: {0.561119, 2666.666666667}, additionally dropping the MinRecursion option speeds it up to the 'naked' value.

Just dropping the MinRecursion but keeping the Method results in longer computational time than I'm willing to wait. So I guess it is a fair guess that the settings of 'Method' are the culprit. Is there any good reason why you don't let Mathematica decide on its own, which method and which strategy to use?