Tag Archives: python

More Buffon’s Needle

In my previous blog entry, I discussed Buffon’s needle (a method of measuring \pi using random numbers). I tried to recreate what would occur in physical reality without using trig functions in the measurement of \pi. The tricky part was trying to make a uniform probability distribution along the angle \theta, which is what you expect with physical needles being dropped. In this entry, I will do something more exact, but less physical.

In reality, the probability density in \theta is uniform due to there being no preferential direction for a needle to point. To simulate this with uniform random numbers without using trigonometric functions (to avoid just measuring python’s \pi) is tricky. The way I did this was assign uniform random numbers (uniform probability density) to a projection of a needle, either perpendicular or parallel to the lines (with a probability of 0.5 assigned to each case). This results in a probability distribution in \theta that isn’t quite uniform, but it’s close enough to approximate \pi after a large number of needles are dropped:
2 \frac{N_T}{N_C} = \pi
where N_T and N_C are the total number of needles dropped and number of needles crossing a line, respectively.

If instead of assigning a uniform probability to either projection fairly, what if we assign it to just one? This was shown in the previous post. As the number of needles dropped increases, the ratio 2 \frac{N_T}{N_C} approaches some number around 2.52.6, shown in the figure below (\pi shown in red). Why is this the case?
DPI1

The reason why has to do with a uniform probability density along a projection (q) does not translate into a uniform probability density along \theta, see figure below. This means needles have a preferred direction. It turns out this can still be used to extract a measurement of \pi!
BuffonS2
The uniform probability density is assigned to the projection q:
\rho_q = 1/L, for 0 \leq q \leq L, and \rho_q=0 elsewhere.
To get this in terms of \theta, some substitutions are made in the probability density integral:
\int_0^L \rho_q dq = 1
The variables q and \theta are related by:
\text{cos}(\theta) = q/L
so by taking the derivative we get the substitution:
dq = -L\text{sin}(\theta) d \theta

So now the probability density integral can be written as:
\int_0^{\pi/2} \rho_q L \text{sin}(\theta) d \theta
which is also:
= \int_0^{\pi/2} \rho_{\theta} d \theta
So it gives the form of the probability density for \theta in terms of the one for q:
\rho_{\theta} = \rho_q L\text{sin}(\theta)

Plug in \rho_q = 1/L to get:
\rho_{\theta} = \text{sin}(\theta)
and it is obviously NOT uniform, it’s just the sine function!
Plotting a histogram of the number of needles per angle \theta will give the shape of the sine function:
DPI2

Following the same steps as in the previous entry, we can find out why when a uniform probability is assigned to one projection, 2 N_T/N_C approaches 2.5 instead of \pi.
For simplicity W=L, and again for a needle to be crossing a line it has to be below the upper bound:
x \leq \frac{L}{2}\text{sin}(\theta)
and the probability of a needle crossing a line is:
P = \int_0^{\frac{L}{2}\text{sin}(\theta)} \int_0^{\pi/2} \rho_x \rho_{\theta} d \theta d x
Plugging in the probability densities: \rho_x = 2/L and \rho_\theta = \text{sin}(\theta):
P = \int_0^{\pi/2} \text{sin}^2(\theta) d \theta = \pi/4

P = \frac{N_C}{N_T} when the number of needles dropped is large, so by combining these equations:
2 \frac{N_T}{N_C} = \frac{8}{\pi} = 2.52
This is the number that is being approached!
But wait, \pi is there, so we can extract it by:
4 \frac{N_C}{N_T} = \pi

By plotting this ratio as a function of number of needles dropped, we see it approaching \pi. See the figure below.
DPI3

Even though this isn’t how needles behave in reality, it’s a simpler way to measure \pi by using random numbers in a simulation. The code used to make these plots is shown below.

By the last drop:
2 \frac{N_T}{N_C} = 2.54
4 \frac{N_C}{N_T} = 3.14

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

CL = 0
NN = 0

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

for ii in range(100000):

    PPy = random.uniform(0.0,1.0) #A-end y position (lines separated by unit distance)

    RRx = random.uniform(0.0,1.0) #x projection
    RRy = numpy.sqrt(1.0 - RRx**2) #y projection

    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)

    PIb.append(4*RAT)

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

print PIa[-1]
print PIb[-1]

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

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

pylab.figure(3)
pylab.xlabel("NT")
pylab.ylabel("4 NC/NT")
pylab.plot(NNa,PIc,'r')
pylab.plot(NNa,PIb,'k')
pylab.savefig('DPI3.png')

pylab.show()
Advertisements

Imaginary Powers

Imaginary numbers (i = \sqrt{-1}) are used to relate Cartesian coordinates to polar coordinates as shown in the figure below. The argument \varphi is the angle from the real axis along the unit circle, and can be expressed as two terms: a real and an imaginary, giving the Cartesian (x,y) position.
e^{i\varphi} = \text{cos}(\varphi) + i \text{sin}(\varphi) = x + i y
(If you want to expand beyond the unit circle just multiply by the desired radius: e^{i\varphi} \rightarrow re^{i\varphi})

euler
This is known as the Euler relation, and it gives us some neat results, like:
e^{i\pi} = -1
e^{i \pi /2} = i
But what happens when you exponentiate i with itself? i^i.

From the Euler relation we know that i=e^{i \pi /2}, so it’s really a matter of substitution:
i^i = e^{i \pi / 2 i} = e^{-\pi / 2} = 0.20788...
So i^i is a real number!

But wait a second, this isn’t the whole answer. We can add or subtract 2\pi from any angle and wind up in the same place on the circle. So i = e^{i \pi /2 + i 2 \pi n}, where n is any integer. So we end up with a set of values:
i^i = e^{- \pi /2 - 2 \pi n}
which are all real!

icurve
The curve above goes through the set of values for i^i, where n is an integer.
i^i = \{...,111.32,0.20788,0.00038820,...\}
But 0.20788… will be looked at as the primary value for now since it is the result when n=0, and having one value keeps things simple.

So what happens when we repeat this exponentiation? i^{i^i}
First, this can be written in a simpler way:
^3i = i^{i^i}
This operation is called “tetration” or “power tower”.
You start from the top and work your way down. So ^3i = i^{i^i} = i^{0.20788}
Now use the Euler relation again (ignoring 2 \pi n terms in the exponential):
i^{0.20788} = e^{i \pi/2 \times 0.20788} = e^{i 0.32653}
So the result is complex:
^3i = 0.94715+0.32076i

We can keep going using results from the last iteration:
^4i = i^{0.94715+0.32076 i} = 0.05009+0.60211 i

Eventually you find that with ^mi as m increases, the results start to converge on one number. This implies there’s a value for ^\infty i. After plugging and chugging we can eventually conclude that:
^\infty i = 0.43828 + 0.36059 i
Convergence can be seen in both the real and imaginary components, around m \approx 70 if we require the difference between iterations to be less than 0.1%. (See figures below).

Refig Imfig

Each result has a different radius associated with it, since the imaginary component of the previous result puts a real term in the exponential of the next. This real term can be expressed as a factor r:
e^{i\varphi} = a + bi \rightarrow i^{a+bi} = e^{i \pi/2 (a + bi)} = e^{-b\pi/2}e^{i a \pi/2} = r e^{i a \pi/2}
The radius can be found by multiplying the result by its complex conjugate and then find the square root (i.e. the Pythagorean theorem).

We can then plot radius of ^mi as a function of m (see plot below).

Imfig

The python code for making these plots is:

import os
import sys
import math
import numpy
import pylab

REAL = []
IMAG = []
RADII = []
ITERA = []

pp = 1j
for n in range(100):
    ppL = pp #save previous value for cut off
    pp = (1j)**pp #interate power tower
    RA = 0.001 #ratio for cut off
    if (abs(pp.real - ppL.real)/(ppL.real + 0.00001) < RA) and (abs(pp.imag - ppL.imag)/(ppL.imag + 0.00001) < RA): #establish a cut off
        print n
        print pp
        break
    else:
        print numpy.sqrt(pp*pp.conjugate()) #radius

        REAL.append(pp.real) #save real component
        IMAG.append(pp.imag) #save imaginary component
        RADII.append(numpy.sqrt(pp*pp.conjugate())) #save radius
        ITERA.append(n+2) #save iteration "m"

RADII = [1]+RADII #add in the initial value for the plot
ITERA = [1]+ITERA
REAL = [0]+REAL
IMAG = [1]+IMAG

pylab.figure(1)
pylab.xlabel("m")
pylab.ylabel("radius")
pylab.plot(ITERA,RADII,'k')

pylab.savefig('Radfig.png')

# pylab.show()

pylab.figure(2)
pylab.xlabel("m")
pylab.ylabel("imaginary component")
pylab.plot(ITERA,IMAG,'k')

pylab.savefig('Imfig.png')

pylab.figure(3)
pylab.xlabel("m")
pylab.ylabel("real component")
pylab.plot(ITERA,REAL,'k')

pylab.savefig('Refig.png')