Buffon's Needle Problem

Author

Eric W. Weisstein
December 13, 2005

This notebook downloaded from http://mathworld.wolfram.com/notebooks/ComputationalGeometry/BuffonsNeedleProblem.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/BuffonsNeedleProblem.html.

©2005 Wolfram Research, Inc. except for portions noted otherwise

Analytic

General

PDF[UniformDistribution[{0,d/2}],x]

PDF[UniformDistribution[{0,Pi/2}],q]

Assuming[{l>0,d>0,0<q<Pi/2},Integrate[Integrate[4/(d Pi), {x, 0, (l*Sin[q])/2}] ,{q,0,Pi/2}]//Simplify
]

Assuming[{l>0,d>0,0<q<Pi/2},Integrate[Integrate[Piecewise[{{2/d, 0 <= x <= d/2}}, 0]2/Pi, {x, 0, (l*Sin[q])/2}] ,{q,0,Pi/2}]//Simplify
]

FullSimplify[%/.l->x d,d>0]

0<x<1

(1/(2π))Integrate[((l Abs[Cos[θ]])/d),{θ,0,2π}]

4(l/(2π d))Integrate[Cos[θ],{θ,0,π/2}]

x==0

%/.d->l

x>1

FullSimplify[(1/(π d))(d(π-2ArcSin[(d/l)])+2l(1-(Sqrt[1-((d^2)/(l^2))])))/.d->l/x,x>1]

FullSimplify[ArcCos[1/x],x>1]

FullSimplify[With[{ϕ0=ArcSec[x]},((2x)/π)(1-Sin[ϕ0])+((2ϕ0)/π)],x>1]

%/.d->l

Probability Plot

Plot[Piecewise[{{((2x)/π),0≤x≤1},{((2 (x-(Sqrt[(x^2)-1])+ ArcSec[x]))/( π)),x>1}}],{x,0,15},PlotStyle->Red,Epilog->{Blue,Dashing[{.02}],Line[{{1,0},{1,1}}]},
AxesLabel->TraditionalForm/@{x,P[x]}]

Estimators

N[Pi^2((1/2)Pi-1),10]

Calculating Pi Using Random Tosses

Code

Throws1[l_,n_]:=Module[{data=RandomReal[{0,1},{n,2}]},
Abs[Subtract@@Function[s,Floor[#[[1]]+s  l/2 Cos[2π#[[2]]]]]/@{1,-1}]&/@data
]

Throws2[l_,n_]:=Module[
{
data=Transpose[{
RandomReal[{0,1},n],
Cos[2π RandomReal[{0,1},n]]}]
},
Abs[Subtract@@Function[s,Floor[#[[1]]+s  l/2 #[[2]]]]/@{1,-1}]&/@data
]

Length[Throws1[1/3,10^5]]//Timing

Length[Throws2[1/3,10^5]]//Timing

Plot

t=Table[Throws2[1/3,10^6],{5}];//Timing

ListPlot[2/3/Take[Accumulate[#]/Range[10^6],{1,-1,1000}]&/@t,
Joined->True,
PlotStyle->{Red,Yellow,Green,Blue,Purple},
Epilog->{Line[{{0,Pi},{10^6,Pi}}]}
]//Quiet

Visualization

BuffonsNeedle[r_,{x_,y_},n_]:=Module[
{
data=Transpose[{Table[{x Random[Real],y Random[Real]},{n}],π Table[Random[Real],{n}]}],
lines,hits,misses
},
lines=(Function[s,First[#]+s  r/2 Through[{Cos,Sin}[Last[#]]]]/@{1,-1})&/@data;
{
Line[{{#,0},{#,y}}]&/@Range[0,Ceiling[x]],
PointSize[.01],
{misses,hits}=Map[Last,Split[Sort[Transpose[{Abs[Subtract@@(Floor/@First/@#)&/@lines],lines}]],(First[#1]==First[#2]&)],{2}];
Print["hits: ",Length[hits],"/",n," = ",N[100Length[hits]/n],"%, misses: ",Length[misses],"/",n," = ",N[100Length[misses]/n],"%,",Overscript[π,"^"]," = ",N[2n r/Length[hits]]];
Transpose[{
{Red,Green},
{Point/@Mean/@#,Line/@#}&/@{hits,misses}
}]
}
]/;0<=r<=1

Show[Graphics[BuffonsNeedle[1/3,{10,5},500]],AspectRatio->Automatic,ImageSize->500]

With[{n=500,x=1/3},((π^2)/(2n))((π/x)-2)]//N

Show[Graphics[BuffonsNeedle[.91,{2,1},10^4]/.{l_Line:>{},p_PointSize:>PointSize[.005]}],AspectRatio->Automatic,ImageSize->500]//Timing