Entwicklung eines Zufallsgenerators für gleichverteilte Mikrozustände (Implementierung in R)

Nach dem Postulat der statistischen Mechanik besitzen alle Mikrozustände, die ein System annehmen kann, die gleiche Wahrscheinlichkeit. Für zahlreiche Simulationen benötigt man einen Zufallsgenerator, der diese gleichverteilten Mikrozustände erzeugt. Dieser Zufallsgenerator wird in der Programmiersprache R entwickelt, die Erklärungen sind aber so allgemein gehalten, dass man sie leicht in eine andere Programmiersprache übersetzen kann.

Einordnung des Artikels

Die Begriffe Mikrozustand und Makrozustand werden so wie in den vorausgegangenen Kapiteln verwendet.

Der Algorithmus für den Zufallsgenerator beruht auf der Beweismethode "Stars and Bars", die in Spezielle Abzählprobleme: Kombinationen mit Wiederholungen und die Beweismethode Stars and Bars besprochen wurde; die dort erzielten Ergebnisse werden hier verwendet und nicht erläutert.

Einführung

In Konzepte der Statistischen Mechanik: Die Gleichwahrscheinlichkeit der Mikrozustände und die Definition der Boltzmann-Entropie wurde das grundlegende Postulat der statistischen Mechanik genannt: Alle Mikrozustände werden mit gleicher Wahrscheinlichkeit angenommen. Durch dieses Wahrscheinlichkeitsmaß auf der Menge der Mikrozustände wird ein Wahrscheinlichkeitsmaß auf der Menge der Makrozustände induziert. Die Definition der Boltzmann-Entropie ergibt sich dann fast zwangsläufig, wenn man nach einer Größe sucht, die die Eigenschaften der thermodynamischen Entropie besitzt.

Für viele Simulationen in diesem Themenkreis benötigt man dann einen Zufallsgenerator, der Mikrozustände gemäß dem Gleichverteilungs-Postulat erzeugt. Für das Modellsystem, das in den vorausgegangenen Kapiteln beschrieben und verwendet wurde, wird hier ein Zufallsgenerator für gleichverteilte Mikrozustände entwickelt. Seine Implementierung erfolgt in der Programmiersprache R, die Erklärungen des Algorithmus sind aber so ausführlich gehalten, dass es leicht fallen sollte, die Implementierung in einer anderen Sprache vorzunehmen.

In Spezielle Abzählprobleme: Kombinationen mit Wiederholungen und die Beweismethode Stars and Bars wurden meherere Zufallsexperimente besprochen, die äquivalent sind und deren Ergebnisse leicht ineinander umgerechnet werden können. Auch das Erzeugen von Mikrozuständen zu gegebener Teilchenzahl N und Gesamtenergie K kann mit diesen Zufallsexperimenten simuliert werden. Es ist daher nur nötig, eines dieser Zufallsexperimente so zu realisieren, dass die Elementarereignisse gleichwahrscheinlich sind. Unten wird dazu das Ziehen ohne Zurücklegen verwendet, das in R mit der Funktion sample() simuliert werden kann und gleichwahrscheinliche Elementarereignisse liefert.

In diesem Kapitel werden noch keine physikalisch relevanten Anwendungen des Zufallsgenerators gezeigt, die Beispiele unten, insbesondere Beispiel 4, sollen aber einen Eindruck vermitteln, welche Fragestellungen man mit seiner Hilfe leicht bearbeiten kann.

Die Anforderungen an den Algorithmus

Für "kleine" Teilchenzahlen N und Energien K·E0, also wenn sich die zugehörigen Binomialkoeffizienten noch berechnen lassen, kann man sämtliche Mikrozustände zu gegebenem N und K berechnen, sie durchnumerieren und mit Hilfe einer Funktion zum Erzeugen von gleichverteilten Zufallszahlen Folgen von Mikrozuständen erzeugen (in R ist runif() die entsprechende Funktion). Die Gleichverteilung sorgt dann dafür, dass jeder Mikrozustand mit gleicher Wahrscheinlichkeit erzeugt werden kann.

Diese Implementierung ist natürlich ein brute force-Algorithmus, der sehr viel mehr Aufgaben erledigt als nötig. Der große Nachteil ist, dass er für große N und K nicht verwendet werden kann. Daher wird im Folgenden ein schlanker Algorithmus mit folgenden Eigenschaften entwickelt:

Allerdings ist man jetzt in der Beweispflicht, dass der Algorithmus eine Gleichverteilung der Mikrozustände erzeugt – beim oben beschriebenen brute force-Algorithmus ist dies offensichtlich.

Implementierung eines Zufallsgenerators für Makrozustände zu gleichverteilten Mikrozuständen

Die Vorgehensweise des Algorithmus

Die Vorgehensweise ist einfach, wenn man die Beweismethode "Stars and Bars" kennt, die in Spezielle Abzählprobleme: Kombinationen mit Wiederholungen und die Beweismethode Stars and Bars beschrieben wurde und hier nicht erklärt wird.

Stattdessen wird nur kurz der Kern des Algorithmus beschrieben:

Die Implementierung

Es wird zunächst die R-Implementierung des Zufallsgenerators gezeigt und erklärt; die Erklärungen sind aber sehr allgemein gehalten, so dass man sie leicht auf eine andere Programmiersprache übertragen kann. In Abbildung 1 wird der Algorithmus unabhängig von einer Programmiersprache für ein spezielles Beispiel dargestellt.

  1. Zu gegebenem N und K wird eine Folge der Länge N-1+K von "Stars and Bars" erzeugt (im Algorithmus sind es dann die Zahlen 0 und 1). Dies geschieht mit Hilfe der Funktion sample(), die eine Ziehung ohne Zurücklegen simuliert. Dies wird unten ausführlich beschrieben, siehe Erklärung zu Zeile 7.
  2. Die 0-1-Folge wird in einen Mikrozustand umgewandelt, der für die N Teilchen beschreibt, in welchem Energiezustand sie sich befinden: (k1, k2, ..., kN) mit k1 + k2 + ... + kN = K. (Der Mikrozustand ist ein Vektor mit N Komponenten, also eine Partition von K der Länge N)
  3. Der Mikrozustand wird in einen Makrozustand verwandelt, der die Besetzungszahlen der Energieniveaus beschreibt: (n0, n1, ..., nK) mit n0 + n1 + ... + nK = N. (Der Makrozustand ist ein Vektor mit K+1 Komponenten.)

Die Funktion zur Berechnung eines Makrozustandes wird mit randomState() bezeichnet, sie besitzt die Eingabewerte:

randomState <- function(n, k, seed = NULL, micro = FALSE, debug = FALSE){
  if(!is.null(seed)){
    set.seed(seed = seed)
  }
  
  # Kombination mit Wdh.: "Stars and Bars" mit sample()
  sb <- sample( x = c( rep(x = 1, times = n - 1), rep(x = 0, times = k) ),
                size = n + k - 1, replace = FALSE )
  
  # Mikrozust. erzeugen
  microstate <- diff( c(0, which(sb == 1), n+k) ) - 1
  
  # Makrozust. erzeugen
  macrostate <- sapply(X = (0:k), FUN = function(i){return( length( which(microstate == i)) ) })
  
  if(debug){
    cat("Kombination mit Wdh.: \n")
    print(sb)
    cat("Mikrozustand: \n")
    print(microstate)
    cat("Makrozustand: \n")
    print(macrostate)
  }
  
  if(micro){
    return( list(microstate = microstate, macrostate = macrostate) )
  } else {
    return(macrostate)
  }
}

Zur Erklärung:

Zeile 7: Mit der Funktion sample() wird die Ziehung aus einer Urne ohne Zurücklegen simuliert; dies wird durch das Argument replace = FALSE bewirkt. Die in der Urne enthaltenen Objekte x sind zu Beginn: n-1 mal die 1 und k mal die 0. Das Argument size gibt an, wie viele Objekte gezogen werden, hier n-1+k, das heißt es wird so lange gezogen bis die Urne leer ist.

Das Resultat dieses Zufallsexperimentes wird im Vektor sb abgespeichert, also eine 0-1-Folge der Länge n-1+k.

In Abbildung 1 unten wird dies durch Bars | für die 1 und Stars * für die 0 symbolisiert (die Symbole wurden dort so gewählt, um die Verbindung mit der Beweismethode "Stars and Bars" herzustellen).

Wenn in sample() keine Wahrscheinlichkeiten vereinbart werden, wird jedes enthaltene Objekt mit gleicher Wahrscheinlichkeit gezogen. Die Menge, aus der das nächste Objekt gezogen wird, wird bei jedem Zug angepasst; die Wahrscheinlichkeiten sind wieder für alle Objekte gleich.

In Abbildung 1 ist ein Beispiel mit N = 6 und K = 4 gezeigt; die Urne enthält zu Beginn die Objekte | | | | | * * * * und es wird | * | * * | | * | gezogen, was in 1 0 1 0 0 1 1 0 1 übersetzt wird.

Zeile 11: Die rechte Seite enthält vier verschachtelte Operationen, die in Abbildung 1 vermutlich leichter zu verstehen sind als im Quelltext:

  1. Mit which(sb == 1) stellt man fest, an welchen Indizes im Vektor sb eine 1 steht, im Beispiel sind dies die Indizes 1 3 6 7 9 .
  2. An diesen Vektor von Indizes wird am Beginn eine 0 und am Ende eine 10 (allgemein: die Zahl n + k) angefügt; dies erledigt c(0, which(sb == 1), n+k) . Im Beispiel entsteht der Vektor 0 1 3 6 7 9 10
  3. Von diesem Vektor bildet man die Differenzen von je zwei aufeinanderfolgenden Komponenten; dadurch verkürzt sich der Vektor um eine Komponente. Im Beispiel entsteht 1 2 3 1 2 1 . Dies geschieht durch den Aufruf von diff().
  4. Jede der Komponenten wird um 1 verkleinert, es entsteht der Vektor 0 1 2 0 1 0 (in Abbildung 1 sind der dritte und vierte Schritt zu einem Schritt zusammengefasst).

Der zuletzt gebildete Vektor hat N = 6 Komponenten und beschreibt einen Mikrozustand, das heißt jede Komponente beschreibt einen Energie-Zustand eines Moleküls. In mehr mathematischer Sprechweise würde man hier sagen: Mit der Folge | * | * * | | * | wird eine Partition der Zahl K = 4 der Länge N = 6 definiert, nämlich 0 + 1 + 2 + 0 + 1 + 0 = 4.

Zeile 14: Um den Mikrozustand in einen Makrozustand zu verwandeln, muss man abzählen, wie oft die Energieniveaus 0, 1, ..., k im Mikrozustand enthalten sind; dies wird hier mit der Funktion sapply() erledigt. An das Argument X werden die k+1 Zahlen von 0 bis k übergeben, also die möglichen Energieniveaus. Mit Hilfe der Funktion FUN wird festgestellt, wie oft im Mikrozustand microstate diese Energieniveaus vorkommen. Der Rückgabewert von sapply() ist dann ein Vektor der Länge k+1: die Besetzungszahlen der Energieniveaus.

(Funktionen wie sapply() sind sehr typisch für die Programmiersprache R und in anderen Sprachen nicht verfügbar. Dort würde man diesen Schritt als Iteration implementieren: Zu jedem möglichen Energieniveau zählt man ab, wie oft es im Mikrozustand vorkommt und erzeugt mit diesen Besetzungszahlen einen Vektor.)

Abbildung 1: Vorgehensweise des Zufallsgenerators, der oben implementiert ist. Es wird durch eine Ziehung ohne Zurücklegen ein (N-1+K)-Tupel von &quot;Stars and Bars&quot; erzeugt, das eindeutig zuerst in einen Mikrozustand und dann in einen Makrozustand umgerechnet wird (genauere Erklärung im Text). Die Vorgehensweise wird an einer speziellen Realisierung der Ziehung erläutert.Abbildung 1: Vorgehensweise des Zufallsgenerators, der oben implementiert ist. Es wird durch eine Ziehung ohne Zurücklegen ein (N-1+K)-Tupel von "Stars and Bars" erzeugt, das eindeutig zuerst in einen Mikrozustand und dann in einen Makrozustand umgerechnet wird (genauere Erklärung im Text). Die Vorgehensweise wird an einer speziellen Realisierung der Ziehung erläutert.

Einfache Beispiele

1. Beispiel: Ist (wie in Abbildung 1) N = 6 und K = 4 und wird durch die Ziehung mit sample() die Folge 1 1 1 1 1 0 0 0 0 erzeugt, so werden in den Schritten des Algorithmus die Folgen erzeugt:

Indizes der 1:
1 2 3 4 5

Ergänzung um 0 und 10:
0 1 2 3 4 5 10

Bildung der Differenzen, anschließend 1 abziehen (Mikrozustand):
0 0 0 0 0 4

Bestimmung der Besetzungszahlen (Makrozustand):
5 0 0 0 1

Der Mikrozustand ist zu interpretieren: Das 6. Molekül befindet sich im 4. angeregten Zustand, alle anderen Moleküle im Grundzustand.

Der Makrozustand zählt die Besetzungszahlen: Der Grundzustand kommt fünfmal vor, der 4. angeregte Zustand einmal, die anderen Energieniveaus sind nicht besetzt.

2. Beispiel: Wird bei der Ziehung ohne Zurücklegen die Folge 1 0 1 0 1 0 1 0 1 erzeugt, so lauten die weiteren Zahlenfolgen:

Indizes der 1:
1 3 5 7 9

Ergänzung um 0 und 10:
0 1 3 5 7 9 10

Bildung der Differenzen, anschließend 1 abziehen (Mikrozustand):
0 1 1 1 1 0

Bestimmung der Besetzungszahlen (Makrozustand):
2 4 0 0 0

Hier befinden sich zwei Moleküle im Grundzustand und 4 Moleküle im ersten angeregten Zustand, die anderen Energieniveaus sind nicht besetzt.

3. Beispiel: Der Algorithmus kann problemlos für sehr große N und K verwendet werden:

randomState(n = 1000, k = 1000, seed = 42)
# 490 256 132  61  33  12  10   3   3   0 ...

In der zweiten Zeile ist der Rückgabewert, also der Makrozustand gezeigt, wobei alle weiteren Besetzungszahlen gleich 0 sind und nicht aufgeführt sind.

randomState(n = 10000, k = 10000, seed = 42)
# 4997 2529 1235  607  323  143   80   38   23   14    8    0    2    1    0 ...

randomState(n = 100000, k = 100000, seed = 42)
# 49943 24988 12535  6324  3103  1545   818   371   181   105    50    23     7     5     0     0     1     1    0 ....

Der letzte Aufruf erfordert etwas Geduld...

4. Beispiel: In den meisten Anwendungen ist man nicht an den Mikrozuständen sondern an den Makrozuständen interessiert – entscheidend ist nur, dass die Mikrozustände gleichverteilt erzeugt werden. Im folgenden Skript werden sechs Mikrozustände und die zugehörigen Makrozustände zu N = 64 und K = 64 erzeugt und ausgegeben. Zusätzlich werden die Besetzungszahlen der Makrozustände als Histogramm dargestellt.

N <- 64; K <- 64;
n.max <- 16
states <- vector(mode = "list", length = 6)

par(mfrow = c(2, 3))
options(width = 80)

for(i in (1:6)){
  states[[i]] <- randomState(n = N, k = K, micro = TRUE)
  
  cat("i = ", i, "\nMikrozustand:\n")
  print(states[[i]]$microstate)
  
  cat("\nMakrozustand:\n")
  print(states[[i]]$macrostate)
  cat("\n")
  
  macr <- paste0(head(states[[i]]$macrostate, n = 8), collapse = " | ")
  main <- paste0("Makrozustand: ", macr)
  plot(x = (1:n.max) - 1, y = states[[i]]$macrostate[1:n.max],
       col = "red", type = "h", 
       xlim = c(0, n.max), ylim = c(0, 36),
       xlab = "n", ylab = "Besetzungszahlen",
       main = main, frame.plot = TRUE, lwd = 5)
  grid()
}

Zur Erklärung:

Zeile 2: Die Besetzungszahlen gehen schnell gegen null und daher werden in den Histogrammen nur die ersten 16 Besetzungszahlen eines Makrozustandes dargestellt. (Verwendung von n.max in Zeile 22.)

Zeile 3: Die Liste states wird vorbereitet, sie wird später durch die Aufrufe von randomState() gefüllt.

Zeile 5 und 6: Graphik- und Ausgabe-Parameter.

Zeile 8 bis 26: In der for-Schleife werden die Zustände erzeugt, ausgegeben und die Diagramme erstellt.

Zeile 9: Aufruf von randomState() ; hier darf man seed nicht setzen, da sonst 6 identische Zustände erzeugt werden. Soll das Skript reproduzierbar sein, muss set.seed() schon vor den Aufrufen von randomState() gesetzt werden.

Zeile 11 bis 16: Ausgaben.

Zeile 18 und 19: Die ersten 8 Besetzungszahlen des Makrozustandes werden für die Bildüberschrift main verwendet.

Zeile 20: Da in R die erste Komponente eines Vektors den Index 1 hat, die Besetzungszahl des Grundzustandes aber k0 lautet, wird für die x-Werte x = (1:n.max) - 1 gesetzt.

Die Ausgaben lauten:

# i =  1 
# Mikrozustand:
#   [1] 0 2 2 1 0 3 1 0 0 0 1 0 2 1 2 1 1 0 0 0 1 5 0 2 1 1 2 0 1 1 0 0 1 1 0 0 1 0
# [39] 4 1 1 2 0 1 1 0 0 0 1 2 1 0 1 1 1 1 7 2 1 1 0 0 1 0
# 
# Makrozustand:
#   [1] 24 27  9  1  1  1  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [26]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [51]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# 
# i =  2 
# Mikrozustand:
#   [1] 3 0 2 1 3 1 1 0 0 0 2 0 1 1 0 0 2 0 0 0 0 0 1 0 3 1 0 0 1 2 1 0 0 1 1 1 1 0
# [39] 0 5 1 0 3 0 3 2 0 1 4 0 0 1 1 0 0 3 3 1 1 3 0 1 1 0
# 
# Makrozustand:
#   [1] 28 21  5  8  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [26]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [51]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# 
# i =  3 
# Mikrozustand:
#   [1] 0 1 1 1 0 1 2 3 2 0 2 0 0 2 1 2 0 0 0 0 0 0 0 1 0 1 0 2 3 1 0 0 4 0 0 0 0 2
# [39] 0 0 2 1 0 2 1 1 0 1 0 1 2 6 4 0 0 0 0 2 0 1 0 3 3 2
# 
# Makrozustand:
#   [1] 31 14 12  4  2  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [26]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [51]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# 
# i =  4 
# Mikrozustand:
#   [1] 1 2 1 0 0 1 0 1 1 0 0 2 0 0 0 2 1 0 0 1 0 0 1 3 5 2 0 0 0 1 2 0 0 0 0 0 0 1
# [39] 0 0 0 2 1 0 1 0 2 0 0 0 0 4 0 8 0 1 3 0 1 3 1 6 1 2
# 
# Makrozustand:
#   [1] 33 16  8  3  1  1  1  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [26]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [51]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# 
# i =  5 
# Mikrozustand:
#   [1] 3 1 1 0 1 3 1 2 0 0 0 4 0 1 0 0 0 0 1 0 0 1 1 1 0 0 0 0 0 0 0 0 0 8 0 5 3 1
# [39] 2 0 1 2 2 0 1 0 0 2 1 3 0 0 0 0 1 0 1 0 4 1 2 0 0 3
# 
# Makrozustand:
#   [1] 33 16  6  5  2  1  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [26]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [51]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# 
# i =  6 
# Mikrozustand:
#   [1] 0 1 2 4 1 1 1 0 3 0 0 2 0 1 0 0 3 0 2 3 1 1 0 0 0 3 0 0 0 0 0 2 1 2 1 0 2 0
# [39] 0 2 0 1 0 2 1 0 3 1 0 3 0 2 0 2 1 1 1 0 4 0 0 3 0 0
# 
# Makrozustand:
#   [1] 30 15 10  7  2  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [26]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
# [51]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

In Abbildung 2 werden die Besetzungszahlen der oben erzeugten Makrozustände dargestellt. Eigentlich hat jeder Makrozustand 65 Besetzungszahlen, da sie aber schnell gegen 0 gehen, werden nur die ersten 16 dargestellt. Um die Diagramme untereinander vergleichen zu können, sind die y-Achsen identisch skaliert, siehe ylim = c(0, 36) in Zeile 22.

Abbildung 2: Die mit dem Zufallsgenerator erzeugten Mikrozustände werden in Makrozustände umgerechnet und die Besetzungszahlen gegen das Energieniveau aufgetragen.Abbildung 2: Die mit dem Zufallsgenerator erzeugten Mikrozustände werden in Makrozustände umgerechnet und die Besetzungszahlen gegen das Energieniveau aufgetragen.

Die Gleichwahrscheinlichkeit der Mikrozustände

Es bleibt noch die für den Einsatz des Zufallsgenerators entscheidende Frage zu klären: Sind die Mikrozustände gleich verteilt?

Die Funktion sample() kann – so wie sie hier konfiguriert ist – zu gegebenem n und k jede Kombination von "Stars and Bars" erzeugen. Und so wie die Wahrscheinlichkeiten in sample() gesetzt werden, ist jede Kombination gleich wahrscheinlich. Und da die Umrechnung in einen Mikrozustand eindeutig ist, sind auch die Mikrozustände gleich verteilt. (Für die Makrozustände gilt dies natürlich nicht, aber das ist gerade der Clou an der Definition der Boltzmann-Entropie.)

Der Beweis, dass alle Mikrozustände mit gleicher Wahrscheinlichkeit erzeugt werden, beruht somit darauf, dass die Funktion sample() für die Ziehung ohne Zurücklegen die Wahrscheinlichkeiten nach jedem Zug neu berechnet. Wenn es eine derartige Funktion in anderen Programmiersprachen nicht gibt, muss man dies erst implementieren.