Igła Buffona

Prezentacja przygotowana przez Tomasza Urbańczyka, studenta V roku fizyki na Wydziale Fizyki, Astronomii i Informatyki Stosowanej UJ, na zaliczenie kursu obsługi programu Mathematica.

Bibliografia:

1) M.Aigner, G.M.Ziegler "Proofs from the book"

2) http://en.wikipedia.org/wiki/Buffon%27s_needle

W rzeczywistości eksperyment polega na rzucaniu igły o długości d na plansze z narysowanymi równoległymi liniami odległymi od siebie o t. Istotne jest aby długość igły była mniejsza od odległości między liniami.  Zliczana jest całkowita ilość rzutów (oznaczana przez n) oraz ilość rzutów, dla których igła przecięła którąś z linii (oznaczana przez a). Znając te wielkości można oszacować wartość liczby π.


W celu przeprowadzenia komputerowej symulacji eksperymentu należy przeprowadzić poniższe rozumowanie.

1) Z punktu widzenia eksperymentu ważne są 2 ,,parametry" określające położenie igły
    -odległość środka igły od najbliższej osi (oznaczana literą x, należy do przedziału [0,t/2] )
    -kąt między igłą a linią (oznaczany literą y, kąt ten należy do przedziału [0,π/2])  

2) Igła przecina oś tylko wtedy gdy 1/2×d×Sin[y] ≥  x

3) Wielokrotne rzucanie igłą można zasymulować za pomocą pętli ,,For" w której ciele znajdują się polecenia RandomReal odpowiedzialne za wylosowanie x i y  z  odpowiednich przedziałów.

4) Zmienne losowe x i y są niezależne ( przynajmniej w teorii - w praktyce mamy do czynienie z kongruencyjnym generatorem liczb pseudolosowych ) zatem podlegają rozkładowi jednostajnemu.

5) Z warunku unormowania funkcji gęstości prawdopodobieństwa i jednostajności rozkładu dostajemy wzory na funkcje gęstości prawdopodobieństwa:
   zmiennej x:    (2/t)dx 
   zmiennej y:     (2/π)dy 
   iloczynu xy (tu korzystamy z niezależności zmiennych losowych):  (4/(π t))dy dx 

6) Prawdopodobieństwo przecięcia linii przez igłę  to całka z funkcji gęstości prawdopodobieństwa:
  
                                               P=(∫_(y=0)^(π/2))(∫_(x=0)^((d/2)Sin[y]))(4/(π t))dxdy=((2d)/(t π))
  
 7) Z drugiej strony prawdopodobieństwo, że igła przetnie linie to stosunek (a/n)
  
 8)  Zatem  ((2d)/(t π))=(a/n)   a stąd wynika że 
  																			π = ((2 d n)/(a t))

Poniższy fragment notebooka rysuje wykres oszacowanej wartości liczby π w zależności od całkowitej ilości rzutów igłą. Wynik n-tego rzutu jest dodawany do wyniku poprzednich n+1 rzutów i na podstawie takich danych liczona jest liczba π dla konkretnego n.   Czerwona linia pokazuje tablicową wartość liczby π

d=5;   (*długość igły*)
t=6;   (*odległość między liniami*)
licznik=1000;  (*ile rzutów*)
ta=List[];
a=0; 

For[n=1,n≤licznik,n++,{x=RandomReal[t/2];y=RandomReal[π/2],
a=If[x≤d/2*Sin[y],a+1,a],If[a≠0,ta=Join[ta,List[N[2*n*d/(t*a)]]],ta]}]
"Po wykonaniu"
 licznik  
"rzutów igłą oszacowana wartość π wynosi"
If[a≠0,N[2*n*d/(t*a)],"nie można oszacować, igła anie razu nie przecięła linii - za mało rzutów"]
graf1=ListPlot[Flatten[ta]];
graf2=Graphics[{Red,Line[{{0,π},{licznik,π}}]}];
Show[graf1,graf2,PlotRange->{{0,licznik},{2,4}},Frame->True,FrameLabel->{"ilość rzutów igłą","wyliczona wartość π"}]

Natomiast dla poniższego wykresu dla każdego konkretnego n liczba π liczona jest "od nowa". Taka procedura znacznie wydłuża czas obliczeń.   Poniższy wykres mówi nam o "niedokładności" z jaką wyliczona została liczba π (dla każdego n liczba π liczona jest niezależnie, dla dużych n π powinna słabo zależeć od n, na poniższym wykresie widać silne fluktuacje wyliczonej liczby π, odzwierciedlające niedokładność metody)

d=5;   (*długość igły*)
t=6;   (*odległość między liniami*)
ile=500;  (*ile rzutów*)

a=0; 
liczbapi[licznik_]:={For[n=1,n≤licznik,n++,{x=RandomReal[t/2];y=RandomReal[π/2];a=If[x≤d/2*Sin[y],a+1,a]}];If[a≠0,wyn=2*licznik*d/(t*a),wyn=0];a=0;N[wyn]}
Tabliczka=List[];
For[g=2,g<ile,g++,Tabliczka=Join[Tabliczka,List[liczbapi[g]]]]
grafika1=ListPlot[Flatten[Tabliczka]];
grafika2=Graphics[{Red,Line[{{0,π},{ile,π}}]}];
Show[grafika1,grafika2,PlotRange->{{0,ile},{2,4}},Frame->True,FrameLabel->{"ilość rzutów igłą","wyliczona wartość π"}]
(*PS Nie sprawdzam czy wpisując nowy element do "Tabliczki" nie wpisuję 0 , bo 0 i tak jest poza zakresem wyświetlania rysunku *)

Poniższy fragment notebooka rysuje wykres obrazujący samo rzucanie igłą.  Zakładam, iż obszar na który może spaść środek igły jest kwadratem o boku (m) równym odległości między liniami (t) pomnożonej przez ilość linii. Aby było możliwe narysowanie igły potrzeba wylosować więcej informacji niż poprzednio.

1) Należy znać położenie środka igły wzdłóż osi poziomej - w tym celu losuję dodatkową zmienną  (z) zakresu  [0,m]

2) Tym razem nie wystarcza znajomość odległości środka igły od najbliższej linii. Trzeba mieć pełną informację o położeniu środka igły wzdłóż osi pionowej. W tym celu najpierw wylosowałem to położenie (p) , następnie (na podstawie wylosowanego położenia) znalazłem odległość środka igły od najbliższej linii (x) za pomocą dwóch następujących poleceń  j = Mod [p , t ]

x = If [ j<t/2 , j , t-j ]

3) Teraz, w przeciwieństwie do poprzednich obliczeń, nie wystarczy znać "nachylenia" igły względem linii (czyli kąta z przedziału [0,π/2], trzeba znać kąt (oznaczony wyrażeniem kat) z przedziału [0,π] (trzeba wiedzieć czy górny koniec igły zwrócony jest w lewo czy w prawo).  Na podstawie wylosowanej zmiennej (kat) obliczam "nachylenie igły" (y) na podstawie wzoru y = If [kat ≤ π/2 , kat , π-kat ]

t=6; (*odległość między liniami*)
ilosclinii=12;  (*ilość linii*)
d=4;   (*długość igły*)
ile=500; (* ile igieł narysować*)



a=0;
ilosclinii=ilosclinii-1;
m=t*ilosclinii;
polowad=d/2;
li={};
linie={};
For[i=0,i<ilosclinii+1,++i,{linie=Join[linie,List[Line[{{-polowad,i*t},{polowad+m,i*t}}]]];}]
For[i=0,i<ile,i++,{
p=RandomReal[m];j=Mod[p,t],x=If[j<t/2,j,t-j],kat=RandomReal[π],y=If[kat≤π/2,kat,π-kat],z=RandomReal[m],
a=If[x≤d/2*Sin[y],a+1,a];
li=Join[li,List[Line[{{polowad*Cos[kat]+z,polowad*Sin[kat]+p},{-polowad*Cos[kat]+z,-polowad*Sin[kat]+p}}]]];}]
"Oszacowana liczba π po"
ile
"rzutach igłą wynosi"
If[a≠0,N[2*ile*d/(t*a)],"igła ani razu nie przecięła linii, nie można oszacować π"]
Graphics[{li,{Red,linie}},PlotRange->{{-polowad,m+polowad},{-polowad,m+polowad}}]

Poniższy fragment Notebooka umożliwia narysowanie planszy z pewnym urozmaiceniem. Gdy igła przecina linię ma kolor czerwony, gdy linii nie przecina ma kolor niebieski.

t=6; (*odległość między liniami*)
ilosclinii=12;  (*ilość linii*)
d=4;   (*długość igły*)
ile=500; (* ile igieł narysować*)



a=0;
ilosclinii=ilosclinii-1;
m=t*ilosclinii;
polowad=d/2;
li={};
linie={};
For[i=0,i<ilosclinii+1,++i,{linie=Join[linie,List[Line[{{-polowad,i*t},{polowad+m,i*t}}]]];}]
For[i=0,i<ile,i++,{
p=RandomReal[m];j=Mod[p,t],x=If[j<t/2,j,t-j],kat=RandomReal[π],y=If[kat≤π/2,kat,π-kat],z=RandomReal[m],
a=If[x≤d/2*Sin[y],a+1,a];
li=Join[li,List[Graphics[{If[x≤d/2*Sin[y],Red,Blue],
Line[{{polowad*Cos[kat]+z,polowad*Sin[kat]+p},{-polowad*Cos[kat]+z,-polowad*Sin[kat]+p}}]}]]]}]

"Oszacowana liczba π po"
ile
"rzutach igłą wynosi"
If[a≠0,N[2*ile*d/(t*a)],"igła ani razu nie przecięła linii, nie można oszacować π"]
ti=Table[Graphics[{li,{Red,linie}},PlotRange->{{-polowad,m+polowad},{-polowad,m+polowad}}],{i,1,5}];
lin=Graphics[{Black,linie},PlotRange->{{-polowad,m+polowad},{-polowad,m+polowad}}];
Show[lin,li]