I'm looking for an elementary derivation of the following formula:
$$\int_0^\infty \frac{\sinh(ax)}{e^{bx}-1}dx=\frac{1}{2a}-\frac{\pi}{2b}\cot\frac{a\pi}{b}$$
...only where $|a|\lt b$, of course, to ensure convergence. Does anyone know of any elementary ways of proving this identity?
I'm fairly sure that a nice, elementary antiderivative does not exist, and I haven't been able to get anywhere with any of the "definite integral tricks" that I have learned thus far.
I give a real method that makes use of series and a bunch of special functions and their associated properties so perhaps this is not quite the elementary method you are after.
Let $$I = \int^\infty_0 \frac{\sinh (ax)}{e^{bx} - 1} \, dx, \quad |a| < b.$$ Rewriting the hyperbolic sine function term appearing in the numerator in terms of exponentials, since the improper integral converges, we may write \begin{align*} I &= \frac{1}{2} \int^\infty_0 \frac{e^{ax}}{e^{bx} - 1} \, dx - \frac{1}{2} \int^\infty_0 \frac{e^{-ax}}{e^{bx} - 1} \, dx\\ &= \frac{1}{2} \int^\infty_0 \frac{e^{-(b - a)x}}{1 - e^{-bx}} - \frac{1}{2} \int^\infty_0 \frac{e^{-(b + a) x}}{1 - e^{-bx}} \, dx. \end{align*}
Recognising the term $\dfrac{1}{1 - e^{-bx}}$ as the sum of a convergent geometric series, namely
$$\frac{1}{1 - e^{-bx}} = \sum^\infty_{n = 0} e^{-bnx}$$
the integral appearing above, after interchanging the summation with the integration, becomes
$$I = \frac{1}{2} \sum^\infty_{n = 0} \int^\infty_0 e^{-(b + nb - a)x} \, dx - \frac{1}{2} \sum^\infty_{n = 0} \int^\infty_0 e^{-(b + nb + a) x} \, dx.$$
Note that as $|a| < b$ the exponents appearing in each of the exponential terms found in each of the integrals will be negative ensuring convergence of the improper integrals. After performing the integrations we are left with
\begin{align*}
I &= \frac{1}{2} \sum^\infty_{n = 0} \frac{1}{(b - a) + nb} - \frac{1}{2} \sum^\infty_{n = 0} \frac{1}{(a + b) + nb}\\
&= \frac{1}{2b} \sum^\infty_{n = 0} \frac{1}{\left (\frac{b - a}{b} \right ) + n} - \frac{1}{2b} \sum^\infty_{n = 0} \frac{1}{\left (\frac{a + b}{b} \right ) + n}.
\end{align*}
Each of these sums can be expressed in terms of the Hurwitz zeta function $\zeta(s,q)$ which is defined as $$\zeta(s,q) = \sum^\infty_{n = 0} \frac{1}{(q + n)^s}.$$ Thus $$I = \frac{1}{2b} \zeta \left (1,1 - \frac{a}{b} \right ) - \frac{1}{2b} \zeta \left (1,1 + \frac{a}{b} \right ).$$
Now the Hurwitz zeta function is related to the polygamma function function of order $m$, which is denoted by $\psi^{(m)} (z)$, by (see here) $$\psi^{(m)} (z) = (-1)^{m + 1} m! \zeta (m + 1, z).$$ Setting $m = 0$ gives $$\zeta (1,z) = - \psi (z),$$ where $\psi (z) = \psi^{(0)} (z)$ is the digamma function. So in terms of digamma functions our integral can be written as $$I = \frac{1}{2b} \left [\psi \left (1 + \frac{a}{b} \right ) - \psi \left (1 - \frac{a}{b} \right ) \right ].$$
Since (see here) $$\psi (z + 1) = \psi (z) + \frac{1}{z},$$ we can write $$\psi \left (1 + \frac{a}{b} \right ) = \psi \left (\frac{a}{b} \right ) + \frac{b}{a},$$ and from the reflection formula for the digamma function, namely $$\psi (1 - z) - \psi (z) = \pi \cot (\pi z),$$ we have $$\psi \left (1 - \frac{a}{b} \right ) = \psi \left (\frac{a}{b} \right ) + \pi \cot \left (\frac{a \pi}{b} \right ),$$ yielding $$\int^\infty_0 \frac{\sinh (ax)}{e^{bx} - 1} \, dx = \frac{1}{2a} - \frac{\pi}{2b} \cot \left (\frac{a \pi}{b} \right ),$$ as claimed.
$$\int_0^\infty \frac{e^{ax}-e^{-ax}}{e^{bx}-1}dx = \sum_{n=1}^\infty \int_0^\infty (e^{-(nb-a)x}-e^{-(nb+a)x})dx\\ = \sum_{n=1}^\infty \frac{1}{nb-a}-\frac1{nb+a}= \frac{1}{2a} - \frac{\pi}{2b} \cot \left (\frac{a \pi}{b} \right )$$
Where the last identity is a consequence of $$\pi \cot(\pi z) - (\frac{1}{z}+\sum_{n=1}^\infty \frac{1}{z-n}+\frac{1}{z+n}) \tag{A}$$ is an everywhere analytic function, $1$-periodic and bounded on $\Re(z) \in [0,2]$ and vanishing at $z=0$.
If you want to avoid Liouville theorem in complex analysis, there are some "elementary proofs" of $(A) = 0$, eg. there, which is historically how Euler proved the product formula for $\sin(z)$ and solved the Basel problem.
Here the $b$ parameter is irrelevant, since it can be removed through a suitable substitution.
It is enough to show that
$$ \int_{0}^{+\infty}\frac{\sinh(ax)}{e^x-1}\,dx = \sum_{n\geq 1}\int_{0}^{+\infty}\sinh(ax)e^{-nx}\,dx = \sum_{n\geq 1}\frac{a}{n^2-a^2}\color{red}{=}\frac{1-\pi a \cot(\pi a)}{2a}$$
holds for any $|a|<1$.
This can be done by considering the Weierstrass product for the sine function
$$ \frac{\sin(\pi z)}{\pi z}=\prod_{n\geq 1}\left(1-\frac{z^2}{n^2}\right)$$
and by applying $-\frac{d}{dz}\log(\cdot)$ to both sides, leading to:
$$ \frac{1}{z}-\pi\cot(\pi z) = \sum_{n\geq 1}\frac{2z}{n^2-z^2}$$
where the evaluation at $z=a$ finishes the job.
For an approach through Fourier series, have a look at page 43 of my notes.
$\newcommand{\bbx}[1]{\,\bbox[15px,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}$
$\ds{\int_{0}^{\infty}{\sinh\pars{ax} \over \expo{bx}-1}\,\dd x = {1 \over 2a} - {\pi \over 2b}\cot\pars{a\pi \over b}:\ {\large ?}}$.
\begin{align} &\bbox[15px,#ffd]{\ds{\int_{0}^{\infty}{\sinh\pars{ax} \over \expo{bx}-1} \,\dd x}} = \int_{0}^{\infty}{\pars{\expo{ax} - \expo{-ax}}/2 \over 1 - \expo{-bx}} \,\expo{-bx}\,\dd x \\[5mm] \stackrel{{\Large t\ =\ \expo{-bx}} \atop {\Large x\ =\ -\ln\pars{t}/b}}{=}\,\,\,& {1 \over 2}\int_{1}^{0}{t^{-a/b} - t^{a/b} \over 1 - t}\,t \,\pars{-\,{\dd t \over bt}} = {1 \over 2b}\bracks{-\int_{0}^{1}{1 - t^{-a/b} \over 1 - t}\,\dd t + \int_{0}^{1}{1 - t^{a/b} \over 1 - t}\,\dd t} \\[5mm] = &\ {1 \over 2b}\bracks{-\Psi\pars{-\,{a \over b} + 1} + \Psi\pars{{a \over b} + 1}} \qquad\pars{~\Psi:\ Digamma\ Function~} \\[5mm] = &\ -\,{1 \over 2b}\bracks{% \Psi\pars{1 - {a \over b}} - \Psi\pars{a \over b} - {1 \over a/b}} \qquad\pars{~\Psi\ Recursive\ Property~} \\[5mm] & = {1 \over 2a} -\,{1 \over 2b}\ \underbrace{\bracks{% \Psi\pars{1 - {a \over b}} - \Psi\pars{a \over b}}} _{\ds{\pi\cot\pars{\pi\,{a \over b}}}}\quad \pars{~Euler\ Reflection\ Formula\ \mbox{for}\ \Psi~} \\[5mm] & = \bbx{{1 \over 2a} - {\pi \over 2b}\cot\pars{\pi a \over b}} \end{align}
