1

Obliczenia transportu

promieniowania metodą

Monte Carlo

Podstawy

2

Podstawy

• Metoda Monte Carlo – dział

matematyki – wykorzystujące

stochastyczne rozwiązywanie równań

• Podejście eksperymentalne (“Metoda

prób statystycznych”)

• Podejście Monte Carlo w celu estymacji

fizycznie mierzalnych wielkości rozbija

się na dwa kroki:

1.

Opracowanie eksperymentu numerycznego

którego wartość oczekiwana odpowiada

fizycznej wartości mierzonej wielkości -

2.

Wykonanie obliczeń – uzyskanie estymacji

mierzonej wielkości.

x

xˆ

3

Podstawy

• Pierwszy krok – trudność zależy od komplikacji

systemu fizycznego

• Jeśli system fizyczny jest z natury

stochastyczny – pierwszy etap jest oczywisty –

symulacja matematyczna winna być odbiciem

zwierciadlanym fizycznej sytuacji. Nazywamy

to symulacją analogową – obliczenia są

idealnym analogiem sytuacji fizycznej

• Transport promieniowania cząstek neutralnych

jest takim typem sytuacji stochastycznej

• CO należy zrobić? – odzwierciedlić w modelu

stochastyczne decyzje podejmowane przez

przyrodę – tak aby można było je „mierzyć”

4

Podstawy

• W procesach które nie są inherentnie

stochastyczne, model eksperymentu

numerycznego jest bardziej złożony i

wymaga:
1.

Wyprowadzenia równania (np.

transportu ciepła, Boltzmana) z rozwiązania

którego można wyciągać wnioski odnośnie

estymowanych wielkości.

2. Rozwinąć i zastosować metodę Monte

Carlo do rozwiązania równania.

5

Pierwszy przykład: obliczyć 

• Zadanie: przeprowadzić numeryczną estymację

 metodą tarczy strzelniczej

• Znamy stosunek powierzchni koła do kwadratu

• Na bazie tego można zaprojektować

eksperyment dostarczający wartości

oczekiwanej liczby 

 

kolo}

w

ie

Pr{trafien

4

2

2

2









r

r

6

przykład

1. Wylosować punkt wewnątrz kwadratu:

A. Wybrać losowo liczbę pomiędzy -1 i 1 dla

położenia we współrzędnej x

B. Wybrać tak samo dla współrzędnej y

2. Sprawdzić czy trafiono koło i zliczać wyniki:
trafienie

s=4 (dlaczego 4?)

pudło: s=0

3. Przeliczyć eksperyment wiele razy. Końcowa

estymacja wyliczona ze wzoru:

1

2

2



 y

x

N

s

N

i

i







1

ˆ



7

Podstawowe spojrzenie na

proces MC

Podstawowe spojrzenie na

proces MC

• Proces Monte Carlo to „czarna skrzynka” –

wchodzi sekwencja liczb losowych

• Wychodzi sekwencja wartości estymowanych

• Czasami wartości estymowane są przybliżone ale

przy wydłużeniu ciągu próba staje się

wystarczająco dobra.

• Potrzebna jest statystyczna ocena próby

8

Formuły statystyczne

Formuły statystyczne

• Trzy najważniejsze formuły

statystyczne:

– Estymowana wartość oczekiwana
– Estymowana wariancja próby
– Estymowana wariancja wartości

oczekiwanej

– Należy je dobrzy rozróżniać

9

Estymator wartości

oczekiwanej

Estymator wartości

oczekiwanej

• Najlepszym estymatorem jest wartość

średnia

• Zasadnicza cecha – jest nieobciążony –

spełnia centralne twierdzenie graniczne (z

N dąży do wartości prawdziwej x)

xˆ

N

x

x

N

i

i







1

ˆ

10

Wartości średnie z rozkładów

ciągłych

Wartości średnie z rozkładów

ciągłych

• Rozkład ciągły f(x), w przedziale (a,b) (tzn:

x=a to x=b).

• Prawdziwa wartość średnia :

• Założyliśmy że f(x) jest prawdziwym

rozkładem prawdopodobieństwa tj.

spełniającym warunki:

x

 

dx

x

f

x

x

b

a





 

1

)

(

b)

(a,

przedziale

w

,

0

0









dx

x

f

x

f

11

Wartość średnia z rozkładów

dyskretnych

Wartość średnia z rozkładów

dyskretnych

• W przypadku rozkładów dyskretnych z M

wartościami x

i

o różnych

prawdopodobieństwach

• Wartość średnia:

• Warunki na prawdopodobieństwa:

i

p







M

i

i

i

x

p

x

1

1

,

0

1









M

i

i

i

p

p

12

przykład: problem 

przykład: problem 

• Zadanie obliczenia , zawiera

rozkład dwumianowy (dwie
odpowiedzi)

• odpowiedź 1 = trafienie w koło:

• odpowiedź 2 = pudło:

• Wartość oczekiwana:



















0

2146

.

0

4

7854

.

0

1

M

i

i

i

x

p

7854

.

0

4

/

4

1

1









p

x

2146

.

0

4

/

1

0

2

2











p

x

13

Estymowana wariancja

próby

Estymowana wariancja

próby

• Wariancja – oczekiwany średni błąd

kwadratowy

• Estymacja wariancji w MC:

 

 





















M

i

i

i

b

a

x

x

p

dx

x

x

x

f

x

1

2

2

2



 

 



























































N

i

N

i

i

N

i

i

i

n

N

x

N

x

N

N

N

x

x

x

S

x

1

2

1

1

2

2

2

2

1

1

ˆ

14

Odchylenie standardowe

próby

Odchylenie standardowe

próby

• Odchylenie standardowe –

pierwiastek z wariancji.

• Analogicznie dla wartości

estymowanych:

• Z własności statystycznych rozkładu

próby wynika możliwość
występowania różnic pomiędzy
indywidualnymi estymacjami.

 

 

n

n

x

S

x

S

2



15

przykład: problem 

przykład: problem 

• Wykorzystując poprzednie dane:

• Jest to rozkład bardzo nie-normalny

 









 

642

.

1

697

.

2

697

.

2

2646

.

0

4

7184

.

0

2

2

1

2

2

































x

x

x

p

x

M

i

i

i

16

Estymowana wariancja wartości

oczekiwanej

Estymowana wariancja wartości

oczekiwanej

• Znaczenie: miara zaufania do

otrzymanej wartości średniej z
pomiaru N próbek

• Można formułować szereg

estymatorów a następnie podawać
je analizie statystycznej

• Przyjmujemy rozwiązania

standardowe z teorii
prawdopodobieństwa

17

wariancja wartości

oczekiwanej

wariancja wartości

oczekiwanej

• wariancja wartości oczekiwanej

wykorzystująca wariancję rozkładu
prawdziwego:

• Musimy jednak wariancję rozkładu

prawdziwego zastąpić jej estymatorem z próby:

 

 

N

x

x

2

2

ˆ







 

 

 

 

 

 

N

x

S

x

S

x

N

x

S

x

S

x

i

i













ˆ

ˆ

ˆ

ˆ

2

2

2

18

przykład: problem 

przykład: problem 

• Nasz przykład dla próby N=10,000 –

wartość odchylenia standardowego

średniej wynosi:

• - WIDZIMY podstawową zasadę

procesu Monte Carlo :

zwiększenie dokładności o jeden rząd

wielkości wymaga zwiększenia 100

krotnego ilości symulowanych historii

(co najmniej)

 

 

 

01642

.

0

000

,

10

642

.

1

ˆ

ˆ











N

x

S

x

S

x

i

2

1



N

19

Liczby pseudolosowe

• Konieczność wykorzystania liczb losowych

w MC

• Co to są liczby pseudolosowe
• Co lepsze liczby losowe czy pseudolosowe
• Liniowe generatory przystające

(kongruentne)

• Wydłużanie okresu generatora – układy

kombinowane generatorów

• Generatory quasi-losowe

20

Ważne reguły

• Nie dysponujemy prawdziwie losowymi

liczbami (i tak naprawdę nie chcemy ich)

• Liczby losowe z komputera nie są prawdziwie

losowe (gdyż są powtarzalne)

• Powtarzalność jest bardzo istotna dla potrzeb

walidacyjnych. Także jest konieczna w

przypadku niektórych metod perturbacyjnych

MC.

• Najbardziej powszechne generatory to liniowe

generatory przystawalne (kongruentne) (LGP)-

generują serie powiązanych liczb całkowitych:

m

i

b

ai

i

n

n

n

n

1

1

1

m

mod

|

)

(















21

możliwe problemy

•

Periodyczność: ciąg liczb ma skończoną długość, a
wiec w końcu nastąpią powtórzenia – grozi złamaniem
podstawowej zasady MC (zwiększania precyzji ze
zwiększaniem próby)

•

Błędny wybór liczb startowych a,b, i m może
prowadzić do powstania korelacji pomiędzy kolejnymi
liczbami (bez uzasadnienia losowego. Np. po bardzo
małej liczbie zawsze następuje inna bardzo mała
liczba)

•

Można pokazać że pełen okres m będzie pokryty przy
wyborze:

–

a(mod 4)=1

–

m = potęga 2

–

b - nieparzyste.

•

Inne problemy: wpływ na pokrycie zagadnienia,
losowość rozkładu cyfr (w zależności od pozycji)

22

Przykład: m=8, a=5, b=3,

i0=1

#

i

ai+b mod

m



0

1

8

0

0

1

0

3

3

.375

2

3

18

2

.25

3

2

13

5

.625

4

5

28

4

.5

5

4

23

7

.875

6

7

38

6

.75

7

6

33

1

.125

23

Wydłużanie okresu generatora

• Nawet przy pełnym pokryciu okresu przez generator

liczb losowych LGP mamy problemy z okresem:

– Okres jest nie większy niż m
– M jest liczba typu INEGER w systemie

komputerowym a:

gdzie N jest liczbą bitów w typie INTEGER (zwykle 16

lub 32)

• Rozwiązanie problemu: wydłużanie okresu generatora –

liniowa kombinacja nieskorelowanych generatorów LGP

• Okres jest iloczynem okresów składowych

1

2





N

m





N



















2

1

FRAC

24

KENO random number

generator

DOUBLE PRECISION FUNCTION FLTRN()
INTEGER X, Y, Z, X0, Y0, Z0, XM, YM, ZM
DOUBLE PRECISION V, X1, Y1, Z1
EXTERNAL RANDNUM
COMMON /QRNDOM/ X, Y, Z
SAVE /QRNDOM/
PARAMETER ( MX = 157, MY = 146, MZ = 142 )
PARAMETER ( X0 = 32363, Y0 = 31727, Z0 = 31657 )
PARAMETER ( X1 = X0, Y1 = Y0, Z1 = Z0 )
X = MOD(MX*X,X0)
Y = MOD(MY*Y,Y0)
Z = MOD(MZ*Z,Z0)
V = X/X1 + Y/Y1 + Z/Z1
FLTRN = V - INT(V)
RETURN
END

25

Generatory quasi-losowe

• Oprócz generatorów pseudolosowych

wykorzystywane też są generatory quasi-
losowe:

–

Nie są one losowe – trzeba uważnie je stosować

–

Nie są stosowane na ogół w Monte Carlo do
transportu promieniowania – (dla nas mają
mniejsze znaczenie) ale mogą być przydatne do
celów specjalnych.

•

Bywają bardzo istotne o ogólnych zastosowaniach
Monte Carlo.

•

-polegają na generowaniu uporządkowanej
sekwencji liczb pokrywającej przedział i
teoretycznie mogącej zastąpić liczby pseuro-losowe
w MC

26

Wybieranie liczb losowych

z rozkładów

• Rozkład jednorodny
• Rozkład dyskretny
• Rozkład ciągły

– Metoda wyboru bezpośredniego
– Metoda odrzucania
– inne

• Przykład – rozkład Gaussa

27

Mapowanie ->x

wykorzystując p(x)

• Podstawowa technika wykorzystuje

odwzorowanie jednoznaczne

y->x

– y= z przedziału (0,1) oraz x z przedziału (a,b)
– Będziemy też wykorzystywać mapowanie odwrotne

28

Mapowanie

mamy:
• (a)=0
• (b)=1
• Funkcja jest niemalejąca w (a,b)

Zagadnienie sprowadza się do:
• Znalezienia funkcji (x)
• Odwrócenia funkcji (x) aby dostać x(),

jest to zamiana wylosowanych liczb

pseudolosowych w liczby podlegające

zadanemu rozkładowi f(x)

29

Mapowanie

• Muszą być spełnione warunki:

f(x)

x

F

dx

x

f

x

dx

x

f

a

x

x

f

dx

x

d

dx

x

f

d

x

x

x

x

a

x

a

ta

dystrybuan

,

)

(

'

)

'

(

)

(

'

)

'

(

)

(

)

(

)

(

)

(

)

(

}

)

d

,

(

Pr{

}

)

d

,

(

Pr{







































30

Rozkład ciągły: metoda wyboru

bezpośredniego

(ogólna procedura)

•

Obliczmy dystrybuantę:

•

Losujemy liczbę pseudolosową 

•

Przyrównujemy dystrybuantę do wylosowanej

liczby pseudolosowej

•

Wykorzystujemy formułę odwrotną by

przekształcić  w x:





x

a

dx

x

f

x

F

'

)

'

(

)

(

)

(x

F





)

(

1







F

x

31

Typowy rozkład gęstości

prawdopodobieństwa

Dystrybuanta

Funkcja odwrotna

dystrybuanty

32

Rozkład jednorodny

• wybrać x jednorodnie z (a,b):

• Tworzymy dystrybuantę.

 

a

b

dx

dx

x

f

x

f

x

f

b

a

x

x

f

b

a

b

a

















1

1

1

)

(

~

)

(

~

)

(

,

1

)

(

~

,

a

b

a

x

dx

a

b

dx

x

f

x

F

x

a

x

a

















'

1

'

)

'

(

)

(

33

Rozkład jednorodny (cd)

•

Przyrównanie:

•

Odwracanie do x():

•

Przykład: wybrać jednorodnie liczbę z

przedziału (-1,1):

a

b

a

x









)

(

a

b

a

x









1

2

))

1

(

1

(

1





















34

Rozkład dyskretny

• Dokonujemy N wyborów stanu i,

każdy opisany
prawdopodobieństwem

• Formujemy dystrybuantę:

i

p

)

(

)

(

)

(

)

(

2

2

1

1

N

N

x

p

x

p

x

p

x

f

































































N

i

N

i

x

x

x

p

x

x

x

p

p

x

x

x

p

x

x

dx

x

f

x

F

1

3

2

2

1

2

1

1

1

,

1

,

,

,

0

'

)

'

(

)

(



35

Rozkład dyskretny

• Przyrównanie:

• Odwracanie do x():

)

(x

F





 

































































1

,

,

,

,

1

1

1

1

1

2

1

1

2

1

1

1

ξ

p

x

p

ξ

p

x

p

p

p

x

p

x

F

x

N-

i

i

N

n

i

i

n-

i

i

n

if

if

if

0

if





36

Rozkład ciągły: metoda wyboru

bezpośredniego

• Formujemy dystrybuantę:

• Przyrównanie:

• Odwracanie do x():

• Przykład:





x

a

dx

x

f

x

F

'

)

'

(

)

(

)

(x

F





)

(

1







F

x











x

e

x

f

x

0

,

)

(

~

37

Rozkład ciągły: metoda

odrzucania

Ogólna procedura:

1. Znajdź wartość maksymalną rozkładu

:

2. Wybierz losowo

3. Losujemy liczbę pseudolosową 

4. Z

atrzymaj jeśli

w innym przypadku odrzuć

)

,

(

gdzie

,

kazdego

dla

)

(

~

max

b

a

x

x

x

f

f





 

a

b

x

p

a,b

x

i







1

)

(

:

max

~

)

(

f

x

f

i





i

x

38

Przeskalowany rozkład gęstości

prawdopodobieństwa

(do metody odrzucania)

Typowy rozkład

gęstości

prawdopodobieństwa

sens metody odrzucania - strzelamy

do tarczy: szare pole odrzucamy,

białe zostawiamy

39

Rozkład ciągły: metoda odrzucania z

funkcją

Podobna zasada ale z wykorzystaniem
funkcji niejednorodnej dystrybuanty p(x)
przeskalowanej stałą C

1. Znajdź taką wartość C że:
2. Wybierz losowo z rozkładem

p(x)

3. Losujemy liczbę pseudolosową 
4. Zatrzymaj jeśli:

Uwaga: sprowadza się do zwykłego odrzucania

gdy p(x) jest stałe oraz:

 

a,b

x

x

f

x

Cp





;

)

(

)

(

 

a,b

x

i



)

(

)

(

i

i

x

Cp

x

f





i

x

max

~

)

(

f

a

b

C







40

Rozkład ciągły: metoda odrzucania - zadania

A) Użyj metody odrzucania

dla rozkładu:

wykorzystując
1. Stałe ograniczenie
2. Funkcję ograniczającą:

B) Inne rozkłady

 

1

,

0

,

)

(

~





x

x

x

f

n

2

0

)

(

~







x

e

x

f

x

,

2

0

,

1

)

(

~









x

kx

x

p

















,

,

)

(

~

2

2

x

e

x

f

x

41

Budujemy dwuwymiarowy analogowy

kod Monte Carlo : 2D MC

• Przegląd wyborów jakie musimy wykonywać

w czasie symulacji MC transportu cząstki.

1.Początkowa pozycja startowa
2.Początkowy kierunek lotu
3.Początkowa energia cząstki
4.Odległość do następnego zderzenia
5.Typ zderzenia
6.Wynik rozpraszania

• Zagadnienie zliczania strumienia
• Zaczynamy od kodu homogenicznego 2D

42

Decyzja 1: Początkowa pozycja

startowa

• Zwykle jest to wybór w oparciu o

wielowymiarowe rozkłady parametrów w
przestrzeni w określonej objętości.

• Matematycznie: należy zdefiniować funkcję

rozkładu cząstek źródłowych w odpowiednim
systemie współrzędnych (wykorzystując
istniejące symetrie) a następnie dokonać
niezależnych wyborów liczb losowych,
niezależnie dla każdego wymiaru obecnego w
rozkładzie.

• Możliwe klasyczne układy współrzędnych:

kartezjański, cylindryczny, sferyczny,
a także złożone (toroidalne)

43

Układ kartezjański

Różniczkowy element objętości:

dz

dy

dx

dV 

44

Układ kartezjański (2)

• Chcemy wybrać punkt z płaskim

rozkładem prawdopodobieństwa
(każdy element objętości jest
równie prawdopodobny)













0

1

0

0

1

0

0

1

0

1

)

(

1

)

(

1

)

(

)

(

)

(

)

(

z

z

z

z

z

f

y

y

y

y

y

f

x

x

x

x

x

f

dz

z

f

dy

y

f

dx

x

f

dV









































45

Układ cylindryczny

(prawoskrętny)

• Różniczkowy element objętości:

dz

dr

rd

dV





46

Układ cylindryczny

(2)

• Płaski rozkład prawdopodobieństwa (każdy

element objętości jest równie prawdopodobny)

• Transformacja do układu kartezjańskiego





0

1

0

0

1

)

(

)

(

2

1

)

(

)

(

)

(

)

(

z

z

z

z

z

f

r

r

r

r

f

f

dz

z

f

dr

r

f

d

f

dV









































z

z

r

y

r

x











sin

cos

47

Układ sferyczny

• Różniczkowy element objętości:



 

dr

rd

d

r

dV







 sin

48

Układ sferyczny (2)

• Płaski rozkład prawdopodobieństwa

(każdy element objętości jest równie
prawdopodobny)

• Transformacja do układu kartezjańskiego























































2

1

)

(

)

(

2

1

cos

sin

)

(

)

(

)

(

)

(

3

0

2

1

f

r

r

r

r

f

f

d

f

dr

r

f

d

f

dV

















cos

sin

sin

cos

sin

r

z

r

y

r

x

49

Wybór ze źródeł złożonych

• W sytuacji złożonego źródła należy

zastosować strategię mieszania
prawdopodobieństwa. Możliwa jest wtedy
kombinacja źródeł z rozkładami o różnej
gęstości, kształcie, rozmiarze.

1.

Cząstkowe źródło wybieramy wykorzystując
całkowitą wydajność źródła
(czasteczki/sekundę)

2.

W wylosowanym źródle losujemy punkt
startowy wykorzystując jego funkcje rozkładu.

50

Źródła z niejednorodnym

rozkładem przestrzennym

• W przypadku niejednorodnym

cząstkowe rozkłady (danej
zmiennej przestrzennej) winny być
przemnożone przez odpowiednie
funkcje rozkładu.

• Przykład:

51

Decyzja 2: Początkowy

kierunek lotu

Wybór kierunku bazuje

na

prawdopodobieństwie

wejścia w różniczkowy

element kąta bryłowego

wypełniającego

pełny kąt bryłowy

(odpowiednik

powierzchni sfery o

jednostkowym

promieniu)



d











d

d

d

sin

52

Decyzja 2: Początkowy kierunek

lotu

(2)

• Uwaga: wybór osi Z jest całkowicie dowolny (w

zależności od potrzeby (istniejących symetrii)

• Na ogół definiujemy  , wtedy

element kąta bryłowego wynosi:

• Znak minus bo opada gdy rośnie





 cos











d

d

d





53

Decyzja 2: Początkowy kierunek

lotu

(3)

• Daje nam to rozkład gęstości

prawdopodobieństwa:

• W ogólności w metodzie Monte Carlo

posługujemy się cosinusami kierunkowymi:









































2

1

)

(

2

1

1

)

(

)

(

)

(

f

f

d

f

d

f

d



























z

y

x

sin

1

cos

1

2

2

54

Decyzja 3: Początkowa

energia cząstki

W ogólności początkowa energia cząstki może

bazować na rozkładzie ciągłym, dyskretnym

bądź wykorzystywać wielogrupowe spektrum

źródła.

• Ciągły: wybór bezpośredni bądź metodą

odrzucania

• Dyskretny : (zwykle kwanty ) energia zwykle

jest sprzężona z wydatkiem

• Wielogrupowe: Multigroup: źródło w każdej

grupie jest wycałkowanym rozkładem p-stwa

na przestrzeni grupy. Zatem wartości

indywidualnych grup mają charakter dyskretny

i tak są traktowane.

i

p

55

Decyzja 4: Oczekiwana odległość do

następnego zderzenia

• Dla ośrodka nieskończonego: p-stwa

zderzenia na dx w odległości x – p-stwo

dotarcia do x bez kolizji razy p-stwo

kolizji na dx

• Zatem funkcja rozkładu gęstości p-stwa:

(jest już znormalizowana w ) 

 

 





dx

e

dx

x

f

t

x

t











 

x

t

t

e

x

f















,

0

56

Oczekiwana odległość do następnego

zderzenia (2)

• dystrybuanta :

• Funkcja odwrotna:

• W ujęciu długości optycznej:
mamy:  

 

x

t

e

x

F







1





t

x











1

ln

x

t



















ln

lub

1

ln











57

Decyzja 5: Typ zderzenia

• Gdy zderzenie następuje w pewnym punkcie musimy

wybrać typ reakcji na podstawie wartości

makroskopowych p. czynnych.

• Stąd mamy prawdopodobieństwa:

• Wybór typu reakcji wykorzystując w/w p-stwa do

budowy rozkładu dyskretnego i odpowiedniej

procedury w tym przypadku  

 

 

 

 

E

E

E

E

c

f

s

t















 

 

 

 

 

 

E

E

p

E

E

p

E

E

p

t

c

c

t

f

f

t

s

s



















58

Decyzja 6: Wynik

rozpraszania

• Wynik rozpraszania cząstki z energią

początkową E dany jest rozkładem:

gdzie M = materiał, „primy” są związane ze stanem po zderzeniu

• losujemy wykorzystując:

• losujemy wykorzystując:























ˆ

ˆ

,

0

0

E

E

M
s

0



 





E

d

E

E

f

M
s



















0

0

0

,

E



 





i

M
s

E

E

E

f

0

,











59

Wynik rozpraszania (2)

• W pewnych przypadkach rozpraszania występuje

jednoznaczna zależność pomiędzy kątem
rozproszenia oraz utratą energii. Wtedy jedna ze
zmiennych jest zależna od drugiej – wystarczy
losować tylko jedną.

• W przypadku wielogrupowym zależności kątowe

w rozpraszaniu pomiędzy grupami są
reprezentowane przez rozwinięcie Legendre’a w
funkcji kąta rozpraszania, lub za pomocą
przedziałów „równo-prawdopodobnych”

60

Zliczanie strumienia:

estymator zderzeń (collision estimator)

•

Występują dwa podstawowe sposoby zliczania estymatora

strumienia (w skończonym obszarze): estymator zderzeń

(collision estimator) oraz estymator długości trasy

(tracklength estimator).

•

Estymator zderzeń bazuje na obliczaniu przyrostu

strumienia na podstawie występujących (symulowanych)

zderzeń.

•

Z gęstości reakcji:

•

mamy po każdym zderzeniu
(infinitezymalny) przyrost strumienia:

•

Można też budować estymatory bazujące na

poszczególnych reakcjach (konsekwencją – zwiększenie

błędu – ale wprowadzamy korelacje)

 

   

r

r

r

RR

g

tg

g













 

 

r

r

tg

g









1



 

.

,

,

,

,

etc

el

s

a

x

r

xg







61

Zliczanie strumienia:

e

stymator długości trasy (tracklength

estimator).

•

definicja

•

Zalety:

–

Możliwa estymacja w próżni

–

Umożliwia określenie co powinno się zdarzyć (np.

ze zbyt małym p-stwem by pojawiło się przy małym

rozmiarze próby)

Oba są powszechnie używane o obliczeniach

reaktorowych, k-eff, - zagadnienia

krytyczności

c

cx

cx

RR







c

V

c

komórce

w

przebyta

droga

całałkowi





c

62

Techniki Redukcji Wariancji

•

Cel: zmniejszyć wariancję estymowanych
parametrów bez zwiększania rozmiaru próby

•

Reguły postępowania (oszukiwanie rozkładów)

•

Najczęściej używane techniki:

–

Unikanie absorpcji (ze zmiana wag w miejscu
absorpcji)

–

Rozszczepianie cząstek (nie mylić z
rozszczepianiem atomów)

–

Wymuszanie zderzeń

–

Rosyjska ruletka

–

Transformacja eksponencjalna

–

Modyfikacja (biasing) źródła (pozycji, energii lub
kierunku)

63

Cel stosowania

•

Aby uzyskać lepszy program Monte Carlo –

lepszy to znaczy generujący wyniki z mniejszą

wariancją

•

Mniejszą wariancja w ujęciu czarnej skrzynki:

oznacza że x-sy są do siebie najbardziej

zbliżone (jak to możliwe)

•

Skuteczny dla lokalnych wielkości – gorzej w

przypadku globalnych parametrów

64

Reguły postępowania (oszukiwanie

rozkładów)

•

Mapowanie prowadzimy wykorzystując (zamiast

naturalnego) zmodyfikowany rozkład

prawdopodobieństwa – do postaci przez nas pożądanej

tak by zagęścić wybór w obszarach naszego

zainteresowania.

•

Zgodność z naturą zachowujemy zmieniając wagę

cząstek

•

Waga cząstki odpowiada części cząstki jaka w

naturalnym rozkładzie pojawiłaby się w danej sytuacji.

•

Jeśli naturalny rozkład daje gęstość p-stwa: a

my chcemy użyć rozkładu o gęstości: to

musimy zastosować korektę wagi:

•

Konieczny wymóg: wszystkie możliwe punkty x

dostępne w naturalnym rozkładzie muszą być

dostępne w nowym – zmodyfikowanym rozkładzie.

 

 

x

p

x

p

w

corr

~



 

x

p

 

x

p

~

65

Przykład: wykorzystanie rozkładu

płaskiego

• Z lenistwa wykorzystujemy do mapowania

najprostszy rozkład płaski

• Nie prowadzi do zniekształcenia wyników

jeśli skorygujemy wagi następująco:

• Chociaż wynik nie będzie zniekształcony,

postępowanie takie na ogół nie jest

rozsądne. (potrzebna jest odpowiednia

motywacja do właściwej zmiany rozkładu)

 

 

 





 



a

b

x

p

a

b

x

p

x

p

x

p

w

corr











1

~

66

Unikanie absorpcji

• Najczęściej używany technika redukcji

wariancji. W wielu kodach transportowych jest
ona wymuszona (nie może być wyłączona
przez użytkownika)

• Pociąga zmianę wag w miejscu absorpcji

• Opis ogólny: Podstawowa zasada polega na

uniemożliwieniu absorpcji cząsteczki; wtedy
cząsteczka będzie żyć dłużej i mieć więcej
szans na bycie zliczaną w innych procesach.

• Losowy wybór wymagający dopasowania:

- typ zderzenia

67

Unikanie absorpcji (2)

• Zarys matematyczny:

–

Załóżmy że są możliwe dwa stany wyjściowe

zderzenia: rozpraszanie i absorpcja. W stanie

naturalnym mamy dystrybucje (tego co powinno

być):

–

rozpraszania z odpowiednim rozkładem p.stwa:

–

absorpcji z odpowiednim rozkładem p.stwa:

•

Nasze decyzje (mapowanie) podejmujemy jednak

tak by zawsze wypadło rozpraszanie (ze 100%

p.stwem)

•

Wymagana korekcja wagi cząstki rozproszonej:

t

s

p







1

 

 

t

s

t

s

corr

x

p

x

p

w















00

.

1

~

t

a

p







2

68

Unikanie absorpcji (3)

• Uwaga: zmieniliśmy naturalnie możliwy wybór

(absorpcja) NIEMOŻLIWYM. Jest to możliwe

tylko dlatego że cząstka w wypadku absorpcji i

tak przestaje żyć (nie wnosi żadnych nowych

wkładów do odpowiedzi systemu)

• Procedura ta oprócz zmniejszenia wariancji (na

historię) gwarantuje zwiększenie czasu obliczeń

na skutek wydłużenia każdej historii.

• Czasami to się nie musi opłacać
• Inne spojrzenie na tę technikę: Dzielimy

cząstkę na dwie części: ułamek który się

rozprasza i ułamek który jest absorbowany.

Kontynuujemy śledzenie tylko pierwszej

podczas gdy druga umiera.

69

Rozszczepianie cząstek (nie

atomów)

• Stanowi podstawę innych technik redukcji

wariancji

• Nie wypełnia wzorca modyfikacji

dystrybucji p.stwa. Polega na unikaniu
dyskretnych decyzji w zamian za śledzenie
wszystkich możliwych opcji decyzyjnych

• Ponieważ probabilistyczne decyzje są

unikane, wpływ tych decyzji na wariancje
też znika.

70

Rozszczepianie cząstek

(2)

• Dwie typowe sytuacje w transporcie cząstek z

wykorzystaniem rozszczepienia (historii)

1. W przypadku zderzenia, historia jednej cząstki

przebiega dalej niezmienionym śladem (tylko z

mniejszą wagą) aż do zakończenia życia a wtedy

wracamy do punktu rozszczepienia by prowadzić

historię cząstki rozproszonej aż do kończ jej życia.

2. Cząstka jest rozszczepiana na dwie lub więcej w

sytuacji wejścia w obszar o większej ważności (tam

gdzie chcemy zmniejszyć wariancje - na ogół by ją

wyrównać do poziomu średniego)

• Matematycznie idea ta polega na śledzeniu

wszystkich możliwych opcji dyskretnych rozkładów

p.stwa z odpowiednim ważeniem ich historii

proporcjonalnie do wartości dyskretnych p.stw.

71

Rozszczepianie cząstek

(3)

• Przykład: zakładamy historię cząstki której

obecna waga wynosi 0.5 napotyka na dyskretną

decyzję z trzema możliwymi stanami:

– stan 1 z 60% prawdopodobieństwem;
– stan 2 z 30% prawdopodobieństwem;
– stan 3 z 10% prawdopodobieństwem.

• Unikamy wyboru, a w zamian śledzimy każdą

opcję.

– Kontynuujemy historię w stanie 1 z wagą 0.3

zapisując w banku (by kontynuować później) historię

w stanie 2 z wagą 0.15 oraz historię 3 z wagą 0.05.

• Po skończeniu każdej historii wracamy do

banku by podjąc pozostające tam cząstki aż

opróżnimy go całkowicie.

72

Wymuszanie zderzeń

• Ogólny opis: Wymuszamy zderzenie w

szczególnym obszarze drogi cząstki.
Wykonamy to w szczególny sposób aby
zapobiec ucieczce cząstki (bez zderzenia)
z tego obszaru.

• Utrzymuje cząstki dłużej przy życiu, zatem

zwiększa ich kontrybucje statystyczne do
wszystkich zliczeń.

•

Losowy wybór wymagający dopasowania:
- dystans do najbliższego zderzenia

73

Wymuszanie zderzeń (2)

• Zarys matematyczny: zasadniczo jest to

technika rozszczepienia cząstki na część

ulegającą

zderzeniu w określonym obszarze

oraz na część oraz

nie ulegającą

zderzeniu

w nim.

• Załóżmy że odległość do granicy wynosi 0

(średnich dróg swobodnych)

• Naturalny rozkład p.stwa ma postać:

–

P.stwo ucieczki:

–

P.stwo uniknięcia ucieczki:

• Nasza decyzja to uniknięcie ucieczki ze 100%

p.stwem

0





e

p

escape

0

1









 e

p

escape

non

74

Wymuszanie zderzeń

(2)

• Korekta wagi wynosi:

• Rozszczepieniowa natura tej techniki powoduje

że musimy wziąć pod uwagę że część cząstki

ucieka poza dany obszar. Ta część powinna być

uwzględniona przy zliczaniu ucieczki poza

najbliższa granicę, w ilości:

, gdzie

w jest wagą cząstki przed przejściem.

• Także tutaj efektywność może się nie poprawić

(tzn. straty z wydłużenia czasu obliczeń

przewyższą zysk z mniejszej wariancji)

–

Faktycznie ta technika nie jest zbyt często używana.

Aczkolwiek czasami występuje konieczność jej użycia

(np. symulacja odpowiedzi małego detektora)

 

 

0

0

1

1

1

~



















e

e

x

p

x

p

w

i

i

corr

0





we

75

Rosyjska ruletka

• Podobnie do rozszczepienia

matematyczne narzędzie wykorzystywane
w technikach zmniejszania wariancji

• Idea: kombinacja kilku cząstek w jedną

przy pomocy selektywnej eliminacji (jak
w rosyjskiej ruletce)

• Matematycznie oznacza to zabijanie

cząstek z dużym p.stwem 1-p (90-99%)
gdzie p oznacza szansę przeżycia.

• Jeśli cząstka przeżyje, jej waga jest

zwiększona o czynnik 1/p

76

Rosyjska ruletka (2)

• Narzędzie to jest wykorzystywane w przypadku

stosowania innych technik redukcji wariancji do

zakańczania historii cząstek których waga jest na

tyle mała że tracą znaczenie.

• Narzędzie to jest konieczne w przypadku

stosowania jednocześnie techniki unikania absorpcji

oraz wymuszania zderzeń (na całym obszarze) gdyż

jest to wtedy jedyna sposób na zakończenie historii

cząstek bo użyte techniki wyeliminowały naturalne

drogi zakończenia historii.

• W praktyce jest wykorzystywane zawsze gdy waga

cząstki spadnie poniżej dolnego progu odcięcia

wagi..

–

Technika ta nie redukuje wariancji wprost, tylko pozwala na

oszczędzenie czasu obliczeniowego poprzez przerywanie

historii o nikłym znaczeniu – słowem pozwala zwiększyć

efektywność.

77

Transformacja

eksponencjalna

• Ogólny opis: idea ta polega na ułatwieniu

cząstce podróży w pożądanym przez nas

kierunku, poprzez rozrzedzenie materiału

w tym kierunku oraz zagęszczenie w

kierunkach niepożądanych.

– Nie powoduje że pożądany kierunek będzie

częściej wybierany a jedynie ułatwi penetrację

materiału w tym kierunku.

• Losowy wybór wymagający dopasowania:

- dystans do najbliższego zderzenia

78

Transformacja eksponencjalna

(2)

• Zarys matematyczny: zmieniamy

całkowity przekrój czynny oraz

uzależniamy go od kierunku lotu:

gdzie p jest ogólnym parametrem

definiowanym przez użytkownika (0<p<1)

oraz 

d

jest wyróżnionym kierunkiem

• Ograniczenia zmienności p do przedziału

(0<p<1) powoduje że zmieniony przekrój

czynny będzie większy od zera i mniejszy

od dwukrotności prawdziwego przekroju

czynnego.

  



d

t

t

p 















ˆ

ˆ

1

ˆ

*

79

Transformacja eksponencjalna

(3)

• Rozważmy przypadek cząstki

lecącej w kierunku zewnętrznej

granicy. Mamy dwa przypadki:

– Cząstka osiągnie granicę
– Cząstka się rozproszy

• Naturalny rozkład p-stw ma postać:

– P.stwo ucieczki:

– P.stwo zderzenia w odległości s:

0





e

p

escape

s

t

t

e

s

p









)

(

80

Transformacja eksponencjalna

(4)

• Wykorzystywany rozkład ma postać:

– P.stwo ucieczki:
– P.stwo zderzenia w odległości s:

• Korekta wagi:

– W przypadku ucieczki:

– W przypadku zderzenia w odległości s:





d

p

escape

e

p













ˆ

ˆ

1

0



 









s

p

d

t

d

t

e

p

s

p

























ˆ

ˆ

1

ˆ

ˆ

1





d

d

p

p

corr

e

e

e

w

























ˆ

ˆ

ˆ

ˆ

1

0

0

0



















d

p

p

d

t

t

corr

p

e

e

p

e

w

d

d













































ˆ

ˆ

1

ˆ

ˆ

1

ˆ

ˆ

ˆ

ˆ

1

0

0

0







81

Modyfikacja źródła

• Jest to najbardziej ogólna technika

redukcji wariancji. Zasada jest
prosta aczkolwiek nie dająca
jasnych kryteriów postępowania:
zamiast wykorzystywać naturalny
rozkład cząstek źródłowych
wykorzystujemy jakiś inny rozkład o
którym myślimy że będzie lepszy.

• Losowy wybór wymagający dopasowania:

- pozycja cząstki, energia i kierunek lotu

82

Modyfikacja źródła (2)

• Zarys matematyczny:

–

Zasadniczo nie ma ścisłych wskazówek jak
postępować – musimy wykorzystywać ogólne
reguły oszukiwania rozkładów.

–

Jeśli fizyka dyktuje rozkład p.stwa:

a my chcemy użyć:

możemy to zrobić stosując korektę wagi:

 

 

x

p

x

p

w

corr

~



 

x

p

 

x

p

~

83

Wybór zmodyfikowanego rozkładu

• W ogólności chcemy zmodyfikować

naturalny rozkład tak by faworyzować
wybory bardziej ważne.

• Pytanie „Co jest bardziej ważne?”
• Odpowiedź:

WAŻNOŚĆ= SPODZIEWANA KONTRYBUCJA
jako kontrybucję rozumiemy zliczenia estymatora
wybranej wielkości fizycznej – np. sygnału
detektora

• Dlatego nasze zadanie polega na takiej

modyfikacji dystrybucji aby kontrybucje do
naszych zliczeń były największe.

84

Wybór zmodyfikowanego rozkładu

(2)

• Jeśli znamy funkcję ważności

(importance) możliwych wyborów I(x),
wtedy optymalną teoretycznie
dystrybucją będzie:

• Możliwe podejścia można podzielić na

kategorie:

–

(1) heurystyczną (bazującą na intuicji)

–

(2) eksperymentalną

–

(3) wykorzystującą sprzężony strumień

 

   

x

I

x

p

x

p

~

~

85

Podejście heurystyczne

•

Modyfikujemy naturalny rozkład tak by

faworyzować wybory cząstek o których myślimy

że są bardziej ważne dla naszego zagadnienia.

–

Przykład 1: położenie cząstki – jeśli interesuje nas

ucieczka przez prawą stronę – z preferencją

wybieramy położenia bliższe prawej strony

–

Przykład 2: kierunek lotu - jeśli interesuje nas

ucieczka przez lewą stronę – z preferencją wybieramy

kierunki w lewo

–

Przykład 3: energia cząstki- jeśli interesuje nas

głęboka penetracja – preferujemy wybory energii dla

których całkowity przekrój czynny jest mały

•

Bazujemy tutaj na naszej intuicji – ocenie jak

bardzo należy faworyzować ważniejsze cząstki.

Jest to metoda prób i błędów – przydatne bywa

doświadczenie.

86

Podejście eksperymentalne

•

Podejście to bazuje na wykonaniu kilku krótkich
przebiegów testowych w celu uzyskania wartości
względnych funkcji ważności dla różnych
możliwych wyborów parametrów źródłowych.

•

Każdy test powinien posiadać ograniczony
obszar (domenę) zmienności poszczególnych
parametrów źródłowych, a po wykonaniu
obliczeń wykorzystujemy uzyskaną funkcji
ważności I(x) do utworzenia skorygowanego
rozkładu:

•

Możemy też zaangażować niepewność obliczonej
funkcji:

 

   

x

I

x

p

x

p

~

~

 

   

 

 

 

x

x

I

x

p

x

I

x

p

x

p

2

2

~

~

~





87

Optymalna modyfikacja -

wyprowadzenie

wzorów

•

Funkcja ważności I(P) jest wykorzystywana do modyfikacji

źródła w następujący sposób:

gdzie:
optymalna (nieznormalizowana dystrybucja

do

wykorzystania w symulacji)
prawdziwa (naturalna) dystrybucja (w

poprzednich

wzorach – p(x) )

Oczekiwana kontrybucja przy wyborze naturalnym - P

 

)

(

)

(

~

~

P

I

P

S

P

p







 

)

(

)

(

~

P

I

P

S

P

p







88

Optymalna modyfikacja 2

Wyprowadzenie:

•

Użycie zamiast daje nam oczekiwaną
wariancję w zmodyfikowanej symulacji optymalnej w
postaci:

=>

gdzie jest prawdziwą kontrybucją
niezmodyfikowaną,

•

Nałożenie warunku optimum (najmniejszej wariancji)

prowadzi do następującego wyniku:

 

 



 

 

 

P

d

I

P

p

P

I

P

S

P

p

I

P

d

I

P

I

w

P

p

I

P

i

P

corr

i

















2

2

2

2

~

)

(

)

(

~

)

(

~





























 

P

p



~

)

(P

I



)

(P

S



 

0

~

~

2









p

p

89

Optymalna modyfikacja - 3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I

P

I

P

S

P

p

P

d

I

P

p

P

I

P

S

I

P

p

P

I

P

S

P

d

P

p

P

I

P

S

I

P

p

P

I

P

S

I

P

p

P

I

P

S

P

d

P

p

P

I

P

S

I

P

p

P

I

P

S

P

p

I

P

p

P

I

P

S

P

d

I

P

p

P

I

P

S

P

p

p

P

d

I

P

p

P

I

P

S

P

p

p

p

p

P

P

P

P

P

)

(

)

(

~

0

~

)

(

)

(

~

)

(

)

(

~

)

(

)

(

2

~

)

(

)

(

~

)

(

)

(

~

)

(

)

(

~

)

(

)

(

~

2

~

)

(

)

(

~

)

(

)

(

~

~

~

)

(

)

(

~

~

~

~

2

2

2

2

2

2











































































































































































































































































if

90

Optymalna modyfikacja (4)

•

Problem pojawić się może gdy nasz wybór

nie gwarantuje kontrybucji do naszego zliczania

•

Natomiast gdy wybór wywołuje kontrybucje

z niepewnością , wtedy optymalna

modyfikacja ma postać:

•

W przypadku poisson’owskim (głęboka

penetracja analogowa) mamy:

Wtedy optimum jest

(oznacza to zanik optymalnej modyfikacji)

)

(

)

(

)

(

P

I

P

I

P













P



)

(P

I



P



)

(P

I



)

(P





 

   

2

2

)

(

)

(

)

(

~

~

P

P

I

P

S

P

p













 

)

(

)

(

~

~

P

I

P

S

P

p







91

Ważenie po celach

(MCNP)

•

Rules:

1.

Gdy cząstka przechodzi z obszaru A do B o

większej ważności wtedy rozszczepiamy w

stosunku:

2.

Gdy cząstka przechodzi z obszaru A do B o

mniejszej ważności wtedy gramy rosyjską
ruletkę z szansą przeżycia:

old

new

I

I

old

new

I

I

92

Okna wagowe w MCNP

Reguły:

Przy przejściu cząstki z obszaru A do nowego B

gdzie następuje zmiana wagi dostosowujemy
wagę tak by mieściła się w wymaganym
przedziale (oknie wagowym)

Stosujemy RR lub rozszczepienie wykorzystując:

desired

current

w

w

high

desired

low

w

w

w



