how to find the asymptotic expansion of the following sum: I need to determine an asymptotic expansion when $q \rightarrow 1$ of the sum
$$S(q)=\sum_{n=0}^{\infty} \frac{q^n}{ (q^n + 1)^2 }.$$
Numerical computations suggest that $S(q)\sim\frac{c}{|q-1|}$ with $c \approx 0.5$.
 A: It is indeed true that $S(q) \sim \frac{1/2}{|1-q|}$.  We'll use the same method as in this answer.
First assume $0 < q < 1$.  The terms of the sum are strictly decreasing in $n$, so
$$
\int_0^\infty \frac{q^x}{(q^x + 1)^2}\,dx \leq \sum_{n=0}^{\infty} \frac{q^n}{(q^n+1)^2} \leq \frac{1}{4} + \int_0^\infty \frac{q^x}{(q^x + 1)^2}\,dx.
$$
To evaluate the integral, make the change of variables $q^x = y$ to get
$$
\begin{align}
\int_0^\infty \frac{q^x}{(q^x + 1)^2}\,dx &= -\frac{1}{\log q} \int_0^1 \frac{dy}{(y+1)^2} \\
&= -\frac{1}{2\log q}.
\end{align}
$$
Thus
$$
\sum_{n=0}^{\infty} \frac{q^n}{(q^n+1)^2} \sim -\frac{1}{2\log q}
$$
as $q \to 1^-$.  Of course $\log q \sim q-1$ as $q \to 1$, so this is equivalent to
$$
\sum_{n=0}^{\infty} \frac{q^n}{(q^n+1)^2} \sim \frac{1}{2(1-q)}
$$
as $q \to 1^-$.
To address the case when $q > 1$ we note that $S(q) = S(1/q)$, so that
$$
S(q) \sim \frac{1}{2(1-\frac{1}{q})} = \frac{q}{2(q-1)} \sim \frac{1}{2(q-1)}
$$
as $q \to 1^+$.  Thus

$$
S(q) \sim \frac{1}{2|1-q|}
$$
  as $q \to 1$.

A: Suppose that $q>1$ and re-write the sum as
$$S(q) = \frac{1}{4} + \sum_{n\ge 1} \frac{q^n}{(q^n+1)^2}$$
The sum term is harmonic and may be re-written as
$$T(x) = \sum_{n\ge 1} \frac{e^{nx}}{(e^{nx}+1)^2}$$
which we will evaluate at $x=\log q.$
We will evaluate $T(x)$ by inverting its Mellin transform. 
Recall the harmonic sum identity
$$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$
where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have
$$\lambda_k = 1, \quad \mu_k = k \quad \text{and} \quad
g(x) = \frac{e^x}{(e^x+1)^2}.$$
We need the Mellin transform $g^*(s)$ of $g(x)$ which is
$$\int_0^\infty \frac{e^x}{(e^x+1)^2} x^{s-1} dx
= \int_0^\infty \frac{e^{-x}}{(1+e^{-x})^2} x^{s-1} dx
= \int_0^\infty \sum_{m\ge 1} (-1)^{m+1} m e^{-mx} x^{s-1} dx\\
\sum_{m\ge 1} (-1)^{m+1} m \int_0^\infty  e^{-mx} x^{s-1}  dx
= \Gamma(s) \sum_{m\ge 1} \frac{(-1)^{m+1} m}{m^s}
= \Gamma(s)  \sum_{m\ge 1} \frac{(-1)^{m+1}}{m^{s-1}}\\
= \Gamma(s)
 \left(1-\frac{1}{2^{s-2}}\right) \zeta(s-1).$$
Therefore the Mellin transform $Q(s)$ of $T(x)$ is given by
$$Q(s) = \Gamma(s)
 \left(1-\frac{1}{2^{s-2}}\right) \zeta(s-1) \zeta(s).$$
The Mellin inversion integral for $Q(s)$ is
$$\frac{1}{2\pi i}\int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds.$$
The two zeta function terms taken together cancel the poles of the gamma function, so we are left with just two poles/residues:
$$\mathrm{Res}(Q(s)/x^s; s=1) = \frac{1}{2x}
\quad\text{and}\quad
\mathrm{Res}(Q(s)/x^s; s=0) = -\frac{1}{8}.$$
Therefore in a neighborhood of zero we have that
$$T(x) \sim \frac{1}{2x}-\frac{1}{8}.$$
Now as $q$ goes to one $\log q$ goes to zero, so this expansion at $x=\log q$ definitely applies. We get that $$S(q) = \frac{1}{4} + \frac{1}{2x}-\frac{1}{8}
= \frac{1}{8} + \frac{1}{2\log q}.$$
This approximation is better than the conjecture from the original question. 
Now use that for $q$ close to one,
$$\frac{1}{2\log q} \sim \frac{1}{2(q-1)} + \frac{1}{4} - \frac{1}{24} (q-1)+\cdots$$
to conclude that
$$S(q) \sim \frac{3}{8} + \frac{1}{2(q-1)}.$$
