Spis treści

Spis treści

1 Kwadratury interpolacyjne w wielu wymiarach 3

1.1 Sformułowanie zadania . . . 3

1.2 Interpolacja na siatkach regularnych . . . 4

1.2.1 Postać wielomianu interpolacyjnego . . . 4

1.2.2 Bład interpolacji_, . . . 6

1.3 Kwadratury interpolacyjne . . . 7

1.3.1 Kwadratury proste . . . 7

1.3.2 Kwadratury złożone . . . 8

1.4 Przekleństwo wymiaru . . . 9

2 Metody Monte Carlo 11 2.1 Wstep, metody niedeterministyczne_, . . . 11

2.2 Klasyczna metoda Monte Carlo . . . 11

2.2.1 Definicja i bład_, . . . 11

2.2.2 Całkowanie z waga_, . . . 13

2.3 Redukcja wariancji . . . 14

2.3.1 Losowanie warstwowe . . . 14

2.3.2 Funkcje kontrolne . . . 16

2.4 Generowanie liczb (pseudo-)losowych . . . 17

2.4.1 Liniowy generator kongruencyjny . . . 17

2.4.2 Odwracanie dystrybuanty i ‘akceptuj albo odrzuć’ . . . 18

2.4.3 Metoda Box-Muller dla rozkładu gaussowskiego. . . 19

3 Metody quasi-Monte Carlo 20 3.1 Co to sa metody quasi-Monte Carlo?_, . . . 20

3.2 Dyskrepancja . . . 21

3.3 Bład quasi-Monte Carlo_, . . . 22

3.3.1 Formuła Zaremby . . . 22

3.3.2 Nierówność Koksmy-Hlawki . . . 24

3.4 Ciagi o niskiej dyskrepancji_, . . . 25

3.4.1 Ciag Van der Corputa_, . . . 25 1

3.4.2 Konstrukcje Haltona i Sobol’a. . . 26 3.4.3 Sieci (t, m, d) i ciagi (t, d)_, . . . 27

Rozdział 1

Kwadratury interpolacyjne w wielu wymiarach

1.1 Sformułowanie zadania

Ostatnie trzy wykłady poświecimy numerycznemu całkowaniu funkcji wielu zmiennych. Dokład-_, niej, dla danej funkcji f : [0, 1]^d→ R chcemy obliczyć (przybliżyć) wartość

I_d(f ) = Z

[0,1]^d

f (~x) d~x = Z 1

0

Z 1 0

· · · Z 1

0

| {z }

d

f (x₁, x₂, . . . , x_d) dx₁dx₂· · · dx_d.

Zakładamy, że powyższa całka istnieje. W ogólniejszym sformułowaniu, chcielibyśmy obliczyć całke z wag_, a ω funkcji f : R_, ^d→ R, która jest postaci

I_d,ω(f ) = Z

R^d

f (~x) ω(~x) d~x.

Waga ω jest tutaj nieujemna i całkowalna.

Zauważmy, że ograniczenie sie w ostatnim przypadku do R_, ^dnie zmniejsza ogólności, gdyż całke_, po dowolnym mierzalnym obszarze D ⊆ R^dmożna wymodelować przyjmujac, że waga ω(~_, x) = 0 dla wszystkich ~x /∈ D.

Zadanie całkowania funkcji wielu zmiennych ma ogromne znaczenie praktyczne i dlatego warto znać skuteczne metody numeryczne jego rozwiazywania._,

Przyk?ad 1.1. Wycena obecnej wartości wielu instrumentów finansowych, w tym tzw. opcji, opiera sie na założeniu, że przyszłe ceny podlegaj_, a losowym zmianom kolejnych odcinkach czaso-_, wych. Obecna wartość opcji obliczana jest jako wartość oczekiwana funkcji wypłaty. Odpowiada to obliczaniu całki oznaczonej funkcji d zmiennych, gdzie d jest liczba odciników czasowych. Jest_, to czesto całka ze standardow_, a d wymiarow_, a wag_, a gaussowsk_, a postaci_,

(2π)^−d/2 Z

R^d

f (ξ1, . . . , ξd) exp^−1

2(ξ₁²+ · · · , ξ_d²)^dξ1· · · dξ_d,

przy czym f jest (zwykle skomplikowana) funkcj_, a wypłaty na końcu okresu, a ξ_, j reprezentuja_, czynniki losowe w kolejnych odcinkach czasu. Wymiar d może wynosić nawet kilka tysiecy._, Z podstawowego wykładu z metod numerycznych każdy z nas wie jak numerycznie całkować funkcje jednej zmiennej. Stosowane metody w znakomitej wiekszości przypadków sprowadzaj_, a_,

3

sie do scałkowania funkcji, która jest kawałkami wielomianem określonego stopnia interpo-_, lujacym funkcj_, e podcałkow_, a. Pomysł ten może być uogólniony na przypadek funkcji wielu_, zmiennych. Aby jednak mówić o kwadraturach interpolacyjnych w wielu wymiarach, musimy najpierw zastanowić sie nad rozwi_, azywalności_, a odpowiedniego zadania interpolacyjnego._,

1.2 Interpolacja na siatkach regularnych

1.2.1 Postać wielomianu interpolacyjnego Niech

a ¬ t₁ < t₂ < · · · < t_r¬ b.

Jeśli f jest funkcja jednej zmiennej, f : [a, b] → R, to wielomian_, p_f(x) =

r

X

j=1

f (t_j) l_j(x),

gdzie l_j jest j-tym wielomianem Lagrange’a, lj(x) =

r

Y

j6=i=1

x − tj

t_i− t_j, 1 ¬ j ¬ r,

(przy czym l₁ ≡ 1 dla r = 1) jest stopnia co najwyżej (r − 1) i interpoluje f w punktach t_j, tzn. przyjmuje w tych punktach te same wartości co f . W przypadku d ­ 2 możemy podobnie zdefiniować ‘wielowymiarowe’ wielomiany Lagrange’a.

W tym celu zakładamy, że na każdej współrzednej dany jest przedział, a w nim układ r punktów_, a^(k) ¬ t^(k)₁ < t^(k)₂ < . . . < t^(k)_r ¬ b^(k), 1 ¬ k ¬ d.

Oznaczajac przez l_, ^(k)_j odpowiednie wielomiany Lagrange’a jednej zmiennej dla k-tego podziału, definiujemy wielomiany Lagrange’a d zmiennych jako

l_j1,...,jd(x₁, . . . , x_d) = l⁽¹⁾_j

1 (x₁) l⁽²⁾_j

2 (x₂) · · · l^(d)_j

d (x_d)

dla wszystkich 1 ¬ j_k ¬ r, 1 ¬ k ¬ d. Dla skrócenia zapisu, bedziemy dalej używać zapisu_, wektorowego ~j = (j₁, . . . , j_d), a 1 ¬ ~j ¬ d bedzie oznaczać, że nierówności zachodz_, a dla każdej_, współrzednej j_{, k}, 1 ¬ k ¬ d. Podobnie, t_~j = (t⁽¹⁾_j

1 , t⁽²⁾_j2 , . . . , t^(d)_j

d ).

Wielomiany l_~j należa do przestrzeni P_{, d}^r wielomianów d zmiennych postaci p(~x) = p(x1, . . . , xd) = ^X

0¬~i¬r−1

a_~i· xⁱ₁¹xⁱ₂²· · · xⁱ_d^d,

gdzie a_~i sa dowolnymi wsoółczynnikami rzeczywistymi. Zauważmy, że p ∈ P_{, d}^r wtedy i tylko wtedy gdy p jest wielomianem stopnia co najwyżej (r − 1) ze wzgledu na każd_, a zmienn_, a x_, k. Lemat 1.1. Jeśli wielomian p ∈ P_d^r zeruje sie we wszystkich r_, ^d punktach t_~j, 1 ¬ ~j ¬ r, to p jest wielomianem zerowym.

Dowód. Dowód przeprowadzimy przez indukcje ze wzgl_, edu na wymiar d. Dla d = 1 lemat jest_, oczywiście prawdziwy, bo na podstawie zasadniczego twierdzenia algebry niezerowy wielomian stopnia co najwyżej (r − 1) nie może mieć r różnych zer.

Korzystaj. Nie kopiuj. Nie rozpowszechniaj. 5/27

Niech d ­ 2. Niech a_~j bed_, a współczynnikami wielomianu p. Dla ustalonej k zdefiniujmy wielo-_, mian p_k∈ P_d−1^r jako

p_k(x₁, . . . , x_d−1) = p(x₁, . . . , x_d−1, t^(d)_k ).

Wielomian ten zeruje sie w r(d − 1) punktach t_, i1,...,id−1. Zapisujac go w postaci_, p_k(x₁, . . . , x_d−1) = ^X

0¬i1,...,id−1¬d−1

b^(k)_i1,...,i

d−1· xⁱ₁¹· · · xⁱ_d−1^d−1, gdzie współczynniki

b^(k)_i1,...,i

d−1 = ^X

0¬i_d¬d−1

a_~i· (t^(d)_k )^id, oraz stosujac założenie indukcyjne mamy, że b_, ^(k)_i

1,...,id−1 = 0. A wiec dla wszystkich wyborów_, indeksów i₁, . . . , i_d−1 wielomian jednej zmiennej ^Pd−1_i

d=0a_~i· t^id zeruje sie w r punktach t = t_, ^(d)s . To zaś wymusza a_~i = 0 dla wszystkich 0 ¬ ~i ¬ d − 1 i w konsekwencji p ≡ 0.

Lemat1.1 wykorzystamy do pokazania nastepuj_, acego twierdzenia._,

Twierdzenie 1.1. Wielomiany l_~j, 1 ¬ ~j ¬ r, tworza baz_, e przestrzeni P_{, d}^r. W szczególności, dim(P_d^r) = r^d.

Dowód. Zauważmy, że podobnie jak w przypadku d = 1, l_~j(t_~i) =

( 1, jeśli ~i = ~j,

0, w przeciwnym przypadku.

Stad, jeśli kombinacja liniowa p =_, ^P_~jα_~jl_~j jest wielomianem zerowym to dla wszystkich ~i 0 = p(t_~i) =^X

~j

α_~jl_~j(t_~i) = α_~i,

czyli układ {l_~i : 1 ¬ ~i ¬ r} jest liniowo niezależny. Z drugiej strony, układ ten rozpina P_d^r, bo dla dowolnego wielomianu p z tej przestrzeni mamy

p = ^X

1¬~j¬r

p(t_~j) l_~j. (1.1)

Rzeczywiście, w przeciwnym przypadku różnica wielomianu p i prawej strony (1.1) byłaby nie- zerowym wielomianem w P_d^r, który zeruje sie we wszystkich r_, ^dpunktach t_~i. To zaś przczyłoby lematowi 1.1.

Stad już jeden krok do nast_, epuj_, acego wniosku podsumowuj_, acego nasze dotychczasowe rozwa-_, żania. Niech

D = [a⁽¹⁾, b⁽¹⁾] × [a⁽²⁾, b⁽²⁾] × · · · × [a^(d), b^(d)] bedzie d wymiarowym prostok_, atem._,

Wniosek 1.1. Dla dowolnej funkcji f : D → R wielomian p_f(~x) = ^X

0¬~j¬r

f (t_~j) l_~j(~j)

jest jedynym wielomianem w P_d^r interpolujacym f w punktach t_{, ~j} tzn. takim, że p_f(t_~j) = f (t_~j)

dla wszystkich 1 ¬ ~j ¬ r.

1.2.2 Bład interpolacji_,

Zastanówmy sie teraz jaki jest bł_, ad otrzymanej interpolacji. Dla uproszczenia b_, edziemy od teraz_, zakładać, że D jest kostka d wymiarow_, a, tzn. wszystkie kraw_, edzie maj_, a t_, a sam_, a długość, któr_, a_, oznaczymy przez H, a wezły na każdej współrz_, ednej_,

t^(k)_j = a^(k)+ u_jH, 1 ¬ j ¬ r, gdzie

0 ¬ u₁ < u2 < · · · < ur ¬ 1 jest pewna ustalon_, a siatk_, a na odcinku jednostkowym._,

W przypadku skalarnym, o ile funkcja f jest r-krotnie różniczkowalna w sposób ciagły, to_, f (x) − p_f(x) = (x − t₁)(x − t₂) · · · (x − t_r)f^(r)(ξ)

r! , przy czym ξ ∈ [a, b] zależy od x. Stad w szczególności mamy_,

|f (x) − p_f(x)| ¬ (b − a)^r

r! kf^(r)k_∞, (1.2)

gdzie kf^(r)k_∞= max_a¬t¬b|f^(r)(t)|. Aby wyprowadzić formułe na bład interpolacji w przypadku_, wielowymiarowym, bedziemy potrzebować pewnego prostego uogólnienia ostatniego wzoru._, Załóżmy, że zamiast dokładnych wartości f (t_i) mamy jedynie wartości przybliżone y_i takie, że bład_,

|y_i− f (t_i)| ¬ δ, 1 ¬ i ¬ r. (1.3)

Niech dalej ˜p_f bedzie wielomianem stopnia co najwyżej (r − 1) interpoluj_, acym dane przybliżone_, yi w punktach t_i. Ponieważ (p_f − ˜p_f) jest wielomianem interpolujacym dane f (t_, j) − y_j, na podstawie wzoru (1.1) mamy

|p_f(x) − ˜pf(x)| ¬ δ ·

r

X

i=1

|l_i(x)| ¬ δ · S_r, gdzie S₁= 1, a dla r ­ 2

Sr= max

0¬z¬1 r

X

i=1 r

Y

i6=j=1

z − uj

ui− u_j      .

Stad i z formuły na bł_, ad interpolacji dla dokładnych danych otrzymujemy_,

|f (x) − ˜p_f(x)| ¬ |f (x) − p_f(x)| + |p_f(x) − ˜p_f(x)| ¬ (b − a)^r

r! kf^(r)k_∞ + δ · S_r. (1.4) Wprowadzimy jeszcze klase F_, r(D) funkcji f : D → R, które w całej swojej dziedzinie sa r-krotnie_, różniczkowalne w sposób ciagły ze wzgl_, edu na każd_, a zmienn_, a. Dla f ∈ F_, r(D) definiujemy

B_r(f ) = max

1¬i¬d

(

∂^rf

∂x^r₁ ∞

, . . . ,

∂^rf

∂x^r_{d ∞}

) .

Twierdzenie 1.2. Niech D = [a⁽¹⁾, a⁽¹⁾+ H] × · · · × [a^(d), a^(d)+ H]. Jeśli f ∈ F_r(D) to dla każdego ~x ∈ D bład interpolacji_,

|f (~x) − pf(~x)| ¬ H^r

r! Cr,dBr(f ), gdzie C_1,d= d, a dla r ­ 2

Cr,d = S_r^d− 1 Sr− 1.

Korzystaj. Nie kopiuj. Nie rozpowszechniaj. 7/27

Dowód. Rozpatrzymy tylko r ­ 2 pozostawiajac przypadek r = 1 jako proste ćwiczenie._, Dla d = 1 nierówność w tezie jest równoważna (1.2). Załóżmy wiec, że d ­ 2. Ponieważ dla_, każdego ustalonego t^(d)_k wielomian (d − 1) zmiennych p_f(x₁, . . . , x_d−1, t^(d)_k ) jest wielomianem interpolacyjnym dla funkcji (d − 1) zmiennych f (x₁, . . . , x_d−1, t^(d)_k ), na podstawie założenia in- dukcyjnego mamy

|f (x₁, . . . , x_d−1, t^(d)_k ) − p_f(x₁, . . . , x_d−1, t^(d)_k )| ¬ H^r

r! B_r(f ) S_r^d−1− 1 S_r− 1

!

. (1.5)

Zauważmy, że dla ustalonych z kolei pierwszych (d − 1) współrzednych x_{, 1}, . . . , x_d−1 wielomian p_f(x₁, . . . , x_d−1, t) jest wielomianem jednej zmiennej t interpoluacym funkcj_, e jednej zmiennej_, f (x1, . . . , xd−1, t) w punktach t^(d)_k na podstawie danych zaburzonych na poziomie δ równym prawej stronie (1.5). Stad i z (1.4) ostatecznie otrzymujemy_,

|f (~x) − p_f(~x)| ¬ H^r

r! Br(f ) + δ · S_r

= H^r

r! Br(f ) 1 +S_r^d−1− 1 Sr− 1 Sr

!

= H^r

r! B_r(f ) S_r^d− 1 Sr− 1

! .

1.3 Kwadratury interpolacyjne

1.3.1 Kwadratury proste

Jesteśmy już dobrze uzbrojeni w mechanizm interpolacyjny i możemy zdefiniować wielowymia- rowe kwadratury interpolacyjne dla całkowania funkcji f : D → R zdefiniowanych na kostce

D = [a⁽¹⁾, a⁽¹⁾+ H] × · · · × [a^(d), a^(d)+ H].

Kwadratury te dane sa równości_, a_,

Q_r,d(f ) = Z

D

p_f(~x) d~x, (1.6)

gdzie p_f ∈ P_d^r jest wielomianem interpolujacym f w punktach t_{, ~j}, 1 ¬ ~j ¬ r.

Chociaż postać (1.6) kwadratury znakomicie nadaje sie do rozważań teoretycznych, nie jest_, jednak praktyczna ze wzgledu na obliczenia. Zauważmy, że_,

Q_r,d(f ) = Z

D

X

~j

f (t_~j)l_~j(~x) d~x = ^X

~j

f (t_~j) Z

D

l_~j(~x) d~x

= H^d·^X

~j

f (t_~j)

d

Y

k=1

Z ₁

0

l_jk(u) du

 ,

gdzie l_j jest j-tym wielomianem Lagrange’a dla punktów u₁, u2, . . . , u_d. Stad, wprowadzaj_, ac_, oznaczenie

a_k = Z ₁

0

l_k(u) du,

kwadrature interpolacyjn_, a można zapisać w postaci_, Q_r,d(f ) = H^d· ^X

1¬j1,...,jd¬r

a_j1a_j2· · · a_jd· f (t⁽¹⁾_j

1 , t⁽²⁾_j

2 , . . . , t^(d)_j

d ).

Zauważmy, że a_k sa współczynnikami jednowymiarowej kwadratury interpolacyjnej Q_, r(f ) = Pr

k=1a_kf (t_k) opartej na punktach u_k, przybliżajacej całk_, e_, ^R₀¹f (u) du. Mówiac inaczej, zdefi-_, niowana przez nas wielowymiarowa kwadratura interpolacyjna jest d-produktem tensorowym wybranej kwadratury jednowymiarowej.

Na koniec tego podrozdziału podamy oszacowanie błedu kwadratury Q_{, r,d}. Ponieważ Z

D

f (~x) d~x − Q_r,d(f ) = Z

D

f (~x) − p_f(~x)
d~x, z twierdzenia 1.2natychmiast otrzymujemy nastepuj_, acy wniosek._,

Wniosek 1.2. Jeśli f ∈ Fr(D) to bład kwadratury interpolacyjnej Q_, r,d jest ograniczony przez 

   Z

D

f (~x) d~x − Qr,d(f )    

¬ H^r+d

r! Cr,dBr(f ).

1.3.2 Kwadratury złożone

Podobnie jak w przypadku funkcji jednej zmiennej, definiujemy kwadratury złożone dla funkcji wielu zmiennych. Dla uproszczenia zakładamy, że całkujemy po kostce jednostkowej [0, 1]^d. Dla danego n wprowadzamy podział kostki na n^d podkostek

i₁− 1 n ,i₁

n



×

i₂− 1 n ,i₂

n



× · · ·

i_d− 1 n ,i_d

n



, 1 ¬ i_k¬ n, 1 ¬ k ¬ d.

Nastepnie na każdej podkostce stosujemy prost_, a kwadratur_, e interpolacyjn_, a opart_, a na siatce_, regularnej składajacej si_, e z r_, ^dpunktów. Skonstruowana w ten sposób kwadratur_, e złożon_, a ozna-_, czymy przez Q⁽ⁿ⁾_r,d.

Przyk?ad 1.2. Jeśli bazowa kwadratur_, a jednowymiarow_, a jest reguła punktu środkowego,_, Q₁(f ) = (b − a) · f

a + b 2

 ,

to wynikowa kwadratur_, a złożon_, a na [0, 1]_, ^d jest po prostu reguła prostokatów_, Q⁽ⁿ⁾_r,d(f ) =

1 n

d

· ^X

1¬i1,,...,id¬n

f

i₁− 1/2

n , · · · ,i_d− 1/2 n

 .

Nasze rozważania wieńczy twierdzenie o błedzie kwadratury złożonej, które natychmiast wynika_, z wniosku1.2oraz sposobu konstrukcji kwadratury.

Twierdzenie 1.3. Kwadratura złożona Q⁽ⁿ⁾_r,d korzysta z co najwyżej N = (r n)^d

wartości funkcji f . Jeśli f ∈ F_r([0, 1]^d) to jej bład_, 

    Z

[0,1]^d

f (~x) d~x − Q⁽ⁿ⁾_r,d(f )     

¬

1 N

r/dr^r r!



C_r,dB_r(f ).

Korzystaj. Nie kopiuj. Nie rozpowszechniaj. 9/27

1.4 Przekleństwo wymiaru

Złożone kwadratury interpolacyjne moga być z powodzeniem stosowane dla niskich wymiarów,_, powiedzmy d = 2, 3. Dla dużych wymiarów d maja one bowiem t_, a niepoż_, adan_, a własność, że_, liczba wezłów rośnie wykładniczo szybko wraz z zag_, eszczaniem siatki. Nawet jeśli weźmiemy_, po 2 punkty na każdej współrzednej to całkowita liczba punktów siatki regularnej wyniesie_, 2^d. Pamietamy, że w wielu praktycznych zastosowaniach d może si_, egać nawet kilku tysi_, ecy. W_, takich przypadkach obliczenie wartości kwadratury jest zadaniem praktycznie niewykonalnym.

To jednak nie koniec złych wiadomości. Przyjrzyjmy sie jeszcze bł_, edowi złożonej kwadratury_, interpolacyjnej. Twierdzenie 1.3 mówi, że bład ten jest ograniczony z góry proporcjonalnie_, do N^−r/d, gdzie N jest liczba wszystkich użytych punktów. To drugi powód do niepokoju,_, uzasadniony poniższym przykładem.

Przyk?ad 1.3. Załóżmy, że chcemy całkować funkcje 360 zmiennych i jako kwadratur_, e bazow_, a_, stosujemy kwadrature Simpsona, dla której r = 4. Górne ograniczenie bł_, edu sugeruje, że aby być_, pewnym wyniku z dokładnościa 10_, ⁻² to musimy obliczyć wartości funkcji w aż 10¹⁸⁰ punktach.

Czy naprawde jest aż tak źle?_,

Rzeczywiście jest tak źle, a nawet gorzej. Okazuje sie, że rz_, edu zbieżości N_, ^−r/d nie da sie_, poprawić w klasie funkcji F_r([0, 1]^d). Mówi o tym nastepuj_, ace twierdzenie._,

Twierdzenie 1.4. Istnieje c = c_r,d > 0 o nastepuj_, acej własności. Dla dowolnej aproksyma-_, cji całki wykorzystujacej N wartości funkcji można znaleźć funkcj_, e f ∈ F_, r([0, 1]^d) dla której Br(f ) = 1, a bład aproksymacji całki wynosi co najmniej c N_, ^−r/d.

Dowód. Załóżmy, że dana aproksymacja całki oblicza wartości funkcji w punktach ~t_j, 1 ¬ j ¬ N . Dowód twierdzenia polega na konstrukcji dwóch funkcji, f₊i f₋, które zeruja si_, e we wszystkich_,

~t_j (a tym samym ich całki sa aproksymowane t_, a sam_, a liczb_, a), dla których B_{, r}(f₊) = 1 = B_r(f₋), ale różnica całek

Z

[0,1]^d

(f₊− f₋)(~t) d~t ­ 2c N^−r/d,

dla pewnej c niezależnej od f i d. Wtedy, przynajmniej dla jednej z tych funkcji bład aproksy-_, macji całki wynosi co najmniej cN^−r/d.

Wybierzmy n taka, że n_, ^d > N i skonstruujmy na [0, 1]^d regularna siatk_, e składaj_, ac_, a si_, e z n_, ^d kostek, każda o krawedzi długości h = 1/n._,

Niech dalej φ : R → R bedzie dowoln_, a funkcj_, a r-krotnie różniczkowaln_, a w sposób ci_, agły speł-_, niajac_, a nast_, epuj_, ace warunki:_,

1. φ(x) = 0 dla x /∈ (0, 1),

2. φ^(j)(0) = 0 = φ^(j)(1) dla 0 ¬ j ¬ r, 3. ^R₀¹φ(t) dt =: a > 0.

Każdej kostce

K_~i := [(i₁− 1)h, i₁h] × · · · × [(id− 1)h, i_dh]

naszej regularnej siatki przyporzadkujemy funkcj_, e_, φ_~i(x₁, . . . , x_d) := h^r

d

Y

k=1

φ(x_k/h − i_k).

Zauważmy, że B_r(φ_~i) = 1 oraz Z

K~i

φ_~i(~t) d~t = a^dh^r+d.

Jasne jest, że istnieje co najmniej n^d− N multi-indeksów ~i (kostek) takich, że żaden z punktów

~t_j nie należy do wnetrza K_{, ~i}. Oznaczmy zbiór takich indeksów przez S i zdefiniujmy funkcje f₊:=^X

~i∈S

φ_~i, f−:= −f₊.

Wtedy obie funkcje zeruja si_, e w ~t_, j, B_r(f₊) = 1 = B_r(f−), oraz Z

[0,1]^d

f₊(~t) d~t = − Z

[0,1]^d

f−(~t) d~t ­ (n^d− N ) a^dh^r+d.

Podstawiajac n = dN_, ^1/d(1 + d/r)^1/de dostajemy ostatecznie Z

[0,1]^d

f₊(~t) d~t ­ c N^−r/d, gdzie c = da^d

r2^r+d(1 + d/r)^1+r/d.

Opisane zjawisko nosi nazwe przekleństwa wymiaru._,

Rozdział 2

Metody Monte Carlo

2.1 Wst ep, metody niedeterministyczne

Poprzedni wykład zakończyliśmy pesymistycznym twierdzeniem 1.4, że nie istnieja efektywne_, metody numerycznego całkowania funkcji wielu zmiennych, ponieważ ma miejsce zjawisko prze- kleństwa wymiaru. Zwróćmy jednak uwage na to, że fakt istnienia przekleństwa wymiaru stwier-_, dziliśmy przy założeniach, że:

(i) model obliczeniowy jest deterministyczny,

(ii) funkcje podcałkowe sa r-krotnie różniczkowalne po każdej zmiennej._,

Można mieć nadzieje, że przekleństwo wymiaru zniknie, albo zostanie złagodzone, gdy przynaj-_, mniej jedno z tych założeń nie bedzie spełnione._,

Ten wykład poświecimy (klasycznej) metodzie Monte Carlo numerycznego całkowania, która_, jest przykładem metody niedeterministycznej, tzn. takiej, która oblicza wynik wykorzystujac_, zjawiska losowe. Chociaż może to brzmieć dziwnie, to właśnie niedeterministyczne zachowanie metody pozwala pokonać przekleństwo wymiaru.

Opisana dalej klasyczna metoda Monte Carlo zwiazana jest ściśle ze Stanisławem Ulamem,_, uczniem Stefana Banacha i reprezentantem Lwowskiej Szkoły Matematycznej. Ulam zastosował metode Monte Carlo do obliczania skomplikowanych całek w ramach ”Manhattan Project”w_, Los Alamos (USA), w czasie II Wojny Światowej.

2.2 Klasyczna metoda Monte Carlo

2.2.1 Definicja i bład_,

Tak jak w poprzednim rozdziale chcemy obliczyć całke_, I_d(f ) :=

Z

[0,1]^d

f (~x) d~x = Z 1

0

Z 1 0

· · · Z 1

0

| {z }

d

f (x₁, x₂, . . . , x_d) dx₁dx₂· · · dx_d.

Zakładamy przy tym, że f : [0, 1]^d→ R jest funkcja, której kwadrat jest całkowalny,_, Z

[0,1]^d

|f (~x)|²d~x < ∞.

11

Definicja 2.1. Klasyczna metoda Monte Carlo polega na przybliżeniu I_d(f ) średnia arytme-_, tyczna wartości funkcji f w losowo wybranych punktach, tzn._,

M Cd,N(f ) = M C_d,N(f ; ~t₁, ~t2, . . . , ~tN) := 1 N ·

N

X

j=1

f (~tj),

gdzie ~t₁, ~t₂, . . . , ~t_N sa punktami wylosowanymi niezależnie od siebie, każdy zgodnie z rozkładem_, jednostajnym na [0, 1]^d.

Konsekwencja zastosowania losowości jest to, że przy różnych realizacjach metody otrzymujemy_, różne wyniki, w zależności od wyboru punktów ~t_j. Wynik M C_d,N(f ) jest wiec zmienn_, a losow_, a,_, której wartość oczekiwana wynosi

E (M C_d,N(f )) = Z

[0,1]^d·N

M C_d,N(f ; ~t₁, . . . , ~tN) d~t₁· · · d~t_N

= 1

N

N

X

j=1

Z

[0,1]^d

f (~t) d~t = I_d(f ).

Ponieważ różnica I_d(f ) − M C_d,N(f ) jest też zmienna losow_, a, za bł_, ad metody Monte Carlo dla_, danej funkcji f przyjmiemy odchylenie standardowe,

e(f ; M C_d,N) :=

q

E (I_d(f ) − M C_d,N(f ))². Twierdzenie 2.1. Dla danej funkcji f bład metody Monte Carlo wynosi_,

e(f ; M C_d,N) = σ(f )

√ N, gdzie

σ(f ) := ^qI_d(f²) − I_d²(f ) jest wariancja funkcji f ._,

Zanim przystapimy do dowodu zauważmy, że σ(f ) jest dobrze zdefiniowan_, a wielkości_, a, bowiem_, nierówność

|I_d(f )| ¬ q

I_d(f²)

jest szczególnym przypadkiem znanej nierówności Schwarza dla całek.

Dowód. Oznaczmy, dla uproszczenia, zmienna losow_, a X = M C_, d,N(f ). Wtedy

E(X − E(X))² = E (X(X − E(X)) − E(X)(X − E(X))) = E(X²) − E²(X). (2.1) Ponadto

E(X²) = E





 1 N





N

X

j=1

f (~tj)





2



 = 1 N²E





N

X

j=1

f²(~t_j) +^X

i6=j

f (~ti)f (~t_j)





= 1

N²

N I_d(f²) + (N²− N )I_d²(f )^ = 1

NI_d(f²) +

 1 − 1

N

 I_d²(f ),

gdzie skorzystaliśmy z niezależności zmiennych losowych f (t_i) i f (t_j) dla i 6= j. Stad i z (2.1)_, dostajemy

e²(f ; M C_d,N) = E(X − E(X))² = 1

NId(f²) +

 1 − 1

N



I_d²(f ) − I_d²(f ) = 1 N

Id(f²) − I_d²(f )^, co kończy dowód.

Korzystaj. Nie kopiuj. Nie rozpowszechniaj. 13/27

Uwaga 2.1. Zauważmy, że w dowodzie pokazaliśmy przy okazji nierówność Schwarza posługujac_, sie narz_, edziami rachunku prawdopodobieństwa._,

Twierdzenie (2.1) mówi, że bład metody Monte Carlo jest proporcjonalny do N_, ^−1/2przy bardzo słabych wstepnych założeniach na funkcj_, e (jedynie całkowalność kwadratu funkcji). Jest to_, istotna poprawa w porównaniu do błedu N_, ^−r/d dla metod deterministycznych. W szczególności ważne jest, że wykładnik 1/2 przy N⁻¹ jest niezależny od wymiaru d, a konsekwencja tego_, pokonanie przekleństwa wymiaru.

Dziwnym może wydawać sie, że przekleństwo wymiaru można zlikwidować używaj_, ac metod nie-_, deterministycznych (losowych). Jednak niczego nie ma za darmo. Należy pamietać, że jest to_, możliwe za cene niepewności wyniku. O ile bowiem metoda deterministyczna produkuje zawsze_, ten sam wynik, metoda niedeterministyczna (taka jak Monte Carlo) produkuje różne wyniki zależnie od konkretnych realizacji zmiennych losowych. Dlatego, mimo iż bład oczekiwany jest_, proporcjonalny do N^−1/2to nie mamy całkowitej pewności, że przy konkretnej realizacji otrzy- many wynik jest tego samego rzedu. Z tego punktu widzenia warto przytoczyć nast_, epuj_, ac_, a_, równość, która wynika z centralnego twierdzenia granicznego; mianowicie, dla dowolnych c₁< c2

mamy

N →∞lim Prob

c1σ(f )

√

N ¬ I_d(f ) − M C_d,N(f ) ¬ c2σ(f )

√ N



= 1

√ 2π

Z c2

c1

e^−t2/2dt,

gdzie Prob oznacza prawdopodobieństwo wzgledem rozkładu jednostajnego na [0, 1]_, ^dN.

2.2.2 Całkowanie z waga_,

Deterministyczne metody interpolacyjne z poprzedniego rozdziału można stosować jedynie do całkowania na d-wymiarowych prostokatach. Metoda Monte Carlo ma oprócz wymienionych_, również i ta zalet_, e, że łatwo j_, a uogólnić na przypadek całkowania z wag_, a. Dla przybliżenia_, wartości

Id,ω(f ) :=

Z

R^d

f (~x)ω(~x) d~x, gdzie Z

R^d

ω(~x) d~x =: W < ∞, możemy bowiem zastosować wzór

M C_d,N^ω (f ) := W N ·

N

X

j=1

f (~tj),

przy czym ~t₁, . . . , ~tN sa tym razem punktami wybranymi losowo i niezależnie od siebie, zgodnie_, z rozkładem na R^d o gestości ω/W ._,

Adaptujac odpowiednio dowód twierdzenia_, 2.1 otrzymujemy nastepuj_, ace wyrażenie na bł_, ad_, uogólnionej metody Monte Carlo.

Twierdzenie 2.2. Niech Id,ω(f²) =^R

R^df²(~x)ω(~x) d~x < ∞. Wtedy e(f ; M C_d,N^ω ) = σ_ω(f )

√N , gdzie

σ_ω(f ) = ^qW I_d,ω(f²) − I_d,ω² (f ).

2.3 Redukcja wariancji

Zauważyliśmy, że zaleta metody Monte Carlo jest nie tylko jej prostota, ale również to, że bł_, ad_, średni wynosi σ_ω(f )N^−1/2. Naturalnym jest teraz pytanie, czy błedu tego nie można poprawić._, Temu celowi służa metody redukcji wariancji, które polegaj_, a w ogólności na redukcji czynnika_, σ_ω(f ). Spośród wielu technik redukcji wariancji skupimy uwage na dwóch: losowaniu warstwo-_, wemu oraz funkcjach kontrolnych. Dla uproszczenia bedziemy zakładać, że całkujemy z wag_, a_, jednostkowa na kostce_,

D = [0, 1]^d. 2.3.1 Losowanie warstwowe

Podzielmy obszar całkowania D na K rozłacznych podzbiorów D_{, i} tak, że D =

K

[

i=1

Di

i zastosujmy Monte Carlo do całkowania po każdym D_i, tzn. całke_, ^R_Df (~x) d~x przybliżymy wielkościa_,

M Cd,N(f ) :=

K

X

i=1

M C_d,N⁽ⁱ⁾

i(f ), gdzie M C_d,N⁽ⁱ⁾

i jest metoda Monte Carlo zastosowan_, a do całki_, I_d⁽ⁱ⁾(f ) :=

Z

Di

f (~x) d~x, 1 ¬ i ¬ K, oraz N =^PK_i=1N_i.

Oznaczmy przez |D_i| objetość d-wymiarow_, a podzbioru D_{, i}. Ponieważ zmienne losowe I_d⁽ⁱ⁾(f ) − M C_d,N⁽ⁱ⁾

i(f ) sa parami niezależne dla 1 ¬ i ¬ K, na podstawie twierdzenia_, 2.2mamy E^(I_d(f ) −M C_d,N(f ))^2 = E





K

X

i=1

I_d⁽ⁱ⁾(f ) − M C_d,Ni(f )

!²



=

K

X

i=1

E^(I_d⁽ⁱ⁾(f ) − M C_d,N⁽ⁱ⁾

i(f ))^2

=

K

X

i=1

1 N_i

|D_i| I_d⁽ⁱ⁾(f²) − (I_d⁽ⁱ⁾(f ))^2.

Przyjmijmy teraz, że

N_i = |D_i| · N, 1 ¬ i ¬ K,

przy czym dla uproszczenia (ale bez utraty ogólności) zakładamy, że wielkości te sa całkowite._, Wtedy otrzymujemy

e(f, M C_d,N) = 1

√N · v u u

tI_d(f²) −

K

X

i=1

1

|D_i|(I_d⁽ⁱ⁾(f ))². (2.2) Bład tak zdefiniowanej metody M C_{, d,N} nie jest wiekszy od bł_, edu klasycznej metody M C_{, d,N} z Twierdzenia 2.1.

Korzystaj. Nie kopiuj. Nie rozpowszechniaj. 15/27

Twierdzenie 2.3. Dla dowolnej funkcji f takiej, że I_d(f²) < ∞ mamy e(f, M C_d,N) ¬ e(f, M C_d,N),

przy czym równość zachodzi tylko wtedy gdy iloraz I_d⁽ⁱ⁾(f )/|D_i| jest stały, niezależnie od i.

Dowód. Rzeczywiście, oznaczajac_,

ai = q

|D_i|, bi= I_d⁽ⁱ⁾(f ) p|D_i|, oraz wykorzystujac nierówność Schwarza dla ci_, agów mamy_,

I_d²(f ) =

K

X

i=1

aibi

!2

¬

K

X

i=1

a²_i

! _K X

i=1

b²_i

!

=

K

X

i=1

b²_i =

K

X

i=1

1

|D_i|(I_d⁽ⁱ⁾(f ))²,

przy czym równość zachodzi tylko wtedy gdy wektory (a₁, . . . , aK)^T i (b₁, . . . , bK)^T sa liniowo_, zależne, co jest równoważne warunkowi w treści twierdzenia.

Prawdziwość tezy pokazuje teraz porównanie wzorów na błedy obu metod._,

Widzimy, że stosujac losowanie warstwowe z ustalonym podziałem na K podzbiory D_{, i} możemy co prawda zmiejszyć bład, ale szybkość zbieżności N_, ^−1/2 pozostaje ta sama. A czy można poprawić zbieżność stosujac różne podziały dla różnych wartości N ? Okazuje si_, e, że tak, o ile_, założymy pewna regularność funkcji f ._,

Aby to uzyskać, najpierw przekształcimy wzór (2.2) na bład metody M C_{, d,N} do postaci

e(f, M C_d,N) = 1

√ N

v u u t

K

X

i=1

I_d⁽ⁱ⁾((f − c_i)²). (2.3)

gdzie

c_i := I_d⁽ⁱ⁾(f )

|D_i| , 1 ¬ i ¬ K.

Załóżmy teraz, że f spełnia warunek Lipschitza ze stała L,_,

|f (~x) − f (~y)| ¬ L · k~x − ~yk∞, ~x, ~y ∈ D.

Wtedy istnieja ~_,x_i∈ D_itakie, że f (~x_i) = c_i, a stad i z lipschitzowskości f mamy, że dla dowolnego_,

~ x ∈ D_i

|f (~x) − ci| ¬ L · k~x − ~xik_∞ ¬ L · diam_∞(D_i), gdzie

diam∞(D_i) := sup {k~x − ~yk∞ : ~x, ~y ∈ Di}

jest średnica zbioru D_, i w normie max. W konsekwencji, ze wzoru (2.3) dostajemy nastepuj_, ace_, oszacowanie błedu:_,

e(f, M C_d,N) ¬ L

√ N

v u u t

K

X

i=1

|D_i| diam²(D_i).

Ustalmy teraz równomierny podział kostki D na K = N podkostek D_i, każda o krawedzi_, długości N^−1/d (zakładamy, bez zmniejszenia ogólności, że N^1/d jest całkowita) tak, że nasza