Using the Euler-Maclaurin formula to approximate Euler's constant, $\gamma := \lim_{n\to\infty}\left(-\ln n+\sum_{k=0}^n\frac1n\right)$

Question

Let $\gamma=\lim_{n\to\infty} F(n)$ where $$F(n)=1+\frac{1}{2}+\frac{1}{3}+\cdots\frac{1}{n}-\ln(n)$$ (This is Euler's constant.) How can I calculate $\gamma$ with $10$ digits of precision using the Euler-Maclaurin Formula?

Answer 1

The Euler-Maclaurin Formula says $$ \scriptsize\sum_{k=1}^n\frac1k=\log(n)+\gamma+\frac1{2n}-\frac1{12n^2}+\frac1{120n^4}-\frac1{252n^6}+\frac1{240n^8}-\frac1{132n^{10}}+R_{12}(n)\tag1 $$ where $R_{12}(n)\ge0$ and $|R_{12}(n)|\le\frac{691}{8192\,n^{12}}$.

Using $n=11$, we get $$ 0.57721566490151\le\gamma\le0.57721566490154 $$ Thus, to $13$ places we get $$ \gamma=0.5772156649015 $$

---

---

The Euler-Maclaurin Formula

A question was raised in comments regarding the Euler-Maclaurin Formula and where $\gamma$ comes from. $$ \begin{align} \sum_{k=1}^n\frac1k &=\int_{1^-}^{n^+}\frac1x\,\mathrm{d}\lfloor x\rfloor\\ &=\int_{1^-}^{n^+}\frac1x\,\mathrm{d}x-\int_{1^-}^{n^+}\frac1x\,\mathrm{d}\left(\{x\}-\tfrac12\right)\tag{2a}\\ &=\log(n)+\frac1{2n}+\frac12-\int_1^n\frac{\{x\}-\tfrac12}{x^2}\,\mathrm{d}x\tag{2b}\\ &=\log(n)+\frac1{2n}+\color{#C00}{\frac12-\int_1^\infty\frac{\{x\}-\tfrac12}{x^2}\,\mathrm{d}x}+\int_n^\infty\frac{\{x\}-\tfrac12}{x^2}\,\mathrm{d}x\tag{2c}\\ &=\log(n)+\color{#C00}{\gamma}+\frac1{2n}-\frac1{12n^2}+2\int_n^\infty\frac{\tfrac12\{x\}^2-\tfrac12\{x\}+\tfrac1{12}}{x^3}\,\mathrm{d}x\tag{2d}\\ &=\log(n)+\gamma+\frac1{2n}-\frac1{12n^2}+6\int_n^\infty\frac{\tfrac16\{x\}^3-\tfrac14\{x\}^2+\tfrac1{12}\{x\}}{x^4}\,\mathrm{d}x\tag{2e}\\ &=\log(n)+\gamma+\frac1{2n}-\frac1{12n^2}+\frac1{120n^4}\\[3pt] &\,+24\int_n^\infty\frac{\tfrac1{24}\{x\}^4-\tfrac1{12}\{x\}^3+\tfrac1{24}\{x\}^2-\tfrac1{720}}{x^5}\,\mathrm{d}x\tag{2f}\\ &\,\ \vdots\\ &=\log(n)+\gamma+\frac1{2n}-\frac1{12n^2}+\frac1{120n^4}-\frac1{252n^6}+\frac1{240n^8}-\frac1{132n^{10}}\\[3pt] &\,+\scriptsize11!\int_n^\infty\frac{6\{x\}^{11}-33\{x\}^{10}+55\{x\}^9-66\{x\}^7+66\{x\}^5-33\{x\}^3+5\{x\}}{239500800\,x^{12}}\,\mathrm{d}x\tag{2g} \end{align} $$ Explanation:
$\text{(2a)}$: $\lfloor x\rfloor=x-\{x\}$
$\text{(2b)}$: evaluate the first integral and integrate the second by parts
$\text{(2c)}$: $\int_1^n=\int_1^\infty-\int_n^\infty$
$\text{(2d)}$: define $\gamma=\frac12-\int_1^\infty\frac{\{x\}-\tfrac12}{x^2}\,\mathrm{d}x$
$\phantom{\text{(2d):}}$ integrate the unbounded integral by parts
$\text{(2e)}$: integrate by parts
$\text{(2f)}$: integrate by parts
$\text{(2g)}$: integrate by parts until we get the formula above

By integrating the tail integrals by parts, we eliminate all of the constant terms. All of the constants are accumulated in $\gamma$.

---

The Remainder

If we denote the degree $k$ polynomial of $\{x\}$ in the remainder integrals above as $p_k(\{x\})$, then the integration by parts gives $p_k=p_{k+1}'$ and we require $\int_0^1p_k(x)\,\mathrm{d}x=0$ so that $p_k(\{x\})$ is continuous for $k\gt1$. Since $p_1(x)=x-\frac12$, it can be shown by induction that $p_{2k-1}\!\left(x+\frac12\right)$ is odd and $p_{2k-1}\in P^-$ when $k$ is odd and $p_{2k-1}\in P^+$ when $k$ is even, where $$ P^{\pm}=\left\{p:\pm\left(\tfrac12-x\right)p(x)\ge0\text{ for }x\in[0,1]\right\}\tag3 $$ Thus, $\frac1{(2k-1)!}$ times the remainder is $$ \begin{align} \sum_{j=n}^\infty\int_0^1\frac{p_{2k-1}(t)}{(j+t)^{2k}}\,\mathrm{d}t &=\sum_{j=n}^\infty\int_0^1\overbrace{\ p_{2k-1}(t)\ \vphantom{\left(\frac1{j^2}\right)}}^{\substack{\in P^+\text{ if $k$ even}\\\in P^-\text{ if $k$ odd}}}\overbrace{\left(\frac1{(j+t)^{2k}}-\frac1{(j+1/2)^{2k}}\right)}^{\in P^+}\,\mathrm{d}t\tag{4a}\\ &\,\,\left\{\begin{array}{}\ge0\text{ if $k$ even}\\\le0\text{ if $k$ odd}\end{array}\right.\tag{4b} \end{align} $$ Furthermore, $$ \begin{align} \left|\sum_{j=n}^\infty\int_0^1\frac{p_{2k-1}(t)}{(j+t)^{2k}}\,\mathrm{d}t\right| &=\left|\sum_{j=n}^\infty\int_0^1p_{2k-1}(t)\left(\frac1{(j+t)^{2k}}-\frac1{(j+1)^{2k}}\right)\,\mathrm{d}t\right|\tag{5a}\\ &\le\sum_{j=n}^\infty\int_0^1|p_{2k-1}(t)|\left(\frac1{j^{2k}}-\frac1{(j+1)^{2k}}\right)\mathrm{d}t\tag{5b}\\ &=\frac1{n^{2k}}\int_0^1|p_{2k-1}(t)|\mathrm{d}t\tag{5c} \end{align} $$

---

Summary

The error is positive when $k$ is even and negative when $k$ is odd. This explains why the asymptotic expansion alternates.

Furthermore, the absolute value of the error is at most $$ \frac{(2k-1)!}{n^{2k}}\int_0^1|p_{2k-1}(t)|\,\mathrm{d}t\tag6 $$ Thus, the error term in $\text{(2g)}$ is positive ($k=6$) and at most $\frac{691}{8192\,n^{12}}$.