Lab 6 – metoda Monte Carlo. Statystyka. Szkic M. Głowacki, L-10, PK
1
Laboratorium komputerowe 6 Metoda Monte Carlo. Statystyka
Zasadniczym celem ćwiczenia jest zastosowanie metody Monte Carlo do całkowania numerycznego (wyznaczania pola powierzchni), a w szczególności do wyznaczenia liczby 𝜋.
Ponadto w ćwiczeniu wykorzystane zostaną podstawowe funkcje statystyczne (wartość średnia, odchylenie standardowe) do zobrazowania działania tej prostej metody statystycznej.
Uwaga:
Wykonanie ćwiczenia nie wymaga żadnych dodatkowych umiejętności z zakresu metod statystycznych. Cała potrzebna wiedza jest przedstawiona poniżej.
1. Wprowadzenie – idea prostej metody Monte Carlo.
Rozważmy kwadrat o boku 𝑎 = 1. Pole takiego kwadratu jest równe 𝑃𝑘𝑤 = 𝑎2 = 1.
Lewy dolny wierzchołek kwadratu znajduje się w początku układu współrzędnych 𝑂(0,0).
W kwadrat wpisujemy ćwiartkę okręgu w taki sposób, że środek okręgu również znajduje się w początku układu współrzędnych.
Pole wycinka koła, ograniczonego dwoma bokami kwadratu i wpisaną w ten kwadrat ćwiartką okręgu (obszar czerwony) jest równe 𝑃 =1
4𝜋𝑟2 =1
4𝜋𝑎2 = 1
4𝜋. Jednak pole to możemy również obliczyć w inny sposób, korzystając z metody Monte Carlo.
Za pomocą generatora liczb pseudolosowych o rozkładzie równomiernym generujemy 𝑁 punktów wewnątrz kwadratu i następnie zliczamy ile punktów znalazło się wewnątrz obszaru ograniczonego okręgiem 𝑥2+ 𝑦2 = 1. W kolejnym kroku korzystamy z zależności:
𝑃 𝑃𝑘𝑤 = 𝑛
𝑁
gdzie 𝑛 to liczba punktów wewnątrz ćwiartki okręgu, 𝑁 to liczba wszystkich punktów, 𝑃 to szukane pole wycinka koła, a 𝑃𝑘𝑤 to pole kwadratu.
1 1
0 x
y
Lab 6 – metoda Monte Carlo. Statystyka. Szkic M. Głowacki, L-10, PK
2
2. Implementacja metody – wyznaczenie liczby 𝜋, ilustracja graficzna.
Korzystając z metody Monte Carlo opisanej we wprowadzeniu (punkt 1) napisać skrypt w Matlabie, który wyznaczy liczbę 𝜋.
Podpowiedź:
Korzystając z zależności:
𝑃 𝑃𝑘𝑤 = 𝑛
𝑁 otrzymujemy:
𝑃 = 𝑛 𝑁𝑃𝑘𝑤 1
4𝜋𝑎2 = 𝑛 𝑁𝑎2
𝜋 = 4𝑛 𝑁
Do generacji liczb pseudolosowych z przedziału (0, 1) można wykorzystać funkcję rand().
Każdy punkt 𝑃 na płaszczyźnie definiujemy za pomocą współrzędnych kartezjańskich 𝑃𝑥 oraz 𝑃𝑦, które możemy wygenerować osobno, przykładowo:
Px = rand(1);
Py = rand(1);
Sprawdzenie, czy punkt należy do obszaru ograniczonego okręgiem możemy wykonać za pomocą instrukcji warunkowej:
if Px*Px + Py*Py < 1 ...
end
Skrypt powinien mieć możliwość łatwej zmiany liczby generowanych punktów 𝑁.
Dodatkowo skrypt ma wykonywać ilustrację graficzną metody: wyświetlać ćwiartkę okręgu oraz wszystkie wygenerowane punkty z podziałem na dwa kolory – jeden kolor dla punktów znajdujących się wewnątrz obszaru ograniczonego okręgiem, a drugi dla pozostałych punktów.
Przykładowy wynik działania skryptu przedstawiono na następnej stronie (ze względu na statystyczny charakter obliczeń wyniki mogą się oczywiście różnić dla kolejnych uruchomień skryptu).
Lab 6 – metoda Monte Carlo. Statystyka. Szkic M. Głowacki, L-10, PK
3 Podpowiedź:
Do narysowania punktu można wykorzystać funkcję plot() z odpowiednim parametrem formatującym wyświetlanie, przykładowo:
plot(Px,Py, 'r*') % kolor punktu: czerwony
% styl punktu: gwiazdka
Każdy punkt może być rysowany osobno – w zmiennych Px i Py trzeba w takiej sytuacji zapisać pojedyncze wartości, tak jak na poprzedniej stronie. Po pierwszym użyciu funkcji plot trzeba pamiętać o instrukcji hold on.
3. Implementacja metody – uśrednianie i rozrzut wyników.
Przerobić skrypt przygotowany w poprzednim punkcie na funkcję obliczającą liczbę 𝜋 w zależności od podanej liczby generowanych punktów 𝑁. W tym przypadku należy pominąć część graficzną, a funkcja ma tylko zwracać liczbę. Przykładowy nagłówek funkcji:
function liczbaPI = metodaMonteCarloPI(N)
Napisać nowy skrypt, który z wykorzystaniem własnej funkcji do wyznaczania liczby 𝜋 wykona wielokrotnie powtórzenia dla kilku wybranych wartości 𝑁, a następnie dla każdego przypadku obliczy z uzyskanych wyników wartość średnią i odchylenie standardowe.
Przyjąć:
𝑁 = 10, 100, 1000, 10000 liczba generowanych punktów
𝑚 = 20 liczba powtórzeń dla każdej wartości 𝑁
Lab 6 – metoda Monte Carlo. Statystyka. Szkic M. Głowacki, L-10, PK
4
Narysować wykres, który będzie prezentował uzyskane wyniki. Na wykresie wyświetlić wartość średnią wraz z odchyleniem standardowym dla każdej wartości 𝑁 (odchylenie standardowe potraktować jako niepewność wyniku). Dla porównania narysować również prostą przedstawiającą wartość ścisłą 𝜋. Oś poziomą narysować w skali logarytmicznej.
Przydatne funkcje:
mean(a) % obliczenie wartości średniej dla
% elementów w wektorze a
std(a) % obliczenie odchylenia standardowego
errorbar(x,y,err) % rys. punktów wraz niepewnościami,
% liczba elementów w wektorze err
% musi być taka sama jak w x i y semilogx(x,y) % rys. wykresu, oś OX układu
% współrzędnych w skali logarytmicznej
Przykładowy wykres znajduje się poniżej.
4. Obliczanie całek oznaczonych – uwagi.
W podobny sposób możemy obliczać całki oznaczone funkcji, korzystając z interpretacji graficznej mówiącej, że pole pod wykresem funkcji jest równe odpowiedniej całce oznaczonej:
𝑃 = ∫|𝑓(𝑥)|𝑑𝑥
𝑏
𝑎
Losowo generujemy punkty w odpowiednim obszarze prostokątnym, a następnie zliczamy ile punktów znalazło się pod krzywą 𝑦 = |𝑓(𝑥)|.
