„A statisztikus fizikai szimulációk alapjai és a Monte Carlo módszer” változatai közötti eltérés

Innen: TételWiki
a
a (A Metropolis algoritmus)
 
(9 közbenső módosítás, amit 2 másik szerkesztő végzett, nincs mutatva)
1. sor: 1. sor:
 
== Statisztikus fizikai szimulációk alapjai ==
 
== Statisztikus fizikai szimulációk alapjai ==
 +
Legyen egy sokaságunk, aminek az energiája <math>E=E(\vec{r})</math>, ahol <math>\vec{r}</math> az összes szabadsái fokot tartalmazza (mechanikai rendszerre <math>\vec{r} = (\vec{q}, \vec{p})</math>). Ha a rendszer követi a Boltzmann-statisztikát, akkor a rendszer <math>[\vec{r}, \vec{r}+d \vec{r}]</math> intervallumban való megtalálási valószínűsége:
 +
 +
<math>P(\vec{r})d\vec{r} = \frac{e^{-\beta E(\vec{r})}}{Z}d\vec{r}</math>
 +
 +
Egy ''A'' mennyiség átlagát a sokaságon a következőképp számoljuk (PS = phase space = fázistér):
 +
 +
<math>\langle A \rangle = \int_{PS} A(\vec{r}) \frac{e^{-\beta E(\vec{r})}}{Z}d\vec{r}</math>
 +
 +
Az integrált csak akkor tudjuk egzaktul kiszámolni, ha a fázistér minden pontjára kiértékeljük a függvényt, és úgy átlagolunk. A Monte-Carlo integrálással ezt a problémát oldhatjuk meg. Az <math>\int d\vec{r}</math> helyére <math>\frac{1}{N} \sum</math> írunk, így egy véges összeget kell elvégeznünk. A pontok számát (N) az határozza meg, hogy mekkora pontossággal akarjuk az integrált meghatározni. Az integrál dimenziójától függetlenül a hiba <math>1/\sqrt{N}</math> szerint csökken, vagyis 10x nagyobb pontossághoz 100x annyi pontot kell kiszámolni.
 +
 
== Molekuladinamika ==
 
== Molekuladinamika ==
 +
A molekuladinamika a részecskék mikroszkopikus dinamikájának követésével foglalkozik. Valódi rendszerben 10<sup>23</sup> nagyságrendű részecske van, ezt a mai számítógépekkel szimulálni lehetetlen. Azonban ennél kevesebb részecskét is elég szimulálnunk ahhoz, hogy a termodinamikai tulajdonságokat vizsgálhassuk. A szimulációkban a párkölcsönhatásokat vesszük figyelembe, ám lehetséges közelítéseket tenni. A részecskék közt a Van der Waals erő hat, ami elég gyorsan lecseng, így távoli részecskék közt elhanyagolható (ezen a pontot különbözik a molekuladinamika és a gravitációs N-test szimuláció, ahol az 1/r<sup>2</sup>-es erő miatt nem hanyagolható el a kölcsönhatás).
 +
 +
A szimulációban van tehát N darab részecske, melyekre a Newton-törvények alapján kiszámítjuk a rájuk ható erőt, majd léptetjük a rendszert. Mivel a kezdőállapotokat általában nem a termodinamikai egyensúlyból indítjuk, meg kell várni, hogy a rendszer beálljon abba. Hogy mikor állt be, azt az ekvipartíció segítségével mutathatjuk ki:
 +
 +
<math>\left \langle \frac{1}{2}mv^2 \right \rangle = \frac{3}{2}kT</math>
 +
 +
Ha a rendszer elérte az egyensúlyi állapotát, mérhetővé válnak a termodinamikai változók (hőmérséklet, nyomás, hőkapacitás, stb.)
 +
 +
===Fizikai mennyiségek mérése a szimulációkban===
 +
====Pillanatnyi hőmérséklet====
 +
Bár azt mondtuk, hogy a rendszert az egyensúlyi állapotában vizsgáljuk, a nem egyensúlyi állapotban is bevezethetünk egy ''pillanatnyi hőmérsékletnek'' nevezett mennyiséget. Termikus egyensúlyban igaz az ekvipartíció:
 +
 +
<math>3(N-1) \cdot \frac{1}{2}k_BT = \left \langle \frac{m}{2} \sum_{i=1}^N v_i^2 \right \rangle</math>,
 +
 +
ahol 3(N-1) a szabadsági fokok száma (a TKP 3 koordinátája van levonva). A nem egyensúlyi helyzetben, bár nem igaz az ekvipartíció, ezt a képletet alkalmazhatjuk a hőmérséklet meghatározására.
 +
 +
====Teljes energia====
 +
<math>E = \frac{m}{2} \sum_{i=1}^N v_i^2 + \sum_{i \neq j} U(|r_i - r_j|)</math>
 +
 +
====Hőkapacitás====
 +
Fluktuáció-disszipáció tétel alapján:
 +
 +
<math>C_V = \frac{1}{k_BT^2}[\langle E^2 \rangle - \langle E \rangle^2]</math>
 +
 +
====Nyomás====
 +
A nyomást a viriál-tétel segítségével fejezhetjük ki:
 +
 +
<math>pV = Nk_BT + \frac{1}{3} \left \langle \sum_{i<j} r_{ij} \cdot F_{ij} \right \rangle</math>
 +
 +
===Verlet-algoritmus===
 +
Molekuladinamikai szimulációkban legtöbbször a Verlet-algoritmust használják a diffegyenletek megoldására. A Runge-Kuttával szemben több előnye van:
 +
* gyorsabb, mert egy lépésben csak egyszer kell a gyorsulásokat számolni
 +
* majdnem olyan pontos, mint a RK4 (<math>O(\tau^4) \text{vs.} O(\tau^5)</math>)
 +
* jól megőrzi az energiát
 +
* időtükrözésre nem változik (ez a részleges egyensúly feltétele miatt fontos)
 +
 +
A Verlet-algoritmus egy lépése ('''R(t)''' a koordináták, '''V(t)''' a sebességek, '''A(t)''' a gyorsulások):
 +
 +
<math>R_{n+1} = 2R_n - R_{n-1} + \tau^2 A_n + O(\tau^4)\,</math>
 +
 +
<math>V_n = \frac{R_{n+1} - R_{n-1}}{2\tau} + O(\tau^2)</math>
 +
 +
Hátrányai:
 +
* két előző lépést használ, így nem indítható 1 kezdeti feltételből
 +
* a sebesség és a pozíció nem egyszerre frissítődik, a sebesség "le van maradva"
 +
 +
Megoldás: velocity-Verlet algoritmus:
 +
 +
<math>R_{n+1} = R_n + \tau V_n + \frac{\tau^2}{2} A_n + O(\tau^3)</math>
 +
 +
<math>V_{n+1} = V_n + \frac{\tau}{2} (A_{n+1} + A_n) + O(\tau^3)</math>
 +
 +
* R-ben már csak <math>O(\tau^3)</math> pontosságú
 +
* ez már nem 2 előző lépést használ
 +
* koordináták és sebességek egyszerre frissülnek
 +
* sebesség is <math>O(\tau^3)</math> pontosságú
 +
 
== A Metropolis algoritmus ==
 
== A Metropolis algoritmus ==
 +
A Metropolis algoritmus egy eljárás, amivel bonyolult valószínűségi eloszlások mintavételezhetőek. Másképp megfogalmazva szeretnénk előállítani egy <math>P</math> eloszlást, a kérdés az, hogy milyen átmeneti valószínűségekkel <math>(W)</math> vezérelt folyamat (egészen pontosan Markov-folyamat) vezet ilyenre? A módszer kihasználja a részletes egyensúly elvét, ez esetünkben azt jelenti, hogy:
 +
 +
:<math>P_i W_{ij} = P_j W_{ji}\,</math> teljesül minden <math> i, j</math>-re.
 +
 +
Ez még kevés feltétel az átmeneti valószínűségekre, ezért a Metropolisz megkötés a következő:
 +
 +
:<math>W_{ij} = \mathrm{min}\left[1, \frac{P_j}{P_i} \right]</math>
 +
 +
Bár általánosan tetszőleges eloszlás előállítható, fizikában általában valamilyen egyensúlyi eloszlást szeretnénk előállítani. Például, ha termális egyensúlyt szeretnénk, akkor <math>P \propto \exp\left( -\beta E \right)</math>, és a Metropolisz megkötésben szereplő hányados is exponenciális alakban írható. A <math>P</math>-t előállító eljárás a következő:
 +
 +
# Induljunk ki egy A konfigurációból.
 +
# Állítsuk elő a megváltoztatott B konfigurációt.
 +
# Számoljuk ki erre a két konfigurációra a megkötésben szereplő hányadost.
 +
# Generáljunk egy 0-1 intervallumon vett egyenletes eloszlású véletlen számot: X.
 +
# Ha <math>W_{AB} > X </math> akkor a B konfigurációt fogadjuk el, ellenkező esetben marad az A.
 +
# Ismételjük 1-től.
 +
 
== A Monte-Carlo módszer ==
 
== A Monte-Carlo módszer ==
 +
A Monte-Carlo módszernek nevezzük az olyan eljárásokat, amelyek a problémákat random számok és valószínűségek felhasználásával oldják meg. Az eljárás során ismétlődően kiértékelünk egy determinisztikus modellt, random számokat használva inputnak. Akkor használják, ha a feladat nagyon összetett, nemlineáris, bonyolultak a határfeltételei, vagy nagyon sok paramétertől függ, illetve sok dimenziós.
  
 +
Használata:
 +
# Állítsuk föl a modellt: y = f(x<sub>1</sub>, x<sub>2</sub>, ..., x<sub>q</sub>)
 +
# Generáljunk random számokat inputnak: x<sub>i1</sub>, x<sub>i2</sub>, ..., x<sub>iq</sub>
 +
# Értékeljük ki a modellt, az eredményt tároljuk el y<sub>i</sub>-ben
 +
# Ismételjük a 2. és 3. lépéseket n-szer
 +
# Elemezzük az eredményeket hisztogram, összesítő statisztikák, stb. segítségével
  
 
{{MSc záróvizsga}}
 
{{MSc záróvizsga}}

A lap jelenlegi, 2011. június 15., 12:50-kori változata

Statisztikus fizikai szimulációk alapjai

Legyen egy sokaságunk, aminek az energiája E=E(\vec{r}), ahol \vec{r} az összes szabadsái fokot tartalmazza (mechanikai rendszerre \vec{r} = (\vec{q}, \vec{p})). Ha a rendszer követi a Boltzmann-statisztikát, akkor a rendszer [\vec{r}, \vec{r}+d \vec{r}] intervallumban való megtalálási valószínűsége:

P(\vec{r})d\vec{r} = \frac{e^{-\beta E(\vec{r})}}{Z}d\vec{r}

Egy A mennyiség átlagát a sokaságon a következőképp számoljuk (PS = phase space = fázistér):

\langle A \rangle = \int_{PS} A(\vec{r}) \frac{e^{-\beta E(\vec{r})}}{Z}d\vec{r}

Az integrált csak akkor tudjuk egzaktul kiszámolni, ha a fázistér minden pontjára kiértékeljük a függvényt, és úgy átlagolunk. A Monte-Carlo integrálással ezt a problémát oldhatjuk meg. Az \int d\vec{r} helyére \frac{1}{N} \sum írunk, így egy véges összeget kell elvégeznünk. A pontok számát (N) az határozza meg, hogy mekkora pontossággal akarjuk az integrált meghatározni. Az integrál dimenziójától függetlenül a hiba 1/\sqrt{N} szerint csökken, vagyis 10x nagyobb pontossághoz 100x annyi pontot kell kiszámolni.

Molekuladinamika

A molekuladinamika a részecskék mikroszkopikus dinamikájának követésével foglalkozik. Valódi rendszerben 1023 nagyságrendű részecske van, ezt a mai számítógépekkel szimulálni lehetetlen. Azonban ennél kevesebb részecskét is elég szimulálnunk ahhoz, hogy a termodinamikai tulajdonságokat vizsgálhassuk. A szimulációkban a párkölcsönhatásokat vesszük figyelembe, ám lehetséges közelítéseket tenni. A részecskék közt a Van der Waals erő hat, ami elég gyorsan lecseng, így távoli részecskék közt elhanyagolható (ezen a pontot különbözik a molekuladinamika és a gravitációs N-test szimuláció, ahol az 1/r2-es erő miatt nem hanyagolható el a kölcsönhatás).

A szimulációban van tehát N darab részecske, melyekre a Newton-törvények alapján kiszámítjuk a rájuk ható erőt, majd léptetjük a rendszert. Mivel a kezdőállapotokat általában nem a termodinamikai egyensúlyból indítjuk, meg kell várni, hogy a rendszer beálljon abba. Hogy mikor állt be, azt az ekvipartíció segítségével mutathatjuk ki:

\left \langle \frac{1}{2}mv^2 \right \rangle = \frac{3}{2}kT

Ha a rendszer elérte az egyensúlyi állapotát, mérhetővé válnak a termodinamikai változók (hőmérséklet, nyomás, hőkapacitás, stb.)

Fizikai mennyiségek mérése a szimulációkban

Pillanatnyi hőmérséklet

Bár azt mondtuk, hogy a rendszert az egyensúlyi állapotában vizsgáljuk, a nem egyensúlyi állapotban is bevezethetünk egy pillanatnyi hőmérsékletnek nevezett mennyiséget. Termikus egyensúlyban igaz az ekvipartíció:

3(N-1) \cdot \frac{1}{2}k_BT = \left \langle \frac{m}{2} \sum_{i=1}^N v_i^2 \right \rangle,

ahol 3(N-1) a szabadsági fokok száma (a TKP 3 koordinátája van levonva). A nem egyensúlyi helyzetben, bár nem igaz az ekvipartíció, ezt a képletet alkalmazhatjuk a hőmérséklet meghatározására.

Teljes energia

E = \frac{m}{2} \sum_{i=1}^N v_i^2 + \sum_{i \neq j} U(|r_i - r_j|)

Hőkapacitás

Fluktuáció-disszipáció tétel alapján:

C_V = \frac{1}{k_BT^2}[\langle E^2 \rangle - \langle E \rangle^2]

Nyomás

A nyomást a viriál-tétel segítségével fejezhetjük ki:

pV = Nk_BT + \frac{1}{3} \left \langle \sum_{i<j} r_{ij} \cdot F_{ij} \right \rangle

Verlet-algoritmus

Molekuladinamikai szimulációkban legtöbbször a Verlet-algoritmust használják a diffegyenletek megoldására. A Runge-Kuttával szemben több előnye van:

  • gyorsabb, mert egy lépésben csak egyszer kell a gyorsulásokat számolni
  • majdnem olyan pontos, mint a RK4 (O(\tau^4) \text{vs.} O(\tau^5))
  • jól megőrzi az energiát
  • időtükrözésre nem változik (ez a részleges egyensúly feltétele miatt fontos)

A Verlet-algoritmus egy lépése (R(t) a koordináták, V(t) a sebességek, A(t) a gyorsulások):

R_{n+1} = 2R_n - R_{n-1} + \tau^2 A_n + O(\tau^4)\,

V_n = \frac{R_{n+1} - R_{n-1}}{2\tau} + O(\tau^2)

Hátrányai:

  • két előző lépést használ, így nem indítható 1 kezdeti feltételből
  • a sebesség és a pozíció nem egyszerre frissítődik, a sebesség "le van maradva"

Megoldás: velocity-Verlet algoritmus:

R_{n+1} = R_n + \tau V_n + \frac{\tau^2}{2} A_n + O(\tau^3)

V_{n+1} = V_n + \frac{\tau}{2} (A_{n+1} + A_n) + O(\tau^3)

  • R-ben már csak O(\tau^3) pontosságú
  • ez már nem 2 előző lépést használ
  • koordináták és sebességek egyszerre frissülnek
  • sebesség is O(\tau^3) pontosságú

A Metropolis algoritmus

A Metropolis algoritmus egy eljárás, amivel bonyolult valószínűségi eloszlások mintavételezhetőek. Másképp megfogalmazva szeretnénk előállítani egy P eloszlást, a kérdés az, hogy milyen átmeneti valószínűségekkel (W) vezérelt folyamat (egészen pontosan Markov-folyamat) vezet ilyenre? A módszer kihasználja a részletes egyensúly elvét, ez esetünkben azt jelenti, hogy:

P_i W_{ij} = P_j W_{ji}\, teljesül minden  i, j-re.

Ez még kevés feltétel az átmeneti valószínűségekre, ezért a Metropolisz megkötés a következő:

W_{ij} = \mathrm{min}\left[1, \frac{P_j}{P_i} \right]

Bár általánosan tetszőleges eloszlás előállítható, fizikában általában valamilyen egyensúlyi eloszlást szeretnénk előállítani. Például, ha termális egyensúlyt szeretnénk, akkor P \propto \exp\left( -\beta E \right), és a Metropolisz megkötésben szereplő hányados is exponenciális alakban írható. A P-t előállító eljárás a következő:

  1. Induljunk ki egy A konfigurációból.
  2. Állítsuk elő a megváltoztatott B konfigurációt.
  3. Számoljuk ki erre a két konfigurációra a megkötésben szereplő hányadost.
  4. Generáljunk egy 0-1 intervallumon vett egyenletes eloszlású véletlen számot: X.
  5. Ha W_{AB} > X akkor a B konfigurációt fogadjuk el, ellenkező esetben marad az A.
  6. Ismételjük 1-től.

A Monte-Carlo módszer

A Monte-Carlo módszernek nevezzük az olyan eljárásokat, amelyek a problémákat random számok és valószínűségek felhasználásával oldják meg. Az eljárás során ismétlődően kiértékelünk egy determinisztikus modellt, random számokat használva inputnak. Akkor használják, ha a feladat nagyon összetett, nemlineáris, bonyolultak a határfeltételei, vagy nagyon sok paramétertől függ, illetve sok dimenziós.

Használata:

  1. Állítsuk föl a modellt: y = f(x1, x2, ..., xq)
  2. Generáljunk random számokat inputnak: xi1, xi2, ..., xiq
  3. Értékeljük ki a modellt, az eredményt tároljuk el yi-ben
  4. Ismételjük a 2. és 3. lépéseket n-szer
  5. Elemezzük az eredményeket hisztogram, összesítő statisztikák, stb. segítségével
MSc záróvizsga tételek
Tételek Soktest rendszerek | Transzportfolyamatok | Véletlen gráfok generálása, tulajdonságai | Elsőrendű és folytonos fázisátalakulások | Válasz- és korrelációs függvények, fluktuáció-disszipáció tétel | Sztochasztikus folyamatok | A statisztikus fizikai szimulációk alapjai és a Monte Carlo módszer | Dinamikai rendszerek, kaotikus viselkedés | Adatelemzés: lineáris és nem lineáris regresszió egy modellen bemutatva | Adatelemzés: bootstrap modellek | TCP hálózat működése | Adatelemzés: ARCH, GARCH folyamatok | Numerikus módszerek | Vizualizációs módszerek