Vergleich numerischer Integrationsschemen am Beispiel des Magnetpendels

Über den Einfluß der Schrittweite auf numerische Simulationen

---

Die Idee zu dem Artikel ist schon etwas älter und kam mir beim schreiben des Artikels zum Magnetpendel. Ich werde hier anhand von zwei Simulationsbeispielen verschiedene numerische Verfahren vergleichen. Im ersten Teil werde ich dazu am Beispiel des gedämpften harmonischen Oszillators genauer auf den Einfluß der Schrittweite auf den Gesamtrechenfehler eingehen. Der zweite Teil beschäftigt sich erneut mit dem Magnetpendel. Das dort beschriebene Problem ist chaotischer Natur und ich fand es interessant einmal an einem solchen Beispiel die verschiedene Integrationsverfahren zu vergleichen. Dieser Artikel soll den Leser für die Problematik der Zeitschrittwahl sensibilisieren und vor allem mit dem häufig anzutreffenden Irrtum aufäumen man müßte die Zeitschritte nur klein genug wählen um letztendlich ein genaues Ergebnis zu bekommen. Auf die mathematischen Hintergründe werde ich hier nur insofern eingehen, als das dies zum unbedingten Verständnis notwendig ist.

Um diesen Artikel zu verstehen sind Grundkenntnisse der in numerischer Integration von Vorteil. Es wäre hilfreich, wenn Leser bereits mit Euler, Runge-Kutta und Adams-Bashforth Verfahren vertraut sind da ich diese Verfahren hier nicht näher erläutern werde. Wer mehr über Runge-Kutte Verfahren lernen möchte, dem empfehle ich folgenden Artikel:

• Das Runge Kutta Verfahren
  Sehr detailierte Einführung in das Runge-Kutta Verfahren 4. Ordnung

Genauigkeit von Integrationsverfahren

Die Genauigkeit eines Integrationsverfahrens wird durch den lokalen Fehler und den globalen Fehler bestimmt. Der lokale Fehler ist der Fehler, der bei einem einzelnen Integrationsschritt auftritt. Der globale Fehler ergibt sich aus der Summe der lokalen Fehler von vielen Integrationsschritten. Normalerweise werden Landau-Symbole verwendet um den lokalen und globalen Fehler anzugeben. Man spricht auch von der Ordnung eines Verfahrens (engl.: big O notation). Der maximale Fehler eines Verfahrens der Ordnung p (für kleine Schrittweiten h) ist eine Funktion proportional zu h^p oder mit Landau-Symbolen: O(h^p). Prinzipiell gilt: Je größer die Ordnung p, desto besser ist das Integrationsschema.

Singlestep methods Method name Abbreviation Local error Global error Euler Euler O(h^2) O(h) Modified Euler Euler mod. O(h^3) O(h^2) Heun Method Heun O(h^3) O(h^2) Runge-Kutta 3th order RK4 O(h^4) O(h^3) Runge-Kutta 4th order RK4 O(h^5) O(h^4) Runge-Kutta 5th order RK5 O(h^6) O(h^5) Multistep methods Method name Abbreviation Local error Global error Adams-Bashforth (2 Steps) ADB2 O(h^3) O(h^2) Adams-Bashforth (3 Steps) ADB3 O(h^4) O(h^3) Adams-Bashforth (4 Steps) ADB4 O(h^5) O(h^4) Adams-Bashforth (5 Steps) ADB5 O(h^6) O(h^5) Adams-Bashforth (6 Steps) ADB6 O(h^7) O(h^6)

Tabelle 1: Liste aller getesteten Integrationsverfahren und ihre Ordnung.

Bild 1: Absoluter Fehler als Summe von Rundungsfehler und Verfahrensfehler.

Das ultimative Ziel einer jeden Simulation ist es, genaue Ergebnisse in einer akzeptablen Zeit zu erlangen. Zunächst steht man also vor der Wahl des Integrationsschemas. Danach muß man eigentlich die Simulation nur noch mit einer akzeptablen Zeitschrittweite durchführen. Das prinzipielle Problem dabei ist allerdings, daß genauere Integrationsschemen mehr Funktionsauswertungen benötigen als einfachere. Dadurch erhöht sich deren Rechenaufwand. Auf der anderen Seite sind aber Verfahren hoher Ordnung selbst bei größeren Zeitschrittweiten genauer als Vergleichbare Verfahren niederer Ordnung. Man hat also die Wahl zwischen einem einfacheren Verfahren mit kleiner Zeitschrittweite und einem aufwendigeren Verfahren mit größerer Zeitschrittweite. Aber Vorsicht: Die Annahme man kann selbst mit Verfahren niederer Ordnung ein "richtiges" Ergebnis errerchnen, wenn man nur die Zeitschrittweite klein genug wählt ist falsch. Um das zu Untermauern sehen wir uns den sogenannten globalen Fehler der Integrationsverfahren mal genauer an. Die folgenden Eigenschaften eines Integrationsschemas sind direkt von der Zeitschrittweite abhängig:

• Genauigkeit der Lösung
  • Rundungsfehler ( ~ 1/h)
  • Globaler Fehler (~ h^p)
• Rechenzeit (~1/h)
• Numerische Stabilität

Bild 1 zeigt den absoluten Fehler als Summe von Rundungsfehler und globalem Fehler des Verfahrens. Die Darstellung erfolgt auf einer logarithmischen Skala. Globaler Fehler und Rundungsfehler Verhalten sich gegenläufig. Während der globale Fehler für kleine Schrittweiten abnimmt, steigt der Rundungsfehler an. Addiert man die beiden Fehler, so wird deutlich, daß es eine optimale Schrittweite gibt, für die der absolute Fehler minimal wird. Daraus folgt, daß eine zu kleine Schrittweite das Ergebnis infolge der dann verstärkt auftretenden Rundungsfehler verschlechtern und nicht verbessern. Soviel zur Theorie, doch kann man das wirklich an einem Beispiel beobachten? Finden wir es raus...

Beispielrechnungen

Der gedämpfte harmonische Oszillator

Wenn man eine Integrationsengine implementiert, ist das erste, was man damit tun sollte sie gegen ein bekanntes Modell, am besten eines dessen analytische Lösung man kennt zu testen. In diesem Falle dient das Modell des gedämpften harmonischen Oszillators als Referenz. Dieses Modell wird durch eine homogene Differentialgleichung 2. Ordnung beschrieben:

Mit: m - Pendelmasse r - Reibungskoeffizient k - Federkonstante Es gelten folgende Anfangsbedingungen: Das Ergebnis dieses Differentialgleichung ist eine Funktion, welche die Position des Oszillators als Funktion der Zeit angibt: wobei:

Bild 2: Schematische Darstellung einer gedämpften Schwingung. (Bildquelle: Wikimedia commons [1] und [2])

Die Differenz zwischen analytischer Lösung (siehe Bild 3) und simulierter Lösung ist der absolute Fehler. Bild 4 zeigt den Zusammenhang zwischen absolutem Fehler und Zeitschrittweite für verschiedene Integrationsverfahren.

Bild 3: Vergleich zwischen analytischer und simulierter Lösung für eine Zeitschrittweite von h=0.005. Alle Verfahren, außer dem von Euler zeigen gute Übereinstimmung mit der exacten Lösung.	Bild 4: Absoluter Fehler in Abhängigkeit der Integrationsschrittweite für verschiedene Integrationsverfahren.

Es wird deutlich, daß der theoretische Zusammenhang zwischen absolutem Fehler und Zeitschrittweite bestätigt wird. Man kann klar sehen, daß es eine optimale Zeitschrittweite für jedes der Verfahren gibt, bei der der absolute Fehler minimal ist. Man sieht auch, daß die Verfahren höherer Ordnung, wie z.B. RK5 eindeutig am genauesten rechnen. Um so höher die Ordnung der Verfahren, um so geringer werden die absoluten Fehler. Es wird auch deutlich, daß die Verfahren höherer Ordnung das Minimimum des absoluten Fehlers erst bei größeren Schrittweiten erreichen. Wenn man sich auf der anderen Seite nun das Eulerschema ansieht, sieht man, daß das Minimum seines absoluten Fehlers bei sehr kleinen Schrittwiten anzutreffen ist. Obwohl diese Schrittweiten um Größenordnungen kleiner sind, als die der anderen Verfahren ist sein Fehler dennoch am größten. Wenn man dazu noch kurz auf Bild 3 schaut, sieht man auch das das Eulerverfahren selbst bei moderaten Schrittweiten versagt. Was lernen wir daraus? Wie wärs damit:

Eulers Methode war zwar ziemlich Visionär für das seine Zeit (1768) aber sie ist einfach nicht mehr das, was man heute noch verwenden sollte. Darüber hinaus wird nun auch am Beispiel klar, daß zu kleine Schrittweiten nicht nur eine ziemliche Verschwendung von Computerresourcen sind, sondern das Ergebnis auch noch schlechter machen.

Die Betrachtung der Genauigkeiten ist ja schön und gut, aber das Verfahren höherer Ordnung genauer sein sollten, liegt ja eigentlich auf der Hand. Was uns wirklich interessiert ist ja, ob sie bei vergleichbarer (oder besserer) Genauigkeit auch schneller sind. Sehen wir uns also mal das Timing genauer an. Um dies zu analysieren, ist in Bild 5 die Abhängigkeit zwischen Zeitschrittweite und Rechenzeit aufgetragen. Wir sehen einen einfachen Zusammenhang, der aussagt, daß die Rechenzeit mit zunehmende Integrationsschrittweite sinkt. Das ganze dargestellt in einem Satz fast perfekter Geraden, die allerdings keinen linearen Zusammenhang beschreiben, da das Diagramm eine doppelt logarithmische Skalierung verwendet. Wenn man die Zeitschrittweite halbiert, verdoppelt man gleichzeitig die Anzahl der zu berechnenden Funktionsauswertungen. Man zahlt also einen ziemlich hohen Preis für kleine Zeitschritte unabhängig davon welche Integrationsmethode man wählt. Da jeder Zeitschrittweite nun eine Rechenzeit für das Modell entspricht kann man abschließend ein Diagram erstellen, welches das Quadrat des absoluten Fehler über der Zeitschrittweite darstellt (Bild 6). Hier wird es jetzt interessant, denn man kann jeder Genauigkeit direkt eine Rechenzeit zuordnen.

Bild 5: Abhängigkeit der Integrationsschrittweite von der Rechenzeit. Wenig überraschend sinkt die Rechenzeit bei Verwendung größere Zeitschrittweiten. Daß Einschritt und Mehrschrittverfahren hier so nah beieinander liegen liegt daran, daß für diese Simulation nicht viel berechnet werden muß, der zusätzliche Rechenaufwand fällt also nicht stark ins Gewicht.

Bild 6: Um zu bestimmen, welches Verfahren bei größtmöglicher Genauigkeit am schnellsten rechnet wird hier der globale Fehler über der Rechenzeit aufgetragen. Fast jede der Kurven hat ein Minimum des globalen Fehlers. Je genauer ein Verfahren ist, desto tiefer liegt sein Minimum. Das beste Verfahren ist dasjenige, das bei kleinstmöglichem Fehler am schnellsten rechnet. Interessanterweise liegen in diesem Fall ADB6 und RK5 gleichauf.

Wenn man jetzt herausfinden möchte, welches Verfahren die Lösung bei vorgegebener Genauigkeit am schnellsten errechnet, so muß man die entsprechenden Daten nur aus obigem Diagramm ablesen. Fodert man zum Beispiel die Quadratsumme des globalen Fehlers soll kleiner als 1e-5 sein, so ergibt sich folgende Tabelle:

Rank	Integration Scheme	Calculation time in ms	Time step width
1	RK 5	0.01	0.02
2	RK 4	0.015	0.01
3	ADB 5 (* see below)	0.03	0.005
4	ADB 4	0.039	0.0035
5	ADB 3	0.09	0.0015
6	Euler mod. Heun	0.5	0.0003
7	RK3 ADB2	0.7	0.0003
9	Euler	36812	4e-9

Der Sieg geht an: Das Runge-Kutta Verfahren 5. Ordnung (RK5)! Es ist in diesem Test überlegen, da es größere Zeitschrittweite erlaubt ohne dabei zu viel an Genauigkeit einzubüßen. Ein starker zweiter Platz geht an das Runge-Kutta Verfahren 4. Ordnung (RK4), welches immer noch besser abschneidet als das Adams-Bashfort Verfahren 5. Ordnung (ADB5). Das ist insofern überraschend, als ADB5 als Mehrschrittverfahren deutlich weniger Funktionsauswertungen benötigt. Bei näherer Betrachtung wird aber klar, daß die Berechnungen für diese Pendelsimulation sehr einfach sind und nur minimalen Aufwand bereiten. Vereinfacht gesagt können Mehrschrittverfahren ihre Vorteile erst dann ausspielen, wenn es etwas zu berechnen gibt. Die Ergebnisse währen also für ein ein anderes Problem (z.B. Dreikörperproblem) anders. Man sollte also die Mehrschrittverfahren nicht auf Basis dieser einfachen Beispielrechnung abschreiben. Betrachtet man nun noch die anderen Verfahren, so fällt noch di für mich enttäuschende Leistung von RK3 auf, das selbst schlechter als das modifizierte Euler Verfahren und das Verfahren von Heun abschneidet. Die Ursache hierfür liegt daran, daß es zwar mehr Funktionsauswertungen benötigt aber keine signifikant größeren Zeitschrittweiten zuläßt.

Das magnetische Pendel

Das physikalische Model

Der nächste Test wird sich nochmals mit dem magnetischen Pendel beschäftigen. Die physikalischen Prinzipien dazu werden nochmals kurz im Video 1 erläutert. Das für die Simulation verwendete Modell ist in Bild 5 dargestellt. Die Simulation wird diesmal jedoch mit verschiedenen Integrationsverfahren durchgeführt. Die Ergebnisse werden dann kurz analysiert und Besprochen.

Video 1: Chaos Theorie erklärt am Beispiel des Magnetischen Pendels ; © 2007 Aldo Cavini Benedetti	Bild 7: Das dem zweiten Test zugrunde liegende Pendelmodell.

Doch bevor wir damit anfangen sollten wir uns mal zwei Punkte der Simulation genauer ansehen. Die Pendelbewegung ist zwar chaotisch aber chaotisch bedeuted ja nicht man kann sie nicht berechnen, es bedeuted lediglich, daß der Rechenfehler exponentiell mit der Zeit anwächst. Man kann also vereinfacht gesagt nicht besonders weit in die Zukunft rechnen, weil früher oder später kleine Variationen in den Anfangsbedingungen Anwachsen werden und letztendlich zu völlig unterschiedliche Ergebnissen führen werden. Der nächste Test wird die Spuren der verschiedenen Integrationsverfahren vergleichen. Alle Verfahren werden mit identischen Anfangsbedingungen initialisiert, so daß sämtliche Abweichungen nur durch die die Ungenauigkeiten der verschiedenen Integrationsverfahren verursacht werden sollten. Wie man in Bild 8 sehen kann, gibt es dabei einige Punkte, in denen es praktisch keinen Unterschied zwischen den verschiedenen Verfahren gibt, wohingegen in anderen Punkten praktisch jedes Integratinsverfahren zu einem anderen Endpunkt konvergiert (Bild 9).

Bild 8: Es gibt Startpunkte von denen aus alle Integrationsverfahren nahezu identische Ergebnisse liefern.	Bild 9: Für andere Startpunkte liefern fast alle Integrationsverfahren unterschiedliche Ergebnisse.

Simulationsergebnisse

Vor diesem Hintergrund ist es interessant nun zu sehen, wie eine vollständige Berechnung aussieht. Kommen wir nun also zu den Ergebnissen. Folgende Bilder wurden mit einer Schrittweite von h=0.015 s berechnet:

Bild 10: Übersicht über die Ergebnisbilder der verschiedenen Integrationsverfahren.

Große Teile des Simulationsgebietes zeigen ein chaotisches Verhalten, der Mittelteil und einige Strukturen sind in allen Simulationen reproduzierbar. Die gute Nachricht ist: Alles sieht ziemlich ähnlich aus! Vielleicht sollten wir mal nachschauen, wie ähnlich genau:

Bild 11: Schwarzweißbild in dem der Wert eines Pixels angibt, wieviele verschiedene Ergebnisse die Integrationsverfahren lieferten. Weiß bedeuted alle Integrationsverfahren lieferten das gleiche Ergebnis, schwarz bedeuted alle lieferten ein anderes Ergebnis.	Bild 12: Ergebnisbild, in dem Pixel nur angezeigt werden, wenn sie in mindestens 9 von 10 Integrationsverfahren gleiche Ergebnisse erzielt haben. Die Pixelfarben sind dem RK5 Verfahren entnommen.

---

Autor: Ingo Berg; Dieser Text und die dazugehörigen Mediendateien stehen, sofern nicht anderweitig angegeben, unter der Creative Commons Namensnennung 3.0 Deutschland Lizenz.