1.2. Einschrittverfahren für Anfangswertprobleme

In diesem Abschnitt wollen wir zunächst eine einfache Möglichkeit zur numerischen Berechnung von Lösungen gewöhnlicher Differentialgleichungen der folgenden Form eines Anfangswertproblems behandeln:

(1.6)\[\begin{split}\begin{split} u'(t) \ &= \ F(t,u(t)), \qquad \forall \ 0 \leq t \leq T,\\ u(0) \: &= \: u_0 \in \R^n, \end{split}\end{split}\]

Die grundlegende Idee ist es sukzessiv (d.h., aufeinander aufbauend) die Lösung zu verschiedenen Zeitschritten zu approximieren. Hierzu führen wir eine sogenannte Zeitdiskretisierung der Differentialgleichung durch indem wir das kontinuierliche Intervall \([0,T] \subset \R^+\) in \(N\in\N\) Teilintervalle mit \(N+1\) Zeitpunkten aufteilen. Der Einfachheit halber wählen wir im Folgenden uniforme Zeitschritte mit einer festgewählten Schrittweite \(\tau > 0\), d.h., wir diskretisieren das Intervall \(\Omega \coloneqq [0,T]\) durch

\[\Omega_\tau \ \coloneqq \ \lbrace t_k \in \Omega \, | \, t_k \coloneqq k \cdot \tau, \ k=0,\ldots,N \rbrace \quad \text{mit} \quad \tau \, \coloneqq \, \frac{T}{N}.\]

Es sei bemerkt, dass ein analoges Vorgehen auch für nicht äquidistante Schrittweiten möglich ist.

Basierend auf der Diskretisierung der Differentialgleichung können wir dann sukzessive numerische Approximationen \(u_\tau(t_k) \in \R^n\) für die unbekannte Lösung \(u\) des Anfangswertproblems zu diesen diskreten Zeitschritten berechnen. Hierbei haben wir die Möglichkeit zeitliche Ableitung \(\frac{d}{dt}\) durch geeignete Differenzenquotienten zu den diskreten Zeitpunkten zu approximieren. Alternativ lässt sich auch eine numerische Quadraturformel (siehe [TR]) für die zugehörige Integralgleichung in diesen Zeitschritten anwenden. Hierbei macht man sich die folgende auf dem Hauptsatz der Differential- und Integralrechnung basierende Beobachtung zu Nutze:

(1.7)\[u(t_{k+1}) \ = \ u(t_k) + \int_{t_k}^{t_{k+1}} u'(t) \,\mathrm{d}t \ = \ u(t_k) + \int_{t_k}^{t_{k+1}} F(t,u(t))\,\mathrm{d}t , \qquad t_k \in \Omega_\tau.\]

Im Fall von Einschrittverfahren, die wir in diesem Abschnitt zunächst behandeln wollen, verwenden wir dabei eine Differenzenformel, die nur einen einzigen Zeitschritt berücksichtigt, d.h. zur Berechnung von \(u_\tau(t_{k+1})\) wird lediglich der Wert \(u_\tau(t_k)\) des vorherigen Zeitschritts verwendet.

Im Folgenden wollen wir zunächst eine allgemeine Form für Einschrittverfahren zur numerischen Lösung von Anfangswertproblemen angeben.

Definition 1.4 (Einschrittverfahren und Verfahrensfunktion)

Sei eine Zeitdiskretisierung \(\Omega_\tau\) des Anfangswertproblems (1.6) für das Zeitinterval \([0,T] \subset \R^+\) mit \(N+1\) äquidistanten Zeitschritten \(t_k = k \cdot \tau\) und \(\tau = T / N\) gegeben.

Dann lassen sich Einschrittverfahren zur numerischen Lösung des Anfangswertproblems in folgender allgemeiner Form angeben:

\[u_\tau(t_{k+1}) \ = \ u_\tau(t_k) + \tau \cdot f_\tau(t, u_\tau(t)),\]

wobei \(f_\tau \colon \Omega_\tau \times \R^n \rightarrow \R^n\) eine numerische Approximation der Funktion \(F\) auf der rechten Seite der gewöhnlichen Differentialgleichung in (1.6) darstellt und allgemein Verfahrensfunktion genannt wird.

Hängt die Verfahrensfunktion \(f_\tau\) nur vom bekannten Wert \(u_\tau(t_k)\) ab, so nennen wir das Einschrittverfahren explizit, da es eine explizite Vorschrift zur Berechnung des nächsten Zeitschritts basierend auf den Werten des aktuellen Zeitschritts liefert. Hängt andererseits die Verfahrensfunktion \(f_\tau\) auch vom unbekannten Wert \(u_\tau(t_{k+1})\) ab, so heißt das Verfahren implizit, da es nur eine implizite Bedingung (im Allgemeinen eine nichtlineare Gleichung) für den neuen Zeitschritt liefert.

Ob ein Einschrittverfahren explizit oder implizit ist, hängt häufig von der Wahl der Approximation der Zeitableitung mittels finiter Differenzen ab, die in folgender Definition eingeführt werden.

Definition 1.5 (Finite Differenzen)

Die Grundidee der sogenannten finiten Differenzen ist die Approximation der Ableitung \(u'(t)\) durch Differenzenquotienten mit Werten auf einem Diskretisierungsgitter \(\Omega_\tau\) für eine feste Schrittweite \(\tau \coloneqq \frac{T}{N}, N \in \mathbb{N}\). Hierbei haben wir verschiedene Möglichkeiten die Zeitableitung zu approximieren und wir unterscheiden folgende Fälle:

\[\begin{split}\begin{aligned} u'(t_k) \ &\approx \ \frac{u(t_{k+1}) - u(t_k)}{\tau} \ = \ \frac{u(t_k + \tau) - u(t_k)}{\tau} \ &=: \ D^+ u(t_k),\\ u'(t_k) \ &\approx \ \frac{u(t_k) - u(t_{k-1})}{\tau} \ = \ \frac{u(t_k) - u(t_k - \tau)}{\tau} \ &=: \ D^- u(t_k),\\ u'(t_k) \ &\approx \ \frac{u(t_{k+1}) - u(t_{k-1})}{2\tau} \ = \ \frac{u(t_k + \tau) - u(t_k - \tau)}{2\tau} \hspace{-1cm} &=: \ D^c u(t_k). \end{aligned}\end{split}\]

Hierbei werden die obigen Differenzenquotienten \(D^+, D^-, D^c\) Vorwärts-, Rückwärts- und zentrale Differenz genannt.

Mittels Taylorformel sieht man ein, dass alle drei finiten Differenzen in Definition 1.5 eine gute numerische Approximation für die erste Ableitung der Funktion \(u\) liefern, wenn die Schrittweite \(\tau > 0\) der Diskretisierung hinreichend klein gewählt wird. Da der Fehler dieser Approximation für \(\tau \rightarrow 0\) bei allen drei Varianten gegen Null konvergiert spricht man auch von konsistenten Verfahren (mehr dazu später). Bei einer quantitativen Analyse des Fehlerglieds in der Taylorentwicklung lässt sich zeigen, dass die Vorwärts- und Rückwärtsdifferenz \(D^+\) und \(D^-\) jeweils die Konsistenzordnung \(1\) besitzen, d.h., der Fehler verschwindet linear bzgl. der Schrittweite \(\tau\), während die zentrale Differenz \(D^c\) sogar Konsistenzordnung \(2\) besitzt, d.h., der Fehler verschwindet quadratisch bzgl. der Schrittweite \(\tau\).

Wir betrachten nun einige bekannte Einschrittverfahren, die häufig zur numerischen Lösung von Anfangswertproblemen genutzt werden.

Beispiel 1.3 (Einschrittverfahren)

Im Folgenden wollen wir drei unterschiedliche Einschrittverfahren diskutieren, die häufig zur numerischen Lösung von Anfangswertproblemen genutzt werden.

1. Das einfachste Beispiel eines Einschrittverfahrens ist das Vorwärts-Euler-Verfahren

   \[u_\tau(t_{k+1}) \ = \ u_\tau(t_k) + \tau \cdot F(t_k,u_\tau(t_k)),\]

   welches auf einer Approximation der Ableitung \(u'(t_k)\) mittels Vorwärtsdifferenz aus Definition 1.5 basiert. Das Vorwärts-Euler Verfahren hat die Verfahrensfunktion

   \[f_\tau(t, u_\tau(t)) \ \coloneqq \ F(t_k,u_\tau(t_k)),\]

   und ist entsprechend nach Definition 1.4 ein explizites Verfahren. Daher wird es häufig auch explizites Euler-Verfahren genannt.

   In der Integralform der Differentialgleichung bedeutet das, dass wir die folgende Approximation benutzen:

   \[\begin{aligned} u(t_{k+1}) - u(t_k) \ = \ \int_{t_k}^{t_{k+1}} F(t,u(t))\, \mathrm{d}t \ \approx \ \tau \cdot F(t_k, u_\tau(t_k)), \end{aligned}\]

   d.h., dass wir das Integral durch die Intervallänge mal dem Wert am linken Intervallrand annähern. Dies entspricht einer numerischen Quadraturformel mit nur einer Stützstelle (siehe [TR]).

2. Ein wenig komplizierter ist das Rückwärts-Euler-Verfahren

   \[u_\tau(t_{k+1}) \ = \ u_\tau(t_k) + \tau \cdot F(t_{k+1},u_\tau(t_{k+1})),\]

   welches auf einer Approximation der Ableitung \(u'(t_k)\) mittels Rückwärtsdifferenz aus Definition 1.5 basiert. In diesem Fall müssen wir zunächst eine (möglicherweise nichtlineare) Gleichung für den neuen Zeitschritt lösen.

   Das Rückwärts-Euler Verfahren hat die Verfahrensfunktion

   \[f_\tau(t, u_\tau(t)) \ \coloneqq \ F(t_{k+1},u_\tau(t_{k+1})),\]

   und ist entsprechend nach Definition 1.4 ein implizites Verfahren. Daher wird es häufig auch implizites Euler-Verfahren genannt.

   In der Integralform der Differentialgleichung bedeutet das, dass wir die folgende Approximation benutzen:

   \[\begin{aligned} u(t_{k+1}) - u(t_k) \ = \ \int_{t_k}^{t_{k+1}} F(t,u(t)) \,\mathrm{d}t \ \approx \ \tau \cdot F(t_{k+1},u_\tau(t_{k+1})), \end{aligned}\]

   d.h., dass wir das Integral durch die Intervallänge mal dem Wert am rechten Intervallrand annähern. Auch dies entspricht einer numerischen Quadraturformel mit nur einer Stützstelle wie beim expliziten Euler-Verfahren.

3. Ein Kompromiss aus den beiden obigen Verfahren ist das sogenannte Crank-Nicolson Verfahren der Form

   \[u_\tau(t_{k+1}) \ = \ u_\tau(t_k) + \frac{\tau}{2} \cdot \left(F(t_{k },u_\tau(t_{k })) + F(t_{k+1},u_\tau(t_{k+1}))\right).\]

   Hierbei ist die Verfahrensfunktion der Mittelwert der Verfahrensfunktionen des Vorwärts- und Rückwärts-Euler-Verfahrens mit

   \[\begin{aligned} f_\tau(t, u_\tau(t)) \ \coloneqq \ \frac{1}{2} \left( F(t_{k },u_\tau(t_{k })) + F(t_{k+1},u_\tau(t_{k+1}))\right) \end{aligned}\]

   und ist entsprechend nach Definition 1.4 ebenfalls ein implizites Verfahren.

   In der Integralform der Differentialgleichung bedeutet das, dass wir die folgende Approximation benutzen:

   \[\begin{split}\begin{aligned} u(t_{k+1}) - u(t_k) \ &= \ \int_{t_k}^{t_{k+1}} F(t,u(t)) \,\mathrm{d}t \\ &\approx \ \frac{\tau}{2} \left( F(t_{k },u_\tau(t_{k })) + F(t_{k+1},u_\tau(t_{k+1}))\right), \end{aligned}\end{split}\]

   Dies entspricht einer Approximation des Integrals mit Hilfe der Trapezregel aus [TR].

Wir sehen, dass ein explizites Verfahren sofort wohldefiniert ist, falls die Verfahrensfunktion \(f_\tau\) eine stetige Funktion auf \(\mathbb{R}^+ \times \mathbb{R}^n\) ist, während bei impliziten Verfahren noch eine Fixpunktgleichung gelöst werden muss. Die Existenz und Eindeutigkeit von Lösungen dieser Gleichung können wir mit dem Banachschen Fixpunktsatz garantieren, wenn wiederum \(\tau\) klein genug ist.

Satz 1.2 (Existenz und Eindeutigkeit von Lösungen)

Sei \(f_\tau:\R\times\R^n\to\R^n\) stetig und Lipschitz-stetig bezüglich dem zweiten Argument mit Lipschitzkonstante \(L > 0\). Dann existiert für \(\tau< \frac{1}{L}\) genau eine Lösung \(u_\tau(t_{k+1})\) der Fixpunktgleichung

\[u(t) = u_\tau(t_k) + \tau\ f_\tau(t,u(t)).\]

Numerisch müssen wir zur Durchführung eines impliziten Verfahrens immer noch ein Gleichungssystem in \(\mathbb{R}^n\) lösen. Ist dieses linear, so können wir die üblichen Verfahren für lineare Gleichungssysteme anwenden. Andernfalls bietet sich die Verwendung eines iterativen Verfahrens wie einer Fixpunktiteration oder des Newton-Verfahrens an (beachte, dass unter der obigen Bedingung \(I_n-\tau f_\tau\) invertierbar ist für \(f_\tau \in C^1\)). Mit dem Wert \(u_\tau(t_k)\) oder einer einfachen Vorhersage in der Zeit (etwa mit dem expliziten Euler-Verfahren) haben wir dafür auch einen sehr guten Startwert.

Nachdem wir die Wohldefiniertheit und numerische Umsetzung von Einschrittverfahren geklärt haben, widmen wir uns nun der Analyse dieser Verfahren. Wir wollen dabei den Fehler

(1.8)\[E_\tau = \max_{t_k\in\Omega_\tau} \Vert u_\tau(t_k) - u(t_k) \Vert\]

abschätzen, wobei \(u\) die exakte Lösung des Anfangswertproblems ist. Wollen wir den Fehler an anderen Stellen \(t \in [0, T]\) abschätzen, so können wir ein Interpolationsverfahren und die entsprechenden Abschätzungen anwenden.

Wir betrachten zunächst die Fehlerdifferenz \(e_\tau \coloneqq u_\tau - u\). Auf Grund es Hauptsatzes der Differential- und Integralrechnung können wir schreiben

\[u(t_{k+1}) \ = \ u(t_k) + \int_{t_k}^{t_{k+1}} F(t, u(t))\,\mathrm{d}t\]

und somit können wir die Fehlerdifferenz entwickeln als:

\[\begin{split}\begin{aligned} e_\tau(t_{k+1}) \ &= \ u_\tau(t_{k+1}) - u(t_{k+1}) \ = \ u_\tau(t_k) + \tau f_\tau(t, u_\tau(t)) - u(t_{k+1})\\ &= \ \underbrace{u_\tau(t_{k}) - u(t_{k})}_{= \, e_\tau(t_k)} \, + \, \tau f_\tau(t, u_\tau(t)) - \int_{t_k}^{t_{k+1}} F(t, u(t)) \, \mathrm{d}t \\ &= \ e_\tau(t_k) + \tau \left( f_\tau(t, u_\tau(t)) - f_\tau(t,u(t)) \right) \\ & \qquad \qquad + \tau \left( f_\tau(t, u(t)) - \frac{1}{\tau}\int_{t_k}^{t_{k+1}} F(t, u(t)) \, \mathrm{d}t \right) \end{aligned}\end{split}\]

Nun benötigen wir zwei zentrale Eigenschaften von Diskretisierungsmethoden:

• Konsistenz: Der Fehler

  \[\begin{aligned} f_\tau(t,u(t)) - \frac{1}{\tau}\int_{t_k}^{t_{k+1}} F(t,u(t))~dt, \end{aligned}\]

  d.h. das Residuum der Lösung des Anfangswertproblems eingesetzt in das numerische Verfahren konvergiert gegen Null für \(\tau \rightarrow 0\).

• Stabilität: Bei der Umsetzung des numerischen Verfahrens mit gegebener rechter Seite werden Störungen in den Anfangsdaten nicht beliebig verstärkt. Insbesondere existiert eine Abschätzung des Fehlers unabhängig von der Zeitschrittweite \(\tau\).

Zusammen ergeben Konsistenz und Stabilität Konvergenz des Verfahrens, d.h. \(E_\tau \rightarrow 0\). Dies halten wir in einer Definition fest:

Definition 1.6 (Konvergenz von Einschrittverfahren)

Sei \(E_\tau\) definiert durch (1.8), dann heißt das Verfahren

1. konvergent, wenn \(E_\tau \rightarrow 0\) für \(\tau \rightarrow 0\),

2. konvergent von der Ordnung \(p\), wenn \(E_\tau = {\cal O}(\tau^p)\) für \(\tau \rightarrow 0\), d.h. es gibt eine Konstante \(C_p\), sodass \(E_\tau \leq C_p \tau^p\) für \(\tau\) hinreichend klein.

1.2.1. Konsistenz von Einschrittverfahren

Gemäß der obigen Motivation definieren wir den globalen Konsistenzfehler als

(1.9)\[K_\tau = \max_{t_k\in\Omega_\tau} \Vert g_\tau(t_k) \Vert,\]

wobei der lokale Konsistenzfehler gegeben ist durch

\[g_\tau(t_k) = f_\tau(t,u(t)) - \frac{1}\tau \int_{t_k}^{t_{k+1}} F(t,u(t))\,\mathrm{d}t.\]

Hierbei stellt \(u\) eine Lösung des Anfangswertproblems (1.1) dar.

Definition 1.7 (Konsistenzfehler)

Sei \(K_\tau\) definiert durch (1.9), dann heißt das Verfahren

1. konsistent, wenn \(K_\tau \rightarrow 0\) für \(\tau \rightarrow 0\),

2. konsistent von der Ordnung \(p\), wenn \(K_\tau = {\cal O}(\tau^p)\) für \(\tau \rightarrow 0\) .

Wir führen die Abschätzung des Konsistenzfehlers im Folgenden an zwei Beispielen durch:

Beispiel 1.4 (Konsistenzfehler des Vorwärts-Euler Verfahrens)

Wir betrachten das Vorwärts-Euler Verfahren unter der Annahme, dass \(F\) bezüglich beider Variablen Lipschitz-stetig ist. Definieren wir \(\varphi(t) = F(t,u(t))\), dann ist \(\varphi\) wegen \(u \in C^1([0,T])\) eine Lipschitz-stetige Funktion auf dem Intervall \([0,T]\) und es gilt

\[\begin{split}\begin{aligned} K_\tau \ = \ \norm{g_\tau(t_k)}_\infty \ &= \ \norm{\varphi(t_k) - \frac{1}\tau \int_{t_k}^{t_{k+1}} \varphi(t)\,\mathrm{d}t}_\infty \ = \ \norm{\frac{1}\tau \int_{t_k}^{t_{k+1}} \varphi(t_k)-\varphi(t)\,\mathrm{d}t}_\infty\\ &\leq\frac{1}\tau \int_{t_k}^{t_{k+1}} \norm{\varphi(t_k)-\varphi(t)}_\infty\,\mathrm{d}t \ \leq \ \frac{1}\tau \int_{t_k}^{t_{k+1}} L \cdot ( t-t_k)\, \mathrm{d}t \\ &= \ \frac{L}{\tau} \cdot \left[\frac{t^2}{2}- t\cdot t_k\right]^{t_{k+1}}_{t_k} \ = \ \frac{L}\tau \cdot (\frac{t_{k+1}^2}{2} - t_{k+1}\cdot t_k - \frac{t_k^2}{2} + t_k^2) \\ &= \ \frac{L}{2\tau}(\underbrace{t_{k+1} - t_k}_{= \, \tau})^2 \ = \ \frac{L}2 \tau. \end{aligned}\end{split}\]

Damit hat das Verfahren also die Konsistenzordnung \(p=1\). Wir sehen im konkreten Beispiel \(F(t,u(t)) = t\) auch sofort, dass man im Allgemeinen nicht Konsistenzordnung zwei erreichen kann.

Häufig verwendet man zur Analyse des Konsistenzfehlers eine Taylorapproximation. Dies wollen wir im folgenden Beispiel demonstrieren.

Beispiel 1.5 (Konsistenzfehler des Crank-Nicholson Verfahrens)

Wir betrachten das Crank-Nicholson Verfahren unter der Annahme, dass \(F\) bezüglich beider Variablen zweimal stetig differenzierbar ist. Dann ist \(u\) ebenfalls zweimal stetig differenzierbar, da

\[u''(t) = \frac{\mathrm{d}}{\mathrm{d}t}F(t,u(t)) = \partial_t F(t,u(t)) + \partial_u F(t,u(t)) \cdot u'(t)\]

stetig ist.

Wir definieren wieder die Hilfsfunktion \(\varphi(t) \coloneqq F(t,u(t))\). Mit Hilfe einer Taylorentwickung zweiter Ordnung von \(\varphi\) in der Stelle \(t_k\) können wir schreiben

\[\begin{split}\begin{aligned} K_\tau \ = \ \norm{g_\tau(t_k)} \ &= \ \norm{\frac{1}2(\varphi(t_k)+\varphi(t_{k+1})) - \frac{1}\tau \int_{t_k}^{t_{k+1}} \varphi(t)\,\mathrm{d}t}_\infty\\ &= \ \Vert \frac{1}{2\tau} \int_{t_k}^{t_{k+1}} \varphi(t_k)+\varphi(t_{k+1}) \, \mathrm{d}t - \frac{1}{2\tau} \int_{t_k}^{t_{k+1}} 2\varphi(t)\,\mathrm{d}t \Vert_\infty \\ &= \ \frac{1}{2\tau} \Vert \int_{t_k}^{t_{k+1}} \varphi(t_k)+\varphi(t_{k+1})-2\varphi(t)\,\mathrm{d}t \Vert_\infty \\ &= \ \frac{1}{2\tau} \Vert \int_{t_k}^{t_{k+1}} \varphi(t_k) + \varphi(t_k) + \varphi'(t_k)(t_{k+1} - t_k) + \frac{\varphi''(\xi_1)}{2}(t_{k+1} - t_k)^2 \\ &\qquad \qquad \qquad - 2 \left[\varphi(t_k) + \varphi'(t_k)(t-t_k) + \frac{\varphi''(\xi_2)}{2}(t-t_k)^2\right] \,\mathrm{d}t \Vert_\infty\\ &= \ \frac{1}{2\tau} \Vert\int_{t_k}^{t_{k+1}} \varphi'(t_k)(t_k+t_{k+1}-2t) + \frac{\varphi''(\xi_1)}{2}\tau^2 - \varphi''(\xi_2)(t - t_k)^2 \,\mathrm{d}t\Vert_\infty\\ &\leq \ \frac{1}{2\tau} \Vert\int_{t_k}^{t_{k+1}} \varphi'(t_k)(t_k+t_{k+1}-2t) \, \mathrm{d}t \Vert_\infty + \frac{1}{2\tau} \Vert \int_{t_k}^{t_{k+1}} \varphi''({\xi})\tau^2 \, \mathrm{d}t \Vert_\infty\\ &= \ \frac{1}{2\tau} \Vert\int_{t_k}^{t_{k+1}} \varphi'(t_k)(t_k+t_{k+1}-2t) \, \mathrm{d}t \Vert_\infty + \frac{\tau^2}{2}\Vert \varphi''({\xi}) \Vert_\infty . \end{aligned}\end{split}\]

Hierbei haben wir \(\xi \coloneqq \operatorname{argmax}_{t\in \lbrace \xi_1, \xi_2\rbrace} ||\varphi''(t)||_\infty\) gesetzt. Darüberhinaus verschwindet der erste Term, da gilt

\[\begin{split}\begin{split} \int_{t_k}^{t_{k+1}} \varphi'(t_k)(t_k+t_{k+1}-2t)\,\mathrm{d}t \ &= \ \varphi'(t_k) \cdot (t_{k+1} + t_k) \cdot (t_{k+1} - t_k) - \varphi'(t_k) \cdot [t^2]^{t_{k+1}}_{t_k} \\ &= \ \varphi'(t_k) \cdot (t^2_{k+1} - t^2_k) - \varphi'(t_k) \cdot (t^2_{k+1} - t^2_k) \ = \ 0. \end{split}\end{split}\]

Somit gilt also für den Konsistenzfehler \(K_\tau \leq C \tau^2\) und damit hat das Crank-Nicholson-Verfahren die Konsistenzordnung \(p=2\).

Wir sehen am konkreten Beispiel \(F(t,u(t)) = t^2\) auch sofort, dass man im Allgemeinen nicht Konsistenzordnung drei erreichen kann.

Um eine Konsistenzordnung \(p \in \mathbb{N}\) zu erhalten, müssen wir fordern, dass die rechte Seite \(F\), ebenso wie die unbekannte Lösung \(u\) der Differentialgleichung \(p\)-mal stetig differenzierbar sind. Aus den Eigenschaften von \(F\) folgt Letzteres aber sofort: wir haben gesehen, dass für \(F\) stetig auch \(u\) stetig differenzierbar folgt. Mit dem Argument aus dem letzten Beispiel sehen wir, dass für \(F\) stetig differenzierbar auch \(u\) zweimal stetig differenzierbar ist. Induktiv können wir durch weiteres differenzieren zeigen, dass \(u\) \(p\)-mal stetig differenzierbar ist, wenn \(F\) \(p-1\)-mal stetig differenzierbar ist.

1.2.2. Stabilität und Konvergenz

Wir widmen uns nun der Frage der Stabilität von Einschrittverfahren. Hierbei verwenden wir eine diskrete Version des Lemmas von Gronwall.

Lemma 1.2 (Diskretes Lemma von Gronwall)

Es sei \(\alpha >0\) und es sei \((\beta_j)_{j \in \N} \subset \R^+_0\) eine Folge nicht-negativer Zahlen. Für die Folge \((u_j)_{j\in\N}\subset\R\) gelte

\[\begin{aligned} u_0 \ \leq \ \alpha\in\R^+, \qquad u_k \ \leq \ \alpha + \sum_{j=0}^{k-1} \beta_j u_j \end{aligned}\]

für \(k\in\N\), dann gilt die Abschätzung

\[\begin{aligned} u_k \ \leq \ \alpha \exp\left(\sum_{j=0}^{k-1} \beta_j\right). \end{aligned}\]

Proof. In der Hausaufgabe zu zeigen. ◻

Lemma 1.3 (Abschätzung für Lipschitz-stetige Verfahrensfunktion)

Sei \(f_\tau \colon \R^+ \times \R^+ \times \R^n \times \R^n \rightarrow \R^n\) die Verfahrensfunktion eines Einschrittverfahrens zum Lösen des Anfangswertproblems (1.6) für eine gegebene Zeitdiskretisierung \(\Omega_\tau\) von \([0,T] \subset \R^+\) mit gewählter Zeitschrittweite \(\tau > 0\). Die Verfahrensfunktion \(f_\tau\) sei Lipschitz-stetig im dritten und vierten Argument, d.h. bezüglich \(u_\tau(t_k)\) und \(u_\tau(t_{k+1})\) mit jeweiligen Lipschitz-Konstanten \(L_1, L_2 > 0\) unabhängig von der gewählten Zeitschrittweite \(\tau\). Dann gilt für alle diskreten Zeitpunkte \(t_k, t_{k+1} \in \Omega_\tau\) für eine Konstante \(L > 0\) die folgende Abschätzung:

\[\begin{split}\begin{split} ||f_\tau(t_k, t_{k+1}, u_\tau(t_k), u_\tau(t_{k+1})) &- f_\tau(t_k, t_{k+1}, u(t_k), u(t_{k+1}))|| \\ &\leq \ L \cdot (||u_\tau(t_k) - u(t_k)|| + ||u_\tau(t_{k+1}) - u(t_{k+1}) ||). \end{split}\end{split}\]

Proof. Da die Verfahrensfunktion \(f_\tau\) Lipschitz-stetig bezüglich des dritten und vierten Arguments ist, gilt

\[\begin{split}\begin{split} ||f_\tau(t_k, t_{k+1}, u_\tau(t_k), u_\tau(t_{k+1})) - f_\tau(t_k, t_{k+1}, u(t_k), u_\tau(t_{k+1}))|| \ &\leq \ L_1 \cdot ||u_\tau(t_k) - u(t_k)||,\\ ||f_\tau(t_k, t_{k+1}, u_\tau(t_k), u_\tau(t_{k+1})) - f_\tau(t_k, t_{k+1}, u_\tau(t_k), u(t_{k+1}))|| \ &\leq \ L_2 \cdot ||u_\tau(t_{k+1}) - u(t_{k+1})|| \end{split}\end{split}\]

Verwenden wir nun die Dreiecksungleichung, so können wir durch Einfügen eines passenden Nullterms wie folgt abschätzen:

\[\begin{split}\begin{split} &\phantom{=} \ \ ||f_\tau(t_k, t_{k+1}, u_\tau(t_k), u_\tau(t_{k+1})) - f_\tau(t_k, t_{k+1}, u(t_k), u(t_{k+1}))|| \\ &= \ || f_\tau(t_k, t_{k+1}, u_\tau(t_k), u_\tau(t_{k+1})) - f_\tau(t_k, t_{k+1}, u(t_k), u(t_{k+1})) \\ & \quad + \underbrace{f_\tau(t_k, t_{k+1}, u_\tau(t_k), u(t_{k+1})) - f_\tau(t_k, t_{k+1}, u_\tau(t_k), u(t_{k+1})) }_{=0}|| \\ &\leq \ || f_\tau(t_k, t_{k+1}, u_\tau(t_k), u_\tau(t_{k+1})) - f_\tau(t_k, t_{k+1}, u_\tau(t_k), u(t_{k+1}))|| \\ & \quad + || f_\tau(t_k, t_{k+1}, u(t_k), u(t_{k+1})) - f_\tau(t_k, t_{k+1}, u_\tau(t_k), u(t_{k+1}))|| \\ &= \ L_1 \cdot || u_\tau(t_k) - u(t_k) || + L_2 \cdot || u_\tau(t_{k+1}) - u(t_{k+1}) || \\ &\leq \ \underbrace{\max(L_1, L_2)}_{=: \, L} \cdot (|| u_\tau(t_k) - u(t_k) || + || u_\tau(t_{k+1}) - u(t_{k+1}) ||) \end{split}\end{split}\]

◻

Mit Hilfe der Abschätzung aus Lemma 1.3 können wir nun zeigen, dass der Fehler des Einschrittverfahrens im Wesentlich durch den Konsistenzfehler bestimmt wird.

Satz 1.3 (Fehlerabschätzung für Einschrittverfahren)

Sei \(u_\tau \colon \Omega_\tau \rightarrow \R^n\) die Lösung eines Einschrittverfahrens mit Lipschitz-stetiger Verfahrensfunktion \(f_\tau\) für eine Zeitdiskretisierung des Zeitintervalls \([0,T]\) gegeben durch \(\Omega_\tau \coloneqq \lbrace t_k \in [0, T] \colon t_k \coloneqq k\cdot \tau, k=0,\ldots,N\rbrace\) mit äquidistanter Zeitschrittweite \(\tau \coloneqq \frac{T}{N} > 0\). Sei außerdem \(u \colon [0, T] \rightarrow \R^n\) die Lösung des Anfangswertproblems (1.1) mit Anfangswert \(u_0 \in \R^n\). Der lokale Konsistenzfehler \(g_\tau(t_k)\) und der globale Konsistenzfehler \(K_\tau\) seien gegeben wie in Definition 1.7.

Dann existiert eine Konstante \(C > 0\), so dass für \(\tau > 0\) hinreichend klein die folgende Abschätzung für den Fehler des Einschrittverfahrens gilt:

\[E_\tau \ = \ \max_{t_k \in \Omega_\tau} \Vert u_\tau(t_k) - u(t_k) \Vert \ \leq \ C \cdot \max_{t_k \in \Omega_\tau} \Vert g_\tau(t_k) \Vert \ = \ C \cdot K_\tau.\]

Proof. Wir definieren uns zunächst die Hilfsfunktion \(v_k \coloneqq \Vert u_\tau(t_{k}) -u(t_k)) \Vert\). Durch Anwendung der Dreiecksungleichung und dem Ausnutzen der Lipschitz-Stetigkeit von \(f_\tau\) im dritten und vierten Argument erhalten wir mit Hilfe von Lemma 1.3 dann

\[\begin{split}\begin{aligned} v_{k+1} \ &= \ \Vert u_\tau(t_{k+1}) - u(t_{k+1})\Vert\\ &= \ \left\Vert u_\tau(t_k) + \tau \cdot f_\tau(t_k, t_{k+1}, u_\tau(t_k), u_\tau(t_{k+1})) - u(t_k) -\int_{t_k}^{t_{k+1}} F(t,u(t)) \, \mathrm{d}t \right\Vert\\ &\leq \ v_k + \tau \cdot \norm{f_\tau(t_k, t_{k+1}, u_\tau(t_k), u_\tau(t_{k+1})) - f_\tau(t_k, t_{k+1}, u(t_k), u(t_{k+1}))}\\ & \hspace{1.01cm} + \tau \cdot \left\Vert f_\tau(t_k, t_{k+1}, u(t_k), u(t_{k+1})) - \frac{1}{\tau} \int_{t_k}^{t_{k+1}} F(t,u(t)) \, \mathrm{d}t \right\Vert\\ &\leq \ v_k + \tau \cdot L \cdot (v_k + v_{k+1}) + \tau \cdot \norm{g_\tau(t_k)}. \end{aligned}\end{split}\]

Somit erhalten wir durch Umstellen also insgesamt

\[\begin{split}\begin{aligned} v_{k+1} \ &\leq \ v_k + \tau \cdot L \cdot (v_k + v_{k+1}) + \tau \cdot \Vert g_\tau(t_k) \Vert\\ \Rightarrow \quad (1 - \tau L) \cdot v_{k+1} \ &\leq \ (1+\tau L) \cdot v_k + \tau \cdot \max_{t_k \in \Omega_\tau}\norm{g_\tau(t_k)} \end{aligned}\end{split}\]

Für hinreichend kleine Schrittweiten \(\tau \leq \frac{1}{2L}\) erhalten wir

\[v_{k+1} \ \leq \ \frac{1+\tau L}{1-\tau L} \cdot v_k + \frac{\tau}{1-\tau L} \cdot \max_{t_k \in \Omega_\tau}\norm{g_\tau(t_k)} \ \leq \! \underbrace{3}_{=: \, \beta} \!\!\! \cdot \, v_k + \underbrace{\frac{1}{L} \cdot \max_{t_k \in \Omega_\tau}\norm{g_\tau(t_k)}}_{=: \, \alpha}.\]

Mittels des diskreten Lemma 1.2 von Gronwall erhalten wir die gewünschte Schranke als

\[v_{k+1} \ \leq \ \alpha \cdot \exp(\beta) \ = \ \frac{e^3}{L} \cdot \max_{t_k \in \Omega_\tau}\norm{g_\tau(t_k)} \ =: \ C \cdot K_\tau.\]

◻

Verwendet man im Beweis von Satz 1.3 die Hilfsfunktion \(v_{k+1} \coloneqq ||u_\tau(t_{k+1}) - u_0||\), so kann man analog uniforme Schranken an das Wachstum der numerischen Lösung \(u_\tau\) beweisen, denn es gilt:

\[||u_\tau(t_{k+1})|| - ||u_0|| \ \leq \ ||u_\tau(t_{k+1}) - u_0|| \ = \ v_{k+1} \ \leq \ C \cdot K_\tau.\]

Da die Abschätzung unabhängig von der Stelle \(t_k \in \Omega_\tau\) gilt, erhalten wir

\[\max_{t_k \in \Omega_\tau} \Vert u_\tau(t_k) \Vert \ \leq \ C \cdot K_\tau + ||u_0|| \ =: \ M.\]

Eine direkte Folgerung der Fehlerabschätzung von Satz 1.3 ist die Äquivalenz von Konsistenz und Konvergenz für allgemeine Einschrittverfahren, wie wir im folgenden Korollar festhalten.

Korollar 1.1 (Äquivalenz von Konsistenz und Konvergenz)

Für ein Einschrittverfahren mit Lipschitz-stetiger Verfahrensfunktion gilt: ist das Verfahren konsistent von der Ordnung \(p\), so ist es auch konvergent von der Ordnung \(p\). Insbesondere folgt also aus der Konsistenz eines Einschrittverfahrens bereits die Konvergenz, da für ein konsistentes Einschrittverfahren gilt

\[E_\tau \ = \ \max_{t_k \in \Omega_\tau} \Vert u_\tau(t_k) - u(t_k) \Vert \ \leq \ C \cdot K_\tau \ \overset{\tau \rightarrow 0}{\longrightarrow} \ 0.\]

Mit Hilfe dieses Korollars und der Abschätzung des Konsistenzfehlers stellen wir fest, dass das Vorwärts- und Rückwärts-Euler Verfahren konvergent von der Ordnung eins sind. Das Crank–Nicholson Verfahren hingegen ist auf Grund der Berechnungen zur Konsistenzordnung in Beispiel 1.5 sogar konvergent von der Ordnung zwei.

1.2.3. Runge–Kutta Verfahren

Bisher haben wir Verfahren kennengelernt, die auf einzelne Funktionsauswertungen von \(F\) an den Zeitschritten \(t_k\) und \(t_{k+1}\) zurückgreifen. Damit haben wir bereits die Konsistenzordnung Eins erreicht für die beiden Euler-Verfahren und als maximale Konsistenzordnung Zwei im Falle des Crank–Nicholson Verfahrens. Eine höhere Konsistenzordnung ist mit so einem Ansatz nicht möglich. Wir widmen uns im Folgenden also der Frage wie wir Einschrittverfahren höherer Ordnung konstruieren können. Aus der Analyse der Konsistenzordnung lässt sich ableiten, dass wir Verfahrensfunktionen der Art konstruieren müssen, so dass die Taylorentwicklung ein Restglied höherer Ordnung liefert.

Eine erste Idee zur Steigerung der Konsistenzordnung ist es Ableitungen von \(F\) in den Zeitschritten \(t_k\) und \(t_{k+1}\) zu berücksichtigen, womit man offensichtlich die Taylor-Entwicklung besser approximieren und somit eine höhere Ordnung erreichen kann. Die Berechnung von Ableitungen von \(F\) ist jedoch potentiell numerisch aufwändig und instabil, deswegen geht alternativ einen anderen Weg und verwendet geschachtelte Funktionsauswertungen.

Die grundlegende Idee ist es das Integral der Anfangswertaufgabe in (1.7) durch eine geeignete Quadraturformel (siehe [TR]) der folgenden Form zu approximieren:

\[\frac{1}{\tau} \int_{t_k}^{t_{k+1}} F(t, u(t)) \,\mathrm{d}t \ \approx \ f_\tau(t_k, u_\tau(t_k)) \ \coloneqq \ \sum_{i=1}^s b_i f_i^k.\]

Hierbei stellen die \(b_i \in \R\) Gewichte der Quadraturformel dar und die Zwischenwerte \(f_i^k \in \R^n\) sollen eine Approximation der Funktionswerte von \(F\) liefern für \(i=1,\ldots,s\) mit

\[f_i^k \ \approx \ F(t_k+c_i\tau,u(t_k + c_i\tau)).\]

Hierbei sind die \(0 \leq c_i \leq 1, i=1,\ldots,s\) aufsteigende Parameter, die Stützstellen für die jeweiligen Funktionsauswertung definieren. Da wir die Werte der unbekannten Funktion \(u(t_k+c_i \tau)\) an den Stützstellen nicht kennen, benötigen wir ebenfalls Approximationen mit Hilfe von Quadraturformeln der folgenden Form:

\[u(t_k + c_i\tau) \ = \ u(t_k) + \int_{t_k}^{t_k + c_i\tau} F(t, u(t)) \, \mathrm{d}t \ \approx \ u_\tau(t_k) + \tau \sum_{j=1}^s a_{ij} f_j^k.\]

Diese Idee führt zur Definition von sogenannten Runge-Kutta Verfahren.

Definition 1.8 (Runge–Kutta Verfahren der Stufe s)

Bei einem Runge–Kutta Verfahren der Stufe \(\mathbf{s}\), \(s \in \N^+\) berechnet man

\[u_\tau(t_{k+1}) \ = \ u_\tau(t_k) + \tau f_\tau(t_k, u_\tau(t_k)) \ = \ u_\tau(t_k) + \tau \sum_{i=1}^s b_i f_i^k.\]

Die Zwischenwerte \(f_i^k\) lassen sich als Funktionsauswertungen von \(F\) wie folgt berechnen

(1.10)\[f_i^k \ \coloneqq \ F\Bigl(t_k + c_i \tau, u_\tau(t_k) + \tau \sum_{j=1}^s a_{ij} f_j^k\Bigr), \qquad i=1,\ldots,s.\]

Um in der Zeit vorwärts zu gehen wählt man die Koeffizienten \(0 \leq c_i \leq 1, i=1,\ldots,s\) als aufsteigende Folge.

Man nennt ein Runge–Kutta Verfahren explizit falls für alle Einträge der Matrix \(A\in \R^{s\times s}\) gilt \(a_{ij}=0, i=1,\ldots,s\) wenn \(j \geq i\), d.h., wenn die Matrix \((a_{ij})_{i,j=1,\ldots,s}\) eine linke, untere Dreiecksmatrix mit Nullen auf der Hauptdiagonale ist.

Die zu bestimmenden Koeffizienten \(b_i,c_i \in \R\) und \(a_{i,j} \in \R\) für \(i,j=1,\ldots,s\) in Definition 1.8 erhält man durch einen Vergleich mit der Taylorentwicklung des zu approximierenden Integrals \(\frac{1}\tau \int_{t_k}^{t_{k+1}} F(t,u(t)) \, \mathrm{d}t\). Das folgende Beispiel soll die Idee zur Herleitung eines Runge–Kutta Verfahren für den einfachen Fall der Stufe \(s=1\) verdeutlichen.

Beispiel 1.6 (Runge–Kutta Verfahren der Stufe 1)

Für ein Runge–Kutta Verfahren der Stufe \(s=1\) ist die Verfahrensfunktion gegeben durch

\[f_\tau(t_k, u_\tau(t_k)) \ \coloneqq \ b_1 f_1^k,\]

und für den Zwischenwert \(f_1^k \in \R^n\) gilt entsprechend:

\[f_1^k \ = \ F(t_k + c_1 \tau, u_\tau(t_k) + \tau a_{11} f_1^k).\]

Wollen wir nun ein explizites Verfahren herleiten, so muss nach Definition 1.8 schon \(a_{11}=0\) gelten. Die Verfahrensfunktion ist also von der Form

\[f_\tau(t_k, u_\tau(t_k)) \ = \ b_1 f_1^k \ = \ b_1 F(t_k+c_1 \tau, u_\tau(t_k)).\]

Die einzige sinnvolle Wahl für den Koeffizienten \(c_1\) ist Null, da wir sonst die Funktion \(F\) zu einer anderen Zeit auswerten als die unbekannte Funktion \(u\). Um ein konsistentes Verfahren zu erhalten sieht man außerdem ein, dass \(b_1=1\) gelten muss. Also erhalten wir schlussendlich das Vorwärts-Euler Verfahren.

Im impliziten Fall können wir interessanterweise eine höhere Ordnung erreichen. Wir berechnen hierzu die Taylor-Entwicklung des Zwischenwerts \(f_1^k\) wie folgt

(1.11)\[\begin{split}\begin{split} f_1^k \ &= \ F(t_k + c_1 \tau, u_\tau(t_k) + \tau a_{11} f_1^k) \\ &= \ F(t_k,u(t_k)) + c_1 \tau \partial_t F(t_k,u(t_k)) + \tau a_{11} D_u F(t_k,u(t_k)) f_1^k + {\cal O}(\tau^2). \end{split}\end{split}\]

Hierbei bezeichnet \(D_u\) die Jacobimatrix bezüglich des zweiten Arguments \(u\) von \(F\). Setzen wir auf der rechten Seite nochmal den ersten Approximationsterm für \(f_1^k\) ein, so folgt

\[f_1^k \ = \ F(t_k,u(t_k)) + c_1 \tau \partial_t F(t_k,u(t_k)) + \tau a_{11} D_u F(t_k,u(t_k)) F(t_k,u(t_k))+ {\cal O}(\tau^2).\]

Andererseits gilt mit einer Taylorentwicklung des unbekannten Funktionswerts \(u(t_{k+1})\) im Punkt \(t_k\):

(1.12)\[\begin{split}\begin{split} \frac{1}\tau \int_{t_k}^{t_{k+1}} F(t,u(t))\,\mathrm{d}t \ &= \ \frac{1}{\tau} \left( u(t_{k+1}) - u(t_k) \right) \\ &= \ \frac{1}{\tau} \left( u(t_k) + \tau u'(t_k) + \frac{\tau^2}{2}u''(t_k) + \mathcal{O}(\tau^3) - u(t_k) \right)\\ &= \ u'(t_k) + \frac{\tau}{2} u''(t_k) + \mathcal{O}(\tau^2)\\ &= \ F(t_k,u(t_k)) + \frac{\tau}2 \partial_t F(t_k,u(t_k))\\ & \hspace{2.58cm} +\frac{\tau}2 D_u F(t_k,u(t_k)) F(t_k,u(t_k)) + {\cal O}(\tau^2). \end{split}\end{split}\]

Damit gilt

\[f_\tau(t_k, u_\tau(t_k)) \ = \ b_1 f_1^k \ \approx \ \frac{1}{\tau}\int_{t_k}^{t_{k+1}} F(t,u(t))\,\mathrm{d}t,\]

vergleichen wir die beiden Formeln (1.11) und (1.12) und sehen ein, dass wir eine Konsistenzordnung von Zwei erreichen können, wenn die Koeffizienten als \(b_1 =1\), \(c_1=\frac{1}2\) und \(a_{11}=\frac{1}2\) gewählt werden.

Die Verfahrensfunktion ist also gegeben durch die Lösung von

\[f_\tau(t_k, u_\tau(t_k)) \ = \ f_1^k \ = \ F(t_k +\frac{\tau}2,u_\tau(t_k) + \frac{\tau}2 f_1^k).\]

Wir können dieses Runge-Kutta Verfahren als Mittelpunktsregel im Intervall \([t_k,t_{k+1}]\) interpretieren, wobei der unbekannte Wert von \(u_\tau\) am Mittelpunkt \(t_k +\frac{\tau}2\) des Intervalls durch das Rückwärts-Euler-Verfahren bestimmt wird.

Im Folgenden betrachten wir noch \(2\)-stufige Runge–Kutta Verfahren.

Beispiel 1.7 (Runge–Kutta Verfahren der Stufe 2)

Für ein Runge–Kutta Verfahren der Stufe \(s=2\) ist die Verfahrensfunktion gegeben durch

\[f_\tau(t_k, u_\tau(t_k)) \ \coloneqq \ b_1 f_1^k +b_2 f_2^k,\]

und für die Zwischenwerte \(f_1^k, f_2^k \in \R^n\) gilt im expliziten Fall mit \(a_{11} = a_{22} = 0\):

\[f_1^k \ = \ F(t_k + c_1 \tau, u_\tau(t_k)), \quad f_2^k \ = \ F(t_k + c_2 \tau, u_\tau(t_k) + \tau a_{21} f_1^k).\]

Analog zur Argumentation im expliziten Fall von Beispiel 1.6 sehen wir wieder, dass nur \(c_1=0\) eine sinnvolle Wahl ist. Eine Taylor-Entwicklung des Zwischenwertes \(f_2^k\) liefert dann

\[\begin{split}\begin{aligned} b_1 f_1^k +b_2 f_2^k \ &= \ (b_1 + b_2) F(t_k , u(t_k)) + b_2 c_2 \tau \partial_t F(t_k , u (t_k))\\ & \hspace{1cm} + b_2 a_{21} \tau D_u F(t_k , u (t_k)) F(t_k , u (t_k)) + {\cal O}(\tau^2). \end{aligned}\end{split}\]

Vergleichen wir dies wieder mit der Taylor-Entwicklung des Integrals in (1.12), so erhalten wir folgendes Gleichungssystem

\[b_1+b_2 \, = \, 1, \qquad b_2 c_2 \, = \, \frac{1}2, \qquad b_2 a_{21} \, = \, \frac{1}2.\]

Eine einfache Lösung ist beispielsweise \(b_1=0\), \(b_2=1\), \(c_2=a_{21} = \frac{1}2\). Diese Wahl der Koeffizienten liefert uns also ein Verfahren der Konsistenzordnung zwei mit der Verfahrensfunktion

\[f_\tau(t_k, u_\tau(t_k)) \ = \ F(t_k +\frac{\tau}2,u_\tau(t_k) + \frac{\tau}2 F(t_k,u(t_k))).\]

Wir können dieses zweistufige Runge-Kutta Verfahren als eine Mittelpunktsregel im Intervall \([t_k,t_{k+1}]\) interpretieren, wobei der unbekannte Wert von \(u_\tau\) am Mittelpunkt \(t_k +\frac{\tau}2\) durch das Vorwärts-Euler-Verfahren bestimmt wird.

Bei der Konstruktion eines Runge-Kutta Verfahrens in Beispiel 1.6 und Beispiel 1.7 haben wir die Koeffizienten \(b_i, c_i \in \R\) und \(a_{i,j} \in \R\) für \(i,j=1,\ldots,s\) basierend auf vernünftigen Überlegungen gewählt. Wie wir speziell im Fall des zweistufigen Runge-Kutta Verfahrens gesehen haben führt dies im Allgemeinen jedoch nicht zu eindeutigen Lösungen. Daher ist es sinnvoll zusätzliche Kriterien an die implizierten Runge-Kutta Verfahren zu stellen, so dass wir eindeutige Koeffizienten herleiten können.

Am einfachsten ist ein Kriterium für die Konsistenz des Runge-Kutta Verfahrens für die unbekannten Koeffizienten abzuleiten, wie das folgende Lemma festhält.

Lemma 1.4 (Konsistenzbedingung von Runge-Kutta-Verfahren)

Ein Runge-Kutta Verfahren der Stufe \(s \in \N^+\) ist genau dann konsistent, wenn die folgende Bedingung erfüllt ist:

\[\sum_{i=1}^s b_i \ = \ 1.\]

Proof. Um zu zeigen, dass ein Runge-Kutta Verfahren konsistent ist müssen wir zeigen, dass der globale Konsistenzfehler

\[K_\tau \ = \ \max_{k=0,\ldots,N} \left\Vert f_\tau(t_k, u(t_k)) - \frac{1}{\tau} \int_{t_k}^{t_{k+1}} F(t, u(t)) \, \mathrm{d}t \right\Vert\]

mindestens von der Fehlerordnung \(\mathcal{O}(\tau)\) ist und für \(\tau \rightarrow 0\) gegen Null geht.

Die Verfahrensfunktion des Runge-Kutta Verfahrens der Stufe \(s\) gegeben ist als

\[f_\tau(t_k, u_\tau(t_k)) \ = \ \sum_{i=1}^s b_i f_i^k\]

und daher können wir die Taylorapproximation erster Ordnung jedes Zwischenwerts \(f_i^k \in \R^n\) für \(i=1,\ldots,s\) analog zu (1.11) betrachten und erhalten somit

\[\sum_{i=1}^s b_i f_i^k \ = \ \sum_{i=1}^s (b_i F(t_k,u(t_k)) + {\cal O}(\tau)) \ = \ \sum_{i=1}^s b_i F(t_k,u(t_k)) + {\cal O}(\tau).\]

Die Taylorapproximation erster Ordnung des Integrals liefert nach (1.12) außerdem

\[\frac{1}\tau \int_{t_k}^{t_{k+1}} F(t,u(t)) \, \mathrm{d}t = F(t_k,u(t_k)) + {\cal O}(\tau).\]

Nun gilt also für den Konsistenzfehler

\[\begin{split}\begin{split} K_\tau \ &= \ \max_{k=0,\ldots,N} \left\Vert \, f_\tau(t_k, u(t_k)) - \frac{1}{\tau} \int_{t_k}^{t_{k+1}} F(t, u(t)) \, \mathrm{d}t \, \right\Vert \\ &= \ \max_{k=0,\ldots,N} \left\Vert \, \sum_{i=1}^s b_i F(t_k,u(t_k)) - F(t_k,u(t_k)) + {\cal O}(\tau) \, \right\Vert. \end{split}\end{split}\]

Wir sehen, dass für den Konsistenzfehler \(K_\tau \in \mathcal{O}(\tau)\) genau dann gilt, wenn \(\sum_{i=1}^s b_i = 1\) ist. ◻

Umgekehrt ist die Konsistenzordnung eines Runge-Kutta Verfahrens nach oben durch dessen Stufe beschränkt, wie folgendes Lemma zeigt.

Lemma 1.5 (Maximale Konsistenzordnung von Runge-Kutta-Verfahren)

Für die Konsistenzordnung \(p \in \N^+\) eines expliziten Runge-Kutta Verfahren der Stufe \(s \in \N^+\) für ein Anfangswertproblem mit einer rechten Seite der Differentialgleichung \(F \in C^\infty([0,T] \times \R^n; \R^n)\) gilt, dass die Konsistenzordnung durch die Stufe des Verfahrens nach oben beschränkt ist, d.h., es gilt \(0 < p \leq s\).

Proof. Wir zeigen die Behauptung für den reellen Fall \(n=1\) und wählen die besonders einfache rechte Seite \(F(t,u) = u\) mit \(u_0=1 \in \R\). Dann ist die Lösung des Anfangswertproblems

\[u'(t) \ = \ F(t, u(t)) \ = \ u(t), \qquad u(0) \ = \ 1,\]

gegeben durch \(u(t)=u_0\cdot e^t = e^t\).

Mit der Taylorentwicklung der Exponentialfunktion \(u(t) = e^t = \sum_{j=0}^\infty \frac{t^j}{j!}\) können wir das Integral als Polynom in \(\tau\) schreiben mit:

\[\begin{split}\begin{aligned} \frac{1}\tau \int_{t_k}^{t_{k+1}} F(t,u(t))\,\mathrm{d}t \ &= \ \frac{1}\tau \int_{t_k}^{t_{k+1}} u(t)\,\mathrm{d}t \\ &= \frac{1}\tau \int_{t_k}^{t_{k+1}} \underbrace{u(t_k - t_k)}_{=\, 1} \cdot \: u(t)\,\mathrm{d}t\\ &= \frac{1}\tau \int_{t_k}^{t_{k+1}} u(t_k) \cdot u(t - t_k)\,\mathrm{d}t\\ &= \ \frac{u(t_k)}\tau \int_{t_k}^{t_{k+1}} \sum_{j=0}^p \frac{(t-t_k)^j}{j!} + {\cal O}(\tau^{p+1}) \,\mathrm{d}t \\ &= \ \frac{u(t_k)}\tau \sum_{j=0}^p \frac{(t_{k+1}-t_k)^{j+1}}{(j+1)!} + {\cal O}(\tau^{p}) \\ &= \ \sum_{j=0}^p \tau^{j} \frac{u(t_k)}{(j+1)!} + {\cal O}(\tau^{p}). \end{aligned}\end{split}\]

Andererseits sehen wir für ein beliebiges explizites Runge-Kutta Verfahren der Stufe \(s \in \N^+\) mit Konsistenzordnung \(p \leq s\), dass die Zwischenwerte gegeben sind durch:

\[\begin{split}\begin{split} f_1^k \ &= \ u(t_k),\\ f_2^k \ &= \ u(t_k) + \tau a_{21} u(t_k),\\ f_3^k \ &= \ u(t_k) + \tau a_{31} u(t_k) + \tau a_{32} u(t_k) + \tau^2 a_{32} a_{21} u(t_k),\\ &\vdots \end{split}\end{split}\]

Wir sehen also, dass jeder Zwischenwert \(f_i^k \in \R^n\) für \(i=1,\ldots, s\) ein Polynom in \(\tau\) vom Grad kleiner gleich \(i-1\) ist und somit ist die Verfahrensfunktion \(f_\tau(t_k, u_\tau(t_k)) \ = \ \sum_{i=1}^s b_i f_i^k\) ein Polynom in \(\tau\) vom Grad kleiner gleich \(s-1\).

Damit können wir keine höhere Ordnung erreichen, da der Fehler ab der Ordnung \(\mathcal{O}(\tau^p)\) im Allgemeinen nicht verschwindet. ◻

Man sieht ein, dass sich mittels der Forderung von Konsistenz in Lemma 1.4 und dem Bestreben nach maximaler Konsistenzordnung \(p = s\) eines \(s\)-stufigen Runge-Kutta Verfahrens, das enstehende Gleichungssystem lösen lässt. Allerdings bleibt hierbei immer noch eine Uneindeutigkeit der Koeffizienten \(c_i \in \R, i=1,\ldots,s\) übrig. Um hier eine systematische Wahl zu treffen, fordert man im Allgemeinheit die Invarianz des Verfahrens gegenüber sogenannter Autonomisierung. Hierzu führen wir zunächst den Begriff einer autonomen Differentialgleichung ein.

Definition 1.9 (Autonome gew. Differentialgleichung)

Sei \(u'(t) = F(t, u(t)), \ u(0) = u_0 \in \R^n\) ein Anfangswertproblem einer gewöhnlichen Differentialgleichung erster Ordnung.

Wir nennen die Differentialgleichung autonom, wenn die Funktion \(F\) auf der rechten Seite nicht explizit von \(t\) abhängt, d.h., es gilt \(F(t,u(t)) = F(u(t))\).

Das allgemeine Anfangswertproblem (1.1) lässt sich autonomisieren, d.h. in ein äquivalentes System autonomer Differentialgleichungen für \(\tilde u = \left(v,u\right)\) mit einer Hilfsfunktion \(v\) umschreiben. Hierzu betrachten wir das Differentialgleichungssystem

\[\begin{split}\begin{split} u'(t) \ &= \ F(v(t),u(t)) \ = \ F(\tilde{u}(t)), \qquad u(0) \, = \, u_0 \in \R^n,\\ v'(t) \ &= \ 1, \hspace{4.7cm} v(0) \, = \, 0. \end{split}\end{split}\]

Man erkennt sofort, dass \(v(t) =t\) gelten muss, was die Äquivalenz zu (1.1) impliziert.

Wir nennen die obige Transformation des Anfangswertproblems Autonomisierung und fordern, dass das Runge-Kutta Verfahren unter der Autonomisierung invariant sein soll. Dies bedeutet, dass wir fordern, dass die Anwendung des numerischen Verfahrens auf das neue System die selbe Lösung \(u_\tau\) liefern soll, wie die Anwendung auf das ursprüngliche Anfangswertproblem (1.1). Darüber hinaus fordern wir, dass die Invarianz unter Autonomisierung für jede geeignete Funktion \(F\) auf der rechten Seite gelten soll.

Schreiben wir die Zwischenwerte für beide Formulierungen, so gilt für das ursprüngliche Anfangswertproblem

\[F(t_k + c_i \tau, u(t_k + c_i \tau)) \ \approx \ F(t_k + c_i \tau,u(t_k) + \tau \sum_{j=1}^s a_{ij} f_j^k) \ =: \ f_i^k,\]

und für die autonome Variante andererseits

\[F(v(t_k + c_i \tau), u(t_k + c_i \tau)) \ \approx \ F(v(t_k) + \tau \sum_{j=1}^s a_{ij},u(t_k) + \tau \sum_{j=1}^s a_{ij} f_j^k) \ =: \ f_i^k.\]

Da wegen der exakten numerischen Integration der linearen Funktion \(t \mapsto t\) immer \(v(t_k) = t_k\) gilt, sehen wir durch Koeffizientenvergleich, dass die Invarianz gegenüber Autonomisierung äquivalent ist zu der Bedingung

(1.13)\[c_i \ = \ \sum_{j=1}^s a_{ij}, \qquad i=1,\ldots,s.\]

Aus diesem Grund bestimmt man die Koeffizienten \(c_i\in \R\) immer aus Gleichung (1.13) und lediglich die Einträge \(a_{ij} \in \R\) der Matrix \(A \in \R^{s\times s}\) werden aus der Konsistenzbedingung mittels Taylorentwicklung bestimmt. Damit ist es möglich ein eindeutiges, explizites \(s\)-stufiges Runge-Kutta Verfahren der Konsistenzordnung \(s \in N^+\) zu konstruieren

1.2.3.1. Butcher-Schema

Allgemein lässt sich ein \(s\)-stufiges Runge-Kutta Verfahren durch die Matrix \(A \in \R^{s\times s}\) und die Vektoren \(b, c \in \R^{s}\) eindeutig repräsentieren. Diese lassen sich kompakt in Form eines Schemas angeben.

Definition 1.10 (Butcher-Schema)

Sei \(s \in \N^+\) die Stufe eines Runge-Kutta Verfahrens mit Koeffizienten, die gegeben sind durch eine Matrix \(A \in \R^{s \times s}\) und den beiden Vektoren \(\vec{c}, \vec{b} \in \R^s\). Dann lässt sich das Runge-Kutta Verfahren kompakt durch das folgende Butcher-Schema repräsentieren:

\[\begin{split}\begin{array}{c|c} \vec{c} & A \\ \hline & \vec{b}^T \end{array} \ = \ \begin{array}{c|cccc} c_1 & a_{11} & a_{12} & \cdots & a_{1s}\\ c_2 & a_{21} & a_{22} & \cdots & a_{2s}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & \cdots & a_{ss}\\ \hline & b_1 & b_2 & \cdots & b_s \end{array}.\end{split}\]

Im folgenden Beispiel geben wir die Butcher-Schemata für vier explizite Runge-Kutta Verfahren im Vergleich an.

Beispiel 1.8 (Explizite Runge–Kutta Verfahren)

\(s = 1\)	$\begin{array}{c|c}
0 & 0\ \hline
& 1
\end{array}$	Vorwärts-Euler Verfahren mit:
\(u_\tau(t_{k+1}) \ = \ u_\tau(t_k) + \tau F(t_k, u_\tau(t_k))\).
Konsistenzordnung \(p=1\).
\(s = 2\)	$\begin{array}{c|cc}
0 & 0 & 0\
1 & 1 & 0\ \hline
& \frac{1}{2} & \frac{1}{2}
\end{array}$	Verbessertes Euler Verfahren mit:
$u_\tau(t_{k+1}) \ = \ u_\tau(t_k) + \frac{\tau}{2} F(t_k, u_\tau(t_k)) \newline
\hphantom{aaaaaaaaaas} + \frac{\tau}{2} F(t_k + \tau, u_\tau(t_k) + \tau F(t_k, u_\tau(t_k)))$.
Konsistenzordnung \(p=2\).
\(s = 3\)	$\begin{array}{c|ccc}
0 & 0 & 0 & 0\
\frac{1}{2} & \frac{1}{2} & 0 & 0\
1 & -1 & 2 & 0\ \hline
& \frac{1}{6} & \frac{4}{6} & \frac{1}{6}
\end{array}$	$u_\tau(t_{k+1}) \ = \ u_\tau(t_k) + \frac{\tau}{6} (f_1^k +4f_2^k + f_3^k), \newline
\hphantom{aaaaa} f_1^k \ = \ F(t_k, u_\tau(t_k)),\newline
\hphantom{aaaaa} f_2^k \ = \ F(t_k + \frac{\tau}{2}, u_\tau(t_k)) + \frac{\tau}{2}f_1^k),\newline
\hphantom{aaaaa} f_3^k \ = \ F(t_k + \tau, u_\tau(t_k)) - \tau f_1^k + 2\tau f_2^k)$.
Konsistenzordnung \(p=3\).
\(s = 4\)	$\begin{array}{c|cccc}
0 & 0 & 0 & 0 & 0\
\frac{1}{2} & \frac{1}{2} & 0 & 0 & 0\
\frac{1}{2} & 0 & \frac{1}{2} & 0 & 0\
1 & 0 & 0 & 1 & 0\ \hline
& \frac{1}{6} & \frac{2}{6} & \frac{2}{6} & \frac{1}{6}
\end{array}$	Standard Runge-Kutta Verfahren mit:
$u_\tau(t_{k+1}) \ = \ u_\tau(t_k) + \frac{\tau}{6} (f_1^k +2f_2^k + 2f_3^k + f_4^k), \newline
\hphantom{aaaaa} f_1^k \ = \ F(t_k, u_\tau(t_k)),\newline
\hphantom{aaaaa} f_2^k \ = \ F(t_k + \frac{\tau}{2}, u_\tau(t_k)) + \frac{\tau}{2}f_1^k),\newline
\hphantom{aaaaa} f_3^k \ = \ F(t_k + \frac{\tau}{2}, u_\tau(t_k)) + \frac{\tau}{2}f_2^k),\newline
\hphantom{aaaaa} f_4^k \ = \ F(t_k + \tau, u_\tau(t_k)) + \tau f_3^k)$.
Konsistenzordnung \(p=4\).