It is well known that an asymptotic expansion of the n-th harmonic number is $$H_{n}= \sum_{k=1}^{n} \frac{1}{k} \sim \ln(n) + \gamma + \frac{1}{2n} -\frac{1}{12n^{2}} + O(n^{-4}).$$
How could we find an asymptotic expansion for the sum $ \displaystyle \sum_{k=0}^n \frac{1}{1+\frac{k}{n}}$ to a similar order?
Let $ \displaystyle f(x) = \frac{1}{1+\frac{x}{n}}$.
Using the Euler-Maclaurin summation formula, we get
$$ \begin{align} \sum_{k=0}^{n} \frac{1}{1+\frac{k}{n}} &\sim \int_{0}^{n} \frac{1}{1+\frac{x}{n}} \, dx + \frac{f(n)+f(0)}{2} + \sum_{m=1}^{\infty} \frac{B_{2m}}{(2m)!} \left(f^{(2m-1)}(n) -f^{(2m-1)}(0) \right) \\ &= n \ln \left( 1+\frac{x}{n} \right) \Bigg|^{n}_{0} + \frac{\frac{1}{2}+1}{2} + \frac{1}{6} \left(\frac{1}{2!} \right) \left(-\frac{1}{n} \frac{1}{(1+\frac{x}{n})^{2}} \right) \Bigg|^{n}_{0} \\ &- \frac{1}{30} \frac{1}{4!} \left( - \frac{6}{n^{3}} \frac{1}{(1+ \frac{x}{n})^{4}}\right) \Bigg|^{n}_{0} + \mathcal{O}(n^{-5}) \\ &= n \ln (2) + \frac{3}{4} + \frac{1}{16n} - \frac{1}{128 n^{3}} + \mathcal{O}(n^{-5}) \end{align}$$
For $n=20$, the above approximation is correct to $8$ digits after the decimal point.
This sum is quite intriguing and invites additional scrutiny. We will now present an alternate derivation of the asymptotic expansion which is totally unorthodox and employs divergent series, yet produces the correct result.
We are studying $$S_n = \sum_{k=0}^n \frac{1}{1+\frac{k}{n}}.$$
Rewrite this as $$S_n = \sum_{k=0}^n \sum_{q\ge 0} (-1)^q \left(\frac{k}{n}\right)^q = \sum_{q\ge 0} \frac{(-1)^q}{n^q} \sum_{k=0}^n k^q.$$
The term for $q=0$ is special and we now extract it from the sum to obtain $$n+1+\sum_{q\ge 1} \frac{(-1)^q}{n^q} \sum_{k=1}^n k^q.$$
We will now work with the remaining sum and apply Faulhaber's formula which we recall says that $$\sum_{k=1}^n k^q = \frac{1}{q+1} \sum_{j=0}^q (-1)^j {q+1\choose j} B_j n^{q+1-j}.$$ Substitution now yields (note that $q+1 > j$ so $q \ge j$) $$n+1 + \sum_{j\ge 0} \sum_{q\ge j} \frac{(-1)^q}{n^q} (-1)^j {q+1\choose j} \times \frac{1}{q+1} \times B_j n^{q+1-j} \\= n+1 + \sum_{j\ge 0} (-1)^j B_j n^{1-j} \sum_{q\ge j} \frac{(-1)^q}{q+1} {q+1\choose j}.$$
There is one more correction to make: we included the term for $q=0$ and $j=0$ which has value $n$ and we need to subtract out this term because we already accounted for it, which finally leaves us with $$ 1 + \sum_{j\ge 0} (-1)^j B_j n^{1-j} \sum_{q\ge j} \frac{(-1)^q}{q+1} {q+1\choose j}.$$
We have now arrived at the moment where we will employ divergent series. The inner sum does not converge but we can assign a value to it evaluating $$Q_j(z) = \sum_{q\ge j} \frac{z^q}{q+1} {q+1\choose j}.$$ and setting $z=-1.$ Note that $$Q_0(z) = \sum_{q\ge 0} \frac{z^q}{q+1} = \frac{1}{z} \sum_{q\ge 0} \frac{z^{q+1}}{q+1} = \frac{1}{z} \log\frac{1}{1-z}.$$ For $j>0$ we get $$Q_j(z) = -\frac{z^{j-1}}{j} + \sum_{q\ge j-1} \frac{z^q}{q+1} {q+1\choose j} = -\frac{z^{j-1}}{j} + \sum_{q\ge j} \frac{z^{q-1}}{q} {q\choose j} \\= -\frac{z^{j-1}}{j} + \sum_{q\ge j} z^{q-1} \frac{(q-1)!}{j! (q-j)!} = -\frac{z^{j-1}}{j} + \frac{1}{j} \sum_{q\ge j} z^{q-1} \frac{(q-1)!}{(j-1)! (q-j)!} \\ = -\frac{z^{j-1}}{j} + \frac{1}{j} \sum_{q\ge j} z^{q-1} {q-1\choose j-1} = -\frac{z^{j-1}}{j} + \frac{1}{j} \sum_{q\ge 0} z^{q+j-1} {q+j-1\choose j-1}.$$ This sum is the Newton binomial and we finally have $$-\frac{z^{j-1}}{j} + \frac{z^{j-1}}{j} \sum_{q\ge 0} z^q {q+j-1\choose j-1} = -\frac{z^{j-1}}{j} + \frac{z^{j-1}}{j} \times \frac{1}{(1-z)^j}.$$
We will take $Q_j(-1)$ as the value of the inner sum that determines the coefficient on $n^{1-j}.$ We get for $j=0$ the value $$Q_0(-1) = \frac{1}{-1} \log \frac{1}{2} = \log 2.$$ For $j\ge 1$ we obtain $$Q_j(-1) = -\frac{(-1)^{j-1}}{j} + \frac{(-1)^{j-1}}{j} \frac{1}{2^j} = \frac{(-1)^{j-1}}{j} \frac{1-2^j}{2^j}.$$
Putting it all together we obtain the asymptotic expansion $$1 + n B_0 \log 2 - B_1 \times \frac{-1}{2} + \sum_{j\ge 2} (-1)^j B_j n^{1-j} \sum_{q\ge j} \frac{(-1)^q}{q+1} {q+1\choose j} \\= n \log 2 + \frac{3}{4} + \sum_{j\ge 2} (-1)^j B_j n^{1-j} \frac{(-1)^{j-1}}{j} \frac{1-2^j}{2^j} \\ = n \log 2 + \frac{3}{4} + \sum_{j\ge 2} B_j \frac{1}{j} \frac{2^j-1}{2^j} n^{1-j}.$$
Actually computing this expansion we get $$n\log 2 +3/4+1/16\,{n}^{-1}-{\frac {1}{128}}\,{ n}^{-3}+{\frac {1}{256}}\,{n}^{-5}-{\frac {17}{4096}}\,{n}^{-7 }+{\frac {31}{4096}}\,{n}^{-9}\\-{\frac {691}{32768}}\,{n}^{-11} +{\frac {5461}{65536}}\,{n}^{-13}-{\frac {929569}{2097152}}\,{ n}^{-15}+{\frac {3202291}{1048576}}\,{n}^{-17}+\ldots$$
This expansion is amazingly precise. Using up to the $n^{-17}$ term we get for $n=20$ the approximation $$14.6160676358538815006576381306$$ while the exact value is $$14.6160676358538815006576332145.$$
Using the terms up to $n^{-23}$ we get for $n=50$ the approximation $$35.4086089655097601622322941347$$ while the exact value is $$35.4086089655097601622322941347$$ i.e. at thirty digits precision we can no longer tell the difference.
I do think this calculation is remarkable in that it makes extensive use of divergent series to obtain a result that appears to be perfectly correct.
$\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\dd}{{\rm d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\dsc}[1]{\displaystyle{\color{red}{#1}}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\imp}{\Longrightarrow} \newcommand{\Li}[1]{\,{\rm Li}_{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$ Since the series has a closed expression in terms of the Digamma Functions, it's interesting to check the expansion by means of the asymptotic behavior of those functions: \begin{align}&\color{#66f}{\large\sum_{k=0}^{n}\frac{1}{1 + k/n}} =n\sum_{k=0}^{n}\frac{1}{k + n} =n\sum_{k=0}^{\infty}\pars{\frac{1}{k + n} - \frac{1}{k + 2n + 1}} \\[5mm]&=n\bracks{\Psi\pars{2n + 1} - \Psi\pars{n}}\tag{1} \end{align} where $\ds{\Psi}$ is the Digamma Function and we used ${\bf 6.3.16}$. With ${\bf 6.3.5}$, expression $\pars{1}$ becomes: \begin{align}&\color{#66f}{\large\sum_{k=0}^{n}\frac{1}{1 + k/n}} =n\bracks{\Psi\pars{2n} - \Psi\pars{n}} + \half \end{align} With the Digamma Asymptotic Expansion ${\bf 6.3.18}$: \begin{align}&\color{#66f}{\large\sum_{k=0}^{n}\frac{1}{1 + k/n}} \sim n\times \\[5mm]&\braces{\!\!\bracks{% \ln\pars{2n} - \frac{1}{4n} - \frac{1}{48n^{2}} + \frac{1}{1920n^{4}} -\frac{1}{16128n^{6}}}\!\! -\!\!\bracks{% \ln\pars{n} - \frac{1}{2n} - \frac{1}{12n^{2}} + \frac{1}{120n^{4}} -\frac{1}{252n^{6}}}} \\[5mm]&-\half \\[1cm]&=n\bracks{\ln\pars{2} + \frac{1}{4n} + \frac{1}{16n^{2}} - \frac{1}{128n^{4}} + \frac{1}{256n^{6}}} + \half \\[5mm]&=\color{#66f}{\large n\ln\pars{2} + \frac{3}{4} + \frac{1}{16n} - \frac{1}{128n^{3}} + \frac{1}{256n^{5}}} \end{align}
