Calculating the asymptotic of integrals For an integral like
$$D_{n}(x) \equiv \int_{0}^{x} \frac{t^{n}}{e^{t}-1} d t$$
The asymptotic values are given as
$$D_{n}(x) \simeq\left\{\begin{array}{ll}
n ! \zeta(n+1)-x^{n} e^{-x}+O\left(x^{n} e^{-2 x}\right), & x \rightarrow \infty \\
x^{n} / n-x^{n+1} / 2(n+1)+O\left(x^{n+2}\right), & x \rightarrow 0
\end{array}\right.$$
For $x\to 0 $, I can just Taylor expand $\frac{1}{e^t-1}$ as $\frac{1}{t}-\frac{1}{2}+\frac{t}{12}-\frac{t^3}{720}+O\left(t^5\right)$ and do term by term integration to to get the given expression.
For $x\to \infty $ the first term is simply gotten by replacing the upper bound of the integral with $\infty$ but how are the correction terms calculated?
 A: The asymptotics you gave is not completely correct. I am going to derive the correct asymptotic form. First, we write
\begin{align*}
D_n (x) & = n!\zeta (n + 1) - \int_x^{ + \infty }  \frac{{t^n }}{{e^t  - 1}}dt 
 \\ & \mathop  = \limits^{t = x(1 + s)} n!\zeta (n + 1) - x^{n + 1} e^{ - x}  \int_0^{ + \infty }  e^{ - xs} \frac{(1 + s)^n }{1 - e^{ - x(1 + s)}}ds .
\end{align*}
Now, by the geometric series,
$$
\int_0^{ + \infty } {e^{ - xs} \frac{{(1 + s)^n }}{{1 - e^{ - x(1 + s)} }}ds}  = \int_0^{ + \infty } {e^{ - xs} (1 + s)^n ds}  + \int_0^{ + \infty } {e^{ - xs} (1 + s)^n \mathcal{O}(e^{ - x(1 + s)} )ds} .
$$
By Watson's lemma,
$$
\int_0^{ + \infty } {e^{ - xs} (1 + s)^n ds}  = \int_0^{ + \infty } {e^{ - xs} (1 + s)^n ds}  = \frac{1}{x} + \mathcal{O}\left( {\frac{1}{{x^2 }}} \right)
$$
and
$$
\int_0^{ + \infty } {e^{ - xs} (1 + s)^n \mathcal{O}(e^{ - x(1 + s)} )ds}  = \mathcal{O}(1)e^{ - x} \int_0^{ + \infty } {e^{ - 2xs} (1 + s)^n ds}  = \mathcal{O}\left( {\frac{{e^{ - x} }}{x}} \right).
$$
Thus,
\begin{align*}
D_n (x) & = n!\zeta (n + 1) - x^n e^{ - x}  + \mathcal{O}(x^{n - 1} e^{ - x} ) + \mathcal{O}(x^n e^{ - 2x} ) \\ & = n!\zeta (n + 1) - x^n e^{ - x}  + \mathcal{O}(x^{n - 1} e^{ - x} ).
\end{align*}
Addendum. It can be shown that for any positive integer $n$,
$$
D_n (x) = n!\zeta (n + 1) - n!x^n \sum\limits_{j = 0}^n {\frac{1}{{(n - j)!}}\frac{1}{x^j}\operatorname{Li}_{j + 1} (e^{ - x} )} 
$$
where $\operatorname{Li}_s$ is the polylogarithm.
A: Not sure, but it seems useful to make the substitution $xu = t$:
\begin{align}
D_n(x) & = \int_0^1 \frac{x^{n+1} u^n}{e^{xu} - 1} \ du \\
& = x^{n+1} \int_0^1 e^{-xu} \frac{u^n}{1 - e^{-xu}} \ du \\
& = x^{n+1} \int_0^1 e^{-xu} u^n \sum_{m=0}^\infty e^{-m x u} \ du \\
& = x^{n+1} \sum_{m=1}^\infty \int_0^1 u^n e^{-xu m} \ du. \\
\end{align}
The first term $\Gamma(n+1) - \Gamma(n+1,x)$ can be written in terms of asymptotic formulae well known for the Gamma function. Though this expansion seems to do better than the first two terms you have there...
A: As Daniel mentioned, using the fact that:
$$\frac{1}{e^t-1}=e^{-t}\frac 1{1-e^{-t}}=e^{-t}\sum_{n=0}^\infty e^{-nt}=\sum_{n=0}^\infty e^{-(n+1)t}=\sum_{n=1}^\infty e^{-nt}$$
and so:
$$D_n(x)=\int_0^x\sum_{k=1}^\infty t^ne^{-kt}dt=\sum_{k=1}^\infty\int_0^xt^ne^{-kt}dt=\sum_{k=1}^\infty\frac{\gamma(n+1,kx)}{k^{n+1}}$$
where $\gamma$ denotes the lower incomplete gamma function, so our series would be:
$$D_n(x)\approx\gamma(n+1,x)+\frac{\gamma(n+1,2x)}{2^{n+1}}+...$$
