# Trouble with PolyLog function evaluation

New poster here, please forgive any formatting issues.

As part of a physics simulation I'm using an analytical expression for the definite integral

$$f(a,b) = \int\limits_{x_0}^{x_1}\frac{x^2}{e^{ax + b}-1}dx$$

in Mathematica. (It's related to the number density of particles in a Bose-Einstein distribution).

For simplicity let's consider the integral from x = 0 to 1. Also a and c are strictly positive. Then Mathematica gives this function containing PolyLogs

(a^2 (-a + 3 I \[Pi] + 3 Log[-1 + E^(a + b)]) +  6 a PolyLog[2, E^(a + b)] + 6 PolyLog[3, E^b] -  6 PolyLog[3, E^(a + b)])/(3 a^3)


which using some identities can be written in manifestly real form as

-((a^2 (a + b - Log[-1 + E^(a + b)]) + 2 a PolyLog[2, E^(-a - b)] +  2 PolyLog[3, E^(-a - b)] - 2 PolyLog[3, E^-b])/a^3)


Now for most values of a and b it works fine, but 3DPlotting the expression as a function of a and b reveals some problems. In particular, when a is small things get pathological.

For probably related reasons, Integrate and NIntegrate give different results for the integral over f(a,b) for e.g. a = 0.00002, b = 2. The analytical expression gives a third result that is even negative for some choices of parameters.

Plotting the integrand reveals a perfectly innocuous looking strictly positive function, btw.

Is there any way I can get the analytical expression to numerically evaluate to correct results for small a?

Thanks!

• Try using WorkingPrecision->30 in your plot. – Carl Woll Mar 23 '17 at 23:39
• I'm really not fond of expressing the Bose-Einstein and Fermi-Dirac functions as compositions of the polylogarithm and $\exp$; so much room for numerical mischief. – J. M.'s discontentment Mar 24 '17 at 0:24
• f[a0_?NumericQ, b0_?NumericQ] := With[{a = SetPrecision[a0, Infinity], b = SetPrecision[b0, Infinity]}, N[(a^2 (-a + 3 I \[Pi] + 3 Log[-1 + E^(a + b)]) + 6 a PolyLog[2, E^(a + b)] + 6 PolyLog[3, E^b] - 6 PolyLog[3, E^(a + b)])/(3 a^3), 3.1]]; <-- change the 3.1 to your desired precision. – Michael E2 Mar 24 '17 at 1:00
• Thank you Michael! That works perfectly. I was trying to change the precision before but evidently I need to invest some more time in learning exactly how that works :) – ButILikeMyStack Mar 24 '17 at 8:44

Following up J.M.'s comment, here's a numeric way to compute somewhat quickly the OP's function via Gaussian quadrature. Some effort has to be put in to determine the number of Gauss-Legendre points to achieve a certain precision. (This could be done somewhat more slowly with NIntegrate of course. NIntegrate is about the same speed as using the adaptive precision approach in my comment, which was used to measure the error below. The approach below is about 50 times faster than either N[] on exact input or NIntegrate[].)

With[{grd41 = NIntegrateGaussRuleData[41, MachinePrecision],
grd21 = NIntegrateGaussRuleData[21, MachinePrecision]},
nf[a_?NumericQ, b_?NumericQ] :=
With[{data = Piecewise[{{grd41, a + 0.05 > 5 b}}, grd21]},
With[{x = data[[1]]},
(x^2/(Exp[a x + b] - 1)).data[[2]]
]];
]


Basic plot:

Plot3D[nf[a, b], {a, 0, 2}, {b, 0, 2}]


Relative error:

Plot3D[(nf[a, b] - #)/# &@ f[a, b] // RealExponent,
{a, 0, 1.95}, {b, 0, 2},
MaxRecursion -> 2, PlotRange -> {-17.6, -5},
AxesLabel -> {a, b, HoldForm@Log10["Rel. err."]}]


It is even faster to use 21 Gauss points throughout, but the error is somewhat greater near b == 0:

With[{data = NIntegrateGaussRuleData[21, MachinePrecision]},
With[{x = data[[1]], weights = data[[2]]},
nf[a_?NumericQ, b_?NumericQ] := (x^2/(Exp[a x + b] - 1)).data[[2]]
]]


• One could use a Gauss-Kronrod pair from NIntegrateGaussKronrodRuleData[] and use one estimate as the answer and the other as a precision check, among other possibilities. It's also a happy coincidence that the OP's desired integration range is the same as the default interval used by the quadrature rules in Mathematica. – J. M.'s discontentment Mar 25 '17 at 3:15