Im Jahre 1967 publizierten J. Orban und A. Bellemans eine etwa 1-seitige Arbeit [1], in der sie auf numerische Instabilitäten bei Molekular-Dynamik-Simulationen (MDS) hinwiesen. Solche Instabilitäten stellen die Gültigkeit und Ausagekraft mancher MDS in Frage, die für entropische Analysen durchgeführt werden. Hat man seine eigene Forschung dem "Zeitreversibilitäts-Paradoxon" verschrieben, so muss sorgfältig auf die Invarianz der Simulation unter Zeitumkehr geachtet werden. Dieses Problem hat man prinzipiell heute im Griff, wenn der richtige Algorithmus für die Simulationen verwendet wird.

Zum Verständnis der angesprochenen Problematik sei zunächst daran erinnert, dass die Grundgleichungen konservativer dynamischer Systeme, also die Hamiltonschen Gleichungen, invariant unter Zeitumkehr sind. Das bedeutet für ein typisches N-Teilchen-Problem, dass keine Zeitrichtung ausgezeichnet ist und daher bei einer Impulsumkehr aller beteiligter Teilchen zu einem bestimmten Zeitpunkt die zurückliegende Trajektorie wieder genau rückwärts durchlaufen wird. Eine beliebige Funktion der Mikrozustände sollte dann naturgemäß nach einer Impulsumkehr ebenfalls rückwärts reproduziert werden können. Orban und Bellemans haben in ihrer Simulationen eines Gases beispielsweise die Boltzmannsche H-Funktion berechnet und erhielten die in der Abbildung gezeigten Ergebnisse:

Abb. aus [1]

Die H-Funktion nimmt im Laufe der Zeit ab und strebt asymptotisch einem Minimum zu - ihrem Geleichgewichtswert. Letzterer wird erreicht, wenn das Gas vollständig relaxiert ist, was i.a. sehr schnell geht. Die Autoren kehrten nun nach 10 bzw. 20 Zeiteinheiten alle Impulse um, und berechneten weiterhin die H-Funktion, die nun eigentlich zu einem symmetrischen Verlauf bezüglich t=10 bzw. t=20 führen müsste. Mit anderen Worten: das System müsste sich aus dem Gleichgewichtszustand heraus in den initialen Nicht-Gleichgewichtszustand zurück bewegen. Tatsächlich ist die Symmetrie bei einer Impulsumkehr bei t=10 in Teil (a) der Abbildung nahezu erfüllt. Bei einer Impulsumkehr zum Zeitpunkt t=20 jedoch, erkennt man eine deutliche Abweichung von der theoretisch zu erwartenden Symmetrie. In den Teilabbildungen (b) und (c) sind auch bei frühzeitiger Impulsumkehr die H-Funktionen kaum zu repoduzieren. Der Unterschied in den Teilabbildungen liegt darin, dass der berechneten Teilchentrajektorie zufällige Fehler überlagert wurden, und zwar in der relativen Größe von (a)10^-8, (b)10^-5 und (c)10^-2. Die Tatsache, dass bei einer Überlagerung von Rauschen eine Entwicklung des Systems zurück zum Nicht-Gleichgewichtszustand unterdrückt wird, selbst dann, wenn man es in entsprechenden Anfangszuständen präpariert (wie im vorliegenden Falle), läßt sich im Lichte einer Arbeit von James Hurley [2] gut verstehen. Wir wollen an dieser Stelle nicht näher darauf eingehen. Das eigentlich erschreckende an Orban und Bellemans Ergebnis ist, dass unter Verwendung eines Standardalgorithmus für Differentialgleichungen, implementiert auf der Basis einer Floating-Point-Arithmetik auf einem Digital-Computer, solches Zufallsrauschen gar nicht vermeidbar ist, da Rundungen nichtlinearer Natur sind und sich als Rauschen bemerkbar machen. Ein konservatives Hamiltonsches System wird durch eine gängige Implementierung auf einem Computer zu einem offenen System, d.h., es wird eine numerische Dissipation überlagert. Die Resultate beziehen sich auf das offene, nicht auf das zu untersuchende konservative System.

Das vorliegende Problem ist ganz eng mit jenen verwandt, die im Zusammenhang mit der Chaosforschung eine zeitlang für Unruhe sorgten. Noch heute gibt es Kritiker an der Chaosforschung, die von computergenierten Artefakten reden. Tatsächlich hat sich auf analytischem und "real-"experimentellem Wege zeigen lassen, dass chaotische Phänomene existieren. Für die Rehabilitierung der "Computerexperimente" sorgten insbesondere die Ergebnisse des Numerikexperten Peter Kloeden und seiner Kollegen [3]. Sie konnten zeigen, dass es eine Klasse von Algorithmen gibt, die bestimmte invariante Strukturen im Phasenraum - z.B. chaotische Attraktoren - bis auf eine beliebig kleine Umgebung reproduzieren können. Mit dem selben Rezept lassen sich analog so genannte symplektische Algorithmen für die Integration Hamiltonscher Systeme ableiten [4]. Diese sind so konstruiert, dass sie gewisse Invarianten bis auf eine beliebig kleine Abweichung erhalten. Die im MDS-Kontext entscheidenden Invarianten sind Energie und Impuls. Ein bedeutsamer Beitrag zur algorithmischen Konsistenz bei MDS stellt die Ableitung eines symplektischen Algorithmus aus dem Prinzip der kleinsten Wirkung dar [5]. Die Anwendung des Extremalverfahrens auf die diskrete Raum-Zeit (des Digital-Computers) liefert eine "diskrete Newtonsche Bewegungsgleichung". Das Resultat ist ein Algorithmus, der den Impuls exakt erhält, keine säkulare Energiedrift erzeugt und exakt-reversibel ist. Die zündende Idee geht auf Gillilan und Wilson [6] zurück, die das Extremalverfahren allerdings nur auf die diskretisierte Zeit anwandten. Dass zudem die Implementierung unter Verwendung der Integer-Arithmetik stattfinden muss, erkannten unabhängig Levesque und Verlet [7] sowie Diebner [8].

Wir können also zusammenfassen, dass konsistenten numerischen Untersuchungen von Hamiltonschen N-Teilchen-Systemen nichts mehr im Wege steht. Die Auflösung des Orban und Bellemans'schen Problems, also die Reproduzierbarkeit der H-Funktion nach Impulsumkehr, wurde in [7, 8] gezeigt. Die im ersten der Beiträge dieser Web-Site beschriebene Entropieuntersuchung [9,10] fand unter Verwendung des symplektischen, exakt-reversiblen Algorithmus statt und findet daher in seiner Aussagekraft nachträglich zusätzliche Rechtfertigung. Eine umfangreiche Bearbeitung dieses und verwandter Themen ist in [11] zu finden.

Literatur:

[1] J. Orban and A. Bellemans:

[6] R.E. Gillilan and K.R. Wilson:

[10] Molekular-Dynamik-Simulationen I

[11] Hans H. Diebner: