Happy Pi Day!

Since it is \pi day, I thought I would show a way to measure \pi. The way I have in mind is known as Buffon’s needle. If you have a flat surface with parallel lines of uniform separation and drop needs of uniform length on it, the ratio between the total number needles dropped and the number of needles crossing a line will be a factor of \pi.

The figure below shows the needles of length L in blue, and the lines in black with separation W. For each needle there are two important values: x and \theta. The x value is the distance between the center of the needle and the nearest line. The \theta value is the acute angle between the needle and the nearest line. The span of available values of x is then 0 \leq x \leq W/2, and the span of \theta is 0 \leq \theta \leq \pi/2.

Buffon3
Considering all possible values of x and \theta, the probability is 1 by definition.
P = \int_{\Delta \theta} \int_{\Delta x} \rho dx d\theta = 1
Where in general, if x and \theta are independent: \rho = \rho(x)\rho(\theta).

A uniform probability densities are assumed for both x and \theta since there is no preferred direction or position for needles. So the probability densities \rho are defined by:
\rho_{\theta} = 2/\pi, for 0 \leq \theta \leq \pi/2, and \rho_{\theta} = 0 elsewhere.
\rho_{x} = 2/W, for 0 \leq x \leq W/2, and \rho_{x} = 0 elsewhere.
A uniform probability density is shown in the figure below for \theta, note that the area is 1.

buffon5
For simplicity, we will only consider the case where the needles are shorter or the same length as the separation distance between the lines: L \leq W.
So for needles crossing a line, the projected half-length of the needle perpendicular to the lines is greater than the distance between the center of the needle to the nearest line:
x \leq \frac{L}{2}\text{sin}(\theta)
This is the upper bound that is then plugged into the probability integral, to get the probability of crossing a line:
P = \int_{0}^{\pi/2} \int_{0}^{\frac{L}{2}\text{sin}(\theta)} \rho_{x} \rho_{\theta} dx d\theta = \frac{2 L}{\pi W}
For the simple case where L = W:
P = \frac{2}{\pi}
After enough drops the ratio between the number of needles crossing a line N_C and the total number of needles N_T will resemble the probability:
P = \frac{N_C}{N_T}
Combining these two equations, we get ratio that can measure \pi:
2 \frac{N_T}{N_C} = \pi

Doing this measurement in real life by dropping needles would make a mess, so I wrote a python code to it instead. Below is a graph showing this ratio as a function of needle drops (shown in black), \pi is shown in red. You can see the black line approaches \pi as the number of needles increases.

BPI

The probability density as a function of \theta is shown below. It is not perfectly flat as in the ideal case. This is because in my code the uniform probability density is along the Cartesian axes (X,Y) rather than along \theta. (The parallel lines are separated along Y). I did this to avoid measuring python’s \pi in its trigonometric functions.

BPI2

A uniform random number is assigned to the projected needle length along X or Y with a range of 0 to 1. The length of a needle is forced to be 1 by using the Pythagorean theorem: X = \sqrt{1-Y^2} or Y = \sqrt{1-X^2}. There is a 0.5 probability that either X or Y is assigned to the uniform random number. This approximates a uniform probability density in \theta. It’s not perfect though.

If however, the uniform random number is assigned to only X, the ratio doesn’t make it to \pi, it stops around 2.6, and the probability density is far from uniform. See the figures below:

BPI1q

BPI2q

The code for these plots is shown below:

import os
import sys
import math
import numpy
import pylab
import random

CL = 0
NN = 0

PIa = []
NNa = []
PIc = []
ANGLES = []

for ii in range(10000):

    PPx = random.uniform(1.0,10.0) #1-10 A-end x position
    PPy = random.uniform(0.0,10.0) #1-10 A-end y position (lines)

    RRx = random.uniform(-1.0,1.0) #B-end x displacement
    RRy = numpy.sqrt(1.0 - RRx**2) #B-end y displacement

    if (random.uniform(0.0,1.0)>0.5): #balance probability density to make it more uniform in terms of acute angle
        SIGN = 1
        if (random.uniform(0.0,1.0)>0.5):
            SIGN = -1
        RRy = random.uniform(0.0,1.0) #B-end y displacement
        RRx = numpy.sqrt(1.0 - RRy**2)*SIGN #B-end x displacement


    QQx = PPx + RRx #B-end x position
    QQy = PPy + RRy #B-end y position
    QQyf = math.floor(QQy)
    PPyf = math.floor(PPy)

    NN = NN + 1
    if (QQyf-PPyf > 0.0): #if crossing line
        CL = CL + 1

    NNa.append(NN)
    PIc.append(math.pi)
    RAT = (1.0*CL)/(1.0*NN)
    if (CL==0):
        PIa.append(0)
    else:
        PIa.append(2.0/RAT)

    ANGLES.append(abs(math.atan(RRy/(RRx+0.0001)))/math.pi)

print PIa[-1]

pylab.figure(1)
pylab.xlabel("NT")
pylab.ylabel("2 NT/NC")
pylab.plot(NNa,PIc,'r')
pylab.plot(NNa,PIa,'k')
pylab.savefig('BPI1q.png')

pylab.figure(2)
pylab.hist(ANGLES, bins=10)
pylab.xlabel("Acute angle in factors of pi")
pylab.ylabel("Needles dropped")
pylab.savefig('BPI2q.png')

pylab.show()
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s