Numerische Mathematik 2 – Anfangswertprobleme für gewöhnliche Differentialgleichungen

Wiederholung: Landau-Notation und Taylor-Entwicklung

Landau-Notation: Seien $f, g :] 0, + \infty [\to R^{n}$ Abbildungen.
Man schreibt $f = O (g)$ , falls $\exists_{c > 0} \exists_{δ > 0} \forall_{x \in] 0, δ [} ∥ f (x) ∥ \leq c ∥ g (x) ∥$ .
Man schreibt $f = o (g)$ , falls $\forall_{ε > 0} \exists_{δ > 0} \forall_{x \in] 0, δ [} ∥ f (x) ∥ \leq ε ∥ g (x) ∥$ .

Beispiel: $f = O (1)$ gilt genau dann, wenn $f$ in einer $δ$ -Umgebung von $0$ beschränkt ist.
$f = o (1)$ ist äquivalent zu $lim_{x \to 0} f (x) = 0$ .
$f = o (x)$ ist äquivalent zu $\tilde{f} = o (1)$ mit $f (x) = x \tilde{f} (x)$ .

Satz (Taylor-Entwicklung):
Seien $U \subset R$ ein Intervall und $f : U \subset R \to R^{n}$ in $x_{0} \in U$ $(m + 1)$ -fach stetig differenzierbar.
Dann gilt $f (x_{0} + h) = f (x_{0}) + \sum_{k = 1}^{m} \frac{1}{k!} f^{(k)} (x_{0}) h^{k} + r_{m} (x_{0}, h)$ mit
$r_{m} (x_{0}, h) = \frac{1}{(m + 1)!} f^{(m + 1)} (y) h^{m + 1}$ für ein $y \in \bar{x_{0}, x_{0} + h}$ , d. h.
$r_{m} (x_{0}, h) = O (h^{m + 1})$ . Es gilt auch $r_{m} (x_{0}, h) = o (h^{m})$ .

Motivation, Beispiele

Bemerkung: Gegeben seien eine Funktion $f : R \times R \to R$ , $t_{0} \in R$ und $u_{0} \in R$ . Gesucht ist eine differenzierbare Funktion $u = u (t) : R \to R$ , sodass $u^{'} (t) = f (t, u (t))$ für $t \geq t_{0}$ . Dieses Problem heißt Anfangswertproblem.

Beispiel: Sei $t$ die Zeit und $P (t)$ eine Population. Für die Zunahme $Δ P := P (t + Δ t) - P (t)$ in der Zeit $Δ t$ soll $Δ P \approx α P (t) Δ t$ mit $α > 0$ gelten. Für $Δ t \to 0$ erhält man die DGL $\frac{d P}{d t} = α P (t)$ . Sie hat die allgemeine Lösung $P (t) = c \cdot e^{α t}$ mit $c$ beliebig (exponentielles Wachstum). Ist ein Anfangswert $P_{0} = P (t_{0})$ gegeben, so bestimmt sich $c$ durch $c = P_{0} e^{- α t_{0}}$ , d. h. die partikuläre Lösung ist $P (t) = P_{0} e^{α (t - t_{0})}$ .

Beispiel: Die DGL $\frac{d P}{d t} = λ P (K - P)$ mit $λ, K > 0$ modelliert logistisches Wachstum. Zum Beispiel gilt für $P \equiv K$ , dass $\frac{d P}{d t} = 0$ , d. h. $P$ ändert sich nicht. Die DGL hat die Lösung $P (t) = \frac{K}{1 + \frac{K}{P_{0} - 1} e^{- λ K t}}$ .

Beispiel: Eine DGL, mit der das aktuelle Bevölkerungswachstum beschrieben werden kann, lautet $\frac{d P}{d t} = α P (t)^{β}$ mit $α > 0$ , $β > 1$ .

Beispiel: Wird die Menge einer radioaktiven Substanz durch $u = u (t)$ beschrieben, so modelliert die DGL $d u = - λ u d t$ , $λ > 0$ den Zerfall der Substanz aufgrund der Radioaktivität. Für $t_{0} = 0$ lautet eine Lösung $u (t) = u_{0} e^{- λ t}$ . Die Halbwertszeit ist die Zeit, in der sich die Menge der Substanz halbiert. Sie ist unabhängig von der aktuellen Menge und beträgt $τ = \frac{\ln (2)}{λ}$ .

Theoretische Grundlagen

Anfangswertproblem: Seien $U \subset R^{n}$ offen (Zustandsraum), $f \in C (R \times U, R^{n})$ , $u_{0} \in U$ , $I \subset R$ und $t_{0} \in I$ . Gesucht ist eine Funktion $u = (u_{1}, \dots, u_{n})^{t} \in C^{1} (I, U)$ mit $u^{'} (t) = f (t, u (t))$ für $t \in I$ und $u (t_{0}) = u_{0}$ . Dieses Problem heißt Anfangswertproblem (AWP).

Beispiel: Im Räuber-Beute-Modell wird mit $y_{1} (t)$ bzw. $y_{2} (t)$ die Population der Beute- bzw. Raubtiere bezeichnet. Die DGLs $y_{1}^{'} (t) = α y_{1} (t) (1 - y_{2} (t))$ und $y_{2}^{'} (t) = β y_{2} (t) (y_{1} (t) - 1)$ modellieren dann den zeitlichen Verlauf der Populationen.

Existenz und Eindeutigkeit der Lösung des Anfangswertproblems

Satz (Peano): Seien $f$ auf einem kompakten Rechteck
$R := {(t, u) \in R \times U | | t - t_{0} | \leq a, ∥ u - u_{0} ∥ \leq b}$ stetig,
$μ := max_{(t, u) \in R} ∥ f (t, u) ∥ < \infty$ und $α := min (a, \frac{b}{μ})$ .
Dann hat das Anfangswertproblem auf $[t_{0} - α, t_{0} + α]$ mindestens eine Lösung.

Satz (Picard-Lindelöf): Sei zusätzlich $f$ in $R$ im zweiten Argument Lipschitz-stetig, d. h. $∥ f (t, w) - f (t, \tilde{w}) ∥ \leq L ∥ w - \tilde{w} ∥$ für alle $(t, w), (t, \tilde{w}) \in R$ .
Dann existiert für $U = R^{n}$ genau eine Lösung $u \in C^{1} ([t_{0} - α, t_{0} + α], R^{n})$ .

Satz (Banachscher Fixpunktsatz): Seien $(X, ∥ \cdot ∥)$ ein Banachraum und $D \subset X$ eine abgeschlossene Teilmenge mit $D \neq \emptyset$ . Sei außerdem $T : D \to X$ eine Abbildung mit $T (D) \subset D$ und $\exists_{0 < c < 1} \forall_{v, \tilde{v} \in D} ∥ T v - T \tilde{v} ∥ \leq c ∥ v - \tilde{v} ∥$ . Dann gibt es genau ein $u \in D$ , sodass $T u = u$ .

Behandlung von Anfangswertproblemen höherer Ordnung

Bemerkung: Ein Anfangswertproblem höherer Ordnung ist ein Anfangswertproblem der Form $y^{(m)} (t) = f (t, y (t), y^{'} (t), \dots, y^{(m - 1)} (t))$ mit $y^{(i)} (t_{0}) = y_{0, i}$ für $i = 0, \dots, m - 1$ .
Es kann in ein System 1. Ordnung umgeformt werden, indem $z_{1} (t) = y (t)$ , $z_{2} (t) = y^{'} (t)$ , …, $z_{m} (t) = y^{(m - 1)} (t)$ gesetzt wird. Damit ist $z^{'} = (z_{1}^{'}, z_{2}^{'}, \dots, z_{m - 1}^{'}, z_{m}^{'})^{t}$
$= (z_{2}, z_{3}, \dots, z_{m}, f (t, z_{1}, z_{2}, \dots, z_{m}))^{t}$ ein System 1. Ordnung mit der Anfangsbedingung
$z (t_{0}) = (y (t_{0}), y^{'} (t_{0}), \dots, y^{(m - 1)} (t_{0}))^{t} = (y_{0, 0}, y_{0, 1}, \dots, y_{0, m - 1})^{t}$ .

Beispiel: Die elastische Schwingung eines fest eingespannten Federpendels, an dem ein Körper mit Masse $m$ hängt, kann durch die DGL $m y^{}'^{'} (t) + r y^{'} (t) + D (y (t) - ℓ) = g (t)$ beschrieben werden, wenn $y (t)$ die Auslenkung darstellt und $y (0)$ und $y^{'} (0)$ gegeben sind. Umgeformt nach $y^{}'^{'}$ ergibt dies $y^{}'^{'} = \frac{1}{m} (g - D (y - ℓ) - r y^{'})$ . Mit $z_{1} = y$ und $z_{2} = y^{'}$ ist $z^{'} = (z_{1}^{'}, z_{2}^{'})^{t} = (z_{2}, \frac{1}{m} (g - D (z_{1} - ℓ) - r z_{2}))^{t}$ ein System 1. Ordnung mit Anfangsbedingung $z (0) = (y_{0, 0}, y_{0, 1})^{t}$ .

Lösung durch Trennung der Variablen

Bemerkung: Eine DGL hat trennbare Veränderliche, falls sie die Form $y^{'} (t) = f (t) g (y)$ mit $y (t_{0}) = y_{0}$ besitzt. In diesem Fall kann sie mit der Gleichung $\frac{1}{g (y)} d y = f (t) d t$ und anschließendem Integrieren, also $\int_{y_{0}}^{y} \frac{1}{g (z)} d z = \int_{t_{0}}^{t} f (s) d s$ , gelöst werden, indem nach $y (t)$ umgeformt und die Integrationskonstante mit der Anfangsbedingung berechnet wird.

Satz (Korrektheit der Trennung der Veränderlichen): Seien $f \in C (I_{t}, R)$ , $g \in C (I_{y}, R)$ und $t_{0}$ bzw. $y_{0}$ seien aus dem Inneren von $I_{t}$ bzw. $I_{y}$ . In diesem Fall ist die obige DGL mit dem eben beschriebenen Algorithmus in einer Umgebung von $t_{0}$ eindeutig lösbar.

Spezielle Typen von DGL 1. Ordnung

autonom: Eine DGL $u^{'} (t) = f (t, u (t))$ heißt autonom, falls $u^{'} (t) = f (u (t))$ .

linear: Eine DGL $u^{'} (t) = f (t, u (t))$ heißt linear, falls $u^{'} (t) = A (t) u (t) + b (t)$ mit
$A \in C (I, R^{n \times n})$ und $b \in C (I, R^{n})$ .
Eine lineare DGL heißt homogen, falls $b \equiv 0$ , sonst heißt sie inhomogen/affin.

Satz (eindeutige Lösbarkeit linearer DGLs):
Sei $u^{'} (t) = A (t) u (t) + b (t)$ eine lineare DGL mit $A \in C (I, R^{n \times n}) \cap L^{\infty} (R, R^{n \times n})$ .
Dann hat das Anfangswertproblem genau eine Lösung in $C^{1} (I, R^{n})$ .

Satz (Lösungen linearer DGLs): Unter den Voraussetzungen von eben gilt:

Die Lösungen der homogenen DGL $u^{'} (t) = A (t) u (t)$ bilden einen $n$ -dimensionalen Unterraum $V \subset C^{1} (R, R^{n})$ mit einer Basis $u_{i} \in C^{1} (R, R^{n})$ , $u_{i} (0) = e_{i}$ , $i = 1, \dots, n$ .
Die normierte Fundamentalmatrix ist $Y_{0} (t) = (u_{1}, \dots, u_{n})$ .
Die Lösungen der inhomogen DGL $u^{'} (t) = A (t) u (t) + b (t)$ bilden einen affinen Unterraum $\tilde{u} + V \subset C^{1} (R, R^{n})$ mit einer speziellen Lösung $\tilde{u}$ . Für die Lösung gilt
$u (t) = Y_{0} (t) u_{0} + \int_{0}^{t} Y_{0} (t) (Y_{0} (s))^{- 1} b (s) d s$ (dabei sei $t_{0} = 0)$ .
Ist die DGL autonom, d. h. ist $u^{'} (t) = A u (t)$ , so gilt $Y_{0} (t) = e^{A t} := \sum_{n = 0}^{\infty} \frac{t^{n}}{n!} A^{n}$ .

Beispiel: $(\begin{matrix} u_{1}^{'} (t) \\ u_{2}^{'} (t) \end{matrix})$ $=$ $(\begin{matrix} u_{2} (t) \\ u_{1} (t) \end{matrix})$ $+$ $(\begin{matrix} 1 \\ 0 \end{matrix})$ , d. h. $A =$ $(\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix})$ . Es gilt $A^{2} = E_{2}$ , d. h.
$Y_{0} (t) = \sum_{n = 0}^{\infty} \frac{t^{n}}{n!} A^{n} = \sum_{k = 0}^{\infty} \frac{t^{2 k + 1}}{(2 k + 1)!} A + \sum_{k = 0}^{\infty} \frac{t^{2 k}}{(2 k)!} E_{2} = \sinh (t) A + \cosh (t) E_{2} =$ $(\begin{matrix} \cosh (t) & \sinh (t) \\ \sinh (t) & \cosh (t) \end{matrix})$ .
Wegen $det Y_{0} (t) = \cosh^{2} (t) - \sinh^{2} (t) = 1$ gilt $Y_{0}^{- 1} (t) =$ $(\begin{matrix} \cosh (t) & - \sinh (t) \\ - \sinh (t) & \cosh (t) \end{matrix})$ und somit ist die Lösung $u (t) = 2$ $(\begin{matrix} \sinh (t) \\ \cosh (t) \end{matrix})$ $-$ $(\begin{matrix} 0 \\ 1 \end{matrix})$ .

Beispiel: Für $u^{'} (t) = t^{3} u + e^{t}$ , $u (0) = 1$ ist die homogene DGL $u_{h}^{'} = t^{3} u_{h}$ , deren Lösung ist $u_{h} (t) = e^{t^{4} / 4} = Y_{0} (t)$ . Die allgemeine Lösungsformel von oben ergibt nun
$u (t) = e^{t^{4} / 4} + e^{t^{4} / 4} \cdot \int_{0}^{t} e^{τ - τ^{4} / 4} d τ$ , jedoch kann das Integral analytisch nicht berechnet werden.
Die bestehenden Möglichkeiten sind nun einerseits das Anwenden einer Quadraturformel für das Integral, zum anderen numerische Verfahren für das Ausgangsproblem.

Einzelschrittverfahren

Einzelschrittverfahren: Angenommen, das Anfangswertproblem besitzt eine eindeutige Lösung $u \in C^{1} (I, R^{n})$ . Seien $t_{0} := 0$ und $I := [0, T]$ mit $T > 0$ .
Ein Schrittweitenvektor ist ein Vektor $h := (h_{0}, \dots, h_{N - 1})^{t} \in [0, T]^{N}$ mit $\sum_{j = 0}^{N - 1} h_{j} = T$ .
Das Gitter $I_{h}$ zu $h$ ist $I_{h} := {0 = t_{0}, t_{1}, \dots, t_{N} = T}$ mit $t_{j} := t_{j - 1} + h_{j - 1}$ .
Das Gitter heißt äquidistant, falls $h_{0} = \dots = h_{N - 1}$ . In diesem Fall sei $h$ skalar ( $h = h_{0}$ ).
Die Gitterweite ist $| h | := max_{j = 0, \dots, N - 1} h_{j}$ .
Das Ziel ist die Bestimmung einer Gitterfunktion $u_{h} : I_{h} \to R^{n}$ . Dabei setzt man $u_{j} := u_{h} (t_{j})$ für $j = 0, \dots, N$ .

Das Eulersche Polygonzugverfahren

Bemerkung: Zur Vereinfachung setzt man $n = 1$ , $I_{h}$ äquidistant und $u_{h} (t_{0}) = u_{0}$ .
Für die exakte Lösung $u$ des Anfangswertproblems gilt $u (t_{1}) = u (t_{0}) + u^{'} (t_{01}) (t_{1} - t_{0})$
$= u_{0} + h f (t_{01}, u (t_{01}))$ mit $t_{01} \in [t_{0}, t_{1}]$ (Taylorformel mit Restglied).
Mittels $t_{01} \approx t_{0}$ erhält man eine Näherung $u_{1} = u_{h} (t_{1})$ für $u (t_{1})$ , wobei $u_{1} = u_{0} + h f (t_{0}, u_{0})$ .

explizites Euler-Verfahren:
Das explizite Euler-Verfahren hat die Iterationsvorschrift $u_{j} := u_{j - 1} + h f (t_{j - 1}, u_{j - 1})$ .

Beispiel: Für $u^{'} (t) = t^{3} u + e^{t}$ , $u (0) = 1$ und $t \in [0, 1]$ erhält man schon für geringe $N$ gute Näherungen. Bei $u^{'} (t) = \sin (t) u (t)$ , $u (0) = 1$ (exakte Lösung $u (t) = e^{1 - \cos (t)}$ ) und $t \in [0, 50]$ benötigt man schon wesentlich größere Werte für $N$ , um sinnvolle Näherungen zu erzeugen.

Allgemeine Definition, Beispiele

explizites Einschrittverfahren: Es seien ein Gitter $I_{h}$ und eine Funktion
$ϕ \in C ([0, T]^{2} \times R^{n}, R^{n})$ gegeben. Dann heißt das Verfahren $u_{j} := u_{j - 1} + h_{j - 1} ϕ (h_{j - 1}, t_{j - 1}, u_{j - 1})$ , $j = 1, \dots, N$ explizites Einschrittverfahren (ESV) und $ϕ$ heißt zugehörige Inkrementfunktion.

Beispiel: Im Euler-Verfahren setzt man $u^{'} (t_{01}) \approx u^{'} (t_{0}) = f (t_{0}, u_{0})$ .
Man kann dies auch anders approximieren: $u^{'} (t_{01}) \approx f (t_{0} + \frac{h}{2}, u (t_{0} + \frac{h}{2}))$ mit
$u (t_{0} + \frac{h}{2}) \approx u (t_{0}) + \frac{h}{2} u^{'} (t_{0}) = u_{0} + \frac{h}{2} f (t_{0}, u_{0})$ . Daraus ergibt sich die neue Iterationsvorschrift $u_{j} := u_{j - 1} + h_{j - 1} f (t_{j - 1} + \frac{h_{j - 1}}{2}, u_{j - 1} + \frac{h_{j - 1}}{2} f (t_{j - 1}, u_{j - 1}))$ , $j = 1, \dots, N$ .
Dieses Verfahren nennt sich modifiziertes explizites Euler-Verfahren.

Beispiel: Ein anderes Verfahren ergibt sich wie folgt: $u^{'} (t_{1}) = u (t_{0}) + (t_{1} - t_{0}) u^{'} (t_{01})$
$= u_{0} + h f (t_{01}, u (t_{01})) = u_{0} + \frac{h}{2} (f (t_{01}, u (t_{01})) + f (t_{01}, u (t_{01}))) \approx u_{0} + \frac{h}{2} (f (t_{0}, u (t_{0})) + f (t_{1}, u (t_{1})))$
$\approx u_{0} + \frac{h}{2} (f (t_{0}, u_{0}) + f (t_{0} + h, u_{0} + h f (t_{0}, u_{0})))$ .
Das sogenannte Verfahren von Heun hat also die Iterationsvorschrift
$u_{j} := u_{j - 1} + \frac{h_{j - 1}}{2} (f (t_{j - 1}, u_{j - 1}) + f (t_{j - 1} + h_{j - 1}, u_{j - 1} + h_{j - 1} f (t_{j - 1}, u_{j - 1})))$ , $j = 1, \dots, N$ .

explizites Euler-Verfahren: Die Inkrementfunktion des expliziten Euler-Verfahrens ist
$ϕ (k, t, w) := f (t, w)$ .

modifiziertes explizites Euler-Verfahren: Die Inkrementfunktion des
modifizierten expliziten Euler-Verfahrens ist $ϕ (k, t, w) := f (t + \frac{k}{2}, w + \frac{k}{2} f (t, w))$ .

Verfahren von Heun: Die Inkrementfunktion des Verfahrens von Heun ist
$ϕ (k, t, w) := \frac{1}{2} (f (t, w) + f (t + k, w + k f (t, w)))$ .

Konsistenz, Konvergenz, Stabilität, numerischer Aufwand

globale Fehlerfunktion/globaler Diskretisierungsfehler:
Die Funktion $e_{h} : I_{h} \to R^{n}$ mit $e_{h} := u |_{I_{h}} - u_{h}$ heißt globale Fehlerfunktion.
Der globale Diskretisierungsfehler ist $\bar{e_{h}} := max_{j = 0, \dots, N} ∥ e_{h} (t_{j}) ∥$ .

lokale Fehlerfunktion/lokaler Diskretisierungsfehler:
Die Funktion $ε_{h} : I_{h} \to R^{n}$ mit $ε_{h} (t_{j}) = \frac{1}{h_{j}} (u (t_{j + 1}) - u (t_{j}) - h_{j} ϕ (h_{j}, t_{j}, u (t_{j})))$ heißt
lokale Fehlerfunktion. Der lokale Diskretisierungsfehler ist $\bar{ε_{h}} := max_{j = 0, \dots, N} ∥ ε_{h} (t_{j}) ∥$ .

Bemerkung: Der lokale Diskretisierungsfehler gibt den Fehler an, der bei einem Schritt gemacht wird. Er kann als Differenz von der Steigung der exakten Lösung $u$ und der Steigung der Approximation $u_{h}$ interpretiert werden.

Konvergenz: Das Einzelschrittverfahren heißt konvergent, falls $\bar{e_{h}} \to 0$ für $| h | \to 0$ .

Konsistenz: Das Einzelschrittverfahren heißt konsistent zu (AWP), falls $\bar{ε_{h}} \to 0$ für $| h | \to 0$ .

Konsistenzordnung:
Das Einzelschrittverfahren heißt konsistent zur Ordnung $p$ zu (AWP), falls $\bar{ε_{h}} = O (| h |^{p})$ .

numerischer Aufwand: Der numerische Aufwand ist die Anzahl der Auswertungen von $f$ .

Beispiel: Konsistenz und numerischer Aufwand der bisher betrachteten Verfahren:
explizites Euler-Verfahren: $p = 1$ und $1$
modifiziertes Euler-Verfahren: $p = 2$ und $2$
Verfahren von Heun: $p = 2$ und $2$

Bemerkung: Der Aufwand pro Zeitschritt ist proportional zu $p$ .

Satz (Konsistenz von Einzelschrittverfahren): Seien $h \in [0, T]^{N}$ ein Schrittweitenvektor, $I_{h}$ ein Gitter und $ϕ \in C ([0, T]^{2} \times R^{n}, R^{n})$ die Inkrementfunktion für ein Einzelschrittverfahren (ESV).

Das Einzelschrittverfahren ist zu (AWP) konsistent genau dann, wenn
$\forall_{t \in I} ϕ (0, t, u (t)) = f (t, u (t))$ .
Seien zusätzlich $f \in C^{p} (I \times R^{n}, R^{n})$ und $ϕ \in C^{p} ([0, T]^{2} \times R^{n}, R^{n})$ .
Dann ist das Einzelschrittverfahren konsistent mit der Ordnung $p$ zu (AWP) genau dann, wenn $\forall_{t \in I} \frac{d^{i}}{d t^{i}} f (t, u (t)) = (i + 1) \frac{\partial^{i}}{\partial k^{i}} ϕ (k, t, u (t)) |_{k = 0}$ für $i = 0, \dots, p - 1$ .

Bemerkung: Was ist der Zusammenhang zwischen dem lokalen Konsistenzfehler $ε_{h}$ und dem globalen Fehler $e_{h}$ ?

Raum der beschränkten Gitterfunktionen: Sei $I_{h}$ ein Gitter zum Schrittweitenvektor $h$ . Die Menge $X_{h} := {v_{h} : I_{h} ∖ {t_{n} = T} \to R^{n} | \exists_{c > 0} \forall_{j = 0, \dots, N - 1} ∥ v_{h} (t_{j}) ∥ \leq c}$ heißt
Raum der beschränkten Gitterfunktionen.
Mit der Norm ${∥ v_{h} ∥}_{\infty} := max_{j = 0, \dots, N - 1} {∥ v_{h} (t_{j}) ∥}_{\infty}$ ist $X_{h}$ ein Banachraum isomorph zu $R^{n N}$ .

diskreter Operator: Seien $I_{h}$ ein Gitter zum Schrittweitenvektor $h$ und $ϕ \in C (I^{2} \times R^{n}, R^{n})$ die Inkrementfunktion für ein Einzelschrittverfahren (ESV). Der Operator $T_{h} : X_{h} \to X_{h}$ ,
$(T_{h} v_{h}) (t_{0}) := v_{h} (t_{0}) - u_{0}$ und $(T_{h} v_{h}) (t_{j}) := \frac{1}{h_{j}} (v_{h} (t_{j + 1}) - v_{h} (t_{j}) - h_{j} ϕ (h_{j}, t_{j}, v_{h} (t_{j})))$ für
$j = 1, \dots, N - 1$ heißt der dem Einzelschrittverfahren (ESV) zugeordnete diskrete Operator.

Bemerkung:
$u_{h}$ ist die Gitterfunktion aus einem Einzelschrittverfahren genau dann, wenn $T_{h} u_{h} = 0$ .
Es gilt $∥ T_{h} (u |_{I_{h}}) ∥ = O (| h |^{p})$ , da $(T_{h} (u |_{I_{h}})) (t_{j}) = ε_{h} (t_{j})$ , falls (ESV) kons. mit Ordn. $p$ ist.

Stabilität: Der Operator $T_{h}$ heißt stabil, falls
$\exists_{c, \bar{h} > 0} \forall_{h \in [0, T]^{N}, | h | < \bar{h}} \forall_{v_{h}^{(1)}, v_{h}^{(2)} \in X_{h}} {∥ v_{h}^{(1)} - v_{h}^{(2)} ∥}_{\infty} \leq c {∥ T_{h} v_{h}^{(1)} - T_{h} v_{h}^{(2)} ∥}_{\infty}$ .

Bemerkung: Sei das Einzelschrittverfahren (ESV) stabil. Dann gilt:
Die Lösung $u_{h}$ von $T_{h} u_{h} = 0$ ist eindeutig, denn ${∥ u_{h} - {\tilde{u}}_{h} ∥}_{\infty} \leq c {∥ T_{h} u_{h} - T_{h} {\tilde{u}}_{h} ∥}_{\infty} = 0$ .
Die Lösung $u_{h}$ von $T_{h} u_{h} = 0$ ist beschränkt, denn
${∥ u_{h} ∥}_{\infty} = {∥ u_{h} - 0 ∥}_{\infty} \leq c {∥ T_{h} u_{h} - T_{h} [0] ∥}_{\infty} = c c_{0}$ für $c_{0} := {∥ T_{h} u_{h} ∥}_{\infty}$ .

Satz (Konvergenz von Einzelschrittverfahren I):
Sei ein ESV mit Inkrementfunktion $ϕ \in C (I^{2} \times R^{n}, R^{n})$ gegeben. Ist das ESV stabil, so gilt:

Ist das ESV konsistent, so ist es auch konvergent.
Ist das ESV konsistent zur Ordnung $p \in N$ , so gilt $\bar{e_{h}} = O (| h |^{p})$ .

Satz (Konvergenz von Einzelschrittverfahren II):
Sei ein ESV mit Inkrementfunktion $ϕ \in C (I^{2} \times R^{n}, R^{n})$ gegeben. Außerdem existiere für $\bar{h}$ fest eine Konstante $M > 0$ mit $\forall_{k \in (0, \bar{h})} \forall_{t \in I} \forall_{w, \tilde{w} \in R^{n}} {∥ ϕ (k, t, w) - ϕ (k, t, \tilde{w}) ∥}_{\infty} \leq M \cdot {∥ w - \tilde{w} ∥}_{\infty}$
(globale Lipschitz-Bedingung an $ϕ$ im dritten Argument).

Das ESV ist stabil.
Ist das ESV konsistent, so ist es auch konvergent.
Ist das ESV konsistent zur Ordnung $p$ , so existiert eine Konstante $c > 0$ , sodass für alle Gitter $I_{h}$ mit $| h | < \bar{h}$ die Abschätzung $\bar{e_{h}} \leq c c_{s} | h |^{p}$ gilt, wobei $c_{s} := e^{M T} (T + 1)$ die Stabilitätskonstante und $I = [0, T]$ ist.

Bemerkung: Die Abschätzung von (iii) ist bzgl. der Stabilitätskonstanten $c_{s}$ bestmöglich, d. h. auf $I = [0, \infty)$ ist nicht mit gleichmäßiger Konvergenz zu rechnen.

Beispiel: Als Beispiel betrachtet man das AWP $u^{'} (t) = a u (t)$ mit $u (0) = 1$ und $a > 0$ . Für die Lösung $u (t) = e^{a t}$ ergibt sich bei Anwendung des expliziten Euler-Verfahrens mit äquidistantem Gitter $u_{j} = (1 + a h)^{j - 1} = (1 + a h)^{t_{j} / h - 1}$ (mit $t_{j} = j h$ ), also $e_{h} (t_{j}) = e^{a t_{j}} - (1 + a h)^{t_{j} / h - 1}$ .

Satz (Konvergenz von Einzelschrittverfahren III):
Sei ein ESV mit Inkrementfunktion $ϕ \in C (I^{2} \times R^{n}, R^{n})$ gegeben. Außerdem existiere für $\bar{h}$ und $ε > 0$ fest eine Konstante $M > 0$ mit
$\forall_{k \in (0, \bar{h})} \forall_{t \in I} \forall_{w, \tilde{w} \in {v \in R^{n} | \exists_{t \in I} {∥ v - u (t) ∥}_{\infty} \leq ε}} {∥ ϕ (k, t, w) - ϕ (k, t, \tilde{w}) ∥}_{\infty} \leq M \cdot {∥ w - \tilde{w} ∥}_{\infty}$
(lokale Lipschitz-Bedingung an $ϕ$ im dritten Argument).
Dann gelten (i), (ii) und (iii) aus obigem Satz:

Das ESV ist stabil.
Ist das ESV konsistent, so ist es auch konvergent.
Ist das ESV konsistent zur Ordnung $p$ , so ist es auch konvergent zur Ordnung $p$ .

Explizite Runge-Kutta-Verfahren

Bemerkung: Seien $p \in N_{0}$ und ein Anfangswertproblem (AWP) mit einer Lösung
$u \in C^{p + 1} (I, R^{n})$ vorgegeben (dies ist z. B. der Fall für $f \in C^{p} (I, R^{n})$ ).
Kann man nun systematisch ein Einzelschrittverfahren mit Konsistenzordnung $p$ konstruieren?

Beispiel: Das Heun-Verfahren $u_{j + 1} = u_{j} + \frac{h_{j}}{2} (f (t_{j}, u_{j}) + f (t_{j} + h_{j}, u_{j} + h_{j} f (t_{j}, u_{j})))$ mit Inkrementfunktion $ϕ (k, t, w) = \frac{1}{2} (f (t, w) + f (t + k, w + k f (t, w)))$ erreicht durch iterative Auswertung von $f$ eine höhere Konsistenzordnung (nämlich $2$ ). Es gehört zu den einfachsten expliziten Runge-Kutta-Verfahren.

explizites Runge-Kutta-Verfahren: Seien $r \in N$ , $α_{2}, \dots, α_{r} \in R$ , $γ_{1}, \dots, γ_{r} \in R$ und $β_{i j}$ für $i = 2, \dots, r$ , $j = 1, \dots, r - 1$ und $i > j$ gegeben.
Das Einzelschrittverfahren (ESV) mit $ϕ (k, t, w) := \sum_{i = 1}^{r} γ_{i} K_{i} (k, t, w)$ und
$K_{1} (k, t, w) := f (t, w)$ ,
$K_{2} (k, t, w) := f (t + α_{2} k, w + k \cdot β_{21} K_{1} (k, t, w))$ ,
…,
$K_{r} (k, t, w) := f (t + α_{r} k, w + k \cdot \sum_{s = 1}^{r - 1} β_{r s} K_{s} (k, t, w))$
heißt allgemeines explizites Runge-Kutta-Verfahren der Stufe $r$ .

Butcher-Tableau: Die Koeffizienten eines allgemeinen Runge-Kutta-Verfahrens können in der Form einer Tabelle (Butcher-Tableau) zusammengefasst werden:

$\begin{aligned} \begin{array}{ccccc} α_{2} & β_{21} \\ ⋮ & ⋮ & ⋱ \\ α_{r} & β_{r 1} & \dots & β_{r, r - 1} \\ γ_{1} & \dots & γ_{r - 1} & γ_{r} \end{array} \end{aligned}$

Beispiel: Das explizite Euler-Verfahren $ϕ (k, t, w) = f (t, w)$ (Stufe $1$ ), das modifizierte Euler-Verfahren $ϕ (k, t, w) = f (t + \frac{k}{2}, w + \frac{k}{2} f (t, w))$ (Stufe $2$ ) und das Verfahren von Heun
$ϕ (k, t, w) = \frac{1}{2} (f (t, w) + f (t + k, w + k f (t, w)))$ (Stufe $2$ ) besitzen folgende Butcher-Tableaus:

$\begin{aligned} \begin{array}{cc} 1 \end{array} \begin{array}{cccc} \frac{1}{2} & \frac{1}{2} \\ 0 & 1 \end{array} \begin{array}{cccc} 1 & 1 \\ \frac{1}{2} & \frac{1}{2} \end{array} \end{aligned}$

Bemerkung: Setzt man $α_{1} := 0$ , so kann man die Koeffizientenfunktionen $K_{i}$ iterativ bestimmen durch die Formel $K_{i} (k, t, w) = f (t + α_{i} k, w + k \cdot \sum_{j = 1}^{i - 1} β_{i j} K_{j} (k, t, w))$ für $i = 1, \dots, r$ .

Beispiel: Für $K_{1} = f (t_{j}, u_{j})$ , $K_{2} = f (t_{j} + \frac{h_{j}}{2}, u_{j} + \frac{1}{2} h_{j} K_{1})$ , $K_{3} = f (t_{j} + \frac{h_{j}}{2}, u_{j} + \frac{1}{2} h_{j} K_{2})$ , $K_{4} = f (t_{j} + h_{j}, u_{j} + h_{j} K_{3})$ ergibt sich ein Runge-Kutta-Verfahren mit der Inkrementfunktion $u_{j + 1} = u_{j} + \frac{h_{j}}{6} (K_{1} + 2 K_{2} + 2 K_{3} + K_{4})$ . Es heißt klassisches Runge-Kutta-Verfahren und besitzt folgendes Butcher-Tableau:

$\begin{aligned} \begin{array}{ccccc} \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & 0 & \frac{1}{2} \\ 1 & 0 & 0 & 1 \\ \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array} \end{aligned}$

Bemerkung: Man kann sich fragen, wieviele Runge-Kutta-Verfahren $r$ -ter Ordnung auch eine Konsistenzordnung von $r$ besitzen.
Für den Fall $r = 2$ ergibt obiger Satz (Konsistenz von Einzelschrittverfahren) die Bedingungen $f (t, w) = ϕ (0, t, w)$ und $\frac{d}{d t} f (t, u (t)) = 2 \cdot \frac{\partial}{\partial k} ϕ (k, t, u (t)) |_{k = 0}$ .
Es gilt $ϕ (0, t, w) = γ_{1} f (t, w) + γ_{2} f (t, w)$ und $\frac{d}{d t} f (t, u (t)) = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial u} \cdot \frac{d u}{d t} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial u} \cdot f (t, u (t))$ .
Für die Ableitung von $ϕ$ gilt $ϕ (k, t, w) = γ_{1} f (t, w) + γ_{2} f (t + α_{2} k, w + k β_{21} f (t, w))$ , also
$2 \cdot {\frac{\partial ϕ}{\partial k} |}_{k = 0} = 2 γ_{2} (α_{2} \frac{\partial f}{\partial t} + β_{21} f (t, u (t)) \frac{\partial f}{\partial u})$ .
Aus Koeffizientenvergleich ergibt sich das nicht-lineare Gleichungssystem $1 = γ_{1} + γ_{2}$ , $2 γ_{2} α_{2} = 1$ , $2 γ_{2} β_{21} = 1$ . Für $γ_{2} \neq 0$ kann man die drei Gleichungen mit vier Unbekannten mit $γ_{2}$ als Parameter auflösen und erhält $γ_{1} = 1 - γ_{2}$ , $α_{2} = \frac{1}{2 γ_{2}}$ und $β_{21} = \frac{1}{2 γ_{2}}$ . Das Butcher-Tableau lautet

$\begin{aligned} \begin{array}{ccc} \frac{1}{2 γ_{2}} & \frac{1}{2 γ_{2}} \\ 1 - γ_{2} & γ_{2} . \end{array} \end{aligned}$ Für $γ_{2} = \frac{1}{2}$ erhält man das Heun-Verfahren und für $γ_{2} = 1$ das modifizierte Euler-Verfahren.

Bemerkung: Allgemein muss man ein nicht-lineares Gleichungssystem lösen. Die Konsistenzordnung eines $r$ -stufigen Runge-Kutta-Verfahrens ist nach oben durch $r$ beschränkt.
Leider gilt i. A. nicht, dass $r$ die maximal erreichbare Konsistenzordnung ist:

$\begin{aligned} \begin{array}{cccccccccc} r & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & \geq 9 \\ p_{max} (r) & 1 & 2 & 3 & 4 & 4 & 5 & 6 & 6 & \leq r - 2 \end{array} \end{aligned}$

Bemerkung: Ein Runge-Kutta-Verfahren der Stufe $r$ ist konsistent genau dann, wenn
$\sum_{i = 1}^{r} γ_{i} = 1$ gilt, denn aufgrund $K_{i} (0, t, w) = f (t, w)$ für $i = 1, \dots, r$ gilt
$ϕ (0, t, w) = \sum_{i = 1}^{r} γ_{i} K_{i} (0, t, w) = f (t, w) \cdot \sum_{i = 1}^{r} γ_{i} \overset{!}{=} f (t, w)$ .

Implizite Runge-Kutta-Verfahren

Bemerkung: Man spricht von einem impliziten Einzelschrittverfahren, falls die Inkrementfunktion $ϕ$ auch von $u_{i + 1} = u_{h} (t_{i + 1})$ abhängt, d. h. $u_{i + 1} = u_{i} + h_{i} ϕ (h_{i} t_{i}, u_{i}, u_{i + 1})$ .
Die Vorteile sind die verbesserte Stabilität und eine höhere mögliche Konsistenzordnung von bis zu $2 r$ . Der Nachteil ist natürlich der höhere numerische Aufwand, da man pro Zeitschritt ein in der Regel nicht-lineares Gleichungssystem lösen muss.

Beispiel: Das implizite Euler-Verfahren ist gegeben durch $u_{i + 1} = u_{i} + h_{i} f (t_{i + 1}, u_{i + 1})$ .

implizites Runge-Kutta-Verfahren:
Seien $r \in N$ , $α_{1}, \dots, α_{r} \in R$ , $γ_{1}, \dots, γ_{r} \in R$ und $b_{i j}$ für $i, j = 1, \dots, r$ gegeben.
Das Einzelschrittverfahren (ESV) mit $ϕ (k, t, w) := \sum_{i = 1}^{r} γ_{i} K_{i} (k, t, w)$ heißt allgemeines
implizites Runge-Kutta-Verfahren der Stufe $r$ , falls das nicht-lineare Gleichungssystem
$K_{1} (k, t, w) := f (t + α_{1} k, w + k \cdot \sum_{s = 1}^{r} b_{1 s} K_{s} (k, t, w))$ ,
…,
$K_{r} (k, t, w) := f (t + α_{r} k, w + k \cdot \sum_{s = 1}^{r} b_{r s} K_{s} (k, t, w))$ ,
erfüllt ist. Die Koeffizienten können analog zum expliziten Fall in einem Butcher-Tableau zusammengefasst werden:

$\begin{aligned} \begin{array}{ccccc} α_{1} & b_{11} & \dots & b_{1 r} \\ ⋮ & ⋮ & ⋮ \\ α_{r} & b_{r 1} & \dots & b_{r r} \\ γ_{1} & \dots & γ_{r} \end{array} \end{aligned}$

Beispiel: Beim impliziten Euler-Verfahren $u_{j + 1} = u_{j} + h_{j} f (t_{j + 1}, u_{j + 1})$ ist $K_{1} = f (t_{j + 1}, u_{j + 1})$ , d. h. $K_{1} = f (t_{j + 1}, u_{j} + h_{j} K_{1})$ . In der Regel ist pro Zeitschritt eine nicht-lineare Gleichung zu lösen. Man kann z. B. eine einfache Iteration $K_{1}^{(ℓ + 1)} = f (t_{j + 1}, u_{j} + h_{h} K_{1}^{(ℓ)})$ bzw.
$u_{j + 1}^{(ℓ + 1)} = u_{j} + h_{j} f (t_{j + 1}, u_{j + 1}^{(ℓ)})$ lösen oder das Newton-Verfahren anwenden.
Das Butcher-Tableau ist folgendes:

$\begin{aligned} \begin{array}{cc} 1 & 1 \\ 1 \end{array} \end{aligned}$

Beispiel: Allgemeine Runge-Kutta-Verfahren der Stufe $r = 1$ besitzen im skalaren Fall $n = 1$ die nicht-lineare Gleichung $K_{1} = f (t + α_{1} k, w + k b_{11} K_{1})$ . Man nimmt an, dass die Gleichung eindeutig lösbar ist mit $K_{1} \in C^{1} (I^{2} \times R, R)$ .
Für die Konsistenz muss $ϕ (0, t, w) = γ_{1} K_{1} (0, t, w) = f (t, w)$ gelten, das stimmt für $γ_{1} = 1$ (für $α_{1} = b_{11} = 1$ erhält man das implizite Euler-Verfahren). Differentiation von obiger Gleichung in $k = 0$ ergibt $\frac{\partial}{\partial k} ϕ (0, t, w) = \frac{\partial f}{\partial t} (t, w) α_{1} + \frac{\partial f}{\partial w} (t, w) b_{11} K_{1} (0, t, w)$ und
$\frac{d}{d t} f (t, u (t)) = \frac{\partial f}{\partial t} (t, u (t)) + \frac{\partial f}{\partial w} (t, u (t)) f (t, u (t))$ . Nach dem Konsistenzsatz muss für $p = 2$ gelten, dass $2 \frac{\partial}{\partial k} ϕ (0, t, u (t)) = \frac{d}{d t} f (t, u (t))$ , also $α_{1} = b_{11} = \frac{1}{2}$ .
Konkret erhält man also $u_{j + 1} = u_{j} + h_{j} K_{1} = u_{j} + h_{j} f (t_{j} + \frac{1}{2} h_{j}, u_{j} + \frac{1}{2} h_{j} K_{1})$
$= u_{j} + h_{j} f (\frac{1}{2} (t_{j} + t_{j + 1}, u_{j} + \frac{1}{2} (u_{j + 1} - u_{j})) = u_{j} + h_{j} f (\frac{1}{2} (t_{j} + t_{j + 1}), \frac{1}{2} (u_{j} + u_{j + 1}))$ ,
da $K_{1} = \frac{1}{h_{j}} (u_{j + 1} - u_{j})$ .

Beispiel: implizites Runge-Kutta-Verfahren der Stufe $r = 2$ und Ordnung $p = 4$ :

$\begin{aligned} \begin{array}{ccc} (3 - \sqrt{3}) / 6 & 1 / 4 & 1 / 4 - \sqrt{3} / 6 \\ (3 + \sqrt{3}) / 6 & 1 / 4 + \sqrt{3} / 6 & 1 / 4 \\ 1 / 2 & 1 / 2 \end{array} \end{aligned}$

Bemerkung: Jedes implizites ESV lässt sich für $| h |$ hinreichend klein als explizites Verfahren darstellen.

Beispiel: Wendet man das implizite Euler-Verfahren auf $u^{'} = a u$ an, so erhält man
$u_{i + 1} = u_{i} + a h u_{i + 1}$ , also $u_{i + 1} = \frac{1}{1 - a h} u_{i} = u_{i} + (\frac{1}{1 - a h} - 1) u_{i} = u_{i} + h \frac{a}{1 - a h} u_{i}$ .

Beispiel: Für die halbimpliziten Runge-Kutta-Verfahren gilt $b_{i s} = 0$ für $i < s$ , also
$K_{1} (k, t, w) := f (t + α_{1} k, w + k b_{11} K_{1} (k, t, w))$ ,
…,
$K_{r} (k, t, w) := f (t + α_{r} k, w + k \cdot \sum_{s = 1}^{r} b_{r s} K_{s} (k, t, w))$ ,
d. h. die einzelnen Gleichungen sind nacheinander lösbar. Das Butcher-Tableau hat dann folgende Form:

$\begin{aligned} \begin{array}{ccccc} α_{1} & b_{11} & 0 \\ ⋮ & ⋮ & ⋱ \\ α_{r} & b_{r 1} & \dots & b_{r r} \\ γ_{1} & \dots & γ_{r} \end{array} \end{aligned}$

Zusammenhang zwischen Runge-Kutta-Verfahren und Quadraturformeln

Bemerkung: Der Zusammenhang zwischen Runge-Kutta-Verfahren und Quadraturformeln gibt einen weiteren Weg zur systematischen Konstruktion von Runge-Kutta-Verfahren zu einer vorgegebenen Konsistenzordnung $p$ .

Bemerkung: Gegeben sei das allgemeine Runge-Kutta-Verfahren $u_{j + 1} = u_{j} + h_{j} \cdot \sum_{i = 1}^{r} γ_{i} K_{i}$ ,
$K_{i} = f (t_{j} + α_{i} h_{j}, u_{j} + h_{j} \cdot \sum_{s = 1}^{r} b_{i s} K_{s})$ .
Im Folgenden wird versucht, eine notwendige Bedingung für die Konsistenzordnung $p$ herzuleiten. Betrachtet man das Anfangswertproblem $u^{'} (t) = g (t)$ , $t \in I$ , $u (0) = u_{0} \in R^{n}$ mit $g : I \to R^{n}$ (Lösung $u (t) = u_{0} + \int_{t_{0}}^{t} g (τ) d τ$ ), so ergibt das Runge-Kutta-Verfahren
$\frac{1}{h_{j}} (u_{j + 1} - u_{j}) = \sum_{i = 1}^{r} γ_{i} g (t_{j} + α_{i} h_{j})$ . Der Konsistenzfehler ist laut Definition
$ε_{h} (t_{j}) = \frac{1}{h_{j}} (u (t_{j + 1}) - u (t_{j}) - h_{j} \sum_{i = 1}^{r} γ_{i} g (t_{j} + α_{i} h_{j})) = \frac{1}{h_{j}} (\int_{t_{j}}^{t_{j + 1}} g (t) d t - h_{j} \sum_{i = 1}^{r} γ_{i} g (t_{j} + α_{i} h_{j}))$ , d. h. um die Konsistenzordnung $p$ zu erreichen, muss
$| \int_{t_{j}}^{t_{j + 1}} g (t) d t - h_{j} \sum_{i = 1}^{r} γ_{i} g (t_{j} + α_{i} h_{j}) | = O (h_{j}^{p + 1})$ gelten.
Dies ist ein Quadraturproblem (Gewichte $γ_{i}$ , Stützstellen $t_{j} + α_{i} h_{j}$ ). Damit können die Koeffizienten $α_{1}, \dots, α_{r}$ und $γ_{1}, \dots, γ_{r}$ bestimmt werden.

Bemerkung: Das Einsetzen der exakten Lösung in das Runge-Kutta-Verfahren ergibt
$\frac{u (t_{j + 1}) - u (t_{j})}{h_{j}} \approx \sum_{i = 1}^{r} γ_{i} K_{i} (h_{j}, t_{j}, u (t_{j}))$ . Daraus folgt $\int_{t_{j}}^{t_{j + 1}} u^{'} (t) d t \approx h_{j} \sum_{i = 1}^{r} γ_{i} K_{i} (h_{j}, t_{j}, u (t_{j}))$ .
Wegen der Gleichung von eben sollte für schon gegebene $α_{i}$ und $γ_{i}$ gelten, dass $K_{i} \approx u^{'} (t_{j} + α_{i} h_{j})$ für $i = 1, \dots, r$ . Aus der Definition der $K_{i}$ folgt damit $K_{i} = f (t_{j} + α_{i} h_{j}, u_{j} + h_{j} \sum_{s = 1}^{r} b_{i s} K_{s}) \approx u^{'} (t_{j} + α_{i} h_{j}) = f (t_{j} + α_{i} h_{j}, u (t_{j} + α_{i} h_{j}))$ . Daraus folgt $u (t_{j} + α_{i} h_{j}) \approx u (t_{j}) + h (t_{j}) \sum_{s = 1}^{r} b_{i s} K_{s}$ , d. h. $\int_{t_{j}}^{t_{j} + α_{i} h_{j}} u^{'} (t) d t \approx \sum_{s = 1}^{r} b_{i s} u^{'} (t_{j} + α_{s} h_{j})$ für $i = 1, \dots, r$ . Somit erhält man ein Quadraturproblem, mit dem sich die $b_{i s}$ bestimmen lassen ( $i, s = 1, \dots, r$ ).
Dies motiviert den folgenden Satz.

Satz (Butcher, Kunzmann, 1969): Es sei ein Runge-Kutta-Verfahren mit $α_{i}, γ_{i} \in R$ für $i = 1, \dots, r$ und $b_{i j} \in R$ für $i, j = 1, \dots, r$ gegeben. Für $p, q \in N$ seien die Koeffizienten so gewählt, dass für alle $g_{1} \in C^{p + 1} (I, R^{n})$ und $g_{2} \in C^{q + 1} (I, R^{n})$ gilt

$| \frac{1}{h_{j}} \int_{t_{j}}^{t_{j + 1}} g_{1} (t) d t - \sum_{s = 1}^{r} γ_{s} g_{1} (t_{j} + α_{s} h_{j}) | = O (h_{j}^{p})$ für $j = 0, \dots, N - 1$ und
$| \frac{1}{h_{j}} \int_{t_{j}}^{t_{j} + α_{i} h_{j}} g_{2} (t) d t - \sum_{s = 1}^{r} β_{i s} g_{2} (t_{j} + α_{s} h_{j}) | = O (h_{j}^{q})$ für $j = 0, \dots, N - 1$ und $i = 1, \dots, r$ .

Dann ist das Runge-Kutta-Verfahren konsistent mit der Ordnung $min {p, q + 1}$ .

Exaktheit einer Quadraturformel: Es seien $g \in C ([0, 1], R)$ und $τ \in (0, 1]$ gegeben.
Sei $Q (g) := \sum_{i = 1}^{r} γ_{i} g (α_{i})$ eine Quadraturformel für das Integral $\int_{0}^{τ} g (t) d t$ , wobei $α_{i} \in [0, 1]$ und $γ_{i} \in R$ für $i = 1, \dots, r$ . $Q$ heißt vom Grad $ℓ$ exakt, falls $Q (p) - \int_{0}^{τ} p (t) d t = 0$ für alle $p \in P_{ℓ}$
( $P_{ℓ}$ Menge der Polynome vom Grad $\leq ℓ$ ).

Satz (Fehler einer Quadraturformel mit Peano-Kern): Seien $ℓ \in N$ , $g \in C^{ℓ + 1} ([0, 1], R)$ ,
$τ \in (0, 1]$ und $Q$ eine Quadraturformel, die vom Grad $ℓ$ exakt ist.
Dann gilt $\int_{0}^{τ} g (t) d t = Q (g) + \int_{0}^{1} π_{ℓ + 1} (t) g^{(ℓ + 1)} (t) d t$ , wobei $π_{ℓ + 1}$ der Peano-Kern
$π_{ℓ + 1} (t) := \frac{1}{(ℓ + 1)!} (((τ - t)_{+})^{ℓ + 1} - (ℓ + 1) \cdot \sum_{i = 1}^{r} γ_{i} ((α_{i} - t)_{+})^{ℓ})$ ist mit
$t \in [0, 1]$ und $α_{+} (t) = max {α (t), 0}$ .

Legendre-Polynom: Für $m \in N_{0}$ ist das Legendre-Polynom $p_{m}$ vom Grad $m$ gegeben durch $p_{m} (t) := \frac{m!}{(2 m)!} \cdot \frac{d^{m}}{d t^{m}} (t^{2} - 1)^{m}$ für $t \in R$ .

Beispiel: Es gilt $p_{0} (t) = 1$ , $p_{1} (t) = t$ , $p_{2} (t) = t^{2} - \frac{1}{3}$ usw.

Lemma (Nullstellen und Orthogonalität der Legendre-Polynome):

Das Legendre-Polynom $p_{m}$ besitzt paarweise verschiedene Nullstellen
$ϱ_{1}, \dots, ϱ_{m}$ mit $- 1 < ϱ_{1} < \dots < ϱ_{m} < 1$ .
Für $m, n \in N$ mit $m \neq n$ gilt $\int_{- 1}^{1} p_{m} (t) p_{n} (t) d t = 0$ .

Satz (Gauß-Quadratur): Seien $g \in C ([- 1, 1], R)$ und $Q (g) := \sum_{i = 1}^{m} ω_{i} g (ϱ_{i})$ die
Gauß-Quadraturformel mit den Stützstellen $ϱ_{i}$ (Nullstellen des Legendre-Polynoms $p_{m}$ ) und den Gewichten $ω_{i} := \int_{- 1}^{1} (\prod_{j = 1, j \neq i}^{m} \frac{t - ϱ_{j}}{ϱ_{i} - ϱ_{j}}) d t$ (Integrale für Lagrange-Polynome) für $m \in N$ .
Dann gilt $Q (p) = \int_{- 1}^{1} p (t) d t$ für alle $p \in P_{2 m - 1}$ ,
d. h. die Gauß-Quadratur ist exakt vom Grad $2 m - 1$ .

Bemerkung: Nun kann man analysieren, wie gut das Runge-Kutta-Verfahren ist, das durch die Gauß-Quadratur bestimmt wird. Dazu wendet man den Satz von Butcher und Kunzmann an.

Mit $t = t_{j} + h_{j} τ$ , $τ \in [0, 1]$ gilt $\frac{1}{h_{j}} \int_{t_{j}}^{t_{j + 1}} g_{1} (t) d t = \int_{0}^{1} g_{1} (t_{j} + h_{j} τ) d τ$
$= \sum_{i = 1}^{r} {\tilde{ω}}_{i} g_{1} (t_{j} + h_{j} {\tilde{ϱ}}_{i}) + \int_{0}^{1} π_{2 r} (τ) g_{1}^{(2 r)} (t_{j} + h_{j} τ) d τ$ , da die Gauß-Quadratur exakt vom Grad $2 r - 1$ ist. Daraus folgt $| \frac{1}{h_{j}} \int_{t_{j}}^{t_{j + 1}} g_{1} (t) d t - \sum_{i = 1}^{r} {\tilde{ω}}_{i} g_{1} (t_{j} + h_{j} {\tilde{ϱ}}_{i}) |$
$\leq max_{τ \in [0, 1]} | π_{2 r} (τ) | \cdot \int_{0}^{1} | g_{1}^{(2 r)} (t_{j} + h_{j} τ) | d τ \leq c | h |^{2 r}$ aufgrund der Beschränktheit von $π$ (bei jeder Ableitung von $g_{1}$ kommt ein Faktor $h_{j}$ hinzu). Dabei ist $γ_{i} := {\tilde{ω}}_{i} = \frac{ω_{i}}{2}$ und $α_{i} := {\tilde{ϱ}}_{i} = \frac{ϱ_{i} + 1}{2}$ für $i = 1, \dots, r$ , weil $\int_{0}^{1} g (τ) d τ = \frac{1}{2} \int_{- 1}^{1} g (\frac{z + 1}{2}) d z \approx \frac{1}{2} \sum_{i = 1}^{r} ω_{i} g (\frac{ϱ_{i} + 1}{2})$ .
Analog wie eben ist $\frac{1}{h_{j}} \int_{t_{j}}^{t_{j} + α_{i} h_{j}} g_{2} (t) d t = \int_{0}^{α_{i}} g_{2} (t_{j} + h_{j} τ) d τ$
$= \sum_{s = 1}^{r} {\hat{ω}}_{s} g_{2} (t_{j} + h_{j} {\hat{ϱ}}_{s}) + \int_{0}^{1} π_{2 r} (τ) g_{2}^{(2 r)} (t_{j} + h_{j} {\hat{ϱ}}_{s}) d τ$ . Daraus folgt wieder
$| \frac{1}{h_{j}} \int_{t_{j}}^{t_{j} + α_{i} h_{j}} g_{2} (t) d t - \sum_{s = 1}^{r} {\hat{ω}}_{s} g_{2} (t_{j} + h_{j} {\hat{ϱ}}_{s}) | \leq c | h |^{2 r}$ mit
$β_{i s} := {\hat{ω}}_{s} = \frac{α_{i} ω_{s}}{2}$ und ${\tilde{ϱ}}_{s} = α_{i} \frac{ϱ_{s} + 1}{2}$ wegen $\int_{0}^{α} g (τ) d τ = \frac{α}{2} \int_{- 1}^{1} g (α \frac{z + 1}{2}) d z$ .

Somit ergibt sich eine Konsistenzordnung von $p = min {2 r, 2 r + 1} = 2 r$ .

Mehrschrittverfahren

Definitionen und Beispiele

Bemerkung: Um die Genauigkeit von Einzelschrittverfahren zu erhöhen, verwendet man nicht nur die letzte, sondern die letzten $k$ Approximationen.

Mehrschrittverfahren: Seien $ψ \in C (I^{k + 2} \times R^{n (k + 1)}, R^{n})$ und $k \in N$ .
Weiter seien $a_{0}, \dots, a_{k} \in R$ und $u_{0} = u (t_{0}), u_{1}, \dots, u_{k - 1} \in R^{n}$ gegeben. Das Verfahren
$\frac{1}{h} (a_{0} u_{j} + a_{1} u_{j + 1} + \dots + a_{k} u_{j + k}) = ψ (h, t_{j}, \dots, t_{j + k}, u_{j}, \dots, u_{j + k})$ mit $j = 0, \dots, N - k$ heißt
$k$ -Mehrschrittverfahren ( $k$ -MSV) mit Verfahrensfunktion $ψ$ . (Das Gitter $I_{h}$ ist also äquidistant.)

Bemerkung: Falls $ψ$ nicht von $u_{j + k}$ abhängt und $a_{k} \neq 0$ gilt, so heißt das $k$ -MSV explizit.
Ein explizites $1$ -MSV ist ein explizites Einzelschrittverfahren.

Bemerkung: Um die Verfahrensgleichung lösen zu können, müssen zunächst die Startwerte $u_{0}, \dots, u_{k - 1}$ bekannt sein. Diese sollte mit einem ESV derselben Konsistenzordnung berechnet werden.

lineares MSV: Falls die Verfahrensfunktion $ψ$ von der Form $ψ (h, t_{j}, \dots, t_{j + k}, u_{j}, \dots, u_{j + k}) = \sum_{i = 0}^{k} b_{i} f (t_{j + i}, u_{j + i})$ mit $b_{0}, \dots, b_{k} \in R$ ist, so heißt das zugehörige $k$ -MSV linear.
Lineare MSV haben also die Form $\frac{1}{h} \sum_{i = 0}^{k} a_{i} u_{j + i} = \sum_{i = 0}^{k} b_{i} f (t_{j + i}, u_{j + i})$

Bemerkung: Auch MSV lassen sich durch Quadraturformeln herleiten. Die Integralgleichung $u (t) = u (s) + \int_{s}^{t} f (r, u (r)) d r$ mit $t > s$ und $t, s \in I$ ist äquivalent zu $u^{'} (t) = f (t, u (t))$ , speziell gilt $u (t_{j + k}) = u (t_{j + k - 1}) + \int_{t_{j + k - 1}}^{t_{j + k}} f (r, u (r)) d r = u (t_{j + k - 1}) + \int_{t_{j + k - 1}}^{t_{j + k}} u^{'} (r) d r$ .

Beispiel: Verwendet man die Trapezregel
$u (t_{j + 1}) = u (t_{j}) + \int_{t_{j}}^{t_{j + 1}} f (r, u (r)) d r \approx u (t_{j}) + \frac{1}{2} h (f (t_{j}, u (t_{j})) + f (t_{j + 1}, u (t_{j + 1})))$ (also $k = 1$ ), so ergibt sich das Trapezverfahren $u_{j + 1} := u_{j} + \frac{1}{2} h (f (t_{j}, u_{j}) + f (t_{j + 1}, u_{j + 1}))$ .

Bemerkung: Eine Idee für weitere Verfahren ist eine bessere Approximation der Integralgleichung durch Ersetzung des Integranden $u^{'} (r)$ durch ein Interpolationspolynom.
Die Interpolationspolynome $p_{j} \in P_{k - 1}$ , $j = 0, \dots, N - k$ sind eindeutig bestimmt durch die $k$ Bedingungen $p_{j} (t_{j + i}) := u^{'} (t_{j + i}) = f (t_{j + i}, u (t_{j + i}))$ für $i = 0, \dots, k - 1$ . Man erhält die veränderte Integralgleichung $u (t_{j + k}) \approx u (t_{j + k - 1}) + \int_{t_{j + k - 1}}^{t_{j + k}} p_{j} (r) d r$ .
Da man dafür allerdings die exakte Lösung $u$ benötigt, kann man auch ${\tilde{p}}_{j} \in P_{k - 1}$ verwenden, die analog definiert sind durch ${\tilde{p}}_{j} (t_{j + i}) := f (t_{j + i}, u_{j + i})$ für $i = 0, \dots, k - 1$ . Verwendet man ${\tilde{p}}_{j}$ statt $p_{j}$ in der Integralgleichung, so erhält man ein explizites $k$ -MSV
$\frac{1}{h} (u_{j + k} - u_{j + k - 1}) = ψ (h, t_{j}, \dots, t_{j + k}, u_{j}, \dots, u_{j + k - 1}) := \frac{1}{h} \int_{t_{j + k - 1}}^{t_{j + k}} {\tilde{p}}_{j} (r) d r$ .

Beispiel: Ein Beispiel für ein so erhaltenes lineares $4$ -Mehrschrittverfahren mit $p = 4$ ist
$\frac{1}{h} (u_{j + 4} - u_{j + 3}) = \frac{1}{24} (55 f (t_{j}, u_{j}) - 59 f (t_{j + 1}, u_{j + 1}) + 37 f (t_{j + 2}, u_{j + 2}) - 9 f (t_{j + 3}, u_{j + 3}))$ .
Es heißt Adams-Bashforth-Verfahren der Stufe $k = 4$ .

Bemerkung: Bei jedem Zeitschritt ist nur eine neue Auswertung von $f$ notwendig
(und zwar in $(t_{j + k - 1}, u_{j + k - 1})$ ).

Beispiel: Ein implizites Verfahren lässt sich analog konstruieren, nur bezieht man dabei
$t_{j + k}, u_{j + k}$ als Stützpunkte für die Interpolation ein.
Diese Verfahren heißen Adams-Moulton-Verfahren.
Ein Beispiel für $k = 4$ und $p = 5$ ist $\frac{1}{h} (u_{j + 4} - u_{j + 3})$
$= \frac{1}{720} (251 f (t_{j + 4}, u_{j + 4}) + 646 f (t_{j + 3}, u_{j + 3}) - 269 f (t_{j + 2}, u_{j + 2}) + 106 f (t_{j + 1}, u_{j + 1}) - 19 f (t_{j}, u_{j})$ .

Bemerkung: Bei den sog. Prädiktor-Korrektor-Verfahren kombiniert man implizite und explizite Verfahren. Seien also $\frac{1}{h} \sum_{i = 0}^{k} a_{i} u_{j + i} = ψ_{1} (h, t_{j}, \dots, t_{j + k}, u_{j}, \dots, u_{j + k - 1})$ ein explizites und $\frac{1}{h} \sum_{i = 0}^{k} α_{i} u_{j + i} = ψ_{2} (h, t_{j}, \dots, t_{j + k}, u_{j}, \dots, u_{j + k})$ ein implizites $k$ -MSV.
Man berechnen nun zuerst den Prädiktor $u_{j + k}^{(p)}$ mit dem expliziten MSV, d. h.
$\frac{1}{h} \sum_{i = 0}^{k - 1} a_{i} u_{j + i} + \frac{1}{h} a_{k} u_{j + k}^{(p)} = ψ_{1} (h, t_{j}, \dots, t_{j + k}, u_{j}, \dots, u_{j + k - 1})$ . Anschließend berechnet man $u_{j + k}$ mit dem impliziten Verfahren, also $\frac{1}{h} \sum_{i = 0}^{k} α_{i} u_{j + i} = ψ_{2} (h, t_{j}, \dots, t_{j + k}, u_{j}, \dots, u_{j + k - 1}, u_{j + k}^{(p)})$ .
Man muss also keine nicht-linearen Gleichungen lösen, sondern man verwendet den Prädiktor als Schätzwert für den wahren Wert $u_{j + k}$ .
Alternativ lässt sich der Prädiktor auch als Startwert für eine Fixpunktiteration verwenden, d. h. $u_{j + k}^{(0)} := u_{j + k}^{(p)}$ und $\frac{1}{h} \sum_{i = 0}^{k - 1} α_{i} u_{j + i} + \frac{1}{h} α_{k} u_{j + k}^{(m + 1)} := ψ_{2} (h, t_{j}, \dots, t_{j + k}, u_{j}, \dots, u_{j + k - 1}, u_{j + k}^{(m)})$ .

Konsistenz und Konvergenz von Mehrschrittverfahren

Fehler von linearen Mehrschrittverfahren:
Es sei $u_{h} : I_{h} \to R^{n}$ durch ein lineares $k$ -MSV gegeben.
$e_{h} := u |_{I_{h}} - u_{h}$ ist die globale Fehlerfunktion.
$\bar{e_{h}} := max_{j = 0, \dots, N} ∥ e_{h} (t_{j}) ∥$ ist der globale Diskretisierungsfehler.
$ε_{h} (t_{j + k}) := \frac{1}{h} \sum_{i = 0}^{k} a_{i} u (t_{j + i}) - \sum_{i = 0}^{k} b_{i} f (t_{j + i}, u (t_{j + i}))$ , $j = 0, \dots, N - k,$ ist die lokale Fehlerfunkt.
$\bar{ε_{h}} := max_{j = 0, \dots, N - k} ∥ ε_{h} (t_{j + k}) ∥$ ist der lokale Diskretisierungsfehler.

Bemerkung: Die Koeffizienten $a_{0}, \dots, a_{k}, b_{0}, \dots, b_{k}$ sollten so bestimmt werden, dass
$\bar{ε_{h}} = O (h^{p})$ . Dafür betrachtet man $ε_{h} (t_{j + k}) = \frac{1}{h} \sum_{i = 0}^{k} a_{i} u (t_{j} + i h) - \sum_{i = 0}^{k} b_{i} u^{'} (t_{j} + i h)$ und setzt für $p_{ℓ} (i) := \frac{1}{ℓ!} i^{ℓ}$ die Taylor-Entwicklungen $u (t_{j} + i h) = \sum_{ℓ = 0}^{p} h^{ℓ} p_{ℓ} (i) u^{(ℓ)} (t_{j}) + O (h^{p + 1})$ bzw.
$u^{'} (t_{j} + i h) = \sum_{ℓ = 0}^{p} h^{ℓ} p_{ℓ} (i) u^{(ℓ + 1)} (t_{j}) + O (h^{p + 1}) = \sum_{ℓ = 1}^{p} h^{ℓ - 1} p_{ℓ - 1} (i) u^{(ℓ)} (t_{j}) + O (h^{p})$
$= \sum_{ℓ = 1}^{p} h^{ℓ - 1} p_{ℓ}^{'} (i) u^{(ℓ)} (t_{j}) + O (h^{p})$ ein. Daraus folgt dann
$ε_{h} (t_{j + k}) = \frac{1}{h} \sum_{i = 0}^{k} a_{i} (\sum_{ℓ = 0}^{p} h^{ℓ} p_{ℓ} (i) u^{(ℓ)} (t_{j})) - \sum_{i = 0}^{k} b_{i} (\sum_{ℓ = 1}^{p} h^{ℓ - 1} p_{ℓ}^{'} (i) u^{(ℓ)} (t_{j})) + O (h^{p})$
$= \sum_{ℓ = 0}^{p} h^{ℓ - 1} u^{(ℓ)} (t_{j}) (\sum_{i = 0}^{k} a_{i} p_{ℓ} (i) - \sum_{i = 0}^{k} b_{i} p_{ℓ}^{'} (i)) + O (h^{p})$ .
Verschwindet der Ausdruck in Klammern, so hat das Verfahren die Konsistenzordnung $p$ . Das beweist folgenden Satz.

Satz (Konsistenz von MSV): Falls die Koeffizienten eines linearen $k$ -MSV
$a_{0}, \dots, a_{k}, b_{0}, \dots, b_{k} \in R$ die Bedingungen $\sum_{i = 0}^{k} a_{i} p_{ℓ} (i) = \sum_{i = 0}^{k} b_{i} p_{ℓ}^{'} (i)$ für $ℓ = 0, \dots, p$ erfüllen, so besitzt das MSV die Konsistenzordnung $p$ .
Dabei ist $p_{ℓ} (i) := \frac{1}{ℓ!} i^{ℓ}$ und $p_{ℓ}^{'} (i) := p_{ℓ - 1} (i) = \frac{1}{(ℓ - 1)!} i^{ℓ - 1}$ für $ℓ \geq 1$ bzw. $p_{0}^{'} (i) := 0$ .

Bemerkung: Diese Bedingungen entsprechen einem LGS mit $p + 1$ Gleichungen und $2 (k + 1)$ Unbekannten. Da die Lösung $a_{0} = \dots = a_{k} = b_{0} = \dots = b_{k} = 0$ keinen Sinn ergibt, ergänzt man manchmal die Normierungsbedingung $\sum_{i = 0}^{k} b_{i} = 1$ .
Damit das Gleichungssystem nicht überbestimmt ist, soll es höchstens so viele Gleichungen wie Variablen geben. Mit der Normierungsbedingung ist dann $p + 2 \leq 2 (k + 1)$ , d. h. die Konsistenzordnung $p$ ist durch $2 k$ nach oben beschränkt.
Bei expliziten Verfahren ist $b_{k} = 0$ , d. h. es gibt eine Variable weniger. Hier ist $p + 2 \leq 2 k + 1$ , also ist die Konsistenzordnung $p$ durch $2 k - 1$ nach oben beschränkt.

Beispiel: Für $k = 1$ soll $p = 2$ erreicht werden, d. h. die Gleichungen $a_{0} + a_{1} = 0$ , $a_{1} = b_{0} + b_{1}$ , $\frac{1}{2} a_{1} = b_{1}$ und $b_{0} + b_{1} = 1$ sollen erfüllt werden. Daraus folgt $a_{0} = - 1$ , $a_{1} = 1$ , $b_{0} = \frac{1}{2}$ und $b_{1} = \frac{1}{2}$ . Man erhält also die Trapezregel $\frac{1}{h} (- u_{j} + u_{j + 1}) = \frac{1}{2} f (t_{j}, u_{j}) + \frac{1}{2} f (t_{j + 1}, u_{j + 1})$ .

Stabilität von Mehrschrittverfahren

erzeugende Polynome:
Sei ein lineares $k$ -Mehrschrittverfahren $\frac{1}{h} \sum_{i = 0}^{k} a_{i} u_{j + i} = \sum_{i = 0}^{k} b_{i} f (t_{j + i}, u_{j + i})$ gegeben.
Dann heißen die Polynome $ϱ (z) := \sum_{i = 0}^{k} a_{i} z^{i}$ und $σ (z) := \sum_{i = 0}^{k} b_{i} z^{i}$ erzeugende Polynome des MSV ( $z \in C$ ).

alternative Schreibweise von linearen MSV: Sei $E$ der Vorwärts-Shift-Operator, d. h. $E y_{j} := y_{j + 1}$ . Dann lässt sich das lineare $k$ -Mehrschrittverfahren
$\frac{1}{h} \sum_{i = 0}^{k} a_{i} u_{j + i} = \sum_{i = 0}^{k} b_{i} f (t_{j + i}, u_{j + i})$ auch durch die erzeugenden Polynome in der Form
$\frac{1}{h} ϱ (E) u_{j} = σ (E) f_{j}$ mit $f_{j} = f (t_{j}, u_{j})$ schreiben, wobei $p (E) y_{j} := \sum_{i = 0}^{k} p_{i} y_{j + i}$ mit einem Polynom $p (z) = \sum_{i = 0}^{k} p_{i} z^{i}$ .

Bemerkung: Nicht jedes konsistente lineare MSV ist konvergent. Es wird eine zusätzliche Stabilitätsbedingung benötigt.

Beispiel: Ein Beispiel für ein instabiles lineares $2$ -MSV mit Konsistenzordnung $p = 3$ ist
$\frac{1}{h} (u_{i + 2} + 4 u_{i + 1} - 5 u_{i}) = 4 f_{i + 1} + 2 f_{i}$ . Die erzeugenden Polynome sind dabei $ϱ (z) = z^{2} + 4 z - 5$ und $σ (z) = 4 z - 2$ . Man wendet das MSV auf das triviale Anfangswertproblem $u^{'} = 0$ , $u (0) = 1$ (d. h. die Lösung ist $u (t) \equiv 1$ ) an.
Sei $u_{1} = 1 + ε h$ leicht gestört. Daraus ergibt sich die Drei-Term-Rekursion $u_{i + 2} + 4 u_{i + 1} - 5 u_{i} = 0$ (rechte Seite verschwindet wegen $f \equiv 0$ ) mit den Startwerten $u_{0} = 1$ und $u_{1} = 1 + ε h$ .
Für spezielle Lösungen betrachtet man die Nullstellen $z_{1} = 1$ und $z_{2} = - 5$ des erzeugenden Polynoms $ϱ (z)$ . Setzt man $u_{i} = z_{1}^{i}$ an, so ist $z_{1}^{i + 2} + 4 z_{1}^{i + 1} - 5 z_{1}^{i} = 0$ genau dann, wenn $z_{1}^{i} ϱ (z_{1}) = 0$ . Wegen $ϱ (z_{1}) = 0$ ist $u_{i} = z_{1}^{i}$ eine spezielle Lösung der Rekursion, analog $u_{i} = z_{2}^{i}$ .
Für die allgemeine Lösung setzt man $u_{i} = A z_{1}^{i} + B z_{2}^{i}$ an, also $u_{i + 2} + 4 u_{i + 1} - 5 u_{i} = 0$ genau dann, wenn $A z_{1}^{i} ϱ (z_{1}) + B z_{2}^{i} ϱ (z_{2}) = 0$ . Die Parameter $A$ und $B$ ergeben sich aus den Startbedingungen $1 = u_{0} = A z_{1}^{0} + B z_{2}^{0} = A + B$ und $1 + ε h = u_{1} = A z_{1}^{1} + B z_{2}^{1} = A - 5 B$ .
Daraus ergibt sich $A = 1 + \frac{ε h}{6}$ und $B = - \frac{ε h}{6}$ . Somit ist die allgemeine Lösung der Rekursion
$u_{i} = A z_{1}^{i} + B z_{2}^{i} = 1 + \frac{ε h}{6} - \frac{ε h}{6} \cdot (- 5)^{i}$ . Für den Fall $ε = 0$ kommt die exakte Lösung heraus. Ist allerdings $u_{1}$ leicht gestört ( $ε > 0$ ), so wird der Fehler durch den Faktor $(- 5)^{i}$ verstärkt, d. h. das MSV ist instabil.

Bemerkung: Die Vorgehensweise lässt sich auf allgemeine $k$ -MSV verallgemeinern. Durch Anwendung der Testgleichung $u^{'} = 0$ , $u (0) = 1$ erhält man die homogene Rekursion bzw. Differenzengleichung $a_{0} u_{j} + \dots + a_{k} u_{j + k} = 0$ für $j = 0, \dots, N - k$ mit Startwerten $u_{0}, \dots, u_{k - 1}$ .

Satz (Lösungen der homogenen Rekursion): Sei $λ \in C$ eine $m$ -fache Nullstelle des erzeugenden Polynoms $ϱ (z)$ , d. h. $ϱ (λ) = ϱ^{'} (λ) = \dots = ϱ^{(m - 1)} (λ) = 0$ . Dann gilt:

$u_{i}^{(1)} := λ^{i}$ , $u_{i}^{(2)} := i λ^{i - 1}$ , …, $u_{i}^{(m)} := D^{m - 1} λ^{i} = i (i - 1) \dots (i - m + 2) λ^{i - m + 1}$
sind spezielle Lösungen der homogenen Rekursion.
Die allgemeine Lösung der homogenen Rekursion ist eine Linearkombination der insgesamt $k$ speziellen Lösungen aus (i).
(Für jede Nullstelle $λ$ von $ϱ (z)$ erhält man entsprechend der Vielfachheit viele spezielle Lösungen, d. h. insgesamt $Grad (ϱ) = k$ viele Lösungen.)

Bemerkung: Sei $λ$ eine Nullstelle von $ϱ (z)$ . Dann gilt für $| λ | > 1$ , dass ${u_{i}} = {λ^{i}}$ exponentiell wächst, und für $| λ | < 1$ , dass ${u_{i}} = {λ^{i}}$ exponentiell fällt.
Für $| λ | = 1$ und Vielfachheit $ℓ$ von $λ$ ist $| u_{i}^{(1)} | = | λ |^{i} = 1$ und $u_{i}^{(ℓ)} = i (i - 1) \dots (i - ℓ + 2) λ^{i - ℓ + 1}$ wächst polynomial für $ℓ \geq 2$ .

stabil: Ein $k$ -Mehrschrittverfahren heißt stabil, falls alle Nullstellen des Polynoms $ϱ (z)$ im abgeschlossenen Einheitskreis liegen und diejenigen auf dem Rand nur einfach sind, d. h.
$ϱ (λ) = 0 \Rightarrow | λ | \leq 1$ und $(ϱ (λ) = 0 \land | λ | = 1) \Rightarrow ϱ^{'} (λ) \neq 0$ .
Wegen der Testgleichung $u^{'} = 0$ spricht man auch von Nullstabilität oder $D$ -Stabilität (nach Dahlquist).

stark/schwach stabil: Das $k$ -MSV heißt stark stabil, falls für alle Nullstellen außer $λ = 1$ gilt, dass $| λ | < 1$ . Ansonsten heißt das $k$ -MSV schwach stabil.

Bemerkung:
Bei konsistenten $k$ -MSV ist $λ = 1$ immer eine Nullstelle von $ϱ (z)$ , denn $ϱ (1) = \sum_{i = 0}^{k} a_{i} = 0$ .
Die Adams-Verfahren (Adams-Bashforth und Adams-Moulton) sind stark stabil, denn hier ist $a_{k} = 1$ , $a_{k - 1} = - 1$ und $a_{k - 2} = \dots = a_{0} = 0$ , d. h. $ϱ (z) = z^{k} - z^{k - 1} = z^{k - 1} \cdot (z - 1)$ .
$λ = 1$ ist einfache Nullstelle, während $λ = 0$ eine $(k - 1)$ -fache Nullstelle ist.

Satz (Dahlquist-Barriere – maximale Konvergenzordnung stabiler linearer MSV):
Ein lineares $k$ -Mehrschrittverfahren $\frac{1}{h} \sum_{i = 0}^{k} a_{i} u_{j + i} = \sum_{i = 0}^{k} b_{i} f (t_{j + i}, u_{j + i})$ , das obige Stabilitätsbedingung erfüllt, hat maximal die Konvergenzordnung

$k + 2$ für $k$ gerade,
$k + 1$ für $k$ ungerade und
$k$ für $\frac{b_{k}}{a_{k}} \leq 0$ (insbesondere für explizite Verfahren).

Die Ordnung $k + 2$ kann nur erzielt werden, wenn alle Nullstellen von $ϱ (z)$ auf dem Rand des Einheitskreises liegen.

Beispiel: Für $k = 1$ wird die maximale Konvergenzordnung $p = 2$ von der Trapezformel erreicht. Für $k = 2$ wird die maximale Konvergenzordnung $p = 4$ vom Milne-Simpson-Verfahren $u_{i + 1} = u_{i - 1} + \frac{h}{3} (f_{i - 1} + 4 f_{i} + f_{i + 1})$ erreicht (schwach stabil).
Das Adams-Bashforth-Verfahren ( $k = 4$ ) ist explizit und erreicht daher nur die Konvergenzordnung $p = 4$ . Das Adams-Moulton-Verfahren ist stark stabil (kann nicht Ordnung $k + 2$ erreichen) und erreicht die Konvergenzordnung $5 = k + 1$ .

Satz (Konvergenz von MSV): Falls ein lineares MSV die Konsistenzordnung $p$ hat und obige Stabilitätsbedingung erfüllt, so ist es auch konvergent mit Ordnung $p$ .

Adaptive Schrittweitensteuerung

Bemerkung: Sei ein Einzelschrittverfahren zur Lösung des Anfangswertproblems (AWP) gegeben. Für ein gegebenes Gitter $I_{h}$ sei $T (I_{h})$ der numerische Aufwand zur Lösung des ESV auf $I_{h}$ (Rechenzeit). Außerdem sei $TOL$ eine gegebene Fehlertoleranz.
Die Aufgabe ist nun, ein Gitter $I_{h}^{opt}$ zu finden mit $\bar{e_{h}} \leq TOL$ und $T (I_{h}^{opt}) \leq T (I_{h})$ für alle Gitter $I_{h}$ mit $\bar{e_{h}} \leq TOL$ .
Man weiß nicht, ob $I_{h}^{opt}$ überhaupt existiert oder ob es eindeutig ist. Durch die sog. adaptive Schrittweitensteuerung versucht man, eine möglichst gute Approximation von $I_{h}^{opt}$ zu finden.

Satz (Fehlerentwicklung): Seien $u \in C^{p + 2} (I, R^{n})$ eine Lösung von (AWP) und $I_{h}$ ein Gitter. Außerdem sei ein stabiles ESV mit Inkrementfunktion $ϕ \in C^{p + 1} (I^{2} \times R^{n}, R^{n})$ , $u_{h} (0) = u_{0}$ und Konsistenzordnung $p$ gegeben.
Dann existiert eine Funktion $e_{0} \in C^{2} (I, R^{n})$ mit $e_{0} (0) = 0$ und
${∥ u_{h} - (u - h^{p} e_{0}) |_{I_{h}} ∥}_{\infty} = O (h^{p + 1})$ .
Es gibt zusätzlich eine Funktion $e_{1} \in C^{3} (I, R^{n})$ mit $e_{1} (0) = 0$ und
${∥ u_{h} - (u - h^{p} e_{0} - h^{p + 1} e_{1}) |_{I_{h}} ∥}_{\infty} = O (h^{p + 2})$ .

Bemerkung: Führt man für ein festes $t \in I_{h}$ das Hilfsproblem $v^{'} (s) = f (s, v (s))$ für $s > t$ und $v (t) = u_{h} (t)$ ein, so gilt $u_{h} (t + h) - v (t + h) = h^{p} e_{0} (t + h) + h^{p + 1} e_{1} (t + h) + O (h^{p + 2}) = h^{p} (e_{0} (t) + h e_{0}^{'} (t)) + h^{p + 1} (e_{1} (t) + h e_{1}^{'} (t)) + O (h^{p + 2}) = h^{p + 1} e_{0}^{'} (t) + O (h^{p + 2})$ wegen $e_{0} (t) = e_{1} (t) = 0$ . Analog gilt $u_{h / 2} (t + h) - v (t + h) = {(\frac{h}{2})}^{p} h e_{0}^{'} (t) + O (h^{p + 2})$ .
Somit ist $u_{h} (t + h) - u_{h / 2} (t + h) = (h^{p + 1} - {(\frac{h}{2})}^{p} h) e_{0}^{'} (t) + O (h^{p + 2})$ ,
d. h. $h e_{0}^{'} (t) = \frac{1}{2^{p} - 1} {(\frac{h}{2})}^{- p} (u_{h} (t + h) - u_{h / 2} (t + h)) + O (h^{2})$ .
Man erhält also für den Fehler der halben Gitterweite die Formel
$u_{h / 2} (t + h) - v (t + h) = Δ_{h} + O (h^{p + 2})$ mit $Δ_{h} := \frac{1}{2^{p} - 1} (u_{h} (t + h) - u_{h / 2} (t + h))$ einem Fehlerschätzer, der nur aus berechenbaren Größen besteht. Man definiert nun den relativen Fehlerschätzer ${\tilde{Δ}}_{h} := \frac{Δ h}{max {1, ∥ u_{h} ∥}}$ und kann daraus einen selbstadaptiven Algorithmus erstellen.

selbstadaptiver Algorithmus mit $(h, h / 2)$ -Gittersteuerung:
Startschrittweite $h_{0} \in [0, T]$ , minimale und maximale Schrittweite $h_{min} < h_{0} < h_{max}$ ,
Fehlertoleranz $TOL > 0$ , Verkleinerungs-/Vergrößerungsfaktoren $k_{min} < 1$ und $k_{max} > 1$ ,
Verfahrensordnung $p \in N$

$\begin{aligned} t := 0; u := u (0); h := h_{0}; \\ w h i l e (t < T) { \\ | {\tilde{Δ}}_{h} | := TOL + 1; \\ w h i l e (| {\tilde{Δ}}_{h} | > TOL) { \\ v := u + h \cdot ϕ (h, t, u); \\ z := u + h / 2 \cdot ϕ (h / 2, t, u); \\ w := z + h / 2 \cdot ϕ (h / 2, t + h / 2, z); \\ | {\tilde{Δ}}_{h} | := \frac{1}{2^{p} - 1} \cdot \frac{| v - w |}{max {1, u}}; \\ h := max {h_{min}, k_{min} \cdot h}; \\ i f h = h_{min} {r e t u r n;} \\ } \\ u := w; h := min {h_{max}, k_{max} \cdot h}; t := t + h; \\ } \end{aligned}$