3
$\begingroup$

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.

trouble

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!

$\endgroup$
  • 1
    $\begingroup$ Try using WorkingPrecision->30 in your plot. $\endgroup$ – Carl Woll Mar 23 '17 at 23:39
  • 1
    $\begingroup$ 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. $\endgroup$ – J. M. will be back soon Mar 24 '17 at 0:24
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Mar 24 '17 at 1:00
  • $\begingroup$ 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 :) $\endgroup$ – ButILikeMyStack Mar 24 '17 at 8:44
2
$\begingroup$

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 = NIntegrate`GaussRuleData[41, MachinePrecision],
  grd21 = NIntegrate`GaussRuleData[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}]

Mathematica graphics

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."]}]

Mathematica graphics

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

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

Mathematica graphics

$\endgroup$
  • $\begingroup$ One could use a Gauss-Kronrod pair from NIntegrate`GaussKronrodRuleData[] 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. $\endgroup$ – J. M. will be back soon Mar 25 '17 at 3:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.