Computation of a general integral of Fermi–Dirac distribution by McDougall–Stoner method

We extended the method of McDougall and Stoner (1938) [15] to the computation of a general integral of the Fermi–Dirac distribution, F(η)F(η). When η>0η>0, the new method splits F(η)F(η) into a sum of three parts, A(η),B(η)A(η),B(η), and C(η)C(η), and integrates them exactly and/or numerically by means of the double-exponential quadrature rules. As η   increases, B(η)/F(η)B(η)/F(η) damps algebraically while C(η)/F(η)C(η)/F(η) decays exponentially. Thus, they can be ignored when η   exceeds certain threshold values depending on the input error tolerance. When A(η)A(η) is exactly computable and η   is sufficiently large such that B(η)B(η) and C(η)C(η) are negligible, the new method was 60–1000 times faster than the numerical integration of F(η)F(η). Even if B(η)B(η) and/or C(η)C(η) are not negligible, their smallness and rapid decay significantly accelerate their numerical quadrature so that the speed-up factor becomes 3–20 unless η   is less than 2.0–4.5. On the other hand, when A(η)A(η) is numerically integrated, the speed-up factor diminishes to 2–5. Still, the superiority of the new method to the numerical integration of F(η)F(η) is unchanged.