Ein numerisches Verfahren

Um Lösungen zu Differentialgleichungen zu finden, Trajektorien und zeitabhängige Graphen zu plotten, sind wir auf die Hilfe von leistungsstarken Computern und Software angewiesen.

Es gibt zahlreiche numerische Verfahren zur Bestimmung von Ableitungen, Integralen, Summen etc. Wir wollen uns im Folgenden einen Auszug aus der Vielfalt der Verfahren anschauen und kurz deren Genauigkeit abschätzen.

Wir betrachten eine Differentialgleichung der Form

y'=f(t,y)

mit Startwert $y_0 = y(t_0)$ und wenden darauf die folgenden Verfahren an, um $y(t)$ über einem Zeitintervall $t_0 < t < T$ numerisch abzuschätzen.

Wir beginnen mit einer simplen Methode, dem Euler-Verfahren, um die grundsätzliche Idee eines Näherungsverfahren zu erfassen und besprechen dann eine weit verbreitete Methode, das Runge-Kutta-Verfahren.

Das Euler-Verfahren

Klar ist, dass ein Computer nicht jeden Punkt einer Kurve berechnen kann, weil es ja unendlich viele davon gibt. Also beschränkt sich ein Näherungsverfahren immer auf eine diskrete Teilmenge. Euler's Methode beschreibt eine sehr einfache Annäherung an die Lösung einer Differentialgleichung für eine endliche Anzahl von Punkten. Die Teilschritte sind

Teile das betrachtete Intervall in $N$ gleich grosse Abschnitte und setze für $n = 0,1,2,\dots,N-1$

t_n = t_0+nh,

wobei $h = \frac{T-t_0}{N}$ die Schrittweite ist.

Ausgehend vom Punkt $(t_0|y_0)$ auf der Kurve approximiert man $y_n = y(t_n)$ . Eine Näherung für $y_1$ bestimmen wir via der Tangente durch $(t_0|y_0)$ bis $t_1$ . Also

y_1 \approx y_0+hf(t_0,y_0)

mit

y' = f(t_0,y_0) \approx \frac{y_1-y_0}{h} = \frac{y_1-y_0}{t_1-t_0}.

Bestimme so $y_2 \approx y_1+f(t_1,y_1)$ und weiter $y_3,\dots,y_n$ .

Allgemein beschrieben hat man das rekursive Schema

y_{n+1} = y_n+hf(t_n,y_n)

mit

t_{n+1} = t_n+h

für $0 \leq n \leq N-1$ .

Das Runge-Kutta-Verfahren

One-step Algorithmen, die durchschnittliche Steigungen einer Funktion $f(t,y)$ in zwei oder mehreren Punkten über einem Intervall $[t_{n-1},t_n]$ verwenden um $y_n$ zu berechnen, heissen Runge-Kutta-Methoden.

Eine weit verbreitete Methode ist das Runge-Kutta-Verfahren vierter Ordnung (RK4). Es verwendet gewichtete durchschnittliche Steigungen um die Mittel- und Endpunkte von Teilintervallen. Algorithmisch formuliert

y_{n+1} = y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)

mit

\begin{align*} k_1 &= f(t_n,y_n)\\ k_2 &= f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1)\\ k_3 &= f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_2)\\ k_4 &= f(t_n+h,y_n+hk_3) \end{align*}

Beweisskizze von Runge-Kutta

Theorem 1: Runge-Kutta-Verfahren (RK4)

Das Ziel ist, die Formeln für das RK4-Verfahren

y_{n+1} = y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)

mit

\begin{align*} k_1 &= f(t_n,y_n)\\ k_2 &= f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1)\\ k_3 &= f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_2)\\ k_4 &= f(t_n+h,y_n+hk_3) \end{align*}

so herzuleiten, dass der Fehler von der Ordnung $O(h^5)$ ist. Die Idee ist, die Taylor-Entwicklung der exakten Lösung mit der Taylor-Entwicklung des RK4-Schritts zu vergleichen und die Koeffizienten so zu wählen, dass sie bis zum Term $h^4$ übereinstimmen.

Proof

Die Herleitung ist algebraisch sehr intensiv. Wir skizzieren hier die wesentlichen Schritte.

Schritt 1: Die "wahre" Lösung – Taylor-Entwicklung von $y(t_{n+1})$

Die exakte Lösung $y(t)$ der Differentialgleichung $y' = f(t, y)$ lässt sich um den Punkt $t_n$ in eine Taylor-Reihe entwickeln: $y(t_n+h) = y(t_n) + h y'(t_n) + \frac{h^2}{2} y''(t_n) + \frac{h^3}{6} y'''(t_n) + \frac{h^4}{24} y''''(t_n) + O(h^5)$ Die Ableitungen von $y$ müssen wir durch die Funktion $f$ und ihre partiellen Ableitungen ausdrücken. Dazu verwenden wir die Kettenregel (wir schreiben $f$ für $f(t_n, y(t_n))$ , $f_t$ für $\frac{\partial f}{\partial t}$ , $f_y$ für $\frac{\partial f}{\partial y}$ usw.):

$y' = f$
$y'' = \frac{d}{dt}f(t, y(t)) = f_t + f_y \cdot y' = f_t + f_y f$
$y''' = f_{tt} + 2f_{ty}f + f_{yy}f^2 + f_y(f_t + f_y f)$
$y'''' = \dots$ (wird extrem kompliziert)

Setzt man dies in die Taylor-Reihe ein, erhält man die exakte Lösung, ausgedrückt durch $f$ und ihre Ableitungen. Dies ist unser "Goldstandard", den wir approximieren wollen.

Schritt 2: Die numerische Lösung – Taylor-Entwicklung von $y_{n+1}$

Nun müssen wir den RK4-Ausdruck $y_{n+1} = y_n + h \cdot \Phi$ ebenfalls als Potenzreihe in $h$ entwickeln. Der Term $\Phi = \frac{1}{6}(k_1+2k_2+2k_3+k_4)$ ist die "effektive Steigung". Dazu müssen wir jeden Term $k_i$ als mehrdimensionale Taylor-Reihe um den Punkt $(t_n, y_n)$ entwickeln.

$k_1 = f(t_n, y_n) = f$ (Das ist einfach.)
Entwicklung von $k_2$ : $k_2 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_1) = f(t_n + \frac{h}{2}, y_n + \frac{h}{2}f)$ Die Taylor-Entwicklung einer Funktion $f(t_0+\Delta t, y_0+\Delta y)$ lautet: $f(t_0+\Delta t, y_0+\Delta y) \approx f(t_0, y_0) + \Delta t \cdot f_t + \Delta y \cdot f_y + \dots$ Mit $\Delta t = h/2$ und $\Delta y = hf/2$ erhalten wir: $k_2 \approx f + \frac{h}{2}f_t + \frac{h}{2}f \cdot f_y + \frac{h^2}{8}f_{tt} + \frac{h^2}{4}f_{ty}f + \frac{h^2}{8}f_{yy}f^2 + O(h^3)$
Entwicklung von $k_3$ und $k_4$ : Dies geschieht analog. Man entwickelt $k_3$ und muss dabei die bereits entwickelte Potenzreihe für $k_2$ einsetzen. Dies führt zu sehr langen und komplizierten Ausdrücken.

Schritt 3: Koeffizientenvergleich

Nun setzen wir die entwickelten Reihen für $k_1, k_2, k_3, k_4$ in die Hauptformel ein:

y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)

Anschliessend sortiert man den gesamten Ausdruck nach Potenzen von $h$ :

y_{n+1} = y_n + h(\dots) + h^2(\dots) + h^3(\dots) + h^4(\dots) + O(h^5)

Der entscheidende Schritt ist nun der Koeffizientenvergleich: Man vergleicht die Terme in den Klammern mit den Koeffizienten der "wahren" Taylor-Reihe aus Schritt 1.

Ordnung $h^1$ : Der Koeffizient von $h$ im RK4 ist $\frac{1}{6}(f + 2f + 2f + f) = f$ . Der Koeffizient von $h$ in der wahren Lösung ist $y' = f$ . Stimmt überein!
Ordnung $h^2$ : Nach dem Einsetzen und Sortieren findet man, dass der Koeffizient von $h^2$ im RK4-Verfahren $\frac{1}{2}(f_t + f_y f)$ ist. Der Koeffizient von $h^2$ in der wahren Lösung ist $\frac{1}{2}y'' = \frac{1}{2}(f_t + f_y f)$ . Stimmt überein!
Ordnung $h^3$ und $h^4$ : Die Magie des RK4-Verfahrens liegt darin, dass die Gewichte $(\frac{1}{6}, \frac{2}{6}, \frac{2}{6}, \frac{1}{6})$ und die Stützpunkte $(0, \frac{1}{2}, \frac{1}{2}, 1)$ genau so gewählt sind, dass sich die extrem komplizierten Terme für die dritte und vierte Ordnung ebenfalls exakt aufheben und mit denen der wahren Taylor-Entwicklung übereinstimmen. Dies führt zu einem System von nichtlinearen Gleichungen für die Koeffizienten, und die "klassische" RK4-Formel ist die eleganteste und am weitesten verbreitete Lösung dieses Systems.

Da die Entwicklungen bis zum Term $h^4$ übereinstimmen, ist der erste Term, in dem sie sich unterscheiden, der Term der Ordnung $h^5$ . Der lokale Fehler eines Schrittes ist also $O(h^5)$ .

Simpson-Regel

Theorem 2: Simpson-Regel für ein Intervall

Das Integral einer Funktion $f(x)$ über das Intervall $[a, b]$ kann durch die Fläche unter einer Parabel approximiert werden, die durch die drei Punkte $(a, f(a))$ , $(\frac{a+b}{2}, f(\frac{a+b}{2}))$ und $(b, f(b))$ verläuft. Die Formel lautet:

\int_a^b f(x)dx \approx \frac{b-a}{6} \left[ f(a) + 4f\left(\frac{a+b}{2}\right) + f(b) \right]

Proof

Die Strategie des Beweises ist wie folgt:

Wir ersetzen die komplizierte Funktion $f(x)$ durch eine einfache Parabel $P(x)$ , die durch die drei gegebenen Punkte verläuft.
Wir integrieren diese Parabel exakt.
Um die Algebra massiv zu vereinfachen, verschieben wir das Koordinatensystem so, dass das Intervall symmetrisch um den Ursprung liegt.

Schritt 1: Vereinfachung des Intervalls

Anstatt über das allgemeine Intervall $[a, b]$ zu integrieren, betrachten wir das symmetrische Intervall $[-h, h]$ . Dies vereinfacht die Berechnungen erheblich, ohne die Allgemeinheit zu verlieren.

Wir setzen den Mittelpunkt des Intervalls in den Ursprung: $x=0$ .
Die Intervallgrenzen sind dann $-h$ und $h$ . Die Breite des Intervalls ist $2h$ .
Die drei Stützpunkte sind:
- $x_0 = -h$ mit Funktionswert $y_0 = f(-h)$
- $x_1 = 0$ mit Funktionswert $y_1 = f(0)$
- $x_2 = h$ mit Funktionswert $y_2 = f(h)$ (Am Ende setzen wir $h = (b-a)/2$ , um die allgemeine Formel zu erhalten.)

Schritt 2: Aufstellen der Parabelgleichung

Eine allgemeine Parabel hat die Form $P(x) = Ax^2 + Bx + C$ . Wir müssen die Koeffizienten $A, B, C$ so bestimmen, dass die Parabel durch unsere drei Punkte $(x_0, y_0), (x_1, y_1), (x_2, y_2)$ geht.

Für $(x_1, y_1) = (0, y_1)$ : $P(0) = A(0)^2 + B(0) + C = y_1 \implies \mathbf{C = y_1}$ .
Für $(x_2, y_2) = (h, y_2)$ : $P(h) = Ah^2 + Bh + C = y_2 \quad (1)$
Für $(x_0, y_0) = (-h, y_0)$ : $P(-h) = A(-h)^2 + B(-h) + C = Ah^2 - Bh + C = y_0 \quad (2)$

Jetzt lösen wir nach $A$ und $B$ . Addieren wir Gleichung (1) und (2): $(Ah^2 + Bh + C) + (Ah^2 - Bh + C) = y_2 + y_0$ $2Ah^2 + 2C = y_2 + y_0$ Setzen wir $C=y_1$ ein: $2Ah^2 + 2y_1 = y_2 + y_0$ $2Ah^2 = y_0 - 2y_1 + y_2 \implies \mathbf{A = \frac{y_0 - 2y_1 + y_2}{2h^2}}$ .

(Den Koeffizienten $B$ benötigen wir für die Integration nicht, da er sich durch die symmetrischen Grenzen aufheben wird, aber zur Vollständigkeit: Subtrahiert man (2) von (1), erhält man $2Bh = y_2 - y_0 \implies B = \frac{y_2 - y_0}{2h}$ .)

Schritt 3: Exakte Integration der Parabel

Wir integrieren nun unsere Parabel $P(x)$ über das vereinfachte Intervall von $-h$ bis $h$ :

\begin{align*} \int_{-h}^{h} P(x)dx &= \int_{-h}^{h} (Ax^2 + Bx + C)\mathrm{d}x \\ &= \left[ A\frac{x^3}{3} + B\frac{x^2}{2} + Cx \right]_{-h}^{h} \\ &= \left( A\frac{h^3}{3} + B\frac{h^2}{2} + Ch \right) - \left( A\frac{(-h)^3}{3} + B\frac{(-h)^2}{2} + C(-h) \right) \\ &= \left( \frac{Ah^3}{3} + \frac{Bh^2}{2} + Ch \right) - \left( -\frac{Ah^3}{3} + \frac{Bh^2}{2} - Ch \right) \\ &= \frac{2Ah^3}{3} + 2Ch \end{align*}

Schritt 4: Einsetzen der Koeffizienten und Vereinfachen

Jetzt setzen wir die gefundenen Ausdrücke für $A$ und $C$ in das Integrationsergebnis ein:

\begin{align*} \text{Fläche} &= \frac{2h^3}{3} \left( \frac{y_0 - 2y_1 + y_2}{2h^2} \right) + 2h(y_1) \\ &= \frac{h}{3} (y_0 - 2y_1 + y_2) + 2hy_1 \\ &= \frac{h(y_0 - 2y_1 + y_2) + 6hy_1}{3} \\ &= \frac{h}{3} (y_0 - 2y_1 + y_2 + 6y_1) \\ &= \frac{h}{3} (y_0 + 4y_1 + y_2) \end{align*}

Schritt 5: Rückkehr zum allgemeinen Intervall

Unser Ergebnis $\frac{h}{3} (y_0 + 4y_1 + y_2)$ gilt für das Intervall $[-h, h]$ der Breite $2h$ . Für das ursprüngliche Intervall $[a, b]$ gilt:

Die Breite ist $b-a$ , also $2h = b-a \implies h = \frac{b-a}{2}$ .
$y_0 = f(a)$ , $y_1 = f(\frac{a+b}{2})$ , $y_2 = f(b)$ .

Wir ersetzen $h$ in unserer Formel:

\text{Fläche} = \frac{(b-a)/2}{3} \left[ f(a) + 4f\left(\frac{a+b}{2}\right) + f(b) \right]

= \frac{b-a}{6} \left[ f(a) + 4f\left(\frac{a+b}{2}\right) + f(b) \right]