Lab 6 – przykładowe rozwiązania M. Głowacki, L-10, PK
1
Laboratorium komputerowe 6 – przykładowe rozwiązania
1. Skrypt obliczający liczbę 𝜋 metodą Monte Carlo z ilustracją graficzną.
n = 0; % liczba punktów wewnątrz ćwiartki okręgu
% narysowanie ćwiartki okręgu x = 0:0.01:1;
y = sqrt(1 - x.*x);
plot(x,y,'r');
hold on;
axis([0 1 0 1]);
grid on;
% generacja i analiza N losowych punktów for i=1:N
Px = rand(1); % współrzędna x punktu
Py = rand(1); % współrzędna y punktu
if Px*Px + Py*Py < 1 % czy punkt leży wewnątrz
% okręgu?
n = n+1;
plot(Px,Py,'r*');
else
plot(Px,Py,'g*');
end
end
liczbaPI = 4*n/N;
disp(liczbaPI);
Lab 6 – przykładowe rozwiązania M. Głowacki, L-10, PK
2 2. Uśrednianie wyników.
Funkcja obliczająca liczbę 𝜋, utworzona na podstawie skryptu z poprzedniej części.
function liczbaPI = metodaMonteCarloPI(N) n = 0;
for i=1:N
Px = rand(1);
Py = rand(1);
if Px*Px + Py*Py < 1
n = n+1;
end
end
liczbaPI = 4*n/N;
end
Skrypt wykonujący odpowiednie obliczenia i wyświetlający wartości średnie wraz z odchyleniem standardowym.
Uwaga:
Przyjmujemy 4 różne wartości dla liczby losowych punktów w metodzie Monte Carlo. Zostają one zapisane w wektorze N. Dla każdej ustalonej liczby punktów wykonywane jest m powtórzeń obliczeń. Wyniki zapisywane są w macierzy PIwyniki o rozmiarach length(N) x m (w tym konkretnym przypadku 4 x m, ale pozostawiamy ogólny zapis, który umożliwia w prosty sposób dodanie kolejnych przypadków). Czyli w każdym i-tym wierszu macierzy otrzymamy w kolumnach wyniki z m powtórzeń obliczeń dla danego N(i).
Lab 6 – przykładowe rozwiązania M. Głowacki, L-10, PK
3
m = 20; % liczba powtórzeń
N = [10 100 1000 10000]; % wektor z liczbami punktów
% powtórzenie m razy obliczeń dla każdego N
for i = 1:length(N) % dla każdej wartości N
for j = 1:m % wykonywane jest m powtórzeń PIwyniki(i,j) = metodaMonteCarloPI(N(i));
end
PIavg(i) = mean(PIwyniki(i,:)); % wartości średnie
PIstd(i) = std(PIwyniki(i,:)); % odchylenia standardowe
end
% narysowanie wartości ścisłej: prosta y = pi,
% wyznaczona na dwóch punktach
% (w tym konkretnym przypadku w zakresie x od 1 do 100000) xsc = [N(1)/10 10*N(end)];
PIsc = [pi pi]; % wartość ścisła – stała pi semilogx(xsc, PIsc, 'b--'); % wykres z osią OX w skali log.
hold on;
% narysowanie wartości średnich z odchyleniem standardowym errorbar(N,PIavg,PIstd,'or');
Uwaga
Nie jest wymagane definiowanie pętli przechodzącej po elementach wektora N. Rozpatrujemy tylko 4 przypadki, więc można rozpisać je po kolei.
% przypadek 1 dla N=10 for j = 1:m
PIwyniki(j) = metodaMonteCarloPI(N(1));
end
PIavg(1) = mean(PIwyniki);
PIstd(1) = std(PIwyniki);
% ...
I analogicznie możemy rozpisać dla pozostałych przypadków, zastępując część skryptu w ramce.
