# Speed up NIntegrate, which takes eigenvalues from another function

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.

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.

• Thank you so much. This was very helpful. One question: when should I use infinite precision (1/2), when to use floating point (0.5) and when it is useful to define exact precision (N [1/2, 30])? I know infinite precision is usually slower but it sometimes yields wrong results. – gurluk Jan 5 '15 at 1:03
• @gurluk What precision to use seems complicated to me. I like to start with exact input (infinite precision, 1/2) even to NIntegrate etc., unless there is a reason not to. Symbolic processing is usually easier for M with exact numbers. Of course, exact answers, if desired, require exact input. Machine floating-point numbers (0.5) are fast and for approximate answers. If they're not accurate, I use arbitrary precision (N[1/2, 30]). When numerical functions suggest higher working precision, I check the stability of the result with that. In some cases, sensitivity analysis is important. Etc. – Michael E2 Jan 5 '15 at 3:53

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?

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?

• Thank you for your help. There is no reason to choose LocalAdaptive, I was just experimenting with methods and LocalAdaptive was the last when I copied my code here. I am trying your suggestions now. – gurluk Jan 5 '15 at 0:58