Całkowanie w czterech wymiarach metodą Monte Carlo
Tomasz Chwiej 20 stycznia 2014
1 Opis problemu
Naszym zadaniem jest numeryczne wyznaczenie wartości poniższej całki:
V =
∫∫ ∞
−∞d2⃗r1d2⃗r2ρ1(⃗r1)ρ2(⃗r2) r12
(1) gdzie:
ρ1(⃗r1) = 1 2πσ2exp
(
−(⃗r1− ⃗R10)2 2σ2
)
= 1
2πσ2exp (
−(x1− X10)2+ (y1− Y10)2 2σ2
)
(2)
ρ2(⃗r2) = 1 2πσ2exp
(
−(⃗r2− ⃗R20)2 2σ2
)
= 1
2πσ2exp (
−(x2− X20)2+ (y2− Y20)2 2σ2
)
(3) oraz
r12=|⃗r1− ⃗r2| =√(x1− x2)2+ (y1− y2)2 (4) Zauważmy że funkcje ρ1i ρ2odpowiadają iloczynom Gaussowskich rozkładów prawdopodobieństwa
ρi(x, y)→ N(Xi0, σ)· N(Yi0, σ)
wobec czego, funkcją podcałkową jest tylko wyraz 1/r12. Całkę V można więc łatwo obliczyć meto- dą Monte Carlo jeśli założymy, że współrzędne położeniowe (x1, y1, x2, y2) mają rozkład normalny N (µ, σ).
x' z' y'
x z y
r0
Rysunek 1: Funkcja ρ2(⃗r) jest przesunięta względem ρ1(⃗r) o wektor ⃗r0. Dokładną wartość całki V można wyznaczyć analitycznie:
Vdok =
√π 2σexp
(
− r02 8σ2
) I0
( r20 8σ2
)
(5) gdzie: I0(x) jest modyfikowaną funkcją Bessel’a pierwszego rodzaju (jej wartość można wyznaczyć stosując funkcję bessi0(float x) z Numerical Recipes), a r0 jest odległością pomiędzy środkami gaussianów r0 =| ⃗R10− ⃗R20|.
1
2 Zadania do wykonania
1. Przyjmujemy σ = 1, wtedy do generowania liczb losowych o rozkładzie N (0, 1) możemy wyko- rzystać schemat Boxa-Millera:
U1, U2 ∈ (0, 1] (6)
Z1 = √−2lnU1cos(2πU2) (7)
Z2 = √−2lnU1sin(2πU2) (8)
Z1, Z2 ∈ N(0, 1) (9)
Uwaga: dla σ̸= 1 zmienne losowe Z1 i Z2 trzeba przeskalować - przemnożyć przez σ.
2. Należy wyznaczyć wartość numeryczną całki V metodą MC dla liczby losowań n = 10i, i = 2, 3, 4, 5 zakładając, że: ⃗R10 = (0, 0), ⃗R20= (x20, 0). Dla każdego n wartość x20 będziemy zmie- niać w zakresie od 0.0 do 6.0 z krokiem ∆x = 0.1. Numeryczną wartość całki liczymy tak:
Vn= 1 n
∑n i=1
ci (10)
gdzie: ci jest wartością funkcji podcałkowej dla współrzędnych losowych ⃗r1 = (x1, y1) i ⃗r2 = (x2, y2) w i-tym losowaniu. Współrzędne losujemy zgodnie z rozkładem N (0, 1). Należy pamię- tać, aby po wylosowaniu liczb pseudolosowych do x2dodać przesunięcie X20, wtedy otrzymamy rozkład przesunięty N (0, 1)→ N(X20, σ) - skalować nie trzeba bo σ = 1 Błąd wartości całki σV¯:
σV¯ =
√ σ2
n (11)
σ2 = 1 n
∑n
i=1
c2i − 1 n
( n
∑
i=1
ci )2
(12)
3. Osobno dla każdej wartości n, na jedym rysunku umieścić wykresy Vdok = f (r0) oraz Vn= f (r0) wraz z zaznaczonym błędem całkowania.
W gnuplocie wykres z błędem rysujemy przy użyciu polecenia:
plot ’100000.dat’ u 1:2:3 w yerrorbars
gdzie: 1- kolumna to współrzedna x20, 2- kolumna to wartość całki, 3 -kolumna to obliczny błąd.
4. W sprawozdaniu proszę dokonać analizy wyników oraz skomentować problem osobliwości funkcji podcałkowej.
