How to solve $\int x/(e^x-1) \, dx$? This might sound very basic, but I tried looking on integral tables and on wolphram math, but I couldn't get the method for solving it.
The steps I went through were:


*

*Part integration with $u=x$.

*That lead to the point I got stuck, $\int \ln(1-e^x)\,dx$. The output was $-\operatorname{Li2}(e^x)$.

*Final result was $\operatorname{Li2}(e^x) - x^2/2 + x\log(1-e^x)$


My questions are:


*

*Is there any shortcut for this integral? Maybe a way to get a cleaner result (i.e. without the Li2 term).

*How can I do it for the interval $0 < x < \infty$?


Here's the result from wolphram
http://www.wolframalpha.com/input/?i=%E2%88%AB+x%2F(e%5Ex+-1) 
Thanks
 A: If I didn't know anything about special functions, I would probably try this:
$$\frac{x}{e^x-1}=\frac{x e^{-x}}{1-e^{-x}}=x \sum_{n=1}^\infty e^{-nx}.$$
This equality holds (and the series converges) for $x>0$; near $x=0$ you may check that the integrand is bounded, so neglecting that point is no loss.
Then you have reason to hope that you can integrate term by term; this is an easy application of integration by parts, giving
$$\int x e^{-nx} dx = x e^{-nx}/(-n) + \int e^{-nx}/n dx = -x e^{-nx}/n - e^{-nx}/n^2.$$
So an antiderivative is
$$\sum_{n=1}^\infty -x e^{-nx}/n - e^{-nx}/n^2$$
and in particular
$$\int_0^\infty \frac{x}{e^x-1} dx = \sum_{n=1}^\infty 1/n^2 = \frac{\pi^2}{6}.$$
There are various convergence theorems that can be used to justify this interchange of infinite summation and integration.
But this is indeed a "special function", so you won't find an answer in terms of a finite combination of elementary functions.
A: To evaluate
\begin{equation}
\int \frac{x}{e^{x}-1} \mathrm{d}x
\tag{1}
\label{eq:161007-1}
\end{equation}
begin with
\begin{equation}
I_{1} = \int \frac{1}{e^{x}-1} \mathrm{d}x = \int \frac{e^{-x}}{1-e^{-x}} \mathrm{d}x = \ln(1-e^{-x})
\tag{2}
\label{eq:161007-2}
\end{equation}
and rewrite our result as
\begin{align}
\ln(1-e^{-x}) &= \ln[-e^{-x}(-e^{-x}+1)] \\
&= \ln(-1) + \ln(e^{-x}) + \ln(1-e^{x}) \\
&= i\pi \,- x + \ln(1-e^{x})
\end{align}
and thus
\begin{equation}
I_{1} = i\pi \,- x + \ln(1-e^{x})
\end{equation}
Now we evaluate equation \eqref{eq:161007-1}, our original integral, via integration by parts. For $\int a\mathrm{d}b = ab - \,\int b\mathrm{d}a$ we have $a = x$, $\mathrm{d}b = (e^{x}-1)^{-1}$, and $b = I_{1}$
\begin{equation}
\int \frac{x}{e^{x}-1} \mathrm{d}x = x I_{1} \,- \int I_{1} \mathrm{d}x
\end{equation}
\begin{equation}
\int I_{1} \mathrm{d}x = \int [i\pi \,- x + \ln(1-e^{x})] \mathrm{d}x = i\pi x \,- \frac{1}{2}x^{2} + \int \ln(1-e^{x}) \mathrm{d}x
\end{equation}
For the integral on the right hand side, we make the substitution $y = e^{x}$,
\begin{equation}
\int \ln(1-e^{x}) \mathrm{d}x = \int \frac{\ln(1-y)}{y} \mathrm{d}y = -\mathrm{Li}_{2}(y) = -\mathrm{Li}_{2}(e^{x})
\end{equation}
Putting all of the pieces together, we have
\begin{align}
\int \frac{x}{e^{x}-1} \mathrm{d}x &= i\pi x \,- x^{2} + x\ln(1-e^{x}) - \left(i\pi x \,- \frac{1}{2}x^{2} -\mathrm{Li}_{2}(e^{x}) \right) \\
&= \mathrm{Li}_{2}(e^{x}) + x\ln(1-e^{x}) \,- \frac{1}{2}x^{2} + \mathrm{const}
\end{align}
A: The integrand $\dfrac{x}{e^x-1}$ is a classical function without a name that generates the very important Bernoulli numbers. 
http://mathworld.wolfram.com/BernoulliNumber.html
Edit: Thus, as these constants (there are an infinite number of them!) are present in tables or can be computed, you can get a good approximated polynomial antiderivative. For example, after truncation above degree 5: $x-x^2/4+x^3/36-x^5/3600$.
But the primitive using the dilogarithm is also interesting in itself. Finding special functions in a computation might be a little disapointing, but it's real life. This said, the dilogarithm is a rather simple (and important) special function. See http://mathworld.wolfram.com/Dilogarithm.html
A: $$\int \frac{x}{e^{x} - 1}\ \text{d}x = \int \frac{x}{e^{x}(1 - e^{-x})}\ \text{d}x$$
The Geometric series always helps us
$$\int xe^{-x}\sum_{k = 0}^{+\infty} (e^{-x})^k\ \text{d}x = \sum_{k = 0}^{+\infty}\int xe^{-x}e^{-kx}\ \text{d}x = \sum_{k = 0}^{+\infty}\int xe^{-x(1+k)}\ \text{d}x $$
That integral, being indefinite, can be done by parts once (just to remove the $x$ from the integrand) and it'll give you
$$\sum_{k = 0}^{+\infty}-\frac{e^{-x(1+k)}(1+x+kx)}{(1+k)^2}$$
This series is convergent to a combinations of functions, one of them is the PolyLog function:
$$\sum_{k = 0}^{+\infty}-\frac{e^{-x(1+k)}(1+x+kx)}{(1+k)^2} = x\ln\left(e^{-x}(e^{x}-1)\right) - \text{PolyLog}[2, e^x]$$
Which is the solution you asked for.
$0$, $+\infty$ range
In this case, starting from the integration by part, you will get an easier result:
$$\int_0^{+\infty} xe^{-x(1+k)}\ \text{d}x = \frac{1}{(1+k)^2}$$
Hence the series would be
$$\sum_{k = 0}^{+\infty}\frac{1}{(1+k)^2} = \sum_{k = 1}^{+\infty}\frac{1}{k^2} = \zeta(2) = \frac{\pi^2}{6}$$
Indeed that series would be nothing than the Riemann Zeta Function evaluated at the point $2$.
