Ilustracja metody
MONTE CARLO
obliczania
całek podwójnych
Często jest tak, iż wiemy, że istnieje całka oznaczona z funkcji f
jednak nie potrafimy jej analitycznie policzyć. Konieczne jest
wtedy zastosowanie jakiejś metody numerycznej.
Metody całkowania oparte o użycie liczb losowych pojawiły się
w latach czterdziestych XX wieku. W Los Alamos w ramach
"Projektu Manhattan" (nad budową bomby jądrowej pracowali
wtedy m.in. John von Neumann, Stanisław Ulam, Richard
Feynman) potrzebne było obliczenie całek dotyczących
rozpraszania i absorpcji neutronów. Wykorzystana została
do tego metoda Monte Carlo.
Nazwa Monte Carlo, kasyna w Monaco, nawiązuje
do losowości gier hazardowych.
Projekt prezentuje zastosowanie metody Monte Carlo do
obliczania całek podwójnych. Zgodnie z Mocnym Prawem
Wielkich Liczb Kołmogorowa, całka dana jest wzorem:
b
n
b - a
n
f (X ) �ł �ł�ł"
"
k
+"f (x)dx
n
k = 1
a
Przybliżenie całki dla całek podwójnych wygląda następująco:
(b - a)2 n
f (xk , yk )
"
+"+"f (x, y)dxdy H"
n
k = 1
D
Gdzie xk i yk to odpowiednie losowe liczby
yk = a + (b - a)vk
xk = a + (b - a)uk
Przykłady
całek
x3 + y3 dxdy
+"+"
D
x3 + y3 dxdy
+"+"
D
D = {(x, y) : 2 e" x e" - 2,0 d" y d" 4 - x2 }
x3 + y3 dxdy
+"+"
D
Obliczenia analityczne
4- x2
2 4- x2 2
1
�ł �ł
(x3 + y3)dy = dx =
+"dx +" +"�ł 4 y4 �ł
�ł łł
- 2 0 - 2
0
2 2
1 1
+"(4 - x2)2 dx = 4 +"(16 - 8x2 + x4)dx =
4
- 2 - 2
2
1 8 1 128
�ł
16x - x3 + x5 �ł =
�ł �ł
4 3 5 15
�ł łł
- 2
Metoda Monte Carlo
D'= [- 2;2]� [- 2;2]
~obszar D zawarty jest w kwadracie
~n liczb ui i vi losujemy z rozkładu jednostajnego U (0;1)
I przekształcamy według wzorów.
~wykorzystując wylosowane liczby możemy przybliżyć całkę:
42 n
3
( (xk 3 + yk 3)
"
+"+"x + y3)dxdy H" n
k = 1
D
Stworzony został wykres przedstawiający rozkład wyników dla wykonanych 50 prób dla każdego n,
wyskalowany został półlogarytmicznie, pominięto ukazanie rozrzutu wartości dla n=1 i n=10 ze
względu na ich duże różnice
Wyniki prób Monte Carlo Wartość całki
10
9
y = 8,53333
8
7
6
10 100 1000 10000 100000 1000000
Ilosc prób Monte Carlo n
Wartość obliczona metodami standardowymi to 8,5(3) . Metodą Monte Carlo, dla 50 prób
przy n=1000000 uzyskujemy wartość średnią 8,535893, co stanowi bardzo dobre
przybliżenie rzeczywistej wartości całki.
Wartosci uzyskane
sin( y)
dxdy
+"+"
y
D
sin(y)
dxdy
+"+"
y
D
D = {(x, y) : 0 d" x d" Ą , x d" y d" Ą }
sin(y)
dxdy
+"+"
y
D
Obliczenia analityczne
x= y
y
Ą Ą
�ł �ł
sin(y) sin(y) sin(y)
dxdy = dx = dy =
+"+" +"dy+" +"�ł x y �ł
�ł �ł
y y
�ł łł
D 0 0 0
x= 0
Ą
=
+"sin( y)dy = - [cos(Ą ) - cos(0)] = 2
0
Metoda Monte Carlo
D'= [0;Ą ]� [0;Ą ]
~obszar D zawarty jest w kwadracie
~n liczb ui i vi losujemy z rozkładu jednostajnego
U (0;1)
I przekształcamy według wzorów na xi i yi
~wykorzystując wylosowane liczby możemy przybliżyć całkę:
2
n
sin(y) 4Ą sin( yk )
dxdy H"
"
+"+"
y n yk
k = 1
D
Wyniki prób Monte Carlo Wartość całki
4
3,5
3
2,5
y = 2
2
1,5
1
10 100 1000 10000 100000 1000000
Ilosc prób Monte Carlo n
Wykres przedstawiający rozkład wyników dla wykonanych
100 prób dla każdego n. Wartość obliczona metodami
standardowymi to 2 . Metodą Monte Carlo, dla 100 prób przy
n=1000000 uzyskujemy wartość średnią 1,99978.
Wartosci uzyskane
2
(
+"+"x + y2)dxdy
D
2
(
+"+"x + y2)dxdy
D
D = {(x, y) : y e" 0, y d" x2 + y2 d" x}
2
(
+"+"x + y2)dxdy
D
Obliczenia analityczne
Aby obliczyć całkę analitycznie,
konieczne jest przejście na współrzędne biegunowe.
Obszar, po którym całkujemy, jest fragmentem koła
y d" x2 + y2 d" x
2
�! � sin� d" � d" � cos�
�! sin� d" � d" cos�
we współrzędnych biegunowych odpowiada on:
Ą
ńł �ł
" = (� ,� ) : 0 d" � d" ,sin� d" � d" cos�
�ł żł
4
ół �ł
Ą Ą
� = cos�
cos�
4 4
1 �ł
�ł
2 2 2 4
( � � d� d� =
�ł
+"+"x + y2)dxdy = +"+" +"d� +"� � d� = +"d� �ł 4 � �ł
�ł
łł
D " 0 sin� 0
� = sin�
Ą
Ą
� =
4
4
1 1 1
= (cos4 � - sin4 � )d� = (cos� sin� ) =
+"
4 4 8
� = 0
0
w obliczenia wykorzystane zostaly calki :
1
4
+"cos � d� = 32 (12� + 8sin(2� ) + sin(4� ))
1
4
+"sin � d� = 32 (12� - 8sin(2� ) + sin(4� ))
1
(cos4 � - sin4 � )� d� = *16sin(2� ) = sin� cos�
+"
32
Metoda Monte Carlo
~obszar D zawarty jest w kwadracie
D'= [0;1]� [0;1]
~n liczb xi i yi losujemy bezpośrednio
z rozkładu jednostajnego
U (0;1)
~wykorzystując wylosowane liczby możemy przybliżyć całkę:
n
1
2
( (xk 2 + yk 2)
"
+"+"x + y2)dxdy H" n
k = 1
D
Wyniki prób Monte Carlo Wartość całki
0,25
0,2
0,15
y = 0,125
0,1
0,05
0
10 100 1000 10000 100000 1000000
Ilosc prób Monte Carlo n
Wartosci uzyskane
e- ( x2 + y2 )dxdy
+"+"
D
e- ( x2 + y2 )dxdy
+"+"
D
D = {(x, y) : x e" 0, y d" x, x2 + y2 d" ln3}
e- ( x2 + y2 )dxdy
+"+"
D
Obliczenia analityczne
e- ( x2 + y2 )dxdy
D = {(x, y) : x e" 0, y d" x, x2 + y2 d" ln3}
+"+"
D
D we współrzędnych biegunowych odpowiada
Ą Ą
" = {( f , p) : - d" f d" ,0 d" p d" ln 3}
2 4
Ą
ln 3
4
- (x2 + y2 ) - p2 - p2
dxdy = df
+"+"e +"+"e dfdp = +" +"pe dp = *
Ą
D " 0
-
2
Ą Ą
�ł - p2
łł
z =
Ą
4 4
ln 3
1 1 Ą
4
�ł śł
- p2
= (- e- p2 0 df = f = * z'= - 2 p =
+"pe dp = �ł śł
+" +"
Ą
-
2 3 4
2
Ą Ą
�ł
dz = - 2 pdpśł
- -
�ł �ł
2 2
1 1
z p2
-
+"e dz = - 2 e-
2
Metoda Monte Carlo
~obszar D zawarty jest w kwadracie D'= [- ln 3;ln 3]� [- ln 3;ln 3]
U (0;1)
~n liczb ui i vi losujemy z rozkładu jednostajnego
i przekształcamy według wzorów.
~wykorzystując wylosowane liczby możemy przybliżyć całkę:
n
2 2
4ln2 3
e- ( x2 + y2 )dxdy H" e- ( xk + yk )
"
+"+"
n
k = 1
D
Wyniki prób Monte Carlo Przewidywana wartość całki
1,75
1,25
y = Ą/4
0,75
0,25
1 10 100 1000 10000 100000 1000000
Ilosc prób Monte Carlo n
Wykres przedstawiający rozkład wyników dla wykonanych 80
prób dla każdego n
Wartosci uzyskane
10xxsin( y)dxdy
+"+"
D
10xxsin( y)dxdy
+"+"
D
D = {(x, y) : 0 d" x d" 1; 0 d" y d" x}
10xxsin( y)dxdy
+"+"
D
Obliczenia analityczne
Całki tej nie da się wyrazić za pomocą funkcji
elementarnych. Tutaj z pomocą przychodzi metoda Monte
Carlo, dzięki której możemy poznać jej przybliżoną wartość.
Metoda Monte Carlo
~obszar D zawarty jest w kwadracie
D'= [0;1]� [0;1]
~n liczb xi i yi losujemy bezpośrednio
z rozkładu jednostajnego
U (0;1)
~wykorzystując wylosowane liczby możemy przybliżyć całkę:
n
10
10xxsin( y)dxdy H" xk xk sin( yk )
"
+"+"
n
k = 1
D
Wyniki prób Monte Carlo Przewidywana wartość całki
4
3
2
y = 1.31947
1
0
1 10 100 1000 10000 100000 1000000
Ilosc prób Monte Carlo n
Wartosci uzyskane
Autorzy projektu:
" Dorota Moskal 171695
" Mateusz Sawicki 171620
" Mariusz Orda 171665
" Michał Piórek 171677
" Mateusz Fuławka 171623
" Miłosz Karolonek 171595