Evaluating the following integral: $\int_{0}^{\infty }\frac{x\cos(ax)}{e^{bx}-1}\ dx$ How can we compute this integral for all $\operatorname{Re}(a)>0$ and $\operatorname{Re}(b)>0$?
$$\int_{0}^{\infty }\frac{x\cos(ax)}{e^{bx}-1}\ dx$$
Is there a way to compute it using methods from real analysis?
 A: Let
$$ I(a)=\int_{0}^{\infty}\frac{\sin(ax)}{e^{bx}-1}\ dx. $$
It is easy to see
$$ I'(a)=\int_{0}^{\infty}\frac{x\cos(ax)}{e^{bx}-1}\ dx. $$
Since
\begin{eqnarray}
I(a)&=&\int_{0}^{\infty}\frac{e^{-bx}\sin(ax)}{1-e^{-bx}}\ dx=\int_{0}^{\infty}\sum_{n=0}^\infty e^{-b(n+1)x}\sin(ax)\ dx\\
&=&\sum_{n=0}^\infty\int_{0}^{\infty} e^{-b(n+1)x}\sin(ax)\ dx=\sum_{n=0}^\infty\frac{a}{a^2+b^2(n+1)^2}.
\end{eqnarray}
Using the result from @projectilemotion, one has
\begin{eqnarray}
I(a)&=&\sum_{n=0}^\infty\frac{a}{a^2+b^2(n+1)^2}=\frac{-b+a\pi \coth(a\pi/b)}{2a b}=\frac12\left(-\frac{1}{a}+\frac\pi b\coth(\frac{a\pi}{b})\right)
\end{eqnarray}
and hence
$$ I'(a)=\frac12\left(\frac1{a^2}-\frac{\pi^2}{b^2}\text{csch}^2(\frac{a\pi}{b})\right).$$
A: We first start by converting our integral into a sum. We have that (by the geometric series):
$$\frac{1}{e^{bx}-1}=\sum_{n=1}^{\infty} e^{-bxn}$$
Additionally, we have that (as can be found by simply finding an antiderivative):
$$\int_0^{\infty} x\cos(ax)e^{-bxn}~dx=\frac{-a^2+b^2 n^2}{(a^2+b^2 n^2)^2}$$
Therefore, the integral in question is equal to the following sum:
$$\int_{0}^{\infty }\frac{x\cos(ax)}{e^{bx}-1}\ dx=\sum_{n=1}^{\infty} \frac{-a^2+b^2n^2}{(a^2+b^2 n^2)^2}=-a^2 S_1+b^2 S_2 \tag{1}$$
Therefore, it suffices to compute the following sums:
$$S_1:=\sum_{n=1}^{\infty} \frac{1}{(a^2+b^2 n^2)^2},\quad S_2:=\sum_{n=1}^{\infty} \frac{n^2}{(a^2+b^2 n^2)^2}$$
To do this, we start from the well-known result that (which can be derived using real analysis, as shown by several answers):
$$\sum_{n=0}^\infty\frac{1}{c^2+n^2}=\frac{1+c\pi\coth (c\pi)}{2c^2}$$
From this, by using the substitution $c:=a/b$ one can easily derive the more general result that (note the difference in the lower limit of the sum):
$$\sum_{n=1}^{\infty} \frac{1}{a^2+b^2 n^2}=\frac{-b+a\pi \coth(a\pi/b)}{2a^2 b} \tag{2}$$
Now, we can compute closed forms for the sums $S_1$ and $S_2$ by differentiating $(2)$ with respect to $a$ and $b$ respectively. We hence have that:
$$S_1=\frac{\pi^2 a^2 \operatorname{csch}^2(a\pi /b)+\pi a b \coth(a\pi /b)-2b^2}{4a^4 b^2}$$
$$S_2=\frac{-\pi^2 a^2 \operatorname{csch}^2(a\pi /b)+\pi a b \coth(a\pi /b)}{4a^2 b^4}$$
Therefore, by $(1)$, the integral is given by:
$$\bbox[5px,border:2px solid #C0A000]{\int_{0}^{\infty }\frac{x\cos(ax)}{e^{bx}-1}\ dx=\frac{1}{2} \left(\frac{1}{a^2}-\frac{\pi ^2    \operatorname{csch}^2(a\pi/b)}{b^2}\right)}$$
A: **my attempt **
$$I=\int_{0}^{\infty }\frac{x\ cos(ax)}{e^{bx}-1}dx=\sum_{n=1}^{\infty }\int_{0}^{\infty }x\ e^{-bnx}cos(ax)dx\\
\\
\\
=\frac{1}{2}\sum_{n=1}^{\infty }\int_{0}^{\infty }x\ e^{-bnx}\ (e^{iax}-e^{-iax})dx=\frac{1}{2}\sum_{n=1}^{\infty }[\int_{0}^{\infty }x\ e^{-(bn-ia)}dx+\int_{0}^{\infty }x\ e^{-(bn+ia)}dx]\\
\\
\\
=\frac{1}{2}\sum_{n=1}^{\infty }(\frac{\Gamma (2)}{(bn-ai)^2}+\frac{\Gamma (2)}{(bn+ia)^2})=\frac{1}{2b^2}\sum_{n=0}^{\infty }\frac{1}{(n-\frac{ai}{b})^2}+\frac{1}{2b^2}\sum_{n=0}^{\infty }\frac{1}{(n+\frac{ai}{b})^2}+\frac{1}{a^2}\\
\\$$
$$=\frac{1}{a^2}+\frac{1}{2b^2}(\Psi ^{1}(\frac{ai}{b})+\Psi ^{1}(\frac{-ai}{b}))\\\\\
\\
but\ we\ know\  \Psi  ^{(1)}(\frac{-ai}{b})=\Psi ^{(1)}(1-\frac{ai}{b})-\frac{b^2}{a^2}\\
\\
\\
\therefore \ I=\frac{1}{a^2}+\frac{1}{2b^2}\left ( \Psi ^{(1)}(1-\frac{ai}{b}) +\Psi ^{(1)}(\frac{ai}{b})-\frac{b^2}{a^2}\right )\\
\\
\\
by\ using\ the\ reflection\ formula\ :\ \Psi ^{(1)}(1-\frac{ai}{b})+\Psi ^{(1)}(\frac{ai}{b})=\frac{\pi ^2}{sin^2(\frac{i\pi a}{b})}\\
\\$$
so we have 
$$\therefore I=\frac{1}{2a^2}+\frac{1}{2b^2}\left ( \frac{-\pi ^2}{sinh^2(\frac{\pi a}{b})} \right )\\
\\
\\
=\frac{1}{2a^2}-\frac{\pi ^2}{2b^2sinh^2(\frac{\pi a}{b})}\ \ \ \ \ \ , b>0$$
note that : 
$$\frac{\pi ^2}{sin^2(\frac{i\pi a}{b})}=-\frac{\pi ^2}{sinh^2(\frac{\pi a}{b})}$$
