13. Kapitel

Die Familie der RUNGE-KUTTA-Verfahren

Alle Verfahren dieser Familie gehören zu den Einschritt-Verfahren. Bei vorgegebener Ordnung p kann die Zahl der erforderlichen Stufen s ≥ p frei gewählt werden. Dadurch ergibt sich eine große Anzahl verschiedener Verfahren.

Um zu verstehen, wie diese Verfahren funktionieren, beschränken wir uns zunächst auf die Lösung einer Differentialgleichung 1. Ordnung. Die Erweiterung auf solche 2. Ordnung, also Bewegungs-Differentialgleichungen, werden wir im Anschluss anhand konkreter Lösungsverfahren besprechen.

Zunächst also das Lösungsprinzip: Innerhalb eines Schrittes werden, um einen neuen Wert aus einem bereits bekannten Wert y_i zu berechnen, s Berechnungsschritte ausgeführt. Dabei werden Zwischenergebnisse erzeugt, die entsprechend a aufsummiert, den neuen Wert y_i+1 ergeben. Für die Summation werden die Zwischenresultate k_j mit den Koeffizienten b_j gewichtet. Die so erhaltene Summe wird nun noch mit der Schrittweite Δt multipliziert und dem alten Wert y_i zu addiert.

(a)

Zeitabhängige Gradientenfunktionen f(t,y) wären zu erwarten, wenn sich die Parameter der Differentialgleichung, z.B. die Gravitationskonstante, der Reibungskoeffizient oder die Federkonstante mit der Zeit verändern würden oder, was viel wahrscheinlicher ist, die Störfunktion einer inhomogenen Differentialgleichung zeitabhängig (z.B. sinusförmig periodisch) ist.

Das besondere an diesem Verfahrenstyp ist nun die Erzeugung der s Zwischenergebnisse. Dafür werden für s unterschiedliche Zeitpunkte t_j (sofern die Gradientenfunktion f(t,y) Zeit-abhängig ist, was bei uns aber sehr selten der Fall ist) und s unterschiedliche Werte der gesuchten Variablen y innerhalb des aktuellen Intervalls berechnet:

(b)

on/off

worin t_j und y_j nach

(c)

und

(d)

berechnet werden.
Die aktuelle Zeit ergibt sich in bekannter Weise zu

(c)

Die Koeffizienten b_j, c_j und a_j,l werden der Stufenzahl s und der gewählten Ordnung p entsprechend berechnet. Diese Koeffizienten sind Konstanten und können der Literatur entnommen werden.
Sehr theoretisch! Um Licht ins Dunkel zu bringen, schauen wir uns gleich ein bekanntes Beispiel an, nämlich

Das klassische RUNGE-KUTTA-Verfahren (RK4)

Das klassische RUNGE-KUTTA-Verfahren ist ein 4-stufiges, explizites Verfahren und damit ein Spezialfall der Verfahrensfamilie der RUNGE-KUTTA-Verfahren. Es zeichnet sich durch eine besonders geringe Anzahl, von 0 verschiedener, Koeffizienten, die zudem noch Brüche einfacher, ganzer Zahlen sind, aus. Dadurch konnten in früheren Zeiten Differentialgleichungen auch ohne Computer numerisch gelöst werde. Das Verfahren hat die Konsistenzordnung p = 4.

Der besseren Anschaulichkeit halber bleiben wir bei der Lösung einer Differentialgleichung 1. Ordnung. Gegeben sei die Gradientenfunktion dy/dt = f(t,y). Dann erfolgt die Berechnung des neuen Wertes y_i+1 entsprechend zu:

()

mit den folgenden Zwischenschritten:

()

Betrachten wir jetzt die Gleichungen . Wir sehen, dass die Gradientenfunktion f(t,y) an 4 verschiedenen Stellen berechnet werden muss, daher handelt es sich um ein 4-stufiges Verfahren! Die Zwischenschritte k_s werden entsprechend der Gewichte, die durch die Koeffizenten laut b) zu berechnen sind, gegeben. Zur besseren Übersichtlichkeit stellen wir die Gewichte zu Spaltenvektoren zusammen. Die b_j nach a):

()

und die c_j nach c):

()

und die a_{j, l} nach d) als Matrix A:

()

Da es, wie schon erwähnt, eine Vielzahl unterschiedlicher RUNGE-KUTTA-Algorithmen gibt, ist eine vereinfachte Schreibweise für die Koeffizienten üblich. Das sog. BUTCHER-Tableau. Es stellt eine schematische Anordnung der Koeffizienten in einer Art Tabelle dar. Die grundlegende Anordnung zeigt :

()

Speziell für das klassische RUNGE-KUTTA-Verfahren nimmt das Butcher-Tableau folgende Gestalt an:

()

Da die rechte Hälfte der durch die Diagonale geteilten A-Matrix stets gleich 0 ist, hat sich die vereinfachte Schreibweise nach durchgesetzt:

()

Nun wollen wir überlegen, wie wir das für Differentialgleichungen 1. Ordnung erläuterte RUNGE-KUTTA-Verfahren auch auf Differentialgleichungen 2. Ordnung anwenden können. Erst diese Lösungsmethode kann uns die mgesuchten Bewegungsgleichungen liefern!
Im Kapitel Verfahren von EULER-CAUCHY für DGLn. 2. Ordnung haben wir die Vorgehensweise ausführlich geschildert. Das Vorgehen besteht darin, dass die Differentialgleichung 2. Ordnung in ein System von zwei Differentialgleichungen 1. Ordnung überführt wird ().

Wir gehen wieder davon aus, dass die Gradientenfunktion f(s, v) nicht zeitabhängig ist, was zu einer Vereinfachung der Lösung des DGl.-Systems führt.
Gemäß haben wir nun zwei Differentialgleichungen zu lösen. Entsprechend der Vorgabe des 4-stufigen RUNGE-KUTTA-Verfahrens sind die Orts-DGl.

(a)

und die Geschwindigkeits-DGl.

(b)

zu lösen.
In jeder der Gleichungen befindet sich jetzt ein Satz von k_j-Koeffizienten, die aber je nach Bestimmung andere sind. Ihrer physikalischen Bedeutung nach sind die k_s,j Geschwindigkeiten und die k_v,j Beschleunigungen.
Es gibt k_sJ-Koeffizienten für die Ortsberechnung und k_vJ-Koeffizienten für die Geschwindigkeitsberechnung! Ihre Berechnung erfolgt in wechselseitiger Abhängigkeit von Ort und Geschwindigkeit, wobei sich diese Abhängigkeiten stets auf die folgende Stufe der Rechnung beziehen:

()

So gesehen sind diese Schritte jeweils einem symplektischen EC-Lösungsschritt vergleichbar.
Die Berechnung der Koeffizienten des Ortes in ist vergleichsweise einfach, da die erste Differentialgleichung (a) naturgemäß keine Ortsabhängigkeit aufweist, es fließt nur die Geschwindigkeit v ein. Die zweite Differentialgleichung (b) erfordert einen höheren Aufwand, da die Gradientenfunktion sowohl vom Ort (z.B. Federauslenkung) und der Geschwindigkeit (z.B. Luftreibung) abhängig sein kann.
So sieht die Implementierung in p5.js aus:

Abb. klassisches 4-stufiges RUNGE-KUTTA-Verfahren

Und so wird das RUNGE-KUTTA-Verfahren angewandt:

Abb. Anwendung des 4-stufiges RUNGE-KUTTA-Verfahren auf ein Feder-Masse-System

Mit der Variablen axis wird der Tatsache vorausschauend Rechnung getragen, dass in der Regel die Bewegungsgleichungen mehrdimensional (2D oder 3D) berechnet werden müssen. Über eine switch/case - Anweisung lassen sich so unterschiedliche Differentialgleichungen leichter implementieren.
Java - und damit auch processing - können nicht mit pointern> arbeiten, so dass die Übergabe einer Funktion (unsere Gradientenfunktion) an einen Algorithmus (den RUNGE-KUTTA-Algorithmus) nicht möglich ist. Damit nun der RK4-Algorithmus nicht für jede neue Bewegungsgleichung neu geschrieben werden muss, habe ich hier zu einem Trick gegriffen: Der RK4-Algorithmus ist nur einmal implementiert und ruft eine Funktion festen Namens auf. Damit nun innerhalb einer Aufgabe beliebige Differentialgleichungen mit dem gleichen Funktionsnamen aufgerufen werden können, werden diese mittels switch/case-Anweisungen unterschieden. Daher werden sowohl in der Funktion als auch beim RK4-Algorithmus die unterscheidende Variable axis verwendet und ausgewertet.
Der Gebrauch erfolgt so:
1. Aufruf der Instanz der Klasse Runge-Kutta 4 Solver
...
ODE.RKsO('y', y1, vy1, dt); // Runge-Kutta 4, 2. Ordnung
y1 = ODE.s;                  // Übergabe des neuen Ortes
vy1 = ODE.v;                 // Übergabe der neuen Geschwindigkeit
...

2. Deklaration der Gradientenfunktion f(). Im Beispiel gibt es nur eine 2D-Bewegung in y-Richtung, also nur eine Bewegungs-DGl. Daher wird der Wert axis = 'y' eingestellt.
...
double f(char axis, double pos, double vel)
{
    switch(axis)
      {
        case 'y':
          return(g - 2*d*vel - omega0*omega0*(pos-l0));
        default:
          return(0);
      }
}
...

3. Die Methode der Klasse Runge-Kutta 4 Solver
...
void RKsO(char axis, double pos, double vel, double dt)
{ // v Variable, delta Schrittweite
    double k1s, k2s, k3s, k4s, k1v, k2v, k3v, k4v;
    s = pos;                   //Orts- und Geschwindigkeitskoordinate
    v = vel;
    k1s = v;                   // Initialisierung
    k1v = f(axis, s, v);
    ...
}
...

Die Implementation des ODE-Solvers nach RUNGE-KUTTA ist in processing leider nicht so unproblematisch wie in p5.js, weil ein Funktionsaufruf innerhalb ein Funktionsdeklaration nicht möglich ist. Wie dennoch eine passable Lösung für processing aussehen kann, zeigt der nebenstehende Geistesblitz:

on/off

Grau ist alle Theorie. Erproben wir den RUNGE-KUTTA-Algorithmus am Beispiel des hängenden Feder-Masse-Systems. Seine Bewegungs-DGl. ist durch gegeben.

Das folgende Beispielprogramm zeigt, wie wirksam die 4-stufige Berechnungsmethode des RUNGE-KUTTA-Verfahrens ist. Dämpfung δ und Eigenfrequenz f₀ können wieder per Schieberegler gewählt werden, ebenso kann die Schrittweite mit dem Button Width verändert werden.

Ergebnisgiskussion: Schon beim ersten Blick auf die Fehlerausschriften Δω und ΔW_ges fallen die im Vergleich zu den bisher besprochenen Verfahren deutlich besseren Werte auf. Aber bei nicht vorhandenem Reibungseinfluss im Feder-Masse-System ist zu bemerken, dass - das negative Vorzeichen beim ΔW_ges zeigt es an - ein geringer Energieverlust im System in Erscheinung tritt. Der sich wegen seines geringen Ausmaßes natürlich erst nach langer Beobachtungsdauer wirklich bemerkbar macht. In der Spielewelt ist dieser Energieverlust natürlich völlig ohne Bedeutung, da unsere Systeme ja generell mit Reibungseinflüssen zu tun haben.
Wohingegen die Parametertreue uneingeschränkt sehr gut ist. Beide Fehlermaße steigen proportional mit der eingestellten Eigenfrequenz, was gleichbedeutend mit einer steiferen Feder ist. Auch der Einfluss der Schrittweite ist deutlich nachweisbar. Je kleiner die Schrittweite desto exakter die Lösung.

Symbol	Bedeutung
p	Konsistenz-Ordnung
s	Anzahl der Stufen (nicht Strecke!)
s_i	Strecke, Weg [m], Rechenergebnis im i-ten Zeitquant
v_i	Geschwindigkeit [m/s], Rechenergebnis im i-ten Zeitquant
a_j,l / c_j / k_j	Zwischenergebnisse
b_j	RUNGE-KUTTA-Koeffizienten
A	BUTCHER-Matrix
t	Zeit [s]
Δt	Zeitquant, Bildwechseldauer [s]