Compute $\sum \frac{1}{k^2}$ using Euler-Maclaurin formula

Question

I read that Euler used the summation formula to calculate the value of the series $\sum_{k =1}^{\infty} \frac{1}{k^2}$ to high precision without too much hassle. The article Dances between continuous and discrete: Euler’s summation formula goes into the calculation, however without too much justification of why it works (especially since the series used to calculate the limit does not converge and one has to truncate it at a certain point). I would be glad if someone could elaborate from a more modern viewpoint on how and why it works.

Answer 1

First, it can be shown, in many ways, that $$ \sum_{k=1}^\infty\frac1{k^2}=\frac{\pi^2}{6}\tag{1} $$ However, the Euler-Maclaurin Summation Formula can be used to numerically sum this series to high precision. Before we attempt to do this, let's discuss a bit about asymptotic series.

In general, asymptotic series, like those arising from the Euler-Maclaurin Summation Formula (EMSF), are divergent. This means that if you tried to sum all the terms arising from the formula, the sum would not converge. For example, the EMSF gives the following asymptotic expansion: $$ \sum_{k=1}^n\frac1k\sim\gamma+\log(n)+\frac1{2n}-\frac1{12n^2}+\frac1{120n^4}-\frac1{252n^6}+\frac1{240n^8}-\frac1{132n^{10}}+\dots\tag{2} $$ where $\gamma=0.5772156649015328606065121$ is the Euler-Mascheroni constant.

This expansion looks well behaved, and it is, up to a point. However, the coefficients grow on the order of $\frac{n!}{(2\pi)^n}$; by the time that we get to the term for $n^{50}$, we have $$ \frac{19802288209643185928499101}{132n^{50}}\tag{3} $$ Due to the growth rate of the coefficients, no matter how large $n$ is, this series cannot converge.

However, if we only use a finite number of terms, the series is a very good approximation for large $n$. As mentioned earlier, the expansion behaves well, up to a point. What this means is that for a given $n$, the terms get smaller, and then, at some point, start blowing up. The point at which the term start blowing up is further along the larger $n$ is. The good part is that if we terminate the sum while the terms are still getting smaller, the approximation is usually as good as the next term.

For example, let's approximate $$ \sum_{k=1}^{1000}\frac1k=7.4854708605503449126565182\tag{4} $$ using the first $4$ terms of $(2)$: $$ \gamma+\log(1000)+\frac1{2000}-\frac1{12000000}=7.4854708605503365793271531\tag{5} $$ The result in $(5)$ is $8.333329365\times10^{-15}$ smaller than the actual value in $(4)$. The next term is $$ \frac1{120000000000000}=8.333333333\times10^{-15}\tag{6} $$

---

Now let's see how the EMSF can be used to approximate $(1)$

The EMSF applied to $\dfrac1{k^2}$ yields $$ \sum_{k=1}^n\frac1{k^2}\sim C-\frac1n+\frac1{2n^2}-\frac1{6n^3}+\frac1{30n^5}-\frac1{42n^7}+\frac1{30n^9}-\frac5{66n^{11}}\tag{7} $$ Note that the EMSF always has a constant that needs to be determined in some other manner. In $(2)$ it was $\gamma$, the Euler-Mascheroni constant. Here, $C$ is the sum of the series; that is, as $n\to\infty$, the left side of $(7)$ tends to the sum, and everything on the right side of $(7)$, except $C$, tends to $0$.

To compute $C$, we will use $n=100$ and truncate the series at the $n^9$ term. The error we get should be less than $\dfrac5{66n^{11}}$, which would give us almost $23$ decimals places.

For $n=100$, the sum on the left of $(7)$ is $$ \sum_{k=1}^{100}\frac1{k^2}=1.6349839001848928650771695\tag{8} $$ For $n=100$, the sum of the terms on the right of $(7)$ other than $C$ is $$ -\frac1n+\frac1{2n^2}-\frac1{6n^3}+\frac1{30n^5}-\frac1{42n^7}+\frac1{30n^9} =-0.0099501666633335713952381\tag{9} $$ Subtracting $(9)$ from $(8)$ gives $$ C\stackrel{.}{=}1.6449340668482264364724076\tag{10} $$ whereas $$ \frac{\pi^2}{6}=1.6449340668482264364724152\tag{11} $$ The value of $(10)$ is $7.6\times10^{-24}$ short of $(11)$ and $\dfrac5{66n^{11}}=7.6\times10^{-24}$.

Using a larger $n$, and possibly more terms of the series, would give more precision.