Numerical integration of $\frac{1}{x} \int_0^x \frac{t}{e^t-1}dt$ I'm trying to compute numerically the non-elementary function below :
$$
\frac{1}{x} \int_0^x \frac{t}{e^t-1}\mathrm{d}t
$$
I tried to evaluate the integral with the quad algorithm, but at values of x < 1, the result is inaccurate compared with the tabled values, maybe because of the indetermination of the integrand at $t=0$.
I also tried to compute the integral with the Simpsons method in which I manually replace the first term of the integrand vector by the value of the limit :
$$
\lim_{t\to 0} \frac{t}{e^t-1} = \lim_{t \to 0} \frac{1}{e^t}=1
$$
It takes however a large number of intervals to attain a precision as low as 3 digits, and it has poor performance.
Any ideas to compute this function more efficiently? Or maybe this function is already computed in a Scipy algorithm?
Thanks for your help
 A: Recall the useful Riemann zeta function:
$$\int_0^\infty\frac t{e^t-1}\ dt=\zeta(2)\Gamma(2)=\frac{\pi^2}6$$
Thus, we know that
$$\frac1x\int_0^x\frac t{e^t-1}\ dt=\frac1x\left(\frac{\pi^2}6-\int_x^\infty\frac t{e^t-1}\ dt\right)$$
This is much easier, since, that's to the geometric series, we can see that
$$\int_x^\infty\frac t{e^t-1}\ dt=\sum_{n=1}^\infty\int_x^\infty te^{-nt}\ dt=\sum_{n=1}^\infty\frac{e^{-nx}(nx+1)}{n^2}=-x\ln(1-e^{-x})+\sum_{n=1}^\infty\frac{e^{-nx}}{n^2}=-x\ln(1-e^{-x})+\operatorname{Li}_2(e^{-x})$$
And for large values of $x$, the first few values of this series suffice.

$$\frac1x\int_0^x\frac t{e^t-1}\ dt=\frac1x\left(\frac{\pi^2}6-\sum_{n=1}^\infty\frac{e^{-nx}}{n^2}\right)+\ln(1-e^{-x})$$

All of this may be found at the end of the Polylogarithm: Relationships to other functions section of Wikipedia.
A: $$\int_{0}^{x}\frac{t\,dt}{e^t-1}=\frac{\pi^2}{6}-\int_{e^x}^{+\infty}\frac{\log(t)\,dt}{t(t-1)}=\frac{\pi^2}{6}+\int_{0}^{e^{-x}}\frac{\log(t)}{1-t}\,dt$$
equals
$$ \frac{\pi^2}{6}+\int_{1-e^{-x}}^{1}\frac{\log(1-t)}{t}\,dt =\text{Li}_2(1-e^{-x})=\sum_{n\geq 1}\frac{(1-e^{-x})^n}{n^2}$$
where the truncated series provides good approximations for large values of $x$.
For any $x\in(0,1)$ it is much better to use the truncated Taylor series:
$$ \int_{0}^{x}\frac{t\,dt}{e^t-1} = x - \frac{x^2}{4} + \frac{x^3}{36}+ O(x^5).$$
A: A clever way is to break the integral into two pieces, and massage the function in the piece having $t<1$ to avoid numerical error accumulation:
$$
\int_0^x\frac{t}{e^t-1}= 
\int_0^1\frac{t}{e^t-1}+ \int_1^x\frac{t}{e^t-1} 
$$
You are having no trouble with the second piece. For the first piece, the trick is to Taylor expand about $t=0$:
$$
\frac{t}{e^t-1}= 1-\frac12 t+\frac1{12}t^2-\frac1{6!}t^4+\frac1{6\cdot 7!}t^6-\frac1{60\cdot 8!}t^8+\frac1{132\cdot9!}t^{10}+ \cdots
$$
Integrating this from zero to one is very easy.
A: $\newcommand{\bbx}[1]{\,\bbox[8px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
In order to perform a numerical integration you must take care of the $\ds{t \over \expo{t} - 1}$ 'integrable singularity'. Note that
$\ds{{t \over \expo{t} - 1}  \sim 1 - {1 \over 2}\,t + {1 \over 12}\,t^{2} - {1 \over 720}\,t^{4} + \,\mrm{O}\pars{t^{6}}}$. 
$$
\mbox{Replace}\quad {t \over \expo{t} - 1}\quad\mbox{by}\quad\mrm{f}\pars{t} \equiv
\left\{\begin{array}{lcl}
\ds{-1} & \mbox{if} & \ds{\phantom{|,}t\phantom{\,|} <\ \texttt{TOL_0}}
\\[2mm]
\ds{\phantom{-}0} & \mbox{if} & \ds{\phantom{|\,}t\phantom{\,|} >\ \texttt{-TOL_0}}
\\[2mm]
\ds{t \over \expo{t} - 1} & \mbox{if} & \ds{\verts{t} > \texttt{TOL_1}}
\\[2mm]
\ds{1.0 - 0.5\,t\pars{1.0 - {t \over 6.0}}} && \mbox{Otherwise}
\end{array}\right.
$$  

$\ds{\texttt{TOL_0}}$ and $\ds{\texttt{TOL_1}}$ are defined in terms of the
  Machine Precision $\texttt{MP}$
  as follows:

$$
\texttt{TOL_0} \equiv \ln\pars{\texttt{MP}}/1.05\,;\qquad
\texttt{TOL_1} \equiv  1.05\pars{720.0\,\,\,\texttt{MP}}^{1/4}
$$
The $\ds{1.05}$-factor is included to avoid rounding errors $\ds{\pars{~namely,\ 5\ \% ~}}$.
$$\bbox[#ffe,10px,border:1px groove navy]{\ds{%
\mbox{With}\ \,\mrm{f}\pars{t}\mbox{, as given above, you can use}\ safely\ \mbox{the}\ Simpson\ Method.}}
$$
For the purpose of the example, I'll include a $\texttt{javascript}$ code
$\ds{\pars{~\mbox{you can run it in a}\ terminal\ \mbox{with}\ \texttt{node}:
\texttt{node integ0.js}~}}$:

// integ0.js

"use strict";
var MACHINEPRECISION = (function() // Machine Precision ( by definition )
{
 var machPrec = 1.0
 var     temp = 1.0;

 while((1.0 + temp) > 1.0) {
      machPrec = temp;
      temp /= 2.0;
 }

 return machPrec;
})();


var f = (function ()
{
 var TOL_0 = Math.log(MACHINEPRECISION)/1.05;
 var TOL_1 = 1.05*Math.pow(720.0*MACHINEPRECISION,0.25);

 console.log("\nTOL_0 = " + TOL_0 + "\nTOL_1 = " + TOL_1);

 return function (t)
 {
  if (t < TOL_0)           return -t;
  if (t > -TOL_0)          return 0;
  if (Math.abs(t) > TOL_1) return t/(Math.exp(t) - 1.0);
                           return 1.0 - 0.5*t*(1.0 - t/6.0);
 };
})();

console.log("\nMACHINE PRECISION = " + MACHINEPRECISION + "\n");


It yields:

$\begin{array}{l}
\ds{\texttt{TOL_0 = -34.32728894201634}}
\\
\ds{\texttt{TOL_1 = 0.0006639455730754197}}
\\
\ds{\texttt{MACHINE PRECISION = 2.220446049250313e-16}}
\\
\end{array}$

Try to implement the Simpson Rule.

