Throwing $\pi$ at a wall

TL;DR

This paper presents a method to approximate pi by counting elastic collisions between two masses and a wall, where the number of collisions encodes the digits of pi, linking physics to mathematical constants.

Contribution

It introduces a novel physical simulation approach to calculate pi using elastic collisions, connecting classical mechanics with number theory.

Findings

Number of collisions equals the first d digits of pi

Method accurately encodes pi in collision counts

Applicable for integer values of d in the collision model

Abstract

We discuss a method for calculating $ \pi $ using elastic collision between two masses $ M $ and $ m $, with $ X = m/M = 10^{2(1-d)} $ where $ d $ is an integer, and a wall. The total number of collisions between $ M $, $ m $ and the wall corresponds to the first $ d $ digits of $ \pi $.

Equations

$$\frac{1}{2}M{v}^2+\frac{1}{2}m{u}^2=K=\frac{1}{2}M{v}_0^2,\text{and}M\mathbf{v}+m\mathbf{u}=\mathbf{P},$$

$$v=\frac{(1-X)v^{\prime }+2X{u}^{\prime }}{1+X},\text{and}u=\frac{2v^{\prime }+(X-1)u^{\prime }}{1+X},$$

$$\frac{v^2}{a^2}+\frac{u^2}{b^2}=1,\text{and}v+Xu=\frac{P}{m}$$

$$\mathbf{w}=\left(\begin{array}{c}v\\ u\end{array}\right),\mathbf{A}=\left(\begin{array}{cc}1-X& 2X\\ 2& X-1\end{array}\right),{\mathbf{w}}^{\prime }=\left(\begin{array}{c}{v}^{\prime }\\ u^{\prime }\end{array}\right).$$

$${\mathbf{w}}_n=\frac{{\mathbf{H}}_{\mathrm{o}}^{(n-1)/2}{\mathbf{w}}_1}{(1+X{)}^{(n-1)/2}},\text{and}{\mathbf{w}}_n=\frac{{\mathbf{H}}_{\mathrm{e}}^{n/2}{\mathbf{w}}_0}{(1+X{)}^{n/2}},$$

$${\mathbf{H}}_{\mathrm{o}}=\mathbf{A}\mathbf{W}=\left(\begin{array}{cc}1-X& -2X\\ 2& 1-X\end{array}\right),{\mathbf{H}}_{\mathrm{e}}=\mathbf{W}\mathbf{A}=\left(\begin{array}{cc}1-X& 2X\\ -2& 1-X\end{array}\right),\text{and}{\mathbf{w}}_1=\frac{\mathbf{A}{\mathbf{w}}_0}{1+X}.$$

$${\mathbf{w}}_{n+2}=\frac{{\mathbf{H}}_{\mathrm{o}}{\mathbf{w}}_n}{1+X},\text{and}{\mathbf{w}}_{n+2}=\frac{{\mathbf{H}}_{\mathrm{e}}{\mathbf{w}}_n}{1+X},$$

$${\mathbf{P}}_{\pm }=\left(\begin{array}{cc}\mp i\sqrt{X}& \pm i\sqrt{X}\\ 1& 1\end{array}\right),{\mathbf{P}}_{\pm }^{-1}=\frac{1}{2i\sqrt{X}}\left(\begin{array}{cc}\mp 1& i\sqrt{X}\\ \pm 1& i\sqrt{X}\end{array}\right),\mathbf{D}=(1-X)\left(\begin{array}{cc}{\lambda }_{-}& 0\\ 0& {\lambda }_{+}\end{array}\right),$$

$${\lambda }_{\pm }=1\pm i\frac{2\sqrt{X}}{1-X}=\frac{(1+X)e^{\pm i{\theta }_0}}{1-X},{\theta }_0={\tan }^{-1}\left(\frac{2\sqrt{X}}{1-X}\right).$$

$${\mathbf{w}}_n={\mathbf{R}}_{+}\left[(n-1)/2\right]{\mathbf{w}}_1,\text{and}{\mathbf{w}}_n={\mathbf{R}}_{-}\left[n/2\right]{\mathbf{w}}_0,$$

$${\mathbf{R}}_{\pm }\left[n\right]=\frac{{\mathbf{P}}_{\pm }{\mathbf{D}}^n{\mathbf{P}}_{\pm }^{-1}}{(1+X{)}^n}=\left(\begin{array}{cc}\cos (n{\theta }_0)& \mp \sqrt{X}\sin (n{\theta }_0)\\ \frac{\pm 1}{\sqrt{X}}\sin (n{\theta }_0)& \cos (n{\theta }_0)\end{array}\right).$$

$$v_n=v_0\cos \left(s{\theta }_0\right),u_n=(-1{)}^{n+1}\frac{v_0}{\sqrt{X}}\sin \left(s{\theta }_0\right),$$

$$A_0=\frac{1}{2}ab{\tan }^{-1}\left(\frac{2\sqrt{X}}{1-X}\right)=\frac{1}{2}ab{\theta }_0.$$

$$A_{n}=\frac{1}{2}ab\lvert\theta_{n+2}-\theta_{n}\rvert=\frac{1}{2}ab\bigg{|}\tan^{-1}\left(\frac{ab(v_{n}u_{n+2}-v_{n+2}u_{n})}{b^{2}v_{n}v_{n+2}+a^{2}u_{n}u_{n+2}}\right)\bigg{|},$$

$$v_n=\frac{(1-X)v_n\mp 2X{u}_n}{1+X},u_n=\frac{\pm 2v_n+(1-X)u_n}{1+X},$$

$$ab(v_n{u}_{n+2}-v_{n+2}{u}_n)$$

$$b^2{v}_n{v}_{n+2}+a^2{u}_n{u}_{n+2}$$

$$A_n=\frac{1}{2}ab{\tan }^{-1}\left(\frac{2\sqrt{X}}{1-X}\right)\equiv A_0$$

$$N=\frac{\pi ab}{A_0}=\frac{2\pi }{{\tan }^{-1}\left(\frac{2\sqrt{X}}{1-X}\right)}=\frac{\pi }{\sqrt{X}}{\left(\sum _{k=0}^{\infty }\frac{(-1{)}^k{X}^k}{2k+1}\right)}^{-1}.$$

$$N\approx \pi 10^{d-1},$$

$$N_{\mathrm{total}}=\mathrm{ceil}(N)-1.$$

Full text

Introduction

Obtaining an accurate value of $\pi$ is one of the oldest problems in mathematics. Early approximations of $\pi$ accurate to one decimal place dating back nearly 4000 years have been discovered in Babylon and Egypt (Beckmann, 1971). The oldest known algorithm for calculating $\pi$ is attributable to Archimedes who placed upper and lower bounds on its value by inscribing and circumscribing n-sided regular polygons on a circle in his treatise Measurement of a Circle. The advent of computers has enabled much faster and more accurate calculations of $\pi$, and the current best algorithm is capable of computing it to 22 trillion decimal places111https://pi2e.ch/blog/2016/10/31/hexadecimal-digits-of-pi/. A recent video published by 3Blue1Brown222https://www.youtube.com/watch?v=HEfHFsfGXjs brought our attention to an unusual method of calculating $\pi$ originally discovered by Galperin (2003). While the efficiency of this method pales in comparison to modern algorithms, it is certainly a novel approach.

An object of mass $M$ moves along a frictionless surface towards a second object of mass $m$ which is stationary, before colliding perfectly elastically with it (i.e. no energy is lost in the collision), propelling the second object towards an immovable wall. The second object undergoes perfectly elastic collision with the wall and is reflected back off the wall and collides again with the first mass as seen in the diagram in Figure 1. This process repeats until both objects are moving away from the wall and the speed of the first object exceeds the speed of the second. Galperin (2003) found that if the ratio of masses, $X=m/M$, is of the form $10^{2(1-d)}$, where $d$ is a positive integer, then the number of collisions is an integer consisting of the first $d$ digits of $\pi$.

Galperin (2003) and Aretxabaleta et al. (2017) consider the position and velocity of both objects and treat the problem as a billiard system where particles collide with each other and immovable boundaries. In this paper we take a different approach to achieve the same result and instead consider the evolution of the velocity phase-space of the system.

Mathematical formulation

Our system consists of masses $M$ and $m$ with velocities ${\bf v}$ and ${\bf u}$, respectively, with magnitudes $v$ and $u$. Collisions between the masses conserve both energy and momentum as we assume perfectly elastic collisions and ignore all external forces. The collision between mass $m$ and the wall switches the sign of the velocity of $m$. Therefore each collision with the wall results in change of momentum of the system while leaving the kinetic energy unaffected. The kinetic energy and momentum of the system may be written, respectively, as

[TABLE]

where $v_{0}$ is the initial speed of $M$. We choose the positive direction to be along the negative horizontal axis (towards the wall). Henceforth we denote ${\bf v}$ and ${\bf u}$ simply as $v$ and $u$ with the direction indicated by their signs. The velocity of $M$ and $m$ after they collide is given by

[TABLE]

where $X=m/M$ and $v^{\prime}$ and $u^{\prime}$ denote the velocity of $M$ and $m$ before the collision, respectively. We may express (1) as

[TABLE]

where $a=v_{0}$ and $b=v_{0}/\sqrt{X}$. Valid solutions for $v$ and $u$ are then given by the points of intersection of a straight line arising from the momentum equation with the kinetic energy ellipse. We may express the collision of $M$ and $m$ using (2) as ${\bf w}={\bf A}{\bf w}^{\prime}/(1+X)$ where

[TABLE]

Similarly the collision with the wall may be expressed as ${\bf w}={\bf W}{\bf w}^{\prime}$ with ${\bf W}={\rm diag}(1,-1)$. Denoting the velocity after $n^{\rm th}$ collision with a subscript $n$ allows us to write for $n$ odd and even, respectively, as

[TABLE]

where

[TABLE]

The odd values of $n$ corresponds to collision between $M$ and $m$ while the even values correspond to collision between $m$ and the wall. It follows from (5) that for odd and even $n$ we have

[TABLE]

respectively. We may write ${\bf H}_{\rm o}={\bf P}_{+}{\bf D}{\bf P}_{+}^{-1}$ and ${\bf H}_{\rm e}={\bf P}_{-}{\bf D}{\bf P}_{-}^{-1}$ with

[TABLE]

where

[TABLE]

This allows us to express (5) as

[TABLE]

for $n$ odd and even, respectively, where

[TABLE]

The matrices ${\bf R}_{+}[n]$ and ${\bf R}_{-}[n]$ are rotation matrices which, respectively, rotate a vector ${\bf w}$ anti-clockwise and clockwise along an ellipse of the form $v^{2}/k+u^{2}/(k/X)=1$ for any real $k$. This implies that for $n$ odd (even) the velocity ${\bf w}_{n}$ is obtained from ${\bf w}_{1}$ (${\bf w}_{0})$ through an anti-clockwise (clockwise) rotation along an ellipse of the form $v^{2}/k+u^{2}/(k/X)=1$ for any real $k$. We can express (10) as

[TABLE]

where $s=(n+1)/2$ for odd values of $n$, and $s=n$ for even values of $n$. This of course denotes the same ellipse as given by (3).

Figure 2 shows the ellipse described by the kinetic energy equation (3) for $X=m/M=0.01$ and $d=2$. The dots lying on the ellipse denote points $(v_{n},u_{n})$ as given by (12) with the arrows trace their development in phase space. The process starts at $(v_{0},0)$, corresponding to $(1,0)$, and continues until the speed of the larger mass, $v_{n}$, exceeds the speed of the smaller mass, $u_{n}$, and the larger mass is moving away from the wall: $|v_{n}|>|u_{n}|$ and $v_{n}<0$. Geometrically this corresponds to ${\rm w}_{n}$ being in the third quadrant of the energy ellipse.

Each solution ${\bf w}_{n}=(v_{n},u_{n})$ as given by (12) is of the form ${\bf w}_{n}=(a\cos\theta_{n},b\sin\theta_{n})$ so that the angle ${\bf w}_{n}$ makes with the horizontal axis satisfies $\tan\theta_{n}=au_{n}/bv_{n}$. In particular for the first collision the angle is $\theta_{0}$ and the area subtended is given by

[TABLE]

The area subtended by successive solutions ${\bf w}_{n}=(v_{n},u_{n})$ and ${\bf w}_{n+2}=(v_{n+2},u_{n+2})$ in the upper/lower portion of the ellipse is given by

[TABLE]

where the absolute value ensures that $A_{n}$ is positive for both positive (odd $n$) and negative (even $n$) angles. Now, using (7) we have

[TABLE]

where upper (lower) sign corresponds to odd (even) values of $n$. This then gives

[TABLE]

Substituting into (14) gives

[TABLE]

Therefore the area subtended by all adjacent solutions are equal.

Now, let $n$ denote the ratio of area of the energy ellipse and $A_{0}$:

[TABLE]

The series is absolutely convergent for $|X|<1$. For $X=1$ (equal mass) we have $N=4$ and for $X=10^{2(1-d)}$, with $d\geq 2$, we may approximate $N$ as

[TABLE]

with maximum error of $(\pi/3)10^{1-d}\ll 1$ for $d\geq 2$. This number approximately corresponds to the number of sectors that the energy ellipse is divided into, $N_{\rm total}$, noting that $N_{\rm total}A_{0}>\pi ab$ for $d\geq 2$, as the final velocity vector $\mathbf{w}_{n}$ is not precisely at $(-v_{0},0)$ as shown in Figure 3.

However, since the total number of collisions must be an integer, we can overcome this discrepancy by rounding, which gives

[TABLE]

noting that for $X=1$, there are 3 collisions but $N$ is precisely equal to 4 which precludes us from using the floor function. Substituting (20) into (21) we find that the integer $N_{\rm total}$ consists of $d$ digits corresponding to first $d$ digits of $\pi$.

Discussion and Conclusions

The framework that we have developed is analogous to Kepler’s Second Law, which states “A line joining a planet and the Sun sweeps out equal areas in equal time” (Kepler, 1609). In fact, this problem has ties to Keplerian motion in another way, in that particle collisions may be treated as an approximation to the dynamics of a satellite performing a slingshot maneuver around a massive object.333http://www.physics.usyd.edu.au/teach_res/mp/doc/mec_slingshot.htm

While the mass required to compute consecutive digits of $\pi$ grows exponentially, rendering physical implementation of this method untenable, it is, nevertheless, possible to determine the first few digits of $\pi$ through physical experiments. The most famous example is Buffon’s Needle (de Buffon, 1777) where a needle is dropped onto floorboards. The fraction of needles that lie across two boards is then dependent on the length of the needles, the distance between floorboards (which are both measurable quantities) and $\pi$. This method has been used to calculate $\pi$ to 5 decimal places using only three thousand needles (Lazzarini, 1901).

We note that while this experiment is entirely non-physical, qualitatively similar phenomena do occur in nature. Explosive astrophysical events such as solar flares (Ellison & Ramaty, 1985), supernovae (Reynolds & Ellison, 1992) produce observable outflows and shock waves. The primary mechanism by which these shock waves are driven is thought to be Fermi acceleration, where charged particles (analogous to the masses in our experiment) are repeatedly reflected by a magnetic mirror (analogous to the wall). This is also thought to be the origin of cosmic rays (Blandford & Eichler, 1987).

We have found a solution to the unusual method of calculating $\pi$ proposed by Galperin (2003). By considering constraints on the velocity phase-space of the system imposed by the conservation of energy and momentum, we take a geometric approach to show that if the ratio of masses of the two objects is of the form $10^{2(1-d)}$ then the number of times they collide is given by the first $d$ digits of $\pi$.

Acknowledgements

We thank A. Zic and M. Wheatland for useful discussions. We thank Grant Sanderson of 3Blue1Brown for bringing our attention to the topic.

MZR is supported by the Australian Research Council through grant DP160102932. DD is supported by an Australian Government Research Training Program Scholarship. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

Software: Matplotlib (Hunter, 2007), Numpy (van der Walt et al., 2011)