# Computing Pi with Colliding Blocks

This notebook attempts to use colliding blocks to calculate the digits of pi, based off a really cool video by 3Blue1Brown


```python
%%html

```





and the paper [Playing Pool with $\pi$](https://www.maths.tcd.ie/~lebed/Galperin.%20Playing%20pool%20with%20pi.pdf) by G. Galperin,

in which a simple problem is described: a block of mass $M$ moves at a constant speed $V$ on a frictionless surface towards a vertical wall. Closer to the wall there is another block of mass $m$ at rest. After the collision, the small block will hit the wall and be reflected back so that the blocks collide again.

If the ratio between the block masses is a power of 100, and the collisions are perfectly elastic, then the total number of collisions $N_\text{collisions}$ (i.e. collisions + reflections) is equal to certain number of digits of $\pi$. Basically, if the masses satisfy the relation $\frac{M}{m} = 100^{d-1}$, then the number of clacks is equal to the first $d$ digits of $\pi$.

## 1. Below is the javascript i attempted to use to get a nice display for the computation simulation

Only problem is that I wanted to try and use a physics engine for the collision calculations. The first one I used is matter.js, but I learned that they don't have continuous collision detection, so when setting the mass of the larger block > 1000kg, the inner mass goes so fast that the engine no longer detects its collisions, and the block just dissapears... I tried again with planck.js, which is a port of box2d and features ccd, however I get a similar problem, where the inner block goes so fast and begins to clip into the larger block at around 1000kg, throwing off the computation.


```python
%%html

```





## 2. Equations and definitions, according to the paper

### Definitions/Conditions
- when the two blocks hit each other we count it as a ***collision***
- when the inner block hits the wall we count it as a ***reflection*** 
- both collisions and reflections are elastic, so energy and momentum are conserved in these processes

We are interested in the total number of "clacks" given by 

\begin{equation} 
N_\text{clacks} = N_\text{collisions} + N_\text{reflections} \tag*{(1)}
\end{equation}

**Pre-condition:** 
- We initialize the simulation with the large block (of mass $M$) moving at some constant speed $V_0$ before colliding with the small block (of mass $m$) initially at rest. 

**Post-condition:** 
- The final state of the simulation will have the two blocks moving independently away from the wall, with the large block moving faster than the small block, guaranteeing no further collisions.

### Equations
**Reflections:** before hitting the wall, the small block moves with positive velocity, which is fully reversed after the collision. This can be proven in general by writting the conservation of total energy and momentum:

\begin{eqnarray}
\frac12 mv^2_\text{before} &=& \frac12 mv^2_\text{after}, \tag*{(2)}\\
mv_\text{before} &=& mv_\text{after}, \tag*{(3)}\\
\end{eqnarray}

and velocity is reversed, so 
\begin{equation} 
v_\text{before} = -v_\text{after} \tag*{(4)}
\end{equation}


**Collisions:** when the blocks collide with one another, their collision is perfectly elastic, so both momentum and energy is conserved, leading to the following:

\begin{align*}
\text{conservation of energy  :} \qquad & \frac12 MV^2_k + \frac12 mv^2_k  = \frac12 MV^2_{k+1} + \frac12 mv^2_{k+1}, \tag*{(5)}\\
\text{conservation of momentum  :} \qquad & MV_k - mv_k = MV_{k+1} + mv_{k+1}. \tag*{(6)}\\
\end{align*}

Doing some work, getting lazy, and using wolfram alpha magic to solve the rest for $V_{k+1}$ and $v_{k+1}$, we get 
\begin{align}
v_{k+1} &= \frac{m-M}{m+M}v_k + \frac{2M}{m+M}V_k \tag*{(7)}, \\
V_{k+1} &= \frac{2m}{m+M}v_k + \frac{M-m}{m+M}V_k \tag*{(8)}
\end{align}

which leads to a pretty nice computation for our simulation

## 3. Calculations


```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
```


```python
# if collision is between two blocks...
def collision(m, M, vec):
    msum = (m + M)
    v_1 = (m - M)/float(msum)
    v_2 = 2*M/float(msum)
    V_1 = 2*m/float(msum)
    V_2 = (M - m)/float(msum)

    v = v_1*vec[0] + v_2*vec[1] # (eq #7)
    V = V_1*vec[0] + V_2*vec[1] # (eq #8)
    return v,V
```


```python
# if collision is a reflection...
def reflection(vec):
    v = (-1)*vec[0] # reverse velocity of small block (eq #4)
    V = vec[1]
    return [v, V]
```


```python
# checks whether large block is moving away faster than smaller block, which is our end condition
def checkspeeds(speed1, speed2):
    if (speed1 <= 0 and speed2 <=0):
        if(abs(speed2) >= abs(speed1)):
            return False
    return True
```


```python
def run_sim(num_digits):
    M = 100**(num_digits-1)
    currentVector = [0,1]
    ctype = 'C' # the large block will always initially collide with the small block
    clackCount = 0
    v, V = [], []
    while(checkspeeds(currentVector[0], currentVector[1])):
        if (ctype == 'C'):
            currentVector = collision(1, M, currentVector)
            v.append(currentVector[0])
            V.append(currentVector[1])
            clackCount += 1 # increment counter on collision
            ctype = 'R'
        elif (ctype == 'R'):
            currentVector = reflection(currentVector)
            v.append(currentVector[0])
            V.append(currentVector[1])
            clackCount += 1 # increment counter on collision
            ctype = 'C'
            
    return clackCount, v, V
```


```python
print(run_sim(6)[0])
```

    314159



```python
print(np.pi)
```

    3.141592653589793


it works!
### now to plot the velocities


```python
def plot_velocities(V, v, num_digits, clack_count, linemarkers=True):
    M = 100**(num_digits-1)
    x = np.arange(0, len(V), 1)
    fig, ax1 = plt.subplots()
    ax1.set_title('Velocity of each block ($m=1, M={}$)'.format(M), fontsize=16)

    colorV = 'g'
    colorv = 'b'

    ax2 = ax1.twinx()
    if linemarkers:
        ax1.plot(x, V, 'o-', color=colorV)
        ax2.plot(x, v, 'o-', color=colorv)
    else:
        ax1.plot(x, V, '-', color=colorV)
        ax2.plot(x, v, '-', color=colorv)
    ax1.grid(False)
    ax2.grid(False)
    ax1.set_xlabel(r'$k$ (number of collisions)', fontsize=16)
    ax1.set_ylabel('velocity large block', fontsize=16, color=colorV)
    ax2.set_ylabel('velocity small block', fontsize=16, color=colorv)
    Vmax = 1.5
    vmax = max(v)*1.1
    xmax = max(x)*1.1
    ax1.plot([0, xmax], [0, 0], 'k-', lw=1)  # y=0 line
    ax1.plot([0, 0], [-2*Vmax, 2*Vmax], 'k-', lw=1)  # x=0 line
    ax1.text(xmax*7/9, 1.3, 'clacks: {}'.format(clack_count), fontsize=16)
    ax1.xaxis.set_ticks(np.arange(0, len(V), np.max([1, int(1/7*len(V))])))
    ax1.set_xlim(-xmax/50, xmax)
    ax1.set_ylim(-Vmax, Vmax)
    ax2.set_ylim(-vmax, vmax);
```


```python
clacks, v, V = run_sim(2)
v_1 = np.asarray(v)
V_1 = np.asarray(V)
```


```python
plot_velocities(V_1, v_1, 
                num_digits=2, 
                clack_count=clacks, 
                linemarkers=True)
```

So you can sort of see the makings of a velocity phase diagram, and I probably could have made one if I figured out a way to stretch the axes enough. 

## 4. Conclusion

Overall, a pretty cool problem. I learned a lot about elastic collisions (I also did some research on relativistic elastic collisions and learned that the equations I used break down when the speed of the colliding bodies is close to the speed of light), physics engines and their pitfalls, the fact that you can model interesting physics phenomenon with geometry, and also the sheer amount of stuff I can cram into a jupyter notebook.

907131c8858a8fd1b4273df939d826dc9d0b59d6
ipynb
Jupyter Notebook
Colliding Pi.ipynb
sirmammingtonham/efficientpicomputation
994e0a6e792ad2c3c92d6729361d16b80353aff2
Unlicense
Colliding Pi.ipynb
sirmammingtonham/efficientpicomputation
994e0a6e792ad2c3c92d6729361d16b80353aff2
Unlicense
Colliding Pi.ipynb
sirmammingtonham/efficientpicomputation
994e0a6e792ad2c3c92d6729361d16b80353aff2
Unlicense
Qwen/Qwen-72B
1. YES
2. YES