Buffon's Needle Simulation In Python

June 06, 2024 Post a Comment

import numpy as np import matplotlib.pylab as plt class Buffon_needle_problem: def __init__(self,x,y,n,m): self.x = x #width of the needle self.y = y #witdh o

Solution 1:

The sampling of the needle's alignment should be a uniform cosine. See the following link for the method: http://pdg.lbl.gov/2012/reviews/rpp2012-rev-monte-carlo-techniques.pdf

Also, there were a few logical problems with the program. Here is a working version.

#!/bin/pythonimport numpy as np

defsample_cosine():
  rr=2.while rr > 1.:
    u1=np.random.uniform(0,1.)
    u2=np.random.uniform(0,1.)
    v1=2*u1-1.
    rr=v1*v1+u2*u2
  cc=(v1*v1-u2*u2)/rr
  return cc

classBuffon_needle_problem:

     def__init__(self,x,y,n,m):
        self.x = float(x)  #width of the needle
        self.y = float(y)  #witdh of the space
        self.r = [] #coordinated of the centre of the needle
        self.z = [] #measure of the alignment of the needle
        self.n = n  #no of throws
        self.m = m  #no of simulations
        self.p = self.x/self.y
        self.pi_approx = []

    defsamples(self):
        # throwing the needlesfor i inrange(self.n):
            self.r.append(np.random.uniform(0,self.y))
            C=sample_cosine()
            self.z.append(C*self.x/2.)
        return [self.r,self.z]

    defsimulation(self):
        # m simulationfor j inrange(self.m):
            self.r=[]
            self.z=[]
            self.samples()
            # n throw
            hits = 0#setting the success to 0for i inrange(self.n):
                #condition for a hitif self.r[i]+self.z[i]>=self.y or self.r[i]-self.z[i]<0.:
                    hits += 1else:
                    continue
            est =self.p*float(self.n)/float(hits)
            self.pi_approx.append(est)
        return self.pi_approx

y = Buffon_needle_problem(1,2,80000,5)

print (y.simulation())

Solution 2:

Buffon's needle work accurately only when the distance between the two lines is double the length of needle. Make sure to cross check it.

I have seen many baffon's online simulation which are doing this mistake. They just take the distance between two adjacent lines to be equal to the needle's length. That's their main logical errors.

Solution 3:

I would say that the problem is that you are defining the alignment of the needle by a simple linear function, when in fact the effective length of the needle from its centre is defined by a sinusoidal function.

You want to calculate the effective length of the needle (at 90° to the lines) by using a function that will calculate it from its angle.

Something like:

self.z.append(np.cos(np.random.uniform(-np.pi/2, np.pi/2))*self.x)

This will give the cosine of a random angle between -90° and +90°, times the length of the needle.

Baca Juga

For reference, cos(+/-90) = 0 and cos(0) = 1, so at 90°, the needle with have effectively zero length, and at 0°, its full length.

I have neither mathplotlib or numpy installed on this machine, so I can't see if this fixes it, but it's definitely necessary.

Solution 4:

Looks like you were committing a simple rounding error. The code below works, though the results are not very close to pi...

import numpy as np
import matplotlib.pylab as plt

classBuffon_needle_problem:
def__init__(self,x,y,n,m):
    self.x = x #width of the needle
    self.y = y #witdh of the space
    self.r = []#coordinated of the centre of the needle
    self.z = []#measure of the alingment of the needle
    self.n = n#no of throws
    self.m = m#no of simulations
    self.pi_approx = []

defsamples(self):
    # throwing the needlesfor i inrange(self.n):
        self.r.append(np.random.uniform(0,self.y))
        self.z.append(np.random.uniform(0,self.x/2.0))
    return [self.r,self.z]

defsimulation(self):
    #self.samples()# m simulationsfor j inrange(self.m):
        self.r=[]
        self.z=[]
        for i inrange(self.n):
            self.r.append(np.random.uniform(0,self.y))
            self.z.append(np.random.uniform(0,self.x/2.0))
        # n throws
        hits = 0# setting the succes to 0for i inrange(self.n):
            # condition for a hitif self.r[i]+self.z[i]>=self.y or self.r[i]-self.z[i] <= 0.0:
                hits += 1else:
                continue
        hits = 2.0*(float(self.x)/self.y)*float(self.n)/float(hits)
        self.pi_approx.append(hits)

    return self.pi_approx

y = Buffon_needle_problem(1,2,40000,5)
print (y.simulation())

Also note that you were using the same sample for all simulations!

Solution 5:

I used Python turtle to approximate the value of Pi:

from turtle import *
from random import *

setworldcoordinates(-100, -200, 200, 200)
ht(); speed(0); color('blue')

drops = 20# increase number of drops for better approximation
hits = 0# hits counter# draw parallel lines with distance 20 between adjacent linesfor i inrange(0, 120, 20): 
    pu(); setpos(0, i); pd()
    fd(100) # length of line# throw needles
color('red')
for j inrange(drops):
    pu()
    goto(randrange(10, 90), randrange(0,100))
    y1 = ycor() # keep ycor of start point
    seth(360*random())
    pd(); fd(20)    # draw needle of length 20
    y2 = ycor() # keep ycor of end pointif y1//20 != y2//20:    # decisive test: if it is a hit then ...
        hits += 1# increase the hits counter by 1print(2 * drops / hits)

Output samples
With 50 drops   3.225806451612903with200 drops  3.3057851239669422with1000 drops 3.1645569620253164

Python Playground