Urban Maths: Computing Pi – Inefficiently!

Urban Maths: Computing Pi – Inefficiently!


At the time of writing, the roads and pavements are icy and people are slipping and sliding all over the place. Seeing this brought to mind a weird way of calculating the digits of \pi by counting the collisions between sliding blocks. More specifically, consider elastic collisions between two blocks of mass M and m that slide without friction, as illustrated in Figure 1. When mass m hits the wall at the right it rebounds elastically.

Figure 1: Sliding blocks
Figure 1: Sliding blocks

All collisions are elastic and there is no friction, and so energy is conserved at all times. When the two blocks collide with each other momentum is conserved during the collision. When mass m collides with the wall momentum is not conserved; the velocity of the mass is simply reversed.

We will assume mass m is initially stationary (v = 0) a short distance from the wall, and mass M is moving towards it from the left with unit velocity (V = 1). If the two masses are equal then when they collide, mass M will stop and transfer its momentum to mass m, which will move to the right with unit velocity. Mass m will now collide with the wall and return with the same speed (velocity will have opposite sign) to collide again with mass M. Now mass m will stop and mass M will move off to negative infinity. There will have been three collisions in total.

Suppose we repeat the previous scenario with mass M now 100 times that of m. In this case m will again move off towards the wall, but M will not come to a complete stop – it will continue moving to the right at a slightly reduced speed. After m bounces off the wall it will collide with M, slowing it a little more, but still not stopping it, while itself shooting off towards the wall again. In fact m will bounce between the wall and M many times before it has not only slowed M to a stop, but also imparted enough momentum to get M moving to the left again at a speed greater than its own. Mass m will undergo 31 collisions before this happens.

Again we repeat the process, but with M now {10\;000} (100^2) times the mass of m. Again there will be multiple collisions before M moves away to negative infinity faster than m. This time m will undergo 314 collisions in total. Repeat with M 100^3 times the mass of m and m will undergo {3141} collisions. If M is 100^4 times the mass of m then m will experience {31\;415} collisions. Note the pattern in the number of collisions: 3, 31, 314, 3141, {31\;415}. These are the digits of \pi, and the pattern continues! For M = 100^Nm, then the number of collisions will be given by the first N + 1 digits of \pi!

If you don’t believe this, I don’t blame you! I didn’t myself when I first came across it, until I had done some numerical calculations to convince myself of the connection. This is a ridiculously inefficient way to calculate \pi of course, but the pattern was first published by G. Galperin [1]. A beautiful visual simulation of the process can be found on 3Blue1Brown’s YouTube channel [2].

We can calculate the velocities of the masses after each collision by making use of conservation of energy, conservation of momentum for the two mass collisions and velocity reversal at the wall, as follows:

Conservation of energy:

(1)   \begin{equation*} \frac{1}{2}MV^2 + \frac{1}{2}mv^2 = E.  \end{equation*}

Conservation of momentum:

(2)   \begin{equation*} MV + mv = P.  \end{equation*}

At the wall

(3)   \begin{equation*} v \to - v.  \end{equation*}

E is the total (kinetic) energy of the system and P is the momentum immediately before and after the two blocks collide (P is different for each such collision). It is convenient to multiply equation (1) by 2/M and divide equation (2) by M to get:

Conservation of energy:

(4)   \begin{equation*} V^2 + rv^2 = \frac{2E}{M}.  \end{equation*}

Conservation of momentum:

(5)   \begin{equation*} V + rv = \frac{P}{M},  \end{equation*}

where

(6)   \begin{equation*} r = \frac{m}{M}. \end{equation*}

Consider the kth collision between the two masses. From equations (4) and (5) we can see that the pre- and post-collision velocities are related by:

    \begin{align*} V^2_k + rv^2_k = V^2_{k-1} + rv^2_{k-1} \\ V_k + rv_k = V_{k-1} + rv_{k-1}. \end{align*}

From these we obtain:

(7)   \begin{equation*} v_k &= \frac{(r-1)v_{k-1} + 2V_{k-1}}{1 + r}  \end{equation*}

(8)   \begin{equation*} V_k &=\frac{2rv_{k-1} + (1-r)V_{k-1}}{1 + r}.  \end{equation*}

Of course we can simply use equation (3) to update the velocity of mass m when it collides with the wall. Below is some pseudocode that calculates the number of collisions, k.

1:  Set a value for r (say, r = 100^{-1})
2:  initialise velocities: v = 0, V = 1
3:  initialise collision counter: k = 0
4:  while V>v do
5:       increment counter k = k + 1
6:       store pre-collision velocities v_\text{old} = v, V_\text{old} = V
7:       update velocities when M and m collide:

8:       if m collides with wall (i.e. if v>0) then
9:            increment counter k = k + 1
10:          update velocity v = -v
11:    end if
12: end while
13: number of collisions = value of k

Why on earth does this process enumerate the digits of \pi? We tend to think of \pi in relation to circles. Is there a circle hidden somewhere here? As it happens, there is! Look more closely at the conservation of energy, equation (4). We can write this as V^2 + (\sqrt{r}v)^2 = 2E/M. Noting that initially V = 1 and v = 0 we have that 2E/M = 1, and if we define a scaled velocity:

(9)   \begin{equation*} vs = \sqrt{r}v, \end{equation*}

we can write the conservation of energy as:

(10)   \begin{equation*} V^2 + vs^2 = 1. \end{equation*}

Figure 2: Collisions with {M = 100m}

Equation (10) is that of a unit circle of course. By suitable modification of the listing in the pseudocode we can obtain the velocities, v and V, at every collision as well as the number of collisions. Calculating the corresponding values of vs from equation (9), and then plotting vs against V, together with a unit circle, we see the collisions in configuration space, as illustrated in Figure 2 (where r = 100^{-1}, or M = 100m). The collisions are marked with the red dots, which we can see all lie on the blue unit circle. The collisions between the two masses are represented by the dots on or below the horizontal axis, those between m and the wall by the dots above the axis. The solid red lines not only connect the collisions in sequence, starting with the rightmost collision at (1, 0), but the slanted lines code the conservation of momentum, and the vertical lines code changes in momentum.

To see how the momentum is coded in Figure 2 we rewrite momentum equation (5), making use of equation (9), and rearrange to get:

(11)   \begin{equation*} vs = \frac{P}{\sqrt{r}M}-\frac{1}{\sqrt{r}}V. \end{equation*}

This represents the slanted lines, all of which have slope -1/\sqrt{r}. These intersect the unit circle at the collision points. The momentum is different for each slanted line of course. It decreases by the equivalent of the current value of 2vs (its coded value in Figure 2) at each wall collision.

Computing-Pi-Inefficiently-figure 3
Figure 3

Well, ok, we can see a circle is involved, but how, exactly, do the digits of \pi appear from the collisions? To see this let’s re-plot the unit circle and the collision dots (this time without the momentum lines) and add radial lines from the origin to the dots, as in Figure 3 (the black dot represents a point where the final momentum line intersects the conservation of energy circle, but it does not represent a collision because at that point, both masses have negative velocities, with the large mass having the more negative velocity, thereby drawing away from the small mass). All the angles between consecutive radial lines, apart from the angle of the sector immediately below the black dot, are the same (I will not prove that here – see [3] for a proof).

Computing-Pi-Inefficiently-figure-4
Figure 4: First sector

To calculate the angle \theta, let’s home in on the first sector (i.e. that involving the first two collisions) as in Figure 4. Prior to the first we have V = 1 and v = 0 so from equations (7) and (8) we find the velocities after the first collision to be v = 2/(1 + r) and V = (1-r)/(1 + r) with vs = 2\sqrt{r}/(1 + r) from equation (9). The length of the straight dotted line between the first two collision points in Figure 4 is now easily found in terms of r from Pythagoras’ theorem to be

    \[ \sqrt{\left( 1 - \frac{(1-r)}{(1 + r)} \right)^2 + \left( \frac{2\sqrt{r}}{(1 + r)} \right)^2} \]

or 2\sqrt{r/(1 + r)}. The angle, \theta, is found from \sin(\theta/2) = 2\sqrt{r/(1 + r)}/2 so that:

(12)   \begin{equation*} \theta = 2\sin^{-1}\left(\sqrt{\frac{r}{1 + r}}\right). \end{equation*}

How many angles \theta (hence, collision sectors) are required to complete a circle? We set s\theta to be 2\pi, so, using equation (12):

(13)   \begin{equation*} s = \frac{\pi}{\sin^{-1}\left(\sqrt{r/(1 + r)}\right)}. \end{equation*}

In general, the value of s from equation (13) will not be an integer, but we must have an integer number of collisions, hence we find the number of collisions, n, from n = \operatorname{floor}(s).1 For example, with r = 100^{-1}, equation (13) results in s = 31.52 so the number of collisions, n = 31. With r = 100^{-2} we get s = 314.17 so n = 314. With r = 100^{-3} we find s = 3141.59 so n = 3141, and so on.

What is so special about the choice of negative powers of 100 for r? Why not choose 99, or 2, or any other number? We can see why by looking carefully at equation (13). Generally, with r = 100^{-N} and N>0, we have

    \[ \sqrt{\frac{r}{(1 + r)}} \approx \sqrt{r} = 10^{-N}. \]

With small values for its argument \sin^{-1}\alpha\approx \alpha, so equation (13) becomes s\approx \pi/10^{-N} or s\approx 10^N\pi. In other words, in our normal base ten notation, s is approximately \pi with the decimal point shifted to the right by N places.

Interesting alternative explanations of the colliding blocks phenomenon are given in [4] (a geometric explanation), and in [5] (a linear algebra explanation).

I guess all of this is cold comfort to those poor souls sliding about on the icy pavements!

Alan Stevens CMath FIMA

Notes

Usually, the \operatorname{floor} function is defined so that \operatorname{floor}(s) is the largest integer less than or equal to s. Here, we slightly modify that to the largest integer strictly less than s. The reason for this is to ensure that when r = 1, n is 3 not 4.

Acknowledgement

‘Urban Maths’ cartoonist: Adrian Metcalfe – www.thisisfruittree.com

References

  1. Galperin G. (2003) Playing Pool with \pi (The number \pi from a billiard point of view), Regul. Chaot. Dyn., vol. 8, no. 4, pp. 375–394.
  2. 3Blue1Brown (2019) The most unexpected answer to a counting puzzlehttps://www.youtube.com/watch?v=HEfHFsfGXjs&t=160s (accessed 14 October 2019).
  3. STEM cell (2019) A solution to 3Blue1Brown’s counting problem, https://www.youtube.com/watch?v=ils7GZqp_iE (accessed 14 October 2019).
  4. 3Blue1Brown (2019) So why do colliding blocks compute pi?, https://www.youtube.com/watch?v=jsYwFizhncE (accessed 14 October 2019).
  5. Dr Peyam (2019) 3b1b colliding block problem with matrices, https://www.youtube.com/watch?v=ODHjcPrQTEY (accessed 14 October 2019).

Reproduced from Mathematics Today, December 2019

Download the article, Urban Maths: Computing Pi – Inefficiently! (pdf)

Image credit: What a game! by Ron Bulovs / Flickr / CC BY 2.0
Published