Buffon’s Needle Problem

This is a simulation of a famous problem in continuous probability.  Drop a needle of length L onto a grid of parallel lines spaced one unit of distance apart.  Find the probability that the needle crosses one of the lines.  If the length of the needle is at most one, the needle can cross at most one line and the probability is predicted to be P =  2*L/Pi.  In particular, if the needle has length 1/2 the probability should be the reciprocal of Pi.

To simulate this we need random real numbers.  We can get them from Mathematica as follows.

Random[Real,{0,1}]

Suppose the lines are drawn in the (X,Y) plane and are horizontal.  Let Y in  [0, 1] be the position of the center of the needle relative to the two nearest lines and let theta in [0,Pi] be the angle of the needle with respect to the horizontal.    A random throw of the needle will be represented by a random point (Y,theta) in the rectangle [0,1] x [0, Pi].  Here is a list of 100 random tosses together with a plot of the corresponding points in the rectangle.

RandomNeedle := {Random[Real,{0,1}],Random[Real,{0,Pi}]}

Table[RandomNeedle,{i,1,100}]
ListPlot[%]

The needle crosses the line below the center if Y <= L*Sin[theta]/2 and crosses the line above if Y>= 1-L*Sin[theta]/2.  We can create random variables with values 1 or 0 according to whether these crossings do or do not occur.   The logical symbol || means “OR” in the function AnyCrossing[  ]

UpCrossing[{y_,theta_},needlelength_]:= If[y≥1-needlelength*Sin[theta]/2,1,0]
DownCrossing[{y_,theta_},needlelength_]:= If[y<=needlelength*Sin[theta]/2,1,0]
AnyCrossing[{y_,theta_},needlelength_]:= If[y≥1-needlelength*Sin[theta]/2 || y<=needlelength*Sin[theta]/2,1,0]

To run the simulation we can just record whether crossing occur for many random trials.  Adding up the corresponding ones and dividing by the total number of trials gives the frequency.

NeedleSimulation[needlelength_,numbertrials_] :=
Table[AnyCrossing[RandomNeedle,needlelength],{i,1,numbertrials}]

CrossingFrequency[needlelength_,numbertrials_]:= Plus@@NeedleSimulation[needlelength, numbertrials]/N[numbertrials]

Now to run the simulations.  Here is what happens for a needle of length 1/2.

NeedleSimulation[1/2,100]
CrossingFrequency[1/2,100]

We can compare the reciprocal of the result with Pi for various numbers of trials.  The convergence to Pi is not very fast.

N[Pi,30]
1/CrossingFrequency[1/2,100]
1/CrossingFrequency[1/2,1000]
1/CrossingFrequency[1/2,1000000]

For needles longer than one unit, we did not compute the exact probability.  It could be done, but here is a simulation for various lengths.

lengthfreqtable = Table[{length, CrossingFrequency[length,10000]},{length,1,10,0.2}]

ListPlot[%,Joined->True,AspectRatio->0.5]