Deriving $\gamma \approx H(n)-\ln(n+1)+\frac{1}{2(n+1)}+\frac{1}{12(n+1)^2}$ The Euler-Mascheroni constant can be represented geometrically by the infinite sum of the areas in blue in the following picture, which is the area between the curve $y=1/x$ and the harmonic numbers.

Thus, the total area can be approximated by taking a finite sum of the first $n$ areas, such that
$$\gamma \approx H(n)-\ln(n+1)$$
This approximation can be improved by noting that, as $n$ increases, these areas approach a triangle with base $\frac{1}{n}-\frac{1}{n+1}$ and height $1$. Therefore, the sum of the remaining areas is
$$A \approx\sum_{k=n+1}^{\infty}\frac{1}{2} \left(\frac{1}{n}-\frac{1}{n+1} \right)=\frac{1}{2(n+1)}$$
Thus,
$$\gamma \approx H(n)-\ln(n+1)+\frac{1}{2(n+1)}$$
Is there also a relatively simple way to derive the following even better approximation?
$$\gamma \approx H(n)-\ln(n+1)+\frac{1}{2(n+1)}+\frac{1}{12(n+1)^2}$$
I saw an almost identical formula in the harmonic numbers Wolfram MathWorld page but it seems that it comes from the rather complex Euler-Maclaurin formula.
 A: Consider one of the intervals from $n$ to $n+1$

After your $\frac{1}{2(n+1)}$ approximation, there is still some area between the $y=1/x$ curve and the hypotenuse of each triangle.
It's possible to approximate the function $y=1/x$ locally by a parabola with the Taylor series at $x=n+1/2$. Let this approximation be $f(x)$, then
$$f(x)=\frac{1}{n+1/2}-\frac{(x-n-1/2)}{(n+1/2)^2}+\frac{(x-n-1/2)^2}{(n+1/2)^3}$$
The hypotenuse of the triangle, represented by the segment AE in the image above, is given by
$$g(x)=\frac{1}{n}+(x-n)\left(\frac{1}{n+1}-\frac{1}{n}\right)$$
Therefore, the missing area is approximately
$$A\approx\int_{n}^{n+1}g(x)-f(x)\text{ }dx=\frac{8n^2+8n+3}{6n(n+1)(2n+1)^3}$$
Note that $(2n+1)^2=4n^2+4n+1$, and the numerator is $2(4n^2+4n+1.5)$. So a new approximation is
$$A\approx\frac{1}{3n(n+1)(2n+1)}=\frac{1}{6n(n+1)(n+1/2)}$$
So the total missing area would be the sum of $\frac{1}{6k(k+1)(k+1/2)}$ from $k=n+1$ to infinity. However, this does not have a "nice" solution due to the $1/2$ term, which is not an integer. Even so, we can get a lower and upper bound by considering the following twos series:
$$S_1=\sum_{k=n+1}^{\infty}\frac{1}{6k(k+1)(k+2)}=\frac{1}{12(n+1)(n+2)}$$
$$S_2=\sum_{k=n+1}^{\infty}\frac{1}{6k(k+1)(k-1)}=\frac{1}{12n(n+1)}$$
In which the $(n+1/2)$ term was changed to $(n+2)$ and $(n-1)$, respectively (-1 and 2 are the closest integers to $1/2$ that are not 0 or 1).
Now, it is interesting to consider the harmonic mean of both results. The reason is that the harmonic mean of $\frac{1}{6k(k+1)(k+2)}$ and $\frac{1}{6k(k+1)(k-1)}$ is exactly $\frac{1}{6k(k+1)(k+1/2)}$, which was the expression we were trying to use initially. So a better approximation for $A$ is to consider the harmonic mean of $S_1$ and $S_2$.
$$A\approx \frac{1}{12(n+1)^2}$$
Therefore,
$$\gamma \approx H(n)-\ln(n+1)+\frac{1}{2(n+1)}+\frac{1}{12(n+1)^2}$$
is a better approximation.
