Igła Buffona

Piotr Sierant
Fizyka teoretyczna I rok, studia pierwszego stopnia
Wydział Fizyki, Astronomii i Informatyki Stosowanej
Uniwersytet Jagielloński
Kraków 2011
Mathematica 7.0

Wstęp

Celem projektu jest przedstawienie problemu tzw. Igły Buffona, zastosowanie go do wyznaczenia stałej π, analiza dokładności tej metody i porównanie jej z innymi sposobami wyznaczania π.

Zagadnienie Igły Buffona został sformułowany przez francuskiego arystokratę, matematyka i encyklopedystę Georges'a-Louisa Leclerca, hrabiego Buffon [1]:

,,Dana jest plansza z zaznaczonymi na niej poziomymi liniami odległymi o d. Upuszczamy na nią igłę o długości l, taką, że l ⩽ d. Jakie jest prawdopodobieństwo tego, że igła przetnie linię?''

Rozwiązanie problemu

Clear["Global`*"]
Graphics[{{Thickness[Large],Black,Line[{{0, -0.5}, {2, -0.5}}], Line[{{0, 0.5}, {2, 0.5}}]}, Arrow[{{0.1, -0.1}, {0.1, -0.5}}], Arrow[{{0.1, 0.1}, {0.1, 0.5}}], {Thickness[Large], Red, Line[{{0.8, 0}, {1.6, 0.6}}],Line[{{0.4, -0.3}, {1.34, -0.6}}]},{Thickness[Large], Gray, Line[{{0.95, -0.1}, {1.9, 0.15}}]}, {Dashed, Line[{{0.9, 0.3}, {1.5, 0.3}}]}, Text[Style[d, Large], {0.1, 0}],Text[Style[α,Large],{1.35,0.35}],Circle[{1.2,0.3},0.2,{0,π/5}]},ImageSize->{500,300}]

Oznaczenia jak w treści zadania. Pomijam dwa przypadki, których prawdopodobieństwo zajścia wynosi 0:
i) w którym igła przetnie obie linie (z założenia l ⩽ d, więc taki przypadek móglby zajść jedynie dla α = 90°, a środek igły leży dokłanie między liniami), 
ii) w którym igła będzie leżeć dokładnie na lini. Przy takich założeniach igła albo w ogóle nie przecina żadnej z linii, albo przecina jedną w dokładnie jednym punkcie.

Widać, że ,,efektywna wysokość" igły wynosi  l Sin[α], wobec czego prawdopodobieństwo, że tak ustawiona igła przetnie linię równa się  ((l Sin[α])/d). Żeby obliczyć prawdopodobieństwo, że dowolnie ustawiona igła przetnie linię należy wynik ten uśrednić po wszystkich możliwych kątach, czyli obliczyć poniższą całkę (przypadki dla 90°< α < 360° są symetryczne do powyższego i dają w wyniku takie samo prawdopodobieństwo):

P=(2/π)(∫_0^(π/2))(l/d)sin(x)\[DifferentialD]x

Z powyższych rozważań wynika praktyczny sposób wyznaczania stałej π. Jeżeli przeprowadzimy n niezależnych rzutów, z których w k rzutach przecięta zostanie linia (przez r oznaczmy (l/d)), to:

P=(k/n)

π=((2 l n)/(d k))=((2 n r)/k)

Wyznaczenie eksperymentalne wartości liczby π

Zadajemy stosunek długości igły do odległości między liniami oraz liczbę n prób w każdym z niezależnych j eksperymentów:

r=Input["Proszę podać r"];
While[r>1 || r<0,r=Input["Niestety podane r nie spełnia warunków zadania, proszę podać r jeszcze raz"];];
n=Input["Proszę podać n"];
j=Input["Proszę podać liczbę eksperymentów jakie należy przeprowadzić (sugeruję n*j < 10000)"];

Losujemy współrzędne środka igły oraz kąt nachylenia do poziomu dla j eksperymentów:

i=0;While[i<j,(Subscript[T, i])=Table[{9RandomReal[],9RandomReal[],2π RandomReal[]},{n}];i++];

Poniższe polecenie tworzy dla każdego eksperymentu dwie listy. W listach (Subscript[e, k]) przechowywane są wyznaczane wartości π, natomiast listy (Subscript[B, k])zawierają położenia kolejnych igieł, które potrzebne będą w dalszej prezentacji danych.

Ce[x_] := If[x>0,IntegerPart[x],IntegerPart[x]-1]
k=0;While[k<j,{
(Subscript[e, k])={};(Subscript[B, k]):={};l1=1;
For[i=0,i<n,{x=Extract[Extract[(Subscript[T, k]),{i}],{1}],y=Extract[Extract[(Subscript[T, k]),{i}],{2}],θ=Extract[Extract[(Subscript[T, k]),{i}],{3}],
If[Ce[y+Abs[(r/2) Sin[θ]]]>Ce[y] || Ce[y-(r/2) Abs[Sin[θ]]]<Ce[y],{l1++,(Subscript[B, k])=Append[(Subscript[B, k]),{Red,Line[{{x- (r/2)Cos[θ],y-(r/2)Sin[θ]},{x+ (r/2)Cos[θ],y+(r/2)Sin[θ]}}]}]},(Subscript[B, k])=Append[(Subscript[B, k]),{Gray,Line[{{x- (r/2)Cos[θ],y-(r/2)Sin[θ]},{x+ (r/2)Cos[θ],y+(r/2)Sin[θ]}}]}]],(Subscript[e, k])=Append[(Subscript[e, k]),{i,2r (i/l1)}]},i++],
k++}]

Podajemy otrzymane wartości π i rysujemy wykresy przedstawiające zależności otrzymanych wartości π w zależności od ilości rzutów dla poszczególnych eksperymentów.

P1[i_]:=ListPlot[(Subscript[e, i]),AxesOrigin->{0,2.5},ImageSize->{900,500},PlotRange->{{((-n)/80),n},{2.2,4.5}}]
P2= Plot[π,{x,0,n},AxesOrigin->{0,2.5},ImageSize->{900,500},PlotStyle->Directive[Red,Thick],PlotRange->{{((-n)/80),n},{2.2,4.5}}];
P3=Graphics[Text[Style["Wyznaczona wartość liczby π na podstawie tego eksperymentu wynosi" ,20],{(n/2.5),4.3}]];
pi[i_]:=Extract[Last[(Subscript[e, i])],{2}]//N
P4[i_]:=Graphics[Text[Style[pi[i] ,27,Green],{(n/1.20),4.3}]];
P7=Graphics[Text[Style["π",20,Red],{((-n)/150),3.14}]];
Manipulate[Show[P1[i],P2,P3,P4[i],P7,Background->LightBrown],{{i,0,"Numer eksperymentu"},0,j-1,1},ControlType->PopupMenu]

Tworzymy ilustrację pozwalającą obserwować proces rzucania igieł w poszczególnych eksperymentach.

P5=Table[Line[{{0,i},{9,i}}],{i,0,9}];
Manipulate[Graphics[{Take[(Subscript[B, a]),b],P5},Background->LightBrown,PlotRegion->{{0.,0.},{0.,0.}},ImageSize->{550,550}],{{a,0,"Numer eksperymentu"},0,j-1,1},{{b,0,"Liczba rzuconych igieł"},0,n,1},ControlPlacement->{Left}]

Porównanie i podsumowanie

Wyniki poprzedniej części pokazują, że powyższa metoda wyznaczania π jest o wiele wolniej zbieżna niż nowoczesne sposoby, dzięki którym superkomputery wyznaczają π z dokładnością do tryliardowych cyfr znaczących. Sensowne wydaje się porównanie z jedną z najstarszych i intuicyjnie najprostszych metod, której autorstwo przypisywane jest Archimedesowi [2], a której ideę przedstawia poniższa wizualizacja:

ngon[p_]:=Polygon[Table[{Sin[2Pi k /p],Cos[2Pi k /p]},{k,p}]]
Manipulate[Graphics[{{EdgeForm[Black],LightBlue,ngon[i]},{Thick,Red,Circle[{0,0},1]},{Thick,Black,Disk[{0,0},0.01]},{Line[{{0,0},{0,1}}]},Text[Style[R,16],{-0.05,0.5}]},ImageSize->470,Background->LightBrown],{{i,5,"Ilość kątów wielokąta [m]"},3,40,1}]

Wyznaczanie π tą metodą polega na zauważeniu, że obwód wielokąta foremnego wpisanego w okrąg (zaznaczony na czerwono) dąży do długości tego okręgu. Długość okręgu związana jest ze stałą π, a na obwód wielokąta dostajemy prosty analityczny wzór. Wobec tego:

π=m Sin[((180°)/m)]

Aby porównać obie metody rysujemy wykres wyznaczanego π od liczby kroków wykonanych przez algorytmy w obu sposobach. Zielona linia dotyczy metody Archimedesa.

P6=Plot[m Sin[((180°)/m)],{m,3,n},AxesOrigin->{0,2.5},ImageSize->{900,500},PlotStyle->Directive[Green],PlotRange->{{0,n},{2.2,4.5}}];
Manipulate[Show[P1[i],P2,P6,P7,Background->LightBrown],{{i,0,"Numer eksperymentu"},0,j-1,1},ControlType->PopupMenu]

Powyższy wykres pokazuje, że wyznaczanie liczby π za pomocą Igły Buffona jest procesem bardzo wolno zbieżnym, który podlega znaczącym fluktuacjom nawet dla dużych wartości n. Metoda ta wypada znacznie słabiej niż znany od starożytności sposób Archimedesa, wobec czego nie ma żadnego znaczenia praktycznego.

Należy więc zadanie o Igle Buffona traktować jako ciekawostkę, która pozwala nam uświadomić sobie, że opisując świat ściśle, przy pomocy matematyki, możemy w wielu zjawiskach natrafić na stałe matematyczne.

Bibliografia

[1] http://en.wikipedia.org/wiki/Buffon's_needle

[2] http://en.wikipedia.org/wiki/Pi