Approximating ζ(2)

Recently, I was watching a series of lectures about Mathematical physics by Dr. Bender on Youtube (link), the series deal a lot with perturbative techniques that are applicable to different areas of applied math and physics. One of the major takeaways of the first lecture was

If you need to find the sum of a series, the worst possible technique is to add the numbers in the series together one at a time. It’s a terrible idea [and] there are so much better ways of summing a series.

- Dr. C. Bender, lecture 1 on MP, minute 43:00.

This reminded me of one of my favorite series:

$$ \zeta(2):=\sum_{n=1}^{\infty}\frac{1}{n^2}=\frac{\pi^2}{6} .$$

There are many articles on the internet (search for the Basel problem) about the history of the problem and showing different solutions for it. Euler’s first solution was very interesting, it involved a factorization of sin(x) that Euler could not justify. However, Euler was able to convince himself and the mathematical community that $ \frac{\pi^2}{6}$ is the right answer by approximating the value of the series up to 15 digits. The question now is

How did Euler manage to approximate the series with 15 digits of accuracy?

This takes me Back to Dr. Bender’s quote: What if we sum the first few terms of the series, that is, we compute

$$ S_m=\sum_{n=1}^m \frac{1}{n^2},$$

For a suitable $ m$. Now, we cannot afford $ m$ to be too big since computing $ S_m$ would be expensive. We also do not want to choose a value of $ m$ that is too small since that would give us a bad approximation. One good approach to choose a suitable value for $ m$ would be to find some bounds on $ |S_m-S|$ where $ S$ is the limit of $ S_m$. Using integrals, one could show that

$$ |S_m-S| = \frac{1}{m}+O(m^{-2})$$

That means that if one wishes to get 6 digits of accuracy, $ m$ has to be at least a million. This was not feasible in the 18th century so Euler had to somehow speed up the convergence of $ S_m$. The first idea was to add the first error term to $ S_m$, that is, we define a new sequence $ S^{(1)}_m$

$$ S^{(1)}_{m}=S_m+\frac{1}{m}$$

This actually works! In hindsight, we can compute the errors for $m=10$

$$ \left|\frac{\pi^2}{6}- S_{10}\right|\approx 9.5\times 10^{-2},\qquad \left|\frac{\pi^2}{6}- S_{10}^{(1)}\right|\approx 4.8\times 10^{-3}$$

We can see that $ S_{10}^{(1)}$ is more than ten times accurate than $ S_{10}$. This may have motivated Euler to obtain a formula for higher order terms and led him to Euler-Maclaurin formula which he established in the same year. When Euler-Maclaurin formula is applied to $ S_m$, we obtain a new sequence $ S^{(q)}_m$

$$ S^{(q)}_{m}=S_m +\sum_{k=0}^{q} \frac{B_{k}}{m^{k+1}},$$

where $ B_k$ are Bernoulli numbers. For example, $$ S^{(10)}_m= S_m+\frac{1}{m}-\frac{1}{2 m^2}+\frac{1}{6 m^3}-\frac{1}{30 m^5}+\frac{1}{42 m^7}-\frac{1}{30 m^9}$$ If we choose $ m=10$, we obtain an approximation that is correct up to 12 decimal places

$$ S^{(10)}_{10}\approx \mathbf{1.64493406684}{7493}$$

By changing $ q$, we could obtain more correct digits. We can also bound the error in terms of $ m$ and $ q$ using the Euler-Maclaurin remainder formulas.