This formula has a splendid history. It was derived in the West in 1671 by James Gregory from the formula for arctan(x) and slightly later and independently by Gottfried Leibniz. However, the same formula (along with many other results involving infinite series) was discovered long before in the 1300s by the great Indian mathematician Madhava. Similar results of equal beauty are the convergent series given by

$$\frac{{\pi }^2}{6}=1+\frac{1}{2^2}+\frac{1}{3^2}+\frac{1}{4^2}+\dots ,\text{(2)}$$ $$\frac{{\pi }^3}{32}=1-\frac{1}{3^3}+\frac{1}{5^3}-\frac{1}{7^3}+\dots ,\text{(3)}$$ and $$\frac{{\pi }^4}{90}=1+\frac{1}{2^4}+\frac{1}{3^4}+\frac{1}{4^4}+\dots \text{(4)}$$ However, in one sense, these formulae are disappointing. If you want to actually calculate $\pi ,{\pi }^2,{\pi }^3$ or ${\pi }^4$ then you probably would not reach for one of these formulae. The reason is that they converge very slowly. If you take formula (1), add up 100 terms and multiply by 4, you get 3.146567747182956, which whilst fairly close to $\pi =3.141592653589793...$ is not a particularly accurate estimate given the effort involved in adding up 100 terms. If you wanted to calculate $\pi$ or ${\pi }^2$ to an accuracy of six decimal places, you would have to take on the order of ${10}^6$ terms of either series (1) or (2), and long before you have added up all of the terms in the series, the rounding errors associated with computer calculations will have accumulated to the point where the accuracy of the answer is severely degraded.

The world needs pi

So what, you (and many pure mathematicians) might say. Surely you don't need to know the value of $\pi$ that accurately, after all the Bible was content to give it to just one significant figure. However, $\pi$ is not any number. It lies at the heart of any technology that involves rotation or waves, and that is much of mechanical and electrical engineering. If rotating parts in, say, a typical jet engine are not manufactured to high tolerance, then the parts simply won't rotate. This typically involves measurements correct to one part in ${10}^4$ and, as these measurements involve $\pi$, we require a value of $\pi$ to at least this order of accuracy to prevent errors. In medical imaging using CAT or MRI scanners, the scanning devices move on a ring which has to be manufactured to a tolerance of one part in ${10}^6$, requiring an even more precise value of $\pi$. However, even this level of accuracy pales into insignificance when we look at modern electrical devices. In high frequency electronics, with frequencies in the order of 1GHz (typical for mobile phones or GPS applications), electrical engineers have to work with functions of the form $u(t)=\cos (2\pi ft)$ where $f\sim {10}^9$ and $t$ is a number close to one. To get the accuracy in the function $u(t)$ needed for GPS to work requires a precision in the value used for $\pi$ in the order of one part in ${10}^{15}$. So, to live in the modern world we really do need to know $\pi$ very accurately. So, what can we do? One possibility is to take a vast number of terms of the series for $\pi$ etc. above, book lots of time on a very expensive computer, sit back and wait (and wait, and wait). Or we can try and accelerate their convergence. So that with only a small number of terms (say 10) we can get 10 significant figures for $\pi$. The nice thing about this method is that the derivation of the formulae is very transparent (well within the reach of a first year undergraduate or even a good A-level student). In principle this method can also be used to find the sum of other slowly convergent series.

Accelerating the convergence of a series

Let's suppose that we have a series $$a_1+a_2+a_3+a_4+\dots +a_n+\dots$$ and we define the sum $S_n$ by $$S_n=a_1+a_2+\dots +a_n.$$ We will assume that this series {\em converges}. This means there is a {\em limiting sum} $S$ so that $$S_n\to S\text{as}n\to \mathrm{\infty }.$$ If we want to work out the value of $S$ then we can simply take the values $S_n$ and let $n$ get very large. However, just adding up a series loses quite a lot of the information contained within it, and is really a very crude thing to do. Maybe we can squeeze more information out of the series and use this information to accelerate the convergence of $S_n$. This means that we add a correction term to $S_n$ so that it approaches $S$ much more rapidly. What is nice about this approach is that it is quite easy to work out the correction terms. To illustrate this idea we will take series (2) for ${\pi }^2/6$. Let's define $$S_n=1+\frac{1}{2^2}+\dots +\frac{1}{n^2}.\text{(5)}$$ As $n$ increases, $S_n$ also increases steadily towards the value of $S={\pi }^2/6$. In figure 1 we plot the values of $S_n$ which you can see are increasing towards the value of ${\pi }^2/6=1.644934066848226$. Suppose that we next look at the difference $E_n=S-S_n$. A plot of is given in figure 2(a) and we can see that this decreases to zero as $n$ increases. But how fast? To estimate this we plot in figure 2(b) the values of $n{E}_n$. It is clear from these figures that while $E\to 0$ as $n\to \mathrm{\infty }$, $n{E}_n$ approaches a constant that looks suspiciously like one. From this, we might guess that $E_n\approx 1/n$ so that to a first approximation $$S=S_n+\frac{1}{n}+\dots$$ We can improve on this guess by assuming that there is a sequence of numbers $B_0,B_1,B_2,B_3,\dots$ such that as $n\to \mathrm{\infty }$ $$S=S_n+\frac{B_0}{n}+\frac{B_1}{n^2}+\frac{B_2}{n^3}+\frac{B_3}{n^4}+\dots \text{(6)}$$ If we can calculate the terms $B_r$ for $r=0,1,2,\dots k$ then we can estimate $S$ from the value of $S_n$ by the expression $$S\approx S_n+\frac{B_0}{n}+\frac{B_1}{n^2}+\frac{B_2}{n^3}+\frac{B_3}{n^4}+\dots +\frac{B_{k-1}}{n^k}.$$

But what are the B_r?

It turns out that the B_r are a very well-known set of numbers. There is a nice and systematic way to calculate their values (but if you want you can skip this calculation and go straight to the result).

If we take (6) and replace $n$ by $n-1$ then we get $$S=S_{n-1}+\frac{B_0}{n-1}+\frac{B_1}{(n-1{)}^2}+\frac{B_2}{(n-1{)}^3}+\frac{B_3}{(n-1{)}^4}+\dots \text{(7)}$$ so that $$S_{n-1}-S_n=B_0(\frac{1}{n}-\frac{1}{(n-1)})+B_1(\frac{1}{n^2}-\frac{1}{(n-1{)}^2})+B_2(\frac{1}{n^3}-\frac{1}{(n-1{)}^3})+\dots \text{(8)}$$ Also, from the definition of $S_n$ we know that $$S_n-S_{n-1}=\frac{1}{n^2}.$$ Combining these two expressions we get $$-\frac{1}{n^2}=B_0(\frac{1}{n}-\frac{1}{(n-1)})+B_1(\frac{1}{n^2}-\frac{1}{(n-1{)}^2})+B_2(\frac{1}{n^3}-\frac{1}{(n-1{)}^3})+\dots \text{(9)}$$ We can now find the values of each of the terms $B_k$ recursively by expanding each of these expressions in powers of $1/n$ and considering $n$ to be very large. For example $$\frac{1}{(n-1)}=\frac{1}{n(1-1/n)}=\frac{1}{n}+\frac{1}{n^2}+\frac{1}{n^3}+\dots ,$$ and $$\frac{1}{(n-1{)}^2}=\frac{1}{n^2(1-1/n{)}^2}=\frac{1}{n^2}+\frac{2}{n^3}+\frac{3}{n^4}+\frac{4}{n^5}+\dots .$$ More generally $$\frac{1}{(n-1{)}^k}=\frac{1}{n^k(1-1/n{)}^k}=\frac{1}{n^k}+\frac{k}{n^{k+1}}+\dots +\frac{\left(\begin{array}{c}k+r-1\\ r\end{array}\right)}{n^{k+r}}+\dots$$ where $$\left(\begin{array}{c}k+r-1\\ r\end{array}\right)=\frac{k(k+1)(k+2)\dots (k+r-1)}{r!}.$$

These are the so-called Taylor series expansions for the above expressions.

Now, we can combine these expressions with the equation (9). This gives $$\begin{array}{rcl}\frac{1}{n^2}& =& B_0(\frac{1}{n^2}+\frac{1}{n^3}+\dots )+B_1(\frac{2}{n^3}+\frac{3}{n^4}+\dots )+\\ & & \dots +B_{k-1}(\frac{k}{n^{k+1}}+\frac{k(k+1)}{2n^{k+2}}+\dots +\left(\begin{array}{c}k+r-1\\ r\end{array}\right)\frac{1}{n^{k+r}}+...)+...\text{(10)}\end{array}$$ To find the value of each of the terms $B_k$ we compare the expressions involving terms of the form $1/n^m$ for $m=2,3,4,\dots$. For example, the coefficient of $\frac{1}{n^2}$ on the left of the above equation is 1 and on the right it is $B_0$. Hence $$B_0=1.$$ If we next look at the terms in $1/n^3$ and $1/n^4$ we get (respectively) $$1/n^3:0=B_0+2B_1,\text{and}1/n^4:0=B_0+3B_1+3B_2.$$ Substituting the value of $B_0$ into the first equation gives $B_1$ and knowing this we can then find $B_2$ from the second equation. This gives $$B_1=-\frac{1}{2}\text{and}{B}_2=\frac{1}{6}.$$ We can then continue inductively, and can work out the terms $B_k$ from the recurrence relations $$B_0=1\text{and}{B}_{k-1}=-\frac{1}{k}(\left(\begin{array}{c}k\\ 2\end{array}\right)B_{k-2}+\left(\begin{array}{c}k\\ 3\end{array}\right)B_{k-3}+\dots +B_0)\text{(11)}$$ Turning the handle we get the following numbers $$B_0=1,B_1=-\frac{1}{2},B_2=\frac{1}{6},B_3=0,B_4=-\frac{1}{30},$$ $$B_5=0,B_6=\frac{1}{42},B_7=0,B_8=-\frac{1}{30},B_9=0,\dots .\text{(12)}$$ We notice that if $k$ is odd and bigger than one then $B_k=0$, and that the values of $B_k$ with $k$ even alternate in sign. The numbers that we have calculated are all quite small, but they get much bigger (rapidly!) as $k$ increases. For example $$B_{50}=\frac{495057205241079648212477525}{66}.$$

These numbers are famous and are called the Bernoulli numbers. They come up everywhere, from number theory to mechanics and beyond. Jakob Bernoulli described them in the book Ars Conjectandi (published posthumously in 1713) in connection with sums of powers of integers and they were discovered almost simultaneously by a number of other mathematicians. Since then, they have played a starring role in mathematical history. For example they are very important in the understanding of both Fermat's last theorem and the Riemann zeta function, and the first computer program written by Ada Lovelace in 1842 was designed to compute them.

Doing a quick sum

So, what's the catch with doing this? Well the problem with this approach is that as the Bernouili numbers $B_k$ increase rapidly as $k$ increases, then ultimately the correction term gets rather large. In fact if $n$ is {\em fixed} then the individual terms in $E_{n,k}$ tend to infinity as $k\to \mathrm{\infty }$. For an {\em optimal error} for a given value of $n$ it generally works best to truncate the series for the correction and to take $k=n$, estimating $S$ by $S_n+E_{n,n}$. This gives an error proportional to $1/n^n.$

We can easily play the same game with the other series mentioned in this article. See this appendix for details.

How well does this work?

Very! To show this we will look at a table of values of $S_n$ for $n=1,5$ and $10$, and then the three corrections given by $$P_n=S_n+\frac{1}{n},$$ $$Q_n=S_n+\frac{1}{n}-\frac{1}{2n^2}+\frac{1}{6n^3},$$ $$R_n=S_n+\frac{1}{n}-\frac{1}{2n^2}+\frac{1}{6n^3}-\frac{1}{30n^5}+\frac{1}{42n^7}-\frac{1}{30n^9}.$$ $$\text{Unknown environment 'center'}$$ We can see that the approximation $R_{10}$ is a very good approximation indeed for ${\pi }^2/6=1.644934066848226$ with an error of only $7\times {10}^{-13}$. This is especially impressive when we see that $S_{10}$ is a very poor approximation to ${\pi }^2/6$. Note that we have only had to add in an extra $6$ terms to do this. To gain the same level of accuracy with $S_n$ we would have to have included another ${10}^{11}$ terms! Readers are invited to work out $S_{20}$ and $R_{20}$ to further test the impressive accuracy of this method.

Is this new?

Not at all. The basic idea of accelerating the convergence of a series certainly goes back to the extraordinary mathematician Leonhard Euler, if not earlier. Faced with an {\em alternating series}, such as the first one (1) that we looked at to calculate $\pi /4$ in which the terms alternate in sign, Euler devised a transformation, which instead of summing the series, summed the various divided differences of the series. This method does not work for series with all positive terms such as the series (2) to calculate ${\pi }^2/6$. However, Euler needed to find the sum of this series as at the time (1735) its sum was unknown. The question of finding it (and finding a proof) was posed in 1644, and many distinguished mathematicians had tried, and failed, to find the sum. Euler's idea was to compare the sum of terms of the form $1/n^2$ with the integral of the function $1/x^2$. As was well known at the time $$\int \frac{dx}{x^2}=-\frac{1}{x}+C,$$ and the sum of the terms $1/n^2$ is an approximation to this integral. If you can find the error in this approximation then you can sum the series by comparing it to the known integral. In a (typical) tour-de-force, Euler was able to calculate this error deriving the (now famous) Euler-Maclaurin (divergent) series for the error. If you take $n$ terms of the original series for the sum and $k$ terms of the Euler-Maclaurin series then you get exactly our formula (13). It is pleasing (for a numerical analyst like me) to recount that Euler used his formula to calculate the sum of the series (2) correct to 20 decimal places (without of course the benefit of any form of caculator). From this calculation he was then able to {\em guess} that the sum of the series was ${\pi }^2/6$. Knowing the answer he was then able to find a proof.

The problem was renamed the Basel problem in honour of the home town of both Euler and the Bernoullis, and Euler's reputation (at the age of only 28) was made for life.

By his method of attack, Euler anticipated a lot of modern mathematics, in which the computer is used as an experimental tool, to gain insight into the solution of a mathematical problem as a first stage in proving the result. Nowadays we use this method all the time, and it continues to rely making highly accurate calculations. But what is nice about the techniques I have described in this article is that we can see that by doing only a little extra mathematics we can do much better than even the most powerful computers.

---

About the author

Chris Budd is Professor of Applied Mathematics at the University of Bath, Vice President of the Institute of Mathematics and its Applications, Chair of Mathematics for the Royal Institution and an honorary fellow of the British Science Association. He is particularly interested in applying mathematics to the real world and promoting the public understanding of mathematics.

He has co-written the popular mathematics book Mathematics Galore!, published by Oxford University Press, with C. Sangwin.