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

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.

-
See this question. – Alex Becker Jul 23 '12 at 18:59
That's the interesting thing about asymptotic series like the one that turns up in Euler-Maclaurin. They diverge, but if you cut them off at the right term, that partial sum is a good approximation... – J. M. Jul 23 '12 at 19:01

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.

-
Thanks for your answer. Something is still bothering me. Let's call $C = C^{(n)}$ since it depends on $n$. How did you justify that the error should be less than 5/(66 n^11)? In particular, eq. (7) holds no matter what term you add on the right hand side (as long as it vanishes for $n \to \infty$), since you can choose $C^{(n)}$ appropriately. What property of EMSF makes this whole thing work? – Jim Jul 25 '12 at 15:20
@Jim: Where did you get the idea that $C$ depends on $n$? Perhaps I used $\stackrel{.}{=}$ where I should have used $\sim$. I have adjusted the answer. The whole point of using these asymptotic expansions is that $C$ doesn't depend on $n$. They would be far less useful if they did. – robjohn Jul 25 '12 at 16:31

$\newcommand{\+}{^{\dagger}} \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{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\down}{\downarrow} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\isdiv}{\,\left.\right\vert\,} \newcommand{\ket}[1]{\left\vert #1\right\rangle} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert} \newcommand{\wt}[1]{\widetilde{#1}}$ Numerically, we must sum up a 'few' terms before we use Euler-Maclaurin ( in the second term ).

For example: \begin{align} \sum_{n = 1}^{\infty}{1 \over n^{2}} &=\sum_{n = 1}^{N}{1 \over n^{2}} + \sum_{n = 1}^{\infty}{1 \over \pars{n + N}^{2}} \approx\sum_{n = 1}^{N}{1 \over n^{2}} +\int_{0}^{\infty}{\dd x \over \pars{x + N}^{2}} =\sum_{n = 1}^{N}{1 \over n^{2}} + {1 \over N} \end{align}

This '$\ds{\tt\mbox{extremely simple formula}}$', for example, yields a relative error of $\ds{1.14\ \%}$ with $\ds{N = 5}$.

Historically; it seems Euler was convinced, by means of the Euler-Maclaurin formula, that the correct value was $\ds{\pi^{2} \over 6}$ before he tried to demonstrate it.

-