\Dan JonesMathematica Project  - Simulation of Buffon's Needle ProblemNovember, 1997\>

Project Description

Use Mathematica to simulate dropping congruent needles randomly on a grid of \parallel, equidistant lines. Calculate an approximate value of π using 	he solution to Buffon's Needle Problem.  For each simulation, show a \graphical representation of the needles and their positions with respect to 	he grid lines.

Buffon's Needle Problem

Statement of problem: what is the probability that a needle dropped on a \grid of parallel lines will land on one of the grid lines? 

The solution \is as follows (see Schroeder, L., "Buffon's needle problem: An exciting \application of many mathematical concepts," Mathematics Teacher, 67 (1974), \183-186):

For p = probability that a needle lands on a grid line, l = 
eedle length, k = distance between grid lines

p = ((2  l)/k)π,  l  ≤  k

p = (((2  l)/k)π  + 1) - (2/π)(((Sqrt[(l^2) - (k^2)])/k) + ArcSin[(k/l)]) ,  l  > k

 Given the value of π, one can calculate the \probability that a needle will land on a grid line.  Conversely, given an \empirical value for the probability that a needle will land on a grid line, \one can calculate an approximate value of π.
				
 For 𝓅 = \empirical probability that a needle lands on a grid line, l = needle length, \k = distance between grid lines
 
𝓅 ≈ ((2  l)/k)π,  l  ≤  k

𝓅 ≈ (((2  l)/k)π  + 1) - (2/π)(((Sqrt[(l^2) - (k^2)])/k) + ArcSin[(k/l)]) ,  l  > k

One can generate an empirical value of 𝓅 by 	ossing n needles on to a grid and counting the number of needles that land \on a line.

For hits = the number of needles that land on a line, and n = 	he number of tosses

𝓅 = (hits/n)

π ≈\(\(\(2  l\ n\)\/\(k\ hits\)
	
\)\),  l  ≤  k
π ≈ ((-2 l (Sqrt[ (l^2) - (k^2)]))/((k l (((hits/n)) - 1))            )) + ((2  l)/(k (((hits/n)) - 1))) -         ((2  ArcSin[(k/l)])/((((hits/n)) - 1))),  l  > k

Mathematica Simulation

In order to get an accurate approximation of π by tossing needles on a \grid, a very large number of needle drops is required (see further discusion \in the Observations section below).

A Mathematica model can be used to simulate the dropping of large numbers of needles on \a grid.  The value of π can be approximated from the results of this \simulation.

Note that for the purposes of the simulation, it is sufficient 	o work with a grid composed of two lines.  Given the premise noted above in 	he Project Description section that the grid lines are parallel and \equidistant and that the needles are randomly dropped, the probability that a 
eedle will land on a grid line is independent of the specific pair of grid \lines that the needle is dropped on.  Therefore, a pattern of needles on a \grid composed of multiple grid lines can be represented by a pattern of 
eedles on a single pair of lines.

Explanation of the Model

\The Mathematica code for the simulation has several different components.  \These components are described in the following subsections.  For each \subsection, code from the simulation is displayed and then an explanation \follows.\>

Module Function

Module[{topframe, bottomline, bottomframe, rtframe, topline, frame,
  esults, θ, centerx, centery, x1, y1, x2, y2 , hitlist, hits, pi,

\The model uses the function "Module" to set up variables as local \variables.  This keeps global variables from "colliding" with variables in the model.  The following example \illustrates how a global variable could collide with a variable from the \model:Assume that the user of the model wants to use the global variable "hits" 	o input the number of needles into the model, as follows: toss[hits, 4, 4].  \The model generates a number of needles equal to the global value of hits, \but because the variable hits is also used inside the model, the value of \hits changes; the original value is lost.  This is a problem if the user 
eeds the original value of hits for further calculations.The "Module" function sets up "hits" as a local variable.  Using the \Module function, the user gets the correct number of needles based on the \global value of hits.  The model uses hits as a local variable, but the \global value of hits remains unchanged after the simulation is run.\>

Options

Sample Function

toss[1000, 4, 4,  Calculation -> True,   GraphicOutput -> True,        TextDescription -> False,  ReportStatistics -> True]

Observations

Accuracy of Results

Using the solution to Buffon's Needle Problem to generate an approximation \of π requires a very large number of needle drops; simulating one million \drops may or may not give a correct result to two decimal places.  The \following are the results from running three simulations with one million 
eedles (needle length equal to distance between grid lines):

3.140
3.143\n3.132
Mean 3.138

Distribution of Results

Given an arbitrary interval containing the mean of the results from a set \of trials, increasing the number of drops increases the probability that a \single trial will fall within the interval.  The following are the results \from running the model 10 times with different numbers of needles (needle \length equal to distance between grid lines):

Number of Drops	Mean of 10 Trials		Standard Deviation of 10 Trials
1,000			       3.144	 			.062
1,000			       3.138	 			9
10,000			       3.135				.022
10,000			       3.141	 			0

The results of trials with 10,000 needles are disbursed in a more 
arrow distribution about the mean than the results of 1,000 needles, as \evidenced by the lower standard deviation values.

Given an arbitrary \interval containing the value π, increasing the number of drops used in \each element of a set of trials does not significantly increase the \probability that the mean of the set of trials will fall within the interval. \ One possible explanation for this is that even with large numbers of 
eedles, the model does not do a good job of evenly distributing the needles \over the grid.  This is apparent from a graphical representation of the 
eedles when running the model with large numbers of needles.  With a run of \10,000 needles, the center of the display is completely filled in, but the \edges are ragged, indicating uneven distribution of needles.  Averaging esults appears to smooth out the unevenness, whether running the model with \sets of trials with 1,000 needles or sets of trials with 10,000 needles.

\Running the model with needle length different than the distance between grid \lines increases the unevenness of the distribution of needles:

Number of Drops	Mean of 10 Trials 		Standard Deviation of 10 Trials
Needle length 40, distance between grid lines 4:
10,000			       \3.136				    .106

Needle length 4, distance between grid lines 40:
\10,000			       3.110				    .158

Random Numbers Generated by Mathematica

With surprising frequency, different trials return the same value of π 	o at least three decimal places.  For instance, for 10,000 needles, k=4, \l=40, the model returned a value of π equal to 3.187 three times out of 	en trials.  This may reflect the manner in which Mathematica generates "random" numbers, and the degree of randomness of the 
umbers.  One hypothesis for this behavior is that Mathematica cycles through a list of random numbers.  When Mathematica reaches the end of the list, it starts over.  If this were the case, then \one would get non-random behavior in a model using a quantity of random 
umbers close to the quantity of random numbers in Mathematica.

Data:

k4l41000a       = {3.185, 3.072, 3.091, 3.195, 3.106, 3.082, 3.150, 3.263, 3.110,         3.185};

k4l41000b       = {3.072, 3.236, 3.026, 3.110, 3.096, 3.012, 3.231, 3.135, 3.241,         3.221};

k4l410000a       = {3.156, 3.114, 3.162, 3.150, 3.136, 3.121, 3.092, 3.152, 3.128,         3.143};

k4l410000b       = {3.139, 3.144, 3.133, 3.178, 3.123, 3.174, 3.118, 3.129, 3.135,         3.141};

k4l41000000 = {3.140, 3.143, 3.132};

k4l4010000       = {3.187, 3.015, 3.187, 2.944, 3.249, 3.177, 3.271, 3.187, 3.089,         3.051};

k40l410000       = {3.454, 3.180, 3.026, 3.008, 3.008, 2.890, 3.175, 3.003, 3.170,         3.185};

<< Statistics`DescriptiveStatistics`

N[Map[        StandardDeviation, {k4l41000a, k4l41000b, k4l410000a, k4l410000b,           k4l41000000, k4l4010000, k40l410000}], 3]

N[Map[        Mean, {k4l41000a, k4l41000b, k4l410000a, k4l410000b, k4l41000000,           k4l4010000, k40l410000}], 4]