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?

**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})}$$
