11th June 2024, 58 min read

Konvergenzresultate für feste Schrittweiten

Original post is here eklausmeier.goip.de/blog/2024/06-11-konvergenzresultate-fuer-feste-schrittweiten.

1. Einführung und grundlegende Begriffe
2. Die Lemmata von Gronwall
3. Notation und Darstellungssatz für Differenzengleichungen
4. Stabilitätsfunktionale für feste Schrittweiten
5. Projektorstabilitätsfunktionale
6. Nichtäquidistante Gitter
7. Die Eigenwerte gewisser tridiagonaler Matrizen
8. Verfahren für parabolische Gleichungen

Es folgt ein recht allgemeiner Konvergenzbeweis für mehrstufige Verfahren, wobei allerdings vorausgesetzt wird, daß mit fester Schrittweite gearbeitet wird. Innerhalb des mehrstufigen Prozesses braucht das verwendete Gitter nicht äquidistant zu sein, wie z.B. bei Runge-Kutta Verfahren. Dabei wird allerdings hier ein etwas längerer Weg eingeschlagen. Zuerst werden in breiter Form Stabilitätsfunktionale vorgestellt und verschiedene, gleichwertige und äquivalente Darstellungen angegeben. Die Beweise für diese Stabilitätsfunktionale enthalten die eigentlichen Konvergenzbeweise, jedoch sind Stabilitätsfunktionale allgemeiner. Sie liefern direkt Stabilitätsungleichungen für die Differenzen zweier Lösungen von Differenzengleichungen, d.h. die Stabilitätsfunktionale liefern direkt Aussagen über das Auseinanderlaufen der Lösungen zweier Differenzengleichungen in Abhängigkeit von Störungen. An Differenzengleichungen werden nur lineare Gleichungen betrachtet, allerdings darf die Inhomogenität beliebig sein, wenn sie nur einer Lipschitzbedingung genügt.

Bevor die eigentlichen Überlegungen bzgl. der Stabilitätsfunktionale angestellt werden, sollen anhand einfacher, vorangestellter Überlegungen, einige grundsätzliche Probleme beleuchtet werden. Danach folgen die sehr wichtigen Aussagen von Gronwall. Das diskrete Lemma von Gronwall spielt eine entscheidende Rolle beim Hauptsatz über Stabilitätsfunktionale. Vielerorts befindet sich das diskrete Lemma von Gronwall versteckt in Konvergenzbeweisen und hier i.d.R. nur in sehr spezialisierter Form. Erst daran anschliessend werden die Stabilitätsfunktionale behandelt und verschiedene Äquivalenzen bewiesen.

1. Einführung und grundlegende Begriffe

Es sei $B$ ein Banachraum und $h \in R$ die Schrittweite. Die Klasse von Verfahren der Form

u_{n + 1} = S u_{n} + h φ (u_{n - 1}), n = 0, 1, \dots, N, u_{i} \in B,

berücksichtigt nicht block-implizite Verfahren, oder überhaupt implizite Verfahren, zumindestens ersteinmal nicht in sofort offenkundiger Weise. Hierbei ist

N := | \frac{b - a}{h} |, also | h | = \frac{| b - a |}{N} .

Subsumiert sind also nicht Verfahren der Vorschrift

A_{1} u_{n + 1} + A_{0} u_{n} = h \cdot (B_{1} F_{n + 1} + B_{0} F_{n}), n = 0, 1, \dots, N .

bzgl. der rein formalen Schreibweise.

Verfahren der leicht allgemeineren Form

\begin{matrix} (*) & u_{n + 1} = S u_{n} + h φ (u_{n + 1}, u_{n}), n = 0, 1, \dots, N, \end{matrix}

berücksichtigen blockimplizite Verfahren und gewöhnliche implizite Verfahren. Die oben angeschriebene Rekurrenz-Vorschrift für $u_{n + 1}$ stellt eine implizite Gleichung für $u_{n + 1}$ dar. An $φ$ muß man daher gewisse Voraussetzungen stellen, um eindeutige Lösbarkeit der impliziten Differenzengleichung zu garantieren. Da die zu integrierende Funktion fast durchweg einer Lipschitzkonstanten genügt, ist es naheliegend dasselbe auch für die Inhomogenität der Differenzengleichung zu fordern. Es möge also gelten

\begin{aligned} | φ ({\hat{u}}_{n + 1}, {\hat{u}}_{n}) - φ (u_{n + 1}, u_{n}) | & \leq K_{1} | {\hat{u}}_{n + 1} - u_{n + 1} |, \\ | φ ({\hat{u}}_{n + 1}, {\hat{u}}_{n}) - φ (u_{n + 1}, u_{n}) | & \leq K_{2} | {\hat{u}}_{n} - u_{n} | . \end{aligned}

1. Satz: Die Differenzengleichung $(*)$ besitzt, für genügend kleine Schrittweiten $| h |$ , eine eindeutige Lösung. Diese eindeutig bestimmte Lösung lässt sich mit Picard-Iteration bestimmen.

Beweis: Die Keplersche Gleichung für $u_{n + 1}$ hat wegen der vorausgesetzten Lipschitz-Stetigkeit in der letzten Komponente von $φ$ , nach dem Banachschen Fixpunktsatz eine eindeutig bestimmte Lösung und lässt sich durch Fxpunktiteration gewinnen.

| h K_{2} | < 1, für geeigntes h .

☐

2. Bemerkung: Für nicht genügend kleines $| h |$ kann die Gleichung in der Tat mehrere oder keine Lösung besitzen.

Es seien betrachtet die beiden Verfahren der Form

\begin{aligned} {\hat{u}}_{n + ℓ} + A_{ℓ - 1} {\hat{u}}_{n + ℓ - 1} + \dots + A_{0} {\hat{u}}_{n} & = h φ ({\hat{u}}_{n + ℓ}, {\hat{u}}_{n + ℓ - 1}, \dots, {\hat{u}}_{n}) + r_{n}, \\ u_{n + ℓ} + A_{ℓ - 1} u_{n + ℓ - 1} + \dots + A_{0} u_{n} & = h φ (u_{n + ℓ}, u_{n + ℓ - 1}, \dots, u_{n}) . \end{aligned}

Das erste Verfahren kann man als gestörtes Verfahren auffassen, während hingegen das zweite Verfahren das eigentliche Verfahren zur Berechnung der numerischen Lösung ist. Es sei

P_{1} := (\begin{matrix} I & 0 & \dots & 0 \end{matrix}), R_{1} := (\begin{matrix} 0 \\ ⋮ \\ 0 \\ I \end{matrix}),

und $δ_{n + ℓ} := {\hat{u}}_{n + ℓ} - u_{n + ℓ}$ und dazu

{\hat{δ}}_{n + ℓ} := φ ({\hat{u}}_{n + ℓ}, \dots, {\hat{u}}_{n}) - φ (u_{n + ℓ}, \dots, u_{n}) .

Es ist also ${\hat{δ}}_{n + ℓ}$ die Differenz der entsprechenden Werte für die Funktion $φ (\cdot)$ , wenn sämtliche Argumente verschieden sind.

Weiter sei

U := (\begin{matrix} u_{0} \\ ⋮ \\ u_{N} \end{matrix}), \hat{U} := (\begin{matrix} {\hat{u}}_{0} \\ ⋮ \\ {\hat{u}}_{N} \end{matrix}), R := (\begin{matrix} r_{0} \\ ⋮ \\ r_{N + ℓ - 1} \end{matrix}) .

Die einzelnen $u_{i}$ und ${\hat{u}}_{i}$ sind aus dem Vektorraum $R$ , nicht notwendig endlichdimensional. Hierbei sind die $A_{i} : B \to B$ stetige, lineare Operatoren zwischen Banach-Räumen. Bei linearen Operatoren ist bekanntlich die Stetigkeit in einem Punkte äquivalent mit der globalen Stetigkeit und dies äquivalent mit der Beschränktheit. $B$ ist hierbei entweder ein reller oder komplexer Banachraum. Die Vektorraumeigenschaften braucht man der Linearität wegen, die Normiertheit für die folgenden Funktionale, und die Vollständigkeit wird benötigt bei der Anwendung des Banachschen Fixpunktsatzes. Beispiele sind $B = C^{k}$ und $B = R^{k}$ , mit $k \geq 1$ . Gelegentlich gelten die Sätze auch in nicht notwendigerweise kommutativen Ringen $B$ . Für $B$ wird im folgenden stets $C^{k}$ gewählt. Die Mengen $C^{k ℓ \times k ℓ}$ wären dann entsprechend zu ersetzen durch $R^{ℓ \times ℓ}$ und andere Mengen entsprechend.

Man beachte, daß die Abschätzung nun abhängig von $h$ ist, die Abschätzung aber ausschliesslich für gewisse sehr stark eingeschränkte Schrittweiten $h$ gilt. Ohne Einschränkung an die Schrittweite $h$ ist der Satz nicht richtig. Die Unabhängigkeit von den $B_{ν}$ , bei

\sum_{ν = 0}^{ℓ} A_{ν} u_{n + ν} = h \sum_{ν = 0}^{ℓ} B_{ν} F_{n + ν}, n = 0, 1, \dots, N,

verlangt eine entsprechende Einschränkung an die Schrittweite $h$ . Für einen praktischen Einsatz ist zusätzlich ein entsprechend großer Stabilitätsbereich erforderlich. In die Größe des Stabilitätsbereiches gehen entscheidend die $B_{ν}$ ein und die Art der Iteration, mit der die impliziten Gleichungen in jedem Zeitschritt gelöst werden. Der Satz verliert ebenfalls seine Gültigkeit bei “langen” Integrationsintervallen. $| b - a |$ wird dann beliebig groß. Der Satz zeigt, daß bei endlichem Integrationsintervall $| b - a |$ , die $| \hat{U} - U |$ -Norm mit der $| {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R |$ -Norm äquivalent ist. Bei unendlich langen Integrationsintervallen, sind diese Normen nicht notwendigerweise mehr äquivalent.

2. Die Lemmata von Gronwall

Das Lemma von Gronwall, Thomas Gronwall, (1877--1932) für den kontinuierlichen Falle lautet

1. Satz: (Lemma von Gronwall) Seien $h$ , $w$ und $k$ stetige, reell-wertige Funktionen auf dem Intervall $[a, b]$ . (Es muß lediglich gelten $(\int_{a}^{x} f)^{'} = f (x)$ , sodaß man mit leicht schwächeren Bedingungen auskäme.) Es gelte auf diesem Intervall die Abschätzung

h (x) \leq w (x) + \int_{a}^{x} k (t) h (t) d t, \forall x \in [a, b] .

Das Integral auf der rechten Ungleichungsseite sei stets nicht-negativ, was beispielsweise für nicht-negative Funktionen $k$ , $w$ , $h$ auf dem Intervall $[a, b]$ , sichergestellt werden kann. Dann gilt die Abschätzung für die Funktion $h$ auf dem gesamten Intervall zu

h (x) \leq w (x) + \int_{a}^{x} \exp (\int_{t}^{x} k (τ) d τ) k (t) w (t) d t, \forall x \in [a, b] .

Beweis: siehe Helmut Werner und Herbert Arndt in Werner/Arndt (1986). Sei

H (x) := \int_{a}^{x} k (t) h (t) d t .

Hiermit gilt dann aufgrund der Stetigkeit von $k$ und $h$

H^{'} (x) = k (x) h (x), H (a) = 0, \forall x \in [a, b] .

Aus dieser Differentialgleichung folgt aufgrund der vorausgesetzten Ungleichung für die Funktion $h$

H^{'} (x) = k (x) h (x) \leq k (x) \cdot (w (x) + H (x)),

also die lineare Differentialungleichung

\begin{matrix} (*) & H^{'} (x) - k (x) H (x) \leq k (x) w (x), H (a) = 0. \end{matrix}

Multiplikation mit

e^{- K (x)} > 0, K (x) := \int_{a}^{x} k (t) d t,

führt zu

e^{- K (x)} [H^{'} (x) - k (x) H (x)] = {(e^{- K (x)} \cdot H (x))}^{'} \overset{\binom{(*)}{↓}}{\leq} e^{- K (x)} \cdot k (x) w (x) .

Integration von $a$ nach $x$ liefert wegen der mittleren Gleichung (Integral ist monotones Funktional)

e^{- K (x)} H (x) - e^{- K (a)} H (a) \leq \int_{a}^{x} e^{- K (t)} \cdot k (t) w (t) d t,

also wegen $H (a) = 0$ somit

H (x) \leq \int_{a}^{x} e^{K (x) - K (t)} \cdot k (t) w (t) d t

und aufgrund der Voraussetzung von $h (x) \leq w (x) + H (x)$ sofort

h (x) \leq w (x) + \int_{a}^{x} (\exp \int_{t}^{x} k (τ) d τ) \cdot k (t) w (t) d t .

☐

2. Folgerung: Gilt $h (x) \leq w + k \int_{a}^{x} h (t) d t$ , mit festen, nicht-negativen Konstanten $w$ und $k$ , so folgt die Abschätzung

h (x) \leq w e^{k \cdot (x - a)}, \forall x \in [a, b] .

Ein völliges Analogon zum kontinuierlichen Lemma von Gronwall macht das diskrete Lemma von Gronwall, welches ebenfalls exponentielles Wachstum anzeigt, wenn eine Funktion geeignet auf beiden Seiten einer Ungleichung vorkommt. Es gilt nun der

3. Satz: (Diskretes Lemma von Gronwall) Es seien $(m + 1)$ positive Zahlen $0 \leq η_{0} \leq η_{1} \leq \dots \leq η_{m}$ vorgegeben. Ferner sei $δ \geq 0$ , $h_{j} \geq 0$ und $x_{j + 1} = x_{j} + h_{j}$ . Es gelte die Ungleichung

ε_{0} \leq η_{0} und ε_{j + 1} \leq η_{j} + δ \sum_{ν = 0}^{j} h_{ν} ε_{ν}, j = 0, \dots, m - 1.

Dann gilt

ε_{j} \leq η_{j} e^{δ \cdot (x_{j} - x_{0})}, j = 0, \dots, m .

Beweis: siehe erneut Helmut Werner und Herbert Arndt in Werner/Arndt (1986). Der Fall $δ = 0$ ist einfach, wegen $e^{0} = 1$ . Sei nun $δ > 0$ . Induktionsverankerung mit $j = 0$ ist klar, ebenfalls einfach, wegen $e^{0} = 1$ . Der eigentliche Beweis reduziert sich jetzt lediglich noch auf den Induktionsschluß von $j$ nach $j + 1$ , wobei $δ > 0$ vorausgesetzt werden kann. Hier gilt nun

\begin{aligned} ε_{j + 1} & \leq η_{j + 1} + δ \sum_{ν = 0}^{j} h_{ν} ε_{ν} \\ \leq η_{j + 1} + δ \sum_{ν = 0}^{j} h_{ν} η_{ν} e^{δ \cdot (x_{ν} - x_{0})} \\ \leq η_{j + 1} \cdot (1 + δ \sum_{ν = 0}^{j} h_{ν} e^{δ \cdot (x_{ν} - x_{0})}) \\ \leq η_{j + 1} e^{δ \cdot (x_{j + 1} - x_{0})} . \end{aligned}

Für die Summe in der Klammer schätzte man ab (Untersumme einer streng monoton steigenden Funktion)

\sum_{ν = 0}^{j} h_{ν} e^{δ \cdot (x_{ν} - x_{0})} \leq \int_{x_{0}}^{x_{j + 1}} e^{δ \cdot (t - x_{0})} d t = \frac{1}{δ} (e^{δ (x_{j + 1} - x_{0})} - 1) .

☐

3. Notation und Darstellungssatz für Differenzengleichungen

Man vgl. auch Matrixpolynome.

Zur multiplen Lipschitzkonstanten von $φ$ :

| φ (u_{ℓ}, \dots, {\hat{u}}_{i}, \dots, u_{0}) - φ (u_{ℓ}, \dots, u_{i}, \dots, u_{0}) | \leq K_{i} \cdot | {\hat{u}}_{i} - u_{i} |, i = 0, \dots, ℓ .

$(ℓ + 1)$ -malige Anwendung der Dreiecksungleichung liefert

| φ ({\hat{u}}_{ℓ}, \dots, {\hat{u}}_{0}) - φ (u_{ℓ}, \dots, u_{0}) | \leq \sum_{i = 0}^{ℓ} K_{i} | {\hat{u}}_{i} - u_{i} | = ⟨ (\begin{matrix} K_{0} \\ ⋮ \\ K_{ℓ} \end{matrix}), (\begin{matrix} | {\hat{u}}_{0} - u_{0} | \\ ⋮ \\ | {\hat{u}}_{ℓ} - u_{ℓ} | \end{matrix}) ⟩

In der Schreibweise von $φ$ seien fortan zahlreiche hier nicht weiter interessierende Argumente der Schreibvereinfachung und der Klarheit wegen weggelassen. Es ist $φ (u_{ℓ}, \dots, u_{1}) = φ (t_{ℓ}, h_{ℓ}, u_{ℓ}, \dots, u_{1})$ .

1. Beispiel: Für die Verfahrensvorschrift der Form

A_{ℓ} u_{n + ℓ} + \dots + A_{0} u_{n} = h (B_{ℓ} F_{n + ℓ} + \dots + B_{0} F_{n}), F_{n + i} := (\begin{matrix} f_{N ℓ + *} \\ ⋮ \\ f_{N ℓ + *} \end{matrix}),

wobei die $f_{k}$ Näherungswerte für $f (t_{k}, y (t_{k}))$ sind. Die Funktion $f$ der Differentialgleichung sei Lipschitz-stetig mit der Lipschitzkonstanten $L$ vorausgesetzt, also

| f (t, \hat{y}) - f (t, y) | \leq L | \hat{y} - y | .

Dann gilt für die obigen Lipschitzkonstanten $K_{i}$ die Verbindung mit der Lipschitzkonstanten der Differentialgleichung zu

K_{i} = | B_{i} | \cdot L, oder ggf. K_{i} = | A_{ℓ}^{- 1} B_{i} | \cdot L .

2. Definition: Es sei $T$ eine gänzlich beliebige Matrix der Größe $k \times k$ . Dann wird der Bidiagonaloperator $[T]$ zur Matrix $T$ der Größe $(N + 1) k \times (N + 1) k$ , wie folgt definiert

[T] := (\begin{matrix} I \\ - T & I \\ ⋱ & ⋱ \\ - T & I \end{matrix}), {[T]}^{- 1} = (\begin{matrix} I \\ T & I \\ ⋮ & ⋮ & ⋱ \\ T^{n} & T^{n - 1} & \dots & I \end{matrix}) .

Rechts daneben steht die Inverse, welche für eine beliebige Matrix $T$ stets existiert.

Die speziellen Operatoren $[\cdot]$ und $[\cdot]^{- 1}$ tauchen im weiteren wiederholt auf. Aufgrund der Häufigkeit, wäre es zweckmässiger, die Rollen von $[\cdot]$ und $[\cdot]^{- 1}$ zu vertauschen, jedoch stände dies dann im Gegensatz zur Schreibweise bei Skeel (1976), Robert David Skeel.

3. Satz: (Eigenschaften von $col$ , $row$ , $diag$ , $[\cdot]$ ) Es gilt

$col A_{ν} B_{ν} = diag A_{ν} col B_{ν}$ .
$col A_{ν} B = (col A_{ν}) B$ ; Rechtsdistributivität des $col$ -Operators.
$row A_{ν} B_{ν} = row A_{ν} diag B_{ν}$ .
$row A B_{ν} = A row B_{ν}$ ; Linksdistributivität des $row$ -Operators.
$diag A_{ν} B_{ν} = diag A_{ν} diag B_{ν}$ ; multiplikative Distributivität des Bidiagonaloperators.
$[S^{- 1} T S] = diag S^{- 1} [T] diag S$ .
${[S^{- 1} T S]}^{- 1} = diag S^{- 1} {[T]}^{- 1} diag S$ .

Beweis: Zu (1):

{col}_{ν = 0}^{n} A_{ν} B_{ν} = (\begin{matrix} A_{0} B_{0} \\ ⋮ \\ A_{n} B_{n} \end{matrix}) = (\begin{matrix} A_{0} \\ ⋱ \\ A_{n} \end{matrix}) (\begin{matrix} B_{0} \\ ⋮ \\ B_{n} \end{matrix}) .

Zu (3):

{row}_{ν = 0}^{n} A_{ν} B_{ν} = (A_{0} B_{0} \dots A_{n} B_{n}) = (A_{0} \dots A_{n}) (\begin{matrix} B_{0} \\ ⋱ \\ B_{n} \end{matrix}) .

Zu (5)

{diag}_{ν = 0}^{N} A_{ν} B_{ν} = (\begin{matrix} A_{0} B_{0} \\ ⋱ \\ A_{n} B_{n} \end{matrix}) = (\begin{matrix} A_{0} \\ ⋱ \\ A_{n} \end{matrix}) (\begin{matrix} B_{0} \\ ⋱ \\ B_{n} \end{matrix}) .

Zu (6): Beachte die Definition von $[T]$ und benutze dann

(\begin{matrix} A_{11} & \dots & A_{1 n} \\ ⋮ & ⋱ & ⋮ \\ A_{m 1} & \dots & A_{m n} \end{matrix}) (\begin{matrix} S_{1} \\ ⋱ \\ S_{n} \end{matrix}) = (\begin{matrix} A_{11} S_{1} & \dots & A_{1 n} S_{n} \\ ⋮ & ⋱ & ⋮ \\ A_{m 1} S_{1} & \dots & A_{m n} S_{n} \end{matrix}),

bzw.

(\begin{matrix} S_{1} \\ ⋱ \\ S_{m} \end{matrix}) (\begin{matrix} A_{11} & \dots & A_{1 n} \\ ⋮ & ⋱ & ⋮ \\ A_{m 1} & \dots & A_{m n} \end{matrix}) = (\begin{matrix} S_{1} A_{11} & \dots & S_{1} A_{1 n} \\ ⋮ & ⋱ & ⋮ \\ S_{m} A_{m 1} & \dots & S_{m} A_{m n} \end{matrix}) .

Zu (7): Folgt aus (4), wegen $(A B)^{- 1} = B^{- 1} A^{- 1}$ , wobei $[T]$ für gänzlich beliebige Matrizen $T$ invertierbar ist. Für $T = 0$ ist $[T]$ die Einheitsmatrix der Größe $(n + 1) k \times (n + 1) k$ . ☐

4. Beispiele: Es gilt

{c o l}_{i = 0}^{ℓ - 1} (X T^{i}) = ({diag}_{i = 0}^{ℓ - 1} X) ({c o l}_{i = 0}^{ℓ - 1} T^{i}) und {r o w}_{i = 0}^{ℓ - 1} (T^{i} Y) = ({r o w}_{i = 0}^{ℓ - 1} T^{i}) ({diag}_{i = 0}^{ℓ - 1} Y) .

Im allgemeinen gilt

\underset{i \in U}{diag} \underset{k \in V}{diag} A_{k} \neq \underset{k \in V}{diag} \underset{i \in U}{diag} A_{i} .

Als nächstes folgt die Darstellung der Differenz der Lösung zweier Differenzengleichungen. Dieser Satz spielt eine wiederholt wichtige Rolle bei den gleich folgenden Hauptsätzen.

5. Satz: (Darstellungssatz) Voraussetzungen: ${\hat{u}}_{n}$ und $u_{n}$ seien die Lösungen der beiden Differenzengleichungen

\begin{aligned} {\hat{u}}_{n + ℓ} + A_{ℓ - 1} {\hat{u}}_{n + ℓ - 1} + \dots + A_{0} {\hat{u}}_{n} & = h φ ({\hat{u}}_{n + ℓ}, \dots, {\hat{u}}_{n}) + r_{n + ℓ} \\ u_{n + ℓ} + A_{ℓ - 1} u_{n + ℓ - 1} + \dots + A_{0} u_{n} & = h φ (u_{n + ℓ}, \dots, u_{n}) \end{aligned}} n = 0, 1, \dots, N .

Die “Störungen” $r_{n + ℓ}$ korrespondieren zum Wert ${\hat{u}}_{n + ℓ}$ . Es seien zur Abkürzung gesetzt

\begin{aligned} δ_{n + ℓ} & := {\hat{u}}_{n + ℓ} - u_{n + ℓ}, \\ {\hat{δ}}_{n + ℓ} & := φ ({\hat{u}}_{n + ℓ}, \dots, {\hat{u}}_{n}) - φ (u_{n + ℓ}, \dots, u_{n}) \end{aligned}} n = 0, \dots, N .

Die Differenzengleichung für ${\hat{u}}_{n}$ habe die Startwerte ${\hat{u}}_{i} := u_{i} + r_{i}$ , für $i = 0, \dots, ℓ - 1$ . Es sei $δ_{i} := r_{i}$ , für $i = 0, \dots, ℓ - 1$ , und $r_{ν} := δ_{ν} := {\hat{δ}}_{ν} := 0$ , für $ν > N$ .

Behauptung:

\begin{aligned} δ_{n} & = P_{1} C_{1}^{n} (\begin{matrix} δ_{0} \\ ⋮ \\ δ_{ℓ - 1} \end{matrix}) + P_{1} \sum_{ν = 0}^{n - ℓ} C_{1}^{n - 1 - ν} R_{1} (r_{ν + ℓ} + h {\hat{δ}}_{ν + ℓ}) \\ = P_{1} C_{1}^{n} (\begin{matrix} δ_{0} \\ ⋮ \\ δ_{ℓ - 1} \end{matrix}) + P_{1} \sum_{ν = 0}^{n - ℓ} C_{1}^{n - 1 - ν} R_{1} r_{ν + ℓ} + P_{1} \sum_{ν = 0}^{n - ℓ} C_{1}^{n - 1 - ν} R_{1} {\hat{δ}}_{ν + ℓ} \end{aligned}

Beweis: Folgt aus dem allgemeinen Satz über die Lösung inhomogener, linearer Matrix-Differenzengleichungen. Die allgemeine Lösung der Differenzengleichung

x_{n + ℓ} + A_{ℓ - 1} x_{n + ℓ - 1} + \dots + A_{0} x_{n} = y_{n}, n = 0, 1, \dots, N

lautet

x_{n} = P_{1} C_{1}^{n} z_{0} + P_{1} \sum_{ν = 0}^{n - 1} C_{1}^{n - 1 - ν} R_{1} y_{ν} .

Mit den obigen Abkürzungen für $δ_{n}$ und ${\hat{δ}}_{n}$ , ergibt sich eine Differenzengleichung für $δ_{n}$ zu

δ_{n + ℓ} + A_{ℓ - 1} δ_{n + ℓ - 1} + \dots + A_{0} δ_{n} = h {\hat{δ}}_{n} + r_{n} .

Diese Gleichung hat die Lösungsdarstellung

δ_{n} = P_{1} C_{1}^{n} (\begin{matrix} δ_{0} \\ ⋮ \\ δ_{ℓ - 1} \end{matrix}) + P_{1} \sum_{ν = 0}^{n - 1} C_{1}^{n - 1 - ν} R_{1} (h {\hat{δ}}_{ν + ℓ} + r_{ν + ℓ}) .

In der ersten Summe verschwinden die letzten $(ℓ - 1)$ Terme, wegen $P_{1} C_{1}^{i} R_{1} = 0$ , für $i = 0, \dots, ℓ - 2$ . Daß hier $P_{1} C_{1}^{ℓ - 1} R_{1} = I$ , braucht man noch nicht. Für $ν > n - ℓ$ ist $n - 1 - ν \leq ℓ - 2$ . Daher folgt genau die behauptete Gleichung, wie oben angegeben. ☐

Die folgende Ungleichung liefert nicht die bestmögliche Abschätzung für $ℓ \geq 2$ , jedoch bleibt sie einfach zu handhaben und wird nachher beim Beweis des Hauptsatzes benötigt.

6. Hilfssatz: (Abschätzungssatz) Voraussetzung: $φ (\cdot)$ sei in jeder Komponente Lipschitz-stetig mit den Lipschitzkonstanten $K_{i}$ . Die Werte $δ_{ν + ℓ}$ und ${\hat{δ}}_{ν + ℓ}$ seien wie oben definiert.

Behauptung:

\begin{aligned} \sum_{ν = 0}^{n - ℓ} | {\hat{δ}}_{ν + ℓ} | & \leq K_{ℓ} | δ_{n} | + (\sum_{i = 0}^{ℓ} K_{i}) (\sum_{ν = 0}^{n - 1} | δ_{ν} |) \\ \leq (\sum_{i = 0}^{ℓ} K_{i}) (\sum_{ν = 0}^{n} | δ_{ν} |) \\ \leq (ℓ + 1) \cdot ({max}_{i = 0}^{ℓ} K_{i}) \sum_{ν = 0}^{n} | δ_{ν} | . \end{aligned}

Beweis: Für $ν = 0, \dots, n - 1$ ist

\begin{aligned} | {\hat{δ}}_{ν + ℓ} | & = | φ ({\hat{u}}_{ν + ℓ}, \dots, {\hat{u}}_{ν}) - φ (u_{ν + ℓ}, \dots, u_{ν}) | \\ \leq K_{0} | δ_{ν} | + K_{1} | δ_{ν + 1} | + \dots + K_{ℓ} | δ_{ν + ℓ} | . \end{aligned}

Sei jetzt, in einer nicht zu Mißverständnissen führenden Doppelbezeichnung, zur Schreibvereinfachung gesetzt $δ_{ν} \leftarrow | δ_{ν} |$ und ${\hat{δ}}_{ν} \leftarrow | {\hat{δ}}_{ν} |$ , d.h. die Betragsstriche werden einfach weggelassen. Nun ist hiermit

\begin{aligned} \sum_{ν = 0}^{n - ℓ} {\hat{δ}}_{ν} & \leq (K_{0} δ_{0} + \dots + K_{ℓ} δ_{ℓ}) + (K_{0} δ_{1} + \dots + K_{ℓ} δ_{ℓ + 1}) + \dots + (K_{0} δ_{n - ℓ + 1} + \dots + K_{ℓ} δ_{n}) \\ = K_{0} (δ_{0} + \dots + δ_{n - ℓ + 1}) \\ + K_{1} (δ_{1} + \dots + δ_{n - ℓ + 2}) \\ + \dots \\ + K_{ℓ} (δ_{ℓ} + \dots + δ_{n}) \end{aligned}

Summation und Abschätzung bzgl. der Spalten im obigen Schema zeigt sofort die erste Abschätzung, wenn man die allerletzte Spalte mit $K_{ℓ}$ und $δ_{n}$ gesondert behandelt. Die weiteren behaupteten Ungleichungen ergeben sich sofort aus der ersten. ☐

Zur handlichen Notation der im folgenden Hauptsatz auftauchenden Stabilitätsfunktionale seien die folgenden abkürzenden Bezeichnungen eingeführt. Es war

P_{1} := (\begin{matrix} I & 0 & \dots & 0 \end{matrix}) \in C^{k \times k ℓ}, R_{1} := (\begin{matrix} 0 \\ ⋮ \\ 0 \\ I \end{matrix}) \in C^{k ℓ \times k} .

und die erste Begleitmatrix lautet $(I = I_{k \times k})$

C_{1} := (\begin{matrix} 0 & I & 0 & \dots & 0 \\ 0 & 0 & I & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ \dots & I \\ - A_{0} & - A_{1} & \dots & - A_{ℓ - 1} \end{matrix}) \in C^{k ℓ \times k ℓ} .

Desweiteren sei

{\overset{―}{P}}_{1} := {diag}_{ν = 0}^{N} P_{1} = (\begin{matrix} I & 0 & \dots & 0 \\ I & 0 & \dots & 0 \\ ⋱ \\ I & 0 & \dots & 0 \end{matrix}) \in C^{(N + 1) k \times (N + 1) k ℓ},

und

{\overset{―}{R}}_{1} := diag (I_{k ℓ \times k ℓ}, {diag}_{ν = 1}^{N} R_{1}) = (\begin{matrix} I_{k ℓ \times k ℓ} \\ 0 \\ ⋮ \\ 0 \\ I \\ ⋱ \\ 0 \\ ⋮ \\ 0 \\ I \end{matrix}) \in C^{(N + 1) k ℓ \times (N + ℓ) k},

wobei

R := {col}_{ν = 0}^{N + ℓ - 1} r_{ν} = (\begin{matrix} r_{0} \\ ⋮ \\ r_{N + ℓ - 1} \end{matrix}) \in C^{(N + ℓ) k} .

Für das Produkt gilt: $[C_{1}]^{- 1} {\overset{―}{R}}_{1} \in C^{(N + 1) k ℓ \times (N + ℓ) k}$ . Es sei $(X, T, Y)$ ein beliebiges Standard-Tripel. Weiter sei

\overset{―}{X} := {diag}_{ν = 0}^{N} X = (\begin{matrix} X \\ X \\ ⋱ \\ X \end{matrix})

und

\overset{―}{Y} := diag [{({c o l}_{i = 0}^{ℓ - 1} X T^{i})}^{- 1}, {diag}_{ν = 1}^{N} Y] = (\begin{matrix} {({col}_{i = 0}^{ℓ - 1} X T^{i})}^{- 1} \\ Y \\ ⋱ \\ Y \end{matrix}) .

Es ist, aufgrund der Biorthogonalitätsbeziehung,

{({col}_{i = 0}^{ℓ - 1} X T^{i})}^{- 1} = ({row}_{i = 0}^{ℓ - 1} T^{i} Y) B,

mit der Block-Hankel-Matrix $B$ zu

B = (\begin{matrix} A_{1} & \dots & A_{ℓ} \\ ⋮ & ⋰ \\ A_{ℓ} \end{matrix}), A_{ℓ} = I .

Die Sonderbehandlung der Blockmatrix bei ${\overset{―}{R}}_{1}$ und $\overset{―}{Y}$ in dem ersten “Diagonalelement” hat seinen Ursprung in der Lösungsdarstellung einer Differenzengleichung für Matrixpolynome der Form

x_{n} = X J^{n} ({col}_{i = 0}^{ℓ - 1} X J^{i}) (\begin{matrix} y_{0} \\ ⋮ \\ y_{ℓ - 1} \end{matrix}) + X \sum_{ν = 0}^{n - 1} J^{n - 1 - ν} Y y_{ν + ℓ} .

Für den Fall $ℓ = 1$ , also $ρ (μ) = I μ - A$ reduzieren sich $P_{1}$ und $R_{1}$ zu Einheitsmatrizen der Größe $n \times n$ . Die Biorthogonalitätsbeziehung schrumpft zu $X = Y^{- 1}$ oder $X^{- 1} = Y$ .

4. Stabilitätsfunktionale für feste Schrittweiten

Man vgl. Peter Albrecht, "Die numerische Behandlung gewöhnlicher Differentialgleichungen: Eine Einführung unter besonderer Berücksichtigung zyklischer Verfahren", 1979. Sowie Peter Albrecht, 1985.

Zuerst sei zur Übersichtlichkeit ein Teil des Beweises des nachfolgenden Hauptsatzes nach vorne gezogen. Später wird dieser Hilfssatz erweitert. Es gibt noch weitere Äquivalenzen zwischen Stabilitätsfunktionalen.

1. Hilfssatz: Voraussetzung: Es sei $C_{1}^{i} := 0 \in C^{k ℓ \times k ℓ}$ , falls $i < 0$ .

Behauptung: Das verkürzte Stabilitätsfunktional ist mit dem ursprünglichen es erzeugenden Stabilitätsfunktional normmässig äquivalent, d.h. es gilt

| [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | \sim | {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | .

Beweis: Der Beweis wird in zwei Teile aufgespalten. Man schätzt beide Stabilitätsfunktionale gegeneinander ab. Die Abschätzung $| {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | \leq | {\overset{―}{P}}_{1} | | [C_{1}]^{- 1} {\overset{―}{R}}_{1} R |$ ist klar, wobei die Zeilensummennorm von ${\overset{―}{P}}_{1}$ unabhängig von $N$ ist. Die andere Abschätzungsrichtung berücksichtigt das Verhalten der Begleitmatrix $C_{1}$ intensiver. Man vergleiche hier auch die beiden nachstehenden Beispiele zur Verdeutlichung des “quasi-nilpotenten” Charakters der Potenzen der Matrizen $C_{1}$ . Man benutzt

| C_{1}^{n} z_{0} | \leq | C_{1}^{ℓ - 1} | | C_{1}^{n - ℓ + 1} z_{0} | = | C_{1}^{ℓ - 1} | {max}_{i = 0}^{ℓ - 1} | P_{1} C_{1}^{n + i - ℓ + 1} z_{0} |,

wegen

| C_{1}^{n} z_{0} | = {max}_{i = 0}^{ℓ - 1} | P_{1} C_{1}^{n + i} z_{0} | .

Diese letzte Identität hat ihre Wurzel in der eben genannten “quasi-nilpotenten” Eigenschaft der Begleitmatrix $C_{1}$ . Das Herausziehen von $C_{1}^{ℓ - 1}$ ist zulässig, da bei der $sup$ -Norm bei $| {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R |$ weiterhin über alle Zeilen das Supremun gebildet wird. Es geht kein Wert bei der Supremunsbildung verloren. Schließlich

| C_{1}^{n - 1 - ν} R_{1} r_{ν + ℓ} | \leq | C_{1}^{ℓ - 1} | | C_{1}^{n - ℓ - ν} R_{1} r_{ν + ℓ} | = | C_{1}^{ℓ - 1} | {max}_{i = 0}^{ℓ - 1} | P_{1} C_{1}^{n - ℓ - ν + i} R_{1} r_{ν + ℓ} |,

wegen $r_{ν} := 0$ , für $ν > N$ . ☐

2. Beispiel: Sei $ℓ = 2$ und sei $N := n := 3$ . Es ist $ρ (μ) = I μ^{2} + A_{1} μ + A_{0} \in C^{k \times k}$ und die Potenzen der ersten Begleitmatrix $C_{1}$ lauten $C_{1}^{ν}$ , für $ν = 1, \dots, N$ :

C_{1} = (\begin{matrix} 0 & I \\ - A_{0} & - A_{1} \end{matrix}), C_{1}^{2} = (\begin{matrix} - A_{0} & - A_{1} \\ A_{1} A_{0} & - A_{0} + A_{1}^{2} \end{matrix}), C_{1}^{3} = (\begin{matrix} A_{1} A_{0} & - A_{0} + A_{1}^{2} \\ A_{0}^{2} - A_{1}^{2} A_{0} & A_{0} A_{1} + A_{1} A_{0} - A_{1}^{3} \end{matrix}) .

Es war

P_{1} = (\begin{matrix} I & 0 \end{matrix}) \in C^{k \times 2 k}, R_{1} = (\begin{matrix} 0 \\ I \end{matrix}) \in C^{2 k \times k} .

Die Matrizen ${\overset{―}{P}}_{1}$ und ${\overset{―}{R}}_{1}$ haben das Aussehen

{\overset{―}{P}}_{1} = (\begin{matrix} I & 0 \\ I & 0 \\ I & 0 \end{matrix}) \in C^{3 k \times 6 k}, {\overset{―}{R}}_{1} = (\begin{matrix} I \\ I \\ 0 \\ I \\ 0 \\ I \end{matrix}) \in C^{6 k \times 4 k} .

Man berechnet

[C_{1}]^{- 1} {\overset{―}{R}}_{1} = (\begin{matrix} I \\ C_{1} & R_{1} \\ C_{1}^{2} & C_{1} R_{1} & R_{1} \\ C_{1}^{3} & C_{1}^{2} R_{1} & C_{1} R_{1} & R_{1} \end{matrix}) \in C^{8 k \times 5 k}

(\begin{matrix} \begin{matrix} I & 0 \\ 0 & I \end{matrix} \\ \begin{matrix} 0 & I \\ - A_{0} & - A_{1} \end{matrix} & \begin{matrix} 0 \\ I \end{matrix} \\ \begin{matrix} - A_{0} & - A_{1} \\ A_{1} A_{0} & - A_{0} + A_{1} \end{matrix} & \begin{matrix} I \\ - A_{1} \end{matrix} & \begin{matrix} 0 \\ I \end{matrix} \\ \begin{matrix} A_{1} A_{0} & - A_{0} + A_{1}^{2} \\ A_{0}^{2} - A_{1}^{2} A_{0} & A_{0} A_{1} + A_{1} A_{0} - A_{1}^{3} \end{matrix} & \begin{matrix} - A_{1} \\ - A_{0} + A_{1}^{2} \end{matrix} & \begin{matrix} I \\ - A_{1} \end{matrix} & \begin{matrix} 0 \\ I \end{matrix} \end{matrix})

An einer weiteren Demonstration ersieht man das sehr schnelle “Großwerden” der überstrichenen Matrizen.

3. Beispiel: Es sei nun $ℓ = 3$ und $N = 3$ . Es ist $ρ (μ) = I μ^{3} + A_{2} μ^{2} + A_{1} μ + A_{0} \in C^{k \times k}$ . Nun berechnet man $C_{1}^{ν}$ für $ν = 1, \dots, N$ :

C_{1} = (\begin{matrix} 0 & I & 0 \\ 0 & 0 & I \\ - A_{0} & - A_{1} & - A_{2} \end{matrix}), C_{1}^{2} = (\begin{matrix} 0 & 0 & I \\ - A_{0} & - A_{1} & - A_{2} \\ A_{2} A_{0} & - A_{0} + A_{2} A_{1} & - A_{1} + A_{2}^{2} \end{matrix})

und

C_{1}^{3} = (\begin{matrix} - A_{0} & - A_{1} & - A_{3} \\ A_{2} A_{0} & - A_{0} + A_{2} A_{1} & - A_{1} + A_{2}^{2} \\ - A_{2} A_{0}^{2} & A_{2} A_{0} + A_{1}^{2} - A_{2} A - 1 & - A_{0} + A_{1} A_{2} + A_{2} A_{1} - A_{2}^{3} \end{matrix}) .

Weiter ist ${\overset{―}{P}}_{1} \in C^{4 k \times 4 k ℓ}$ und ${\overset{―}{R}}_{1} \in C^{4 k ℓ \times 6 k}$ mit

{\overset{―}{P}}_{1} = (\begin{matrix} I & 0 & 0 \\ I & 0 & 0 \\ I & 0 & 0 \\ I & 0 & 0 \end{matrix}), {\overset{―}{R}}_{1} = (\begin{matrix} I \\ I \\ I \\ 0 \\ 0 \\ I \\ 0 \\ 0 \\ I \\ 0 \\ 0 \\ I \end{matrix}) .

Nun ist $[C_{1}]^{- 1} {\overset{―}{R}}_{1} \in C^{12 k \times 6 k}$ mit

[C_{1}]^{- 1} {\overset{―}{R}}_{1} = (\begin{matrix} I \\ C_{1} & R_{1} \\ C_{1}^{2} & C_{1} R_{1} & R_{1} \\ C_{1}^{3} & C_{1}^{2} R_{1} & C_{1} R_{1} & R_{1} \end{matrix}) \in C^{4 k ℓ \times 6 k},

also

(\begin{matrix} \begin{matrix} I \\ I \\ I \end{matrix} \\ \begin{matrix} 0 & I & 0 \\ 0 & 0 & I \\ - A_{0} & - A_{1} & - A_{2} \end{matrix} & \begin{matrix} 0 \\ 0 \\ I \end{matrix} \\ \begin{matrix} 0 & 0 & I \\ - A_{0} & - A_{1} & - A_{2} \\ A_{2} A_{0} & - A_{0} + A_{2} A_{1} & - A_{1} + A_{2}^{2} \end{matrix} & \begin{matrix} 0 \\ I \\ - A_{2} \end{matrix} & \begin{matrix} 0 \\ 0 \\ I \end{matrix} \\ \begin{matrix} - A_{0} & - A_{1} & - A_{3} \\ A_{2} A_{0} & - A_{0} + A_{2} A_{1} & - A_{1} + A_{2}^{2} \\ - A_{2} A_{0}^{2} & A_{2} A_{0} + A_{1}^{2} - A_{2} A - 1 & - A_{0} + A_{1} A_{2} + A_{2} A_{1} - A_{2}^{3} \end{matrix} & \begin{matrix} I \\ - A_{2} \\ - A_{1} + A_{2}^{2} \end{matrix} & \begin{matrix} 0 \\ I \\ - A_{2} \end{matrix} & \begin{matrix} 0 \\ 0 \\ I \end{matrix} \end{matrix})

Das zugrunde liegende Schema ist hier

\begin{matrix} \begin{matrix} 1 & 1 & 1 \\ 2 & 2 & 2 \\ 3 & 3 & 3 \end{matrix} \\ \begin{matrix} 2 & 2 & 2 \\ 3 & 3 & 3 \\ 4 & 4 & 4 \end{matrix} \\ \begin{matrix} ⋮ & ⋮ & ⋮ \end{matrix} \end{matrix}

Es folgt nun der angekündigte Hauptsatz, aus dem sich leicht ein entsprechender Konvergenzsatz für sehr allgemeine Diskretisierungsverfahren ableiten lässt. Desweiteren zeigt der Satz mehrere Querverbindungen zwischen verschiedenen Stabilitätsfunktionalen auf. In gewissen Situationen hat jedes der vorkommenden Funktionale seine spezifischen Vor- und Nachteile, und es lohnt sich mehrere Darstellungen, oder äquivalente Funktionale zur Verfügung zu haben. Insbesondere sollte jede der Darstellungen in gegenseitiger Befruchtung gepflegt werden. Später werden noch zwei andere Darstellungen hinzukommen, die bei gewissen Untersuchungen abermals vereinfachend wirken.

4. Hauptsatz: Voraussetzungen: $(P_{1}, C_{1}, R_{1})$ sei das erste Begleiter-Tripel zum Matrixpolynom

ρ (μ) := I μ^{ℓ} + A_{ℓ - 1} μ^{ℓ - 1} + \dots + A_{1} μ + A_{0},

vom Grade $ℓ \geq 1$ . Die Funktion $φ$ sei Lipschitz-stetig in jeder Komponente mit den Lipschitzkonstanten $K_{i}$ , also

| φ (u_{ℓ}, \dots, {\hat{u}}_{i}, \dots, u_{0}) - φ (u_{ℓ}, \dots, u_{i}, \dots, u_{0}) | \leq K_{i} \cdot | {\hat{u}}_{i} - u_{i} |, für i = 0, \dots, ℓ .

Die Potenzen der Matrix $C_{1}$ seien beschränkt durch die obere Schranke $D$ , also $| C_{1}^{ν} | \leq D$ , $\forall ν \in N$ . Seien $ξ$ und $\hat{ξ}$ definiert durch

ξ := | P_{1} | D | R_{1} | K_{ℓ}, \hat{ξ} := | P_{1} | D | R_{1} | (\sum_{i = 0}^{ℓ} K_{i}) .

Die Größe $\hat{ξ}$ ist eine Funktion von mehreren Veränderlichen, es ist $\hat{ξ} = \hat{ξ} (P_{1}, D, R_{1}, K_{0}, \dots, K_{ℓ})$ . Es sei $| b - a | \neq 0$ . Die Schrittweite $h$ sei so gewählt, daß erstens $| b - a | / | h |$ natürlich ist und zweitens gleichzeitig gilt

| h | < {\begin{cases} 1 / ξ, & falls ξ > 0; \\ \infty, & falls ξ = 0 . \end{cases}

und $N$ sei implizit definiert durch $N | h | = | b - a |$ .

Behauptungen: (1) Beide (möglicherweise) impliziten Differenzengleichungen

\begin{aligned} {\hat{u}}_{n + ℓ} + A_{ℓ - 1} {\hat{u}}_{n + ℓ - 1} + \dots + A_{0} {\hat{u}}_{n} & = h φ ({\hat{u}}_{n + ℓ}, \dots, {\hat{u}}_{n}) + r_{n + ℓ} \\ u_{n + ℓ} + A_{ℓ - 1} u_{n + ℓ - 1} + \dots + A_{0} u_{n} & = h φ (u_{n + ℓ}, \dots, u_{n}) \end{aligned}

besitzen für jedes $n$ , eine eindeutig bestimmte Lösung $u_{n + ℓ}$ bzw. ${\hat{u}}_{n + ℓ}$ , die man mit Picard-Iteration berechnen kann.

(2) Für die maximale normmässige Abweichung $| {\hat{u}}_{n} - u_{n} |$ gilt die beidseitige Abschätzung bzgl. der additiven Störglieder $r_{n}$ , wie folgt

c_{1} | {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | \leq | \hat{U} - U | \leq c_{2} | {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | \leq c_{3} N | R | .

(3) Die positiven Konstanten $c_{i}$ , für $i = 1, 2, 3$ , sind gegeben durch

c_{1} = \frac{1}{1 + \hat{ξ} | b - a |}, c_{2} = \frac{1}{1 - | h | ξ} \exp \frac{\hat{ξ} | b - a |}{1 - | h | ξ}, c_{3} = c_{2} | P_{1} | D | R_{1} | .

(4) Die Abschätzung bei (3) ist unabhängig von der Wahl des Standard-Tripels, d.h. es gilt

{\overset{―}{X}}_{1} [T_{1}]^{- 1} {\overset{―}{Y}}_{1} R = {\overset{―}{X}}_{2} [T_{2}]^{- 1} {\overset{―}{Y}}_{2} R,

für zwei beliebige Standard-Tripel $(X_{1}, T_{1}, Y_{1})$ und $(X_{2}, T_{2}, Y_{2})$ zum Matrixpolynom $ρ$ .

(5) Das verkürzte Funktional $| [C_{1}]^{- 1} {\overset{―}{R}}_{1} R |$ ist ebenfalls Stabilitätsfunktional und zum unverkürzten Funktional äquivalent, unabhängig von $N$ , d.h. es gilt

| {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | \sim | [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | .

(6) Verkürzte Stabilitätsfunktionale sind bei Wechsel des Standard-Tripels untereinander äquivalent, jedoch nicht notwendig mehr gleich. Es gilt

| [T_{1}]^{- 1} {\overset{―}{Y}}_{1} R | \sim | [T_{2}]^{- 1} {\overset{―}{Y}}_{2} R | .

Beweis: Zur Abkürzung werde wieder benutzt

δ_{n + ℓ} := {\hat{u}}_{n + ℓ} - u_{n + ℓ}, {\hat{δ}}_{n + ℓ} := φ ({\hat{u}}_{n + ℓ}, \dots, {\hat{u}}_{n}) - φ (u_{n + ℓ}, \dots, u_{n}) .

Zu (1): Beide Differenzengleichungen stellen für jedes $n$ eine Lipschitz-stetige Keplersche Gleichung in ${\hat{u}}_{n + ℓ}$ bzw. $u_{n + ℓ}$ dar. Die Fixpunktgleichungen bzgl. $\hat{F}$ und $F$ , mit

{\hat{u}}_{n + ℓ} = \hat{F} ({\hat{u}}_{n + ℓ} := h φ ({\hat{u}}_{n + ℓ}, \dots) + \hat{ψ}, bzw. u_{n + ℓ} = F (u_{n + ℓ}) := h φ (u_{n + ℓ}, \dots) + ψ,

sind kontrahierend, falls $| h | K_{ℓ} < 1$ . Durch die oben vorausgesetzte Einschränkung an $h$ , nämlich $| h | ξ < 1$ , ist diese hinreichende Bedingung für Kontraktion erfüllt. Auf einem geeigneten vollständigen Teilraum, lässt sich dann Existenz und Eindeutigkeit eines Fixpunktes deduzieren.

Zu (2): a) Nach dem Hilfssatz über die Darstellung der Differenz der Lösung zweier Differenzengleichungen (siehe Darstellungssatz), folgt sofort durch Umstellung, die Abschätzungskette

\begin{aligned} | P_{1} C_{1}^{n} (\begin{matrix} r_{0} \\ ⋮ \\ r_{ℓ - 1} \end{matrix}) + P_{1} \sum_{ν = 0}^{n - 1} C_{1}^{n - 1 - ν} R_{1} r_{ν + ℓ} | & \leq | δ_{n} | + | h | | P_{1} | D | R_{1} | \sum_{ν = 0}^{n - ℓ} | {\hat{δ}}_{ν + ℓ} | \\ \leq | δ_{n} | + | h | | P_{1} | D | R_{1} | (\sum_{i = 0}^{ℓ} K_{i}) \sum_{ν = 0}^{n - 1} | δ_{ν} | \\ \leq | δ_{n} | + | b - a | | P_{1} | D | R_{1} | (\sum_{i = 0}^{ℓ} K_{i}) {sup}_{ν = 0}^{n - 1} | δ_{ν} | \\ \leq {sup}_{ν = 0}^{n} | δ_{ν} | + | b - a | | P_{1} | D | R_{1} | (\sum_{i = 0}^{ℓ} K_{i}) {sup}_{ν = 0}^{n} | δ_{ν} | \\ = (1 + \hat{ξ} | b - a |) {sup}_{ν = 0}^{n} | δ_{ν} | . \end{aligned}

Hierbei wurde die Abschätzung

\sum_{ν = 0}^{n - ℓ} | {\hat{δ}}_{ν + ℓ} | \leq K_{ℓ} | δ_{n} | + (\sum_{i = 0}^{ℓ - 1} K_{i}) (\sum_{ν = 0}^{n - 1} | δ_{ν} |)

des Abschätzungssatzes benutzt, desweiteren die Gleichung $N | h | = | b - a |$ und schließlich die Abschätzung $\sum_{ν = 0}^{n} | δ_{ν} | \leq N {sup}_{ν = 0}^{n - 1} | δ_{ν} |$ . Durchmultiplikation mit

\frac{1}{1 + \hat{ξ} | b - a |}

liefert die erste Ungleichung der Behauptung (2), wobei sich entsprechend die Konstante $c_{1}$ ergibt.

b) Wiederum nach dem Satz über die Darstellung der Differenz der Lösungen zweier Differenzengleichungen (siehe Darstellungssatz), folgt sofort beim Übergang zu Normen

\begin{aligned} | δ_{n} | & \leq | h | | P_{1} | D | R_{1} | \sum_{ν = 0}^{n - ℓ} | {\hat{δ}}_{ν + ℓ} | + | P_{1} C_{1}^{n} (\begin{matrix} r_{0} \\ ⋮ \\ r_{ℓ - 1} \end{matrix}) + P_{1} \sum_{ν = 0}^{n - 1} C_{1}^{n - 1 - ν} R_{1} r_{ν + ℓ} | \\ \leq | h | \underset{= \hat{ξ}}{\underset{⏟}{| P_{1} | D | R_{1} | (\sum_{i = 0}^{ℓ} K_{i})}} \sum_{ν = 0}^{n - 1} | δ_{ν} | + | h | \underset{= ξ}{\underset{⏟}{| P_{1} | D | R_{1} | K_{ℓ}}} | δ_{n} | + | {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R |, \end{aligned}

wobei wieder der Abschätzungssatz benutzt wurde, also durch Umformung

(1 - | h | ξ) | δ_{n} | \leq | h | \hat{ξ} \sum_{ν = 0}^{n - 1} | δ_{ν} | + | {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | .

Wegen der Voraussetzung an $h$ , nämlich $| h | < 1 / ξ$ , ist $1 - | h | ξ > 0$ . Mit Hilfe des diskreten Lemmas von Gronwall, wobei man in der dort angegebenen Bezeichnung setzt

ε_{j + 1} \leftarrow | δ_{n} |, η_{j} \leftarrow \frac{| {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R |}{1 - | h | ξ}, δ \leftarrow \frac{\hat{ξ}}{1 - | h | ξ}, h_{ν} \leftarrow | h |,

erhält man jetzt die Abschätzung

| δ_{n} | \leq \frac{| {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R |}{1 - | h | ξ} \exp \frac{\hat{ξ} | b - a |}{1 - | h | ξ} .

Anhand dieser Darstellung ersieht man auch das Zustandekommen der Konstanten $c_{2}$ . Die Konstante $c_{3}$ ergibt sich sofort durch typische Abschätzungen.

Zu (4): Das Standard-Tripel $(X_{1}, T_{1}, Y_{1})$ ist ähnlich zum Standard-Tripel $(X_{2}, T_{2}, Y_{2})$ genau dann, wenn

X_{2} = X_{1} S, T_{2} = S^{- 1} T_{1} S, Y_{2} = S^{- 1} Y_{1},

oder

X_{1} = X_{2} S^{- 1}, T_{1} = S T_{2} S^{- 1}, Y_{1} = S Y_{2} .

Nun ist

\begin{aligned} {\overset{―}{X}}_{1} [T_{1}]^{- 1} {\overset{―}{Y}}_{1} R & = ({diag}_{ν = 0}^{N} X_{1}) [T_{1}]^{- 1} diag [({r o w}_{i = 0}^{ℓ - 1} T_{1}^{i} Y) B, {diag}_{ν = 1}^{N} Y_{1}] R \\ = ({diag}_{ν = 0}^{N} (X_{2} S^{- 1})) [S T_{2} S^{- 1}]^{- 1} diag {{r o w}_{i = 0}^{ℓ - 1} [{(S T_{2} S^{- 1})}^{i} S Y_{2}] B, {diag}_{ν = 1}^{N} (S Y_{2})} R \\ = ({diag}_{ν = 0}^{N} X_{2}) [T_{2}]^{- 1} diag [({r o w}_{i = 0}^{ℓ - 1} T_{2}^{i} Y) B, {diag}_{ν = 1}^{N} Y_{2}] R \\ = {\overset{―}{X}}_{2} [T_{2}]^{- 1} {\overset{―}{Y}}_{2} R . \end{aligned}

Hierbei wurden die Recheneigenschaften der Operatoren $diag$ , $row$ und $[\cdot]$ benutzt.

Zu (5): Dies wurde im vorausgeschickten Hilfssatz bewiesen.

Zu (6): Mit der gleichen Notation wie beim Beweis zu (4) rechnet man

\begin{aligned} [T_{1}]^{- 1} {\overset{―}{Y}}_{1} R & = [S T_{2} S^{- 1}]^{- 1} diag {{r o w}_{i = 0}^{ℓ - 1} [{(S T_{2} S^{- 1})}^{i} S Y_{2}] B, {diag}_{ν = 1}^{N} S Y_{2}} R \\ = ({diag}_{ν = 0}^{N} S) [T_{2}]^{- 1} diag [{r o w}_{i = 0}^{ℓ - 1} (T_{2}^{i} Y_{2}) B, {diag}_{ν = 1}^{N} Y_{2}] {\overset{―}{Y}}_{2} R \\ = ({diag}_{ν = 0}^{N} S) [T_{2}]^{- 1} {\overset{―}{Y}}_{2} R . \end{aligned}

Durch Multiplikation von links mit ${diag}_{ν = 0}^{N} S^{- 1}$ folgt sofort

[T_{2}]^{- 1} {\overset{―}{Y}}_{2} R = ({diag}_{ν = 0}^{N} S^{- 1}) [T_{1}]^{- 1} {\overset{―}{Y}}_{1} R .

Damit sind beide verkürzten Funktionale äquivalent. ☐

5. Bemerkungen: Zur Voraussetzung: Aufgrund der Einschränkung der Schrittweite $h$ in Abhängigkeit der Konstanten $ξ$ , ist das Ergebnis nur von praktischer Bedeutung bei kurzen Integrationsintervallen und nicht-steifen Differentialgleichungsproblemen, also Problemen mit kleiner Lipschitzkonstante. Bei steifen Problemen werden die Konstanten $ξ$ , $\hat{ξ}$ , $c_{1}$ , $c_{2}$ , $c_{3}$ schnell unangemessen groß. Die Konstanten $ξ$ und $\hat{ξ}$ enthalten direkt als multiplikativen Faktor die obere Schranke $D$ für die Matrixpotenzen. $ξ$ seinerseits geht exponentiell in die Abschätzung ein. Die Aufspaltung in zwei sehr ähnliche Konstanten $ξ$ und $\hat{ξ}$ geschieht nur, weil $ξ$ i.a. kleiner ist als $\hat{ξ}$ und damit schärfere Schranken liefert. Man könnte mit $\hat{ξ}$ alleine auskommen. Dabei würde man $ξ$ vollständig durch $\hat{ξ}$ ersetzen. $\hat{ξ}$ lässt sich wiederum ersetzen durch $| P_{1} | D | R_{1} | (ℓ + 1) ({max}_{i = 0}^{ℓ} K_{i})$ , man vergl. hier den entsprechenden Hilfssatz mit den diesbezüglichen Abschätzungen, siehe Abschätzungssatz. Erfüllt die erste Begleitmatrix $C_{1}$ die Bedingung $| C_{1}^{ν} | \leq D$ , $\forall ν \in N$ , so auch jede zu $C_{1}$ ähnliche Matrix, allerdings mit u.U. verändertem $D$ .

Numerische Ergebnisse von Tischer/Sacks-Davis (1983)3 und Tischer/Gupta (1985)2 zeigen, daß selbst bei steifen Differentialgleichungen das Stabilitätsfunktional die richtige Konvergenzordnung anzeigt und dies obwohl $| h | ξ < 1$ , verletzt ist. Autoren Peter E. Tischer, Ron Sacks-Davis und Gopal K. Gupta. Dies deutet darauf hin, daß die Ergebnisse des Hauptsatzes allgemeiner gelten als der den Hauptsatz tragende Beweis.

Kerninhalt der Beweise im Hauptsatz sind der Darstellungssatz, das diskrete Lemma von Gronwall und die Abschätzungen im Hilfssatz 6. Die obere Schranke $D$ für die Matrixpotenzen $| C_{1}^{i} |$ hängt u.U. ab von der Dimension der zugrundeliegenden Differentialgleichung $\dot{y} = f (t, y)$ , aufgrund der Beziehung $| C_{1} \otimes I | = | C_{1} |$ , falls die zur Maximumnorm kompatible Zeilensummennorm verwendet wird. Die Lipschitzkonstanten $K_{i}$ sind abhängig von der Lipschitzkonstante von $f$ .

Zu (1): Die Lösungen der möglicherweise impliziten Differenzengleichungen müssen nicht mit Picard-Iteration berechnet werden. Ebenso gut kann das Newton-Raphson-Iterationsverfahren oder das Newton-Kantorovich-Iterationsverfahren benutzt werden. Der Startfehlersatz von Liniger liefert eine obere Schranke für die Anzahl der nötigen Iterationen. Es zeigt sich, daß eine einzige Iteration vielfach vollkommen ausreicht. Weitere Iterationen schaffen keinerlei Verbesserungen an denen man interessiert ist. Für nicht genügend kleine Schrittweiten $| h |$ können in der Tat die beiden angegebenen Differenzengleichungen und damit die entsprechenden Keplerschen Gleichungen keine oder mehr als eine Lösung besitzen.

Zu (2): Die Stabilitätsfunktionale $| {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{P}}_{1} R |$ , $| [C_{1}]^{- 1} {\overset{―}{R}}_{1} R |$ und $| R |$ sind unabhängig von $φ (\cdot)$ , und unabhängig von den Lipschitzkonstanten $K_{i}$ , jedoch abhängig von $N$ und damit letztlich abhängig von $h$ und/oder der Länge des Integrationsintervalles $| b - a |$ . Die Konstanten $c_{1}$ , $c_{2}$ , $c_{3}$ hängen von den Lipschitzkonstanten $K_{i}$ ab. Da bei dem Hauptsatz allerdings als zentrale Voraussetzung die Lipschitzkonstanten eingehen und beim obigen Beweis auch benötigt werden, hängen in soweit auch die Funktionale hiervon ab. Das durch die Konstante $c_{2}$ induzierte exponentielle Wachstum kann bei den gegebenen Voraussetzungen des Hauptsatzes nicht so ohne weiteres verbessert werden, wie z.B. die beiden Differentialgleichungen $\dot{y} = 0$ und $\dot{y} = y$ mit $y (0) = 1$ zeigen, wenn man das explizite Eulerverfahren $y_{n + 1} = y_{n} + h f_{n}$ anwendet. Daß hierdurch auch das qualitative Verhalten gänzlich überschätzt werden kann, zeigen die beiden Differentialgleichungen $\dot{y} = 0$ und $\dot{y} = - y$ , mit $y (0) = 1$ , wenn man das implizite Eulerverfahren $y_{n + 1} = y_{n} + h f_{n + 1}$ verwendet. Dieses Verhalten ist schon beim kontinuierlichen Lemma von Gronwall und den hieraus sich ableitenden Sätzen wohl bekannt. Dort sind allerdings auch Sätze bekannt, die dieses falsche Voraussagen des qualitativen Verhaltens vermeiden, siehe Hairer/Wanner/N\o rsett (1987). {Hairer, Ernst}{Wanner, Gerhard}_{N\o rsett, Syvert Paul}% Hier benutzt man u.a. die logarithmische Norm $μ$ definiert zu

μ (A) := lim_{ε \to 0, ε > 0} \frac{| I + ε A | - 1}{ε}, A \in C^{k \times k} .

Für die euklidische-, Maximum- und die 1-Norm ergeben sich

\begin{aligned} μ (A) & = λ_{m a x}, λ_{m a x} größter Eigenwert von \frac{1}{2} (A^{⊤} + A), \\ μ (A) & = {max}_{k = 1}^{n} (Re a_{i i} + \sum_{i \neq k} | a_{i k} |), \\ μ (A) & = {max}_{i = 1}^{n} (Re a_{i} i + \sum_{k \neq i} | a_{k i} |) . \end{aligned}

Zu (4): Die Aussage (4) zeigt, daß das Stabilitätsfunktional unabhängig von der Basisdarstellung und Basiswahl ist. Die Abschätzungen sind invariant unter der Wahl des Standard-Tripels. Vielfach geeignet ist das Stabilitätsfunktional zum Jordan-Tripel, zum ersten Begleiter-Tripel $(P_{1}, C_{1}, R_{1})$ oder zum zweiten Begleiter-Tripel $(P_{1} B^{- 1}, C_{2}, B R_{1})$ , mit der Block-Hankel-Matrix zu

B := (\begin{matrix} A_{1} & A_{2} & \dots & A_{ℓ - 1} & I \\ A_{2} & ⋰ & ⋰ & ⋰ \\ ⋮ & ⋰ & ⋰ \\ A_{ℓ - 1} & I \\ I \end{matrix}) .

Zu (5) und (6): Gleichheit beider verkürzter Stabilitätsfunktionale ist nicht mehr zu erwarten. Desweiteren erkennt man, daß die Rechtseigenvektoren, also die Spalten in $X_{1}$ bzw. $X_{2}$ (falls einer der beiden zu einem Jordan-Tripel gehört) keinen “kalkülmässigen Einfluß” haben, entgegen den Linkseigenvektoren. Natürlich haben die Rechtseigenvektoren Einfluß auf das Gesamtverhalten, denn ändern sich die Rechtseigenvektoren, repräsentiert durch $X$ , so ändern sich sich nach der Biorthogonalitätsbeziehung

(row_{i = 0}^{ℓ - 1} T^{i} Y) B ({col}_{i = 0}^{ℓ - 1} X T^{i}) = I_{k ℓ \times k ℓ},

auch die Linkseigenvektoren, repräsentiert durch $Y$ . Jedoch brauchen die Rechtseigenvektoren oder gar die Inverse von $col (X T^{i})$ nicht berechnet zu werden. Dies ist hier mit “kalkülmässig” unabhängig gemeint.

6. Corollar: Voraussetzung: $ξ$ , $c_{1}$ , $c_{2}$ wie beim Hauptsatz.

Behauptung: $ξ \to 0$ , $c_{1} \to 1$ , $c_{2} \to 1$ , falls ${max}_{i = 0}^{ℓ} K_{i} \to 0$ .

D.h. die beidseitige Ungleichungskette entartet zu einer Gleichungskette, falls alle Lipschitzkonstanten $K_{ρ}$ gegen Null gehen, was insbesondere bei Quadraturproblemen auftritt. Eine Anfangswertaufgabe für Differentialgleichungen enthält mit $\dot{y} = f (t)$ , $y (a) = 0$ , $y (b) => ?$ , das Quadraturproblem $\int_{a}^{b} f (t) d t$ als Spezialfall.

In Komponentenschreibweise liest sich der Hauptsatz wie folgt.

7. Hauptsatz: Voraussetzung: Es sei

| h | < \frac{1}{| P_{1} | D | R_{1} | K_{ℓ}}

Behauptung: (2) Für die maximale normmässige Abweichung $| {\hat{u}}_{n} - u_{n} |$ gilt die beidseitige Abschätzung bzgl. der additiven Störglieder $r_{n}$ , wie folgt

\begin{aligned} c_{1} {sup}_{n = 0}^{N} | P_{1} C_{1}^{n} (\begin{matrix} r_{0} \\ ⋮ \\ r_{ℓ - 1} \end{matrix}) + P_{1} \sum_{ν = 0}^{n - 1} C_{1}^{n - 1 - ν} R_{1} r_{ν + ℓ} | & \leq {sup}_{n = 0}^{N} | {\hat{u}}_{n} - u_{n} | \\ \leq c_{2} {sup}_{n = 0}^{N} | P_{1} C_{1}^{n} (\begin{matrix} r_{0} \\ ⋮ \\ r_{ℓ - 1} \end{matrix}) + P_{1} \sum_{ν = 0}^{n - 1} C_{1}^{n - 1 - ν} R_{1} r_{ν + ℓ} | \\ \leq c_{3} N {sup}_{n = 0}^{N} | r_{n} | . \end{aligned}

(4) Die Abschätzung bei (3) ist invariant unter der Wahl eines Standard-Tripels, d.h. es gilt

\begin{matrix} {sup}_{n = 0}^{N} | \hat{X} {\hat{T}}^{n} (row_{i = 0}^{ℓ - 1} {\hat{T}}^{i} \hat{Y}) B (\begin{matrix} r_{0} \\ ⋮ \\ r_{ℓ - 1} \end{matrix}) + \hat{X} \sum_{ν = 0}^{n - 1} {\hat{T}}^{n - 1 - ν} \hat{Y} r_{ν + ℓ} | \\ = {sup}_{n = 0}^{N} | X T^{n} (row_{i = 0}^{ℓ - 1} T^{i} Y) B (\begin{matrix} r_{0} \\ ⋮ \\ r_{ℓ - 1} \end{matrix}) + X \sum_{ν = 0}^{n - 1} T^{n - 1 - ν} Y r_{ν + ℓ} |, \end{matrix}

für zwei Standard-Tripel $(X, T, Y)$ und $(\hat{X}, \hat{T}, \hat{Y})$ .

8. Bemerkung: Zu (2): Man erkennt an der Komponentendarstellung, daß $r_{n}$ eingeht in das Stabilitätsfunktional, ohne von rechts “echt” mit $R_{1}$ (bzw. $Y$ ) multipliziert zu werden; $ν = n - ℓ$ , für $ν > n - ℓ$ , also $n - 1 - ν \leq ℓ - 2$ , man denke an $P_{1} C_{1}^{n - 1 - ν} R_{1} = 0$ . M.a.W. für $ν = n - ℓ$ kann $r_{n}$ nicht in den Kernbereich von $P_{1} C_{1}^{n - 1 - ν} R_{1}$ gelangen. Im Falle von Diskretisierungsverfahren, wo die $r_{ν}$ die lokalen Diskretisierungsfehlervektoren darstellen, hat dies zur Konsequenz: Die Ordnung in $h$ der $r_{ν}$ kann nie überschritten werden. Durch Summation kann selbstverständlich eine Reduktion der Ordnung anfallen. Ist beispielsweise $r_{ν} = O (h^{p + 1})$ , so ist $| {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | = O (h^{p + 1 + ε})$ mit $ε > 0$ unmöglich. Für ein Diskretisierungsverfahren ist dennoch ein Ordnungssprung größer 1 möglich, falls gewisse Komponenten von $r_{ν} \in C^{k}$ bei der Ordnungsfindung unberücksichtigt bleiben. Dies ist z.B. bei Runge-Kutta-Verfahren der Fall.

5. Projektorstabilitätsfunktionale

Im weiteren sei vorausgesetzt, daß die Eigenwerte von $C_{1}$ auf dem Einheitskreis nur aus der 1 bestehen, also nicht von der Form $e^{i φ}$ [ $φ \neq 0 (\mod 2 π)$ ] sind, da andernfalls die typischen Projektoreigenschaften verloren gehen. Sei $E$ diejenige Matrix, die lediglich die spektralen Eigenschaften von $C_{1}$ zu dem (den) dominanten Eigenwert(en) $μ = 1$ trägt. $T$ sei die Transformationsmatrix von $C_{1}$ auf Jordannormalform, also

C_{1} = T J T^{- 1}, J = diag (1, \dots, 1, weitere Jordanblöcke) .

Die Matrix $E$ filtert aus dieser Darstellung nur den (die) Eigenwert(e) $μ = 1$ heraus, also

E := T \hat{J} T^{- 1}, \hat{J} := diag (1, \dots, 1, 0, \dots, 0) .

Es zeigt sich, daß $E$ nicht speziell von der Matrix $C_{1}$ abhängt.

Die Eigenwerte $μ = 1$ sind ja gerade diejenigen Eigenwerte, die für den dominanten lokalen Fehler verantwortlich sind. Es gilt $\sum_{ν = 0}^{N} 1 \to \infty$ , falls $N \to \infty$ , aber $\sum_{ν = 0}^{\infty} μ^{ν} < C$ , falls $| μ | < 1$ . Bei $\sum_{ν = 0}^{N} 1$ $(N \to \infty)$ verliert man gerade eine $h$ -Potenz. Es war $N | h | = | b - a |$ , und $N \to \infty ⟺ | h | \to 0$ .

Die Matrix $E$ hat nun eine Reihe von Recheneigenschaften, die in nachstehendem Satz zusammengefaßt sind. $N$ ist hier wie üblich die Menge der natürlichen Zahlen von eins an.

1. Satz: (Projektorsatz) $E$ sei wie oben definiert. $S$ sei eine beliebige Matrix, ähnlich zu $C_{1}$ , also $S = H^{- 1} C_{1} H$ . Dann gelten

$E$ ist idempotent, d.h. $E^{2} = E$ , also allgemein $E^{i} = E$ , für alle $i \in N$ . $E$ ist damit ein Projektor.
$E$ ist unabhängig von $C_{1}$ . $E$ hängt nur ab von $V$ , beim Standard-Tripel $(X, V, Y)$ zu $C_{1}$ .
$S E = E = E S$ . $S^{ν} E = E = E S^{ν}$ , $\forall ν \in N$ .
$S^{n} = E + (S - E)^{n}$ , $\forall n \in N$ .
$[E]^{- 1} [S] = [S - E]$ , also $| [E]^{- 1} [S] | = 1 + | S - E |$ .
$[S]^{- 1} [E] = [S - E]^{- 1}$ , also $| [S]^{- 1} [E] | = 1 + | S - E | + \dots + | (S - E)^{N} | .$
Es gilt

[S - E]^{- 1} {diag}_{ν = 0}^{N} E = {diag}_{ν = 0}^{N} E = ({diag}_{ν = 0}^{N} E) [S - E]^{- 1}

und

[S - E] {diag}_{ν = 0}^{N} E = {diag}_{ν = 0}^{N} E = ({diag}_{ν = 0}^{N} E) [S - E], [S]^{- 1} = [S - E]^{- 1} + [E]^{- 1} - I .

Beweis: Zu (1): $E^{2} = (T \hat{J} T^{- 1}) (T \hat{J} T^{- 1}) = T {\hat{J}}^{2} T^{- 1} = T \hat{J} T^{- 1} = E$ .

Zu (2): Liegt an der Ähnlichkeit von Standard-Tripeln.

Zu (3): $S E = (T J T^{- 1}) (T \hat{J} T^{- 1}) = T J \hat{J} T^{- 1} = T \hat{J} T^{- 1} = E$ . Für $E = E S$ verfährt man analog.

Zu (4): Aufgrund der Vertauschbarkeit von $S$ und $E$ nach (2) und der Projektoreigenschaft nach (1) rechnet man

(S - E)^{n} = \sum_{ν = 0}^{n} (\binom{n}{ν}) (- 1)^{ν} S^{n - ν} E^{ν} = S^{n} + \sum_{ν = 1}^{n} (\binom{n}{ν}) (- 1)^{ν} E = S^{n} - E,

wegen

0 = (1 - 1)^{n} = \sum_{ν = 0}^{n} (\binom{n}{ν}) (- 1)^{ν} und (\binom{n}{0}) = 1.

Zu (5): Berechne direkt anhand der Definition von $[X]$ das Matrixprodukt $[E]^{- 1} [S]$ aus. Auf der ersten Subdiagonalen erscheint stets $E - S$ und auf allen darunterliegenden Diagonalen steht $E - E S,$ was nach Behauptung (3) gleich der Nullmatrix ist.

Zu (6): Folgt sofort aus (5). Invertierung von $[E]^{- 1} [S]$ liefert $[S]^{- 1} [E]$ . Die Invertierbarkeit ist aufgrund der Definition von $[X]$ gesichert. Die angegebenen Normen (Zeilensummennorm oder Spaltensummenorm) ergeben sich unmittelbar.

Zu (7): Ist klar. ☐

2. Beispiel: Zu (5): Für das Produkt $[E]^{- 1} [S]$ im Falle $N = 3$ berechnet man

(\begin{matrix} I & 0 & 0 & 0 \\ E & I & 0 & 0 \\ E^{2} & E & I & 0 \\ E^{3} & E^{2} & E & I \end{matrix}) (\begin{matrix} I & 0 & 0 & 0 \\ - S & I & 0 & 0 \\ 0 & - S & I & 0 \\ 0 & 0 & - S & I \end{matrix}) = (\begin{matrix} I & 0 & 0 & 0 \\ E - S & I & 0 & 0 \\ E - E S & E - S & I & 0 \\ E - E S & E - E S & E - S & I \end{matrix})

Es ist $E - E S = 0$ .

Bei den Ausdrücken hinter $sup (\cdot)$ steht immer eine endliche Menge, für welches natürlich stets das Maximum existiert. Warum schreibt man $sup (\cdot)$ und nicht $max (\cdot)$ ? Weil man später beim Konvergenzsatz zu $N \to \infty$ übergehen will und man dann zu $lim sup (\cdot)$ gelangt.

3. Satz: Voraussetzungen: $| S^{ν} | \leq D$ , $\forall ν \in N$ , $S = X J Y$ , wobei $J$ die (bis auf Permutation eindeutige) Jordanmatrix ist. $X$ ist die Matrix, die die Rechtseigenvektoren trägt und $Y$ enthält in entsprechend umgekehrter Reihenfolge die Linkseigenvektoren, d.h. es gilt $S X = X J$ und $Y S = J Y$ . Weiter gilt $S = X J x^{- 1} = Y^{- 1} J Y = X J Y$ .

Behauptungen: (1) ${\hat{c}}_{1} | [S]^{- 1} R | \leq | \hat{U} - U | \leq {\hat{c}}_{2} | [S]^{- 1} R | \leq {\hat{c}}_{3} N | R |$ .

(2) $| [S]^{- 1} R | \sim | [E]^{- 1} R | \sim | [J]^{- 1} \overset{―}{Y} R |$ .

(3) Die Konstanten ${\hat{c}}_{1}$ , ${\hat{c}}_{2}$ und ${\hat{c}}_{3}$ sind gegeben durch

{\hat{c}}_{1} = \frac{c_{1}}{1 + | S - E |}, {\hat{c}}_{2} = c_{2} \sum_{ν = 0}^{\infty} | (S - E)^{ν} |, {\hat{c}}_{3} = {\hat{c}}_{2} N max (1, | E |) .

Die Werte ${\hat{c}}_{1}$ und ${\hat{c}}_{2}$ sind unabhängig von $N$ , ${\hat{c}}_{3}$ hingegen nicht. $c_{1}, c_{2}$ wie beim Hauptsatz.

Beweis: Man schätzt $| [E]^{- 1} R |$ und $| [S]^{- 1} R |$ gegeneinander ab. Es ist

[E]^{- 1} R = [E]^{- 1} [S] [S]^{- 1} R = [S - E] \cdot [S]^{- 1} R .

Durchmultiplikation mit $[S - E]^{- 1}$ würde jetzt liefern $[S - E]^{- 1} [E]^{- 1} R = [S]^{- 1} R$ . Alternativ könnte man rechnen

[S]^{- 1} R = [S]^{- 1} [E] [E]^{- 1} R = [S - E]^{- 1} \cdot [E]^{- 1} R .

Nach dem Projektorsatz sind $| [S - E] |$ und $| [S - E]^{- 1} |$ für alle $N$ beschränkt und damit sind beide Normen äquivalent. Die letzte Äquivalenz hat ihre Ursache in

[S]^{- 1} R = \overset{―}{X} [J]^{- 1} \overset{―}{Y} R .

☐

In Komponentenschreibweise liest sich der obige Satz wie folgend.

4. Satz: Voraussetzungen: $| S^{ν} | \leq D$ , $\forall ν \in N$ . $E$ ist wie oben und es gelten die restlichen Voraussetzungen des Hauptsatzes.

Behauptungen: (1) Die Stabilitätsnormen

{sup}_{n = 0}^{N} | \sum_{ν = 0}^{n} S^{n - ν} r_{ν} | und {sup}_{n = 0}^{N} | r_{n} + E \sum_{ν = 0}^{n - 1} r_{ν} |

sind zueinander äquivalent.

(2) Es gilt die zweiseitige Fehlerabschätzung

{\hat{c}}_{1} {sup}_{n = 0}^{N} | r_{n} + E \sum_{ν = 0}^{n - 1} r_{ν} | \leq {sup}_{n = 0}^{N} | {\hat{u}}_{n} - u_{n} | \leq {\hat{c}}_{2} {sup}_{n = 0}^{N} | r_{n} + E \sum_{ν = 0}^{n - 1} r_{ν} | \leq {\hat{c}}_{3} {sup}_{n = 0}^{N} | r_{n} | .

5. Satz: Voraussetzungen: $(P_{1}, C_{1}, R_{1})$ sei das erste Begleitertripel, $(X, J, Y)$ sei das Jordan-Tripel zum Matrixpolynom

ρ (μ) = A_{ℓ} μ^{ℓ} + A_{ℓ - 1} μ^{ℓ - 1} + \dots + A_{0} \in C^{k \times k},

$C_{1}^{ν} \leq D$ , $\forall ν \in N$ . Die Matrix $\hat{J}$ filtere aus $J$ nur die Jordanblöcke zum Eigenwert $μ = 1$ heraus. $J$ selber enthalte keine weiteren dominanten Eigenwerte, $J$ ist also stark $D$ -stabil. Zwischen $C_{1}$ und $J$ besteht grundsätzlich der Zusammenhang

C_{1} {col}_{i = 0}^{ℓ - 1} X J^{i} = ({col}_{i = 0}^{ℓ - 1} X J^{i}) J .

Nun sei ${\hat{C}}_{1}$ die entsprechend zu $\hat{J}$ ähnliche Matrix. ${\hat{C}}_{1}$ ist wie $\hat{J}$ Projektor.

Behauptung: Das verkürzte Projektorstabilitätsfunktional ist äquivalent zum Projektorstabilitätsfunktional, welches seinerseits äquivalent ist zum ursprünglichen Stabilitätsfunktional, d.h. es gilt unabhängig von $N$ , daß

| [{\hat{C}}_{1}]^{- 1} {\overset{―}{R}}_{1} R | \sim | {\overset{―}{P}}_{1} [{\hat{C}}_{1}]^{- 1} {\overset{―}{R}}_{1} R | \sim | [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | \sim | {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | .

Diese Äquivalenzen sind unabhängig von der Wahl der Standard-Tripel, d.h. es gilt genauso

| [\hat{J}]^{- 1} \overset{―}{Y} R | \sim | \overset{―}{X} [\hat{J}]^{- 1} \overset{―}{Y} R | \sim | [J]^{- 1} \overset{―}{Y} R | \sim | \overset{―}{X} [J]^{- 1} \overset{―}{Y} R | .

Beweis: Die dritte Äquivalenz wurde schon im Hauptsatz postuliert und bewiesen, desgleichen die Invarianz vom Standard-Tripel. Für die weiteren gegenseitigen Äbschätzungen rechnet man

| [{\hat{C}}_{1}]^{- 1} {\overset{―}{R}}_{1} R | = | [{\hat{C}}_{1}]^{- 1} [C_{1}] [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | = | [C_{1} - {\hat{C}}_{1}] [C_{1}]^{- 1} {\overset{―}{R}}_{1} R |

Für die Rückabschätzung rechnet man

| {\overset{―}{P}}_{1} [C_{1}]^{- 1} {\overset{―}{R}}_{1} R | = | {\overset{―}{P}}_{1} [C_{1}]^{- 1} [{\hat{C}}_{1}] [{\hat{C}}_{1}]^{- 1} {\overset{―}{R}}_{1} R | = | {\overset{―}{P}}_{1} [C_{1} - {\hat{C}}_{1}]^{- 1} [{\hat{C}}_{1}]^{- 1} {\overset{―}{R}}_{1} R | .

Die Beschränktheit der Normen von $[C_{1} - {\hat{C}}_{1}]$ und $[C_{1} - {\hat{C}}_{1}]^{- 1}$ , und zwar gänzlich unabhängig von $N$ , wurde schon vorher gezeigt. An dieser Stelle benutzt man dann $| C_{1}^{ν} | \leq D$ , $\forall ν \in N$ . Die Beschränktheit von $| {\overset{―}{P}}_{1} |$ , unabhängig von $N$ , ist ebenfalls klar. ☐

Wünscht man lediglich ein Konvergenzresultat, so beschränke man sich auf das diskrete Lemma von Gronwall, den Darstellungssatz und beim Abschätzungssatz genügt völlig die letzte, sehr leicht einzusehende Abschätzung. Schließlich beim Hauptsatz genügt lediglich (1), (2) und (3). Weitere Vereinfachungen ergeben sich, falls man sich auf lineare Matrixpolynome der Form $ρ (μ) = I μ - S$ beschränkt, also überall lediglich den Fall $ℓ = 1$ betrachtet. Die untersuchten Verfahrenstypen bleiben dabei die gleichen, man verliert also letztlich nichts an Allgemeinheit. Die Werte $| P_{1} |$ , $| R_{1} |$ entfallen dann völlig. Die Beweise werden kürzer, aber ggf. muß $S$ in den Anwendungen auf Diskretisierungsprobleme unnötig groß gewählt werden. Allerdings kann häufig $S$ kleiner als $C_{1}$ ausfallen, jedoch spielt $C_{1}$ für praktische Rechnungen nicht die entscheidende Rolle, vielmehr ist es $Y$ .

6. Nichtäquidistante Gitter

1. Aufgrund der Konsistenzbedingung $ρ (1) = 0$ für jede Stufe eines zusammengesetzten Verfahrens, gilt $S w = w$ ( $S \in C^{k \times k}$ ), mit $w = (1, \dots, 1)^{⊤} \in C^{k}$ , wenn man das Verfahren in der Form $u_{n + 1} = S u_{n} + h φ (u_{n}, u_{n + 1})$ notiert. Bei zyklischer oder auch nicht-zyklischer Kombination mehrerer Verfahren gelangt man zu $u_{n + 1} = S_{ν} S_{ν - 1} \dots S_{2} S_{1} u_{n} + \dots$ . Als notwendige und hinreichende Bedingung für Stabilität erhält man

‖ S_{i} S_{i + 1} \dots S_{j - 1} S_{j} ‖ beschränkt, \forall i > j .

Ein typischer Fall ist die Benutzung eines Grundverfahrens $u_{n + 1} = S u_{n} + h φ$ und Variation der Schrittweite $h$ . Gilt die obige Stabilitätsbedingung, so spricht man auch von schrittwechsel-stabil. Die Stabilitätsbedingung in der obigen Form ist nicht einfach zu verifizieren. Hinreichende Kriterien sind allerdings viel einfacher zu handhaben, man fordert also mehr und zwar:

Bei einem Schrittweitenwechsel sind sämtliche Matrizen $S_{i}$ identisch, oder/und
nach einem Schrittweitenwechsel wird ein Sonderverfahren mit einer Matrix $T$ nachgeschaltet, mit der Eigenschaft, daß $S_{i} T = T$ gilt, "Matrix-fressende Eigenschaft".

Den ersten Fall kann man so erreichen, daß man die $α_{i j}$ -Koeffizienten eines Verfahrens vorgibt und anschließend die $β_{i j}$ -Koeffizienten berechnet. Das Verfahren ist offensichtlich stabil (die zu $α_{i j}$ gehörende Matrix $S$ war ja stabil) und es konvergiert mit der gleichen Konvergenzordnung wie $S$ , wenn man die $β_{i j}$ so wählt, daß eine ausreichend hohe Konsistenzordnung erreicht wird. Das ist aber stets möglich nach dem Dimensionssatz für die Konsistenzmatrix $C_{p + 1, k}$ .

2. Beim zweiten Fall absorbiert bildlich gesprochen die Matrix $T$ sämtliche vorhergehenden Matrizen. Dies lässt sich bewerkstelligen, wenn die Spaltenvektoren von $T$ Rechsteigenvektoren von $S_{i}$ sind, zum Eigenwert 1. Da aber alle $S_{i}$ konsistent sind, gilt $S_{i} w = w$ $\forall i$ . Damit hat $T$ die Gestalt

T = (ε_{1} w, ε_{2} w, \dots, ε_{k} w), ε_{1}, \dots, ε_{k} \in C

Aufgrund der Konsistenz von $T$ muß gelten $T w = w$ , also

ε w = 1, ε = (ε_{1}, \dots, ε_{k}),

im Falle von $w = (1, \dots, 1)$ also $ε_{1} + \dots + ε_{k} = 1$ .

Egal wie $w$ aussieht, $T$ ist ein Projektor, also $T^{2} = T$ , oder was dasselbe ist: $\ker T$ und $Im T$ sind zueinander komplementäre % Unterräume des $C^{k}$ ( $A_{1} \cap A_{2} = \emptyset$ , $A_{1} + A_{2} = C^{k}$ ). Offensichtlich ist aber nicht jeder konsistente Projektor von der Gestalt $T = (ε_{1} w, \dots, ε_{k} w)$ . Dennoch zeigt sich, daß man bei einem stark stabilen, konsistenten Projektor nur einige Zeit warten muß, bis er die gewünschte Form annimmt, also man eine gewisse Potenz dieser Matrix zu bilden hat. Eine äquivalente Charakterisierung liefert

3. Satz: Spektraldarstellung schrittwechsel-stabiler Matrizen.

Voraussetzung: Es seien $T_{ν} \in C^{k \times k}$ ( $ν = 1, \dots, k - 1$ ) mit

T_{1} = (\begin{matrix} 0 \\ 0 \\ ⋱ \\ 0 \\ 1 \end{matrix}), T_{2} = (\begin{matrix} 0 & 1 \\ 0 \\ ⋱ \\ 0 \\ 1 \end{matrix}), \dots, T_{k - 1} = (\begin{matrix} 0 & 1 \\ 0 & 1 \\ 0 \\ ⋱ \\ 1 \end{matrix}) .

Behauptung: Es gilt

\begin{aligned} T = (ε_{1} w, \dots, ε_{k} w) \\ ε w = 1 \end{aligned}} ⟺ T \sim T_{1} = T_{ν}^{ν}

Beweis: “ $\Rightarrow$ ”: Die Matrix $T = (ε_{1} w, \dots, ε_{k} w)$ hat den Rang genau 1: Jeder Minor mit 2 oder mehr Zeilen verschwindet (Spalten Vielfaches voneinander oder Nullspalte); aus $w \neq 0$ folgt $T \neq 0$ . Damit ist $T$ ähnlich zu einer der Matrizen $T_{ν}$ , $ν = 1, \dots, k - 1$ . Da eine Projektoreigenschaft invariant unter einem Basiswechsel ist % ( $P^{2} = P$ $\Rightarrow$ $S^{- 1} P S S^{- 1} P S = S^{- 1} P S$ ) muß $T$ ähnlich zu $T_{1}$ sein.

“ $\Leftarrow$ ”: Es sei $w$ Rechtseigenvektor von $T$ und $X$ sei die Matrix der Rechtsjordanvektoren und $Y$ sei die Matrix der Linksjordanvektoren, also $T = X T_{1} Y$ , mit

X = (*, \dots, *, w), Y = (\begin{matrix} * \\ ⋮ \\ * \\ v \end{matrix}) .

Multiplikation von rechts mit $T_{1}$ filtert aus $X$ gerade $w$ heraus, Multiplikation von links mit $T_{1}$ filtert aus $Y$ gerade $v$ heraus, also

X T_{1} = (0, \dots, 0, w), T_{1} Y = (\begin{matrix} 0 \\ ⋮ \\ 0 \\ v \end{matrix}) .

Offensichtlich hat $T = (X T_{1}) (T_{1} Y)$ dann die verlangte Gestalt $T = (v_{1} w, \dots, v_{k} w)$ (dyadisches Produkt), mit $v w = 1$ aufgrund der Biorthogonalitätsbeziehung $X Y = Y X = I$ . ☐

Der Beweis zeigt gleichzeitig, daß $ε T = ε$ , also $ε$ Linkseigenvektor von $T$ ist, was man natürlich auch so gesehen hätte. Die Beschränkung beim zweiten Fall auf ein einziges Sonderverfahren, ist nach dem ersten Fall unerheblich, wenn man z.B. immer das gleiche Sonderverfahren mehrmals anwendet. Wie man sieht, muß man $T_{ν}$ nur sooft wiederholen, wie der Nilpotenzgrad von 0 angibt, also $ν$ -mal. Bei einem $k$ -Schritt Adams-Moulton-Verfahren also $(k - 1)$ -mal und bei einem Runge-Kutta-Verfahren einmal.

4. Satz: Die stark stabile Matrix $A \in C^{n \times n}$ habe die Eigenwerte $λ_{1} = 1$ und $| λ_{i} | < 1$ ( $i = 2, \dots, n$ ). Es gelte $T w = w$ , $v_{1} T = v_{1} T$ , $v_{1} w = 1$ , $v_{1} = row (v_{1 i})$ und es sei $T^{\infty} := (v_{11} w, \dots, v_{1 n} w)$ . Dann gilt: $T^{ν} \to T^{\infty}$ ( $ν \to \infty$ ).

Beweis: Für $S := T - T^{\infty}$ gilt wegen $T T^{\infty} = T^{\infty} = (T^{\infty})^{ν} = T^{\infty} T$ offensichtlich $S^{ν} = T^{ν} - T^{\infty}$ ( $ν \in N$ ). Mit den Linkseigenvektoren $v_{2}, \dots, v_{n}$ zu $λ_{2}, \dots, λ_{n}$ ergibt sich $v_{i} T T^{\infty} = v_{i} T^{\infty} = λ_{i} v_{i} T^{\infty}$ , also $v_{i} T^{\infty} = 0$ , somit $v_{i} (T - T^{\infty}) = v_{i} S = λ_{i} v_{i}$ , für $i = 2, \dots, n$ . Weiter ist $v_{1} T^{\infty} = v_{1} = v_{1} T$ , daher $v_{1} (T - T^{\infty}) = v_{1} S = 0$ , folglich hat $S$ die Eigenwerte $0, λ_{2}, \dots, λ_{n}$ , ergo $ρ (S) < 1$ . ☐

Erneut muß $w$ nicht gleich $(1, \dots, 1)^{⊤}$ sein. $v_{1} w = 1$ lässt sich immer erreichen. Bei stark stabilen konsistenten Matrizen ist entweder $T \sim T_{1} ν$ , was äquivalent ist mit $T^{ν} = (ε_{1} w, \dots, ε_{n} w)$ , oder aber zumindestens konvergiert eine Potenz von $T$ gegen diese Gestalt. Das Wiederholen eines stark stabilen Zykluses hat hierin seine Erklärung.

7. Die Eigenwerte gewisser tridiagonaler Matrizen

Sei

A = tridiag (c, a, b) := (\begin{matrix} a & b \\ c & ⋱ & ⋱ \\ ⋱ & ⋱ & b \\ c & a \end{matrix}) \in R^{n \times n},

und es sei $c \cdot b > 0$ . Dann ist $A$ diagonalisierbar mit den $n$ Eigenwerten

λ_{i} = a + 2 \sqrt{b c} \cos \frac{i π}{n + 1}, i = 1, \dots n .

Sei

E = tridiag (1, 0, 1) = (\begin{matrix} 0 & 1 \\ 1 & ⋱ & ⋱ \\ ⋱ & ⋱ & 1 \\ 1 & 0 \end{matrix}) \in R^{n \times n} .

$E$ ist diagonalisierbar mit den Eigenwerten

λ_{i} = 2 \cos \frac{i π}{n + 1}, i = 1, \dots n .

Weiter gelten

B = tridiag (- 1, a, - 1) = (\begin{matrix} a & - 1 \\ - 1 & ⋱ & ⋱ \\ ⋱ & ⋱ & - 1 \\ - 1 & a \end{matrix}) \in C^{n \times n}, λ_{i} = a - 2 \cos \frac{i π}{n + 1}, i = 1, \dots n .

und

T = tridiag (- 1, 2, - 1) = (\begin{matrix} 2 & - 1 \\ - 1 & ⋱ & ⋱ \\ ⋱ & ⋱ & - 1 \\ - 1 & 2 \end{matrix}) \in R^{n \times n}, λ_{i} = 4 \sin \frac{i π}{2 (n + 1)}, i = 1, \dots n .

Mit $T$ gilt dann für $C := 1 - α T = tridiag (α, 1 - 2 α, α)$ , also

C = (\begin{matrix} 1 - 2 α & α \\ α & ⋱ & ⋱ \\ ⋱ & ⋱ & α \\ α & 1 - 2 α \end{matrix}) \in R^{n \times n}, λ_{i} = 1 - 4 α \sin \frac{i π}{2 (n + 1)}, i = 1, \dots n .

8. Verfahren für parabolische Gleichungen

Parabolische, partielle Differentialgleichugen kann man durch Semidiskretisierung der Ortsvariablen auf ein i.a. vergleichsweise großes gewöhnliches Differentialgleichungssystem umformen. Dieses kann man dann mit den üblichen Verfahren numerisch lösen. Ein anderer Weg ist, vollständig zu diskretisieren. Dies bietet u.U. die Möglichkeit Verbindungen zwischen Orts- und Zeitdiskretisierungen zu nutzen. Dies soll hier kurz dargestellt werden.

Betrachtet werde die inhomogene Wärmeleitungsgleichung

u_{t} = σ u_{x x} + f (t, x, u), σ > 0 (Materialkonstante)

mit den Rand- und Anfangsdaten

\begin{aligned} u (0, x) & = η (x), \\ u (t, 0) & = u (t, x_{e}), \end{aligned} \begin{aligned} für alle x \in [0, x_{e}], \\ für alle t \in [0, t_{e}] . \end{aligned}

Die Differentialausdrücke für $u_{t}$ und $u_{x x}$ werden jetzt durch Differenzenausdrücke ersetzt und zwar 1)

\begin{aligned} u_{t} & = \frac{u (t + τ) - u (t)}{τ} + O (τ), Vorwärtsdifferenz \\ u_{x x} & = \frac{u (x + h) - 2 u (x) + u (x - h)}{h^{2}} + O (h), zentrale Differenz \end{aligned}

\begin{aligned} u_{t} & = \frac{u (t + τ) - u (t - τ)}{2 τ} + O (τ^{2}), zentrale Differenz \\ u_{x x} & = \frac{u (x + h) - 2 u (x) + u (x - h)}{h^{2}} + O (h), zentrale Differenz \end{aligned}

\begin{aligned} u_{t} & = \frac{u (t + τ) - u (t)}{τ} + O (τ), Vorwärtsdifferenz \\ u_{x x} & = \frac{u (x + h) - (u (t + τ, x) + u (t - τ, x)) + u (x - h)}{h^{2}} + O (h^{2}) . \end{aligned}

\begin{aligned} u_{t} & = \frac{u (t + τ) - u (t)}{τ} + O (τ), Vorwärtsdifferenz \\ u_{x x} & = \frac{u (t + τ, x + h) - 2 u (t + τ, x) + u (t + τ, x - h)}{h^{2}} + O (h), zentrale Differenz bei (t + τ, x) . \end{aligned}

5) Man wende die Trapezregel ( $ϑ = \frac{1}{2}$ -Verfahren) an, wobei jedoch, wie oben dauernd geschehen, $f$ nicht mit in die Implizitheit mit hineinbezogen wird.

\begin{aligned} u_{t} & = \frac{u (t + τ) - u (t)}{τ} + O (τ), Vorwärtsdifferenz \\ u_{x x} & = (\frac{u (x + h) - 2 u + u (x - h)}{h^{2}} + \frac{u (t + τ, x + h) - 2 u (t + τ, x) + u (t + τ, x - h)}{h^{2}}) / 2 + O (h^{2}), \\ Mittelwert zweier zentraler Differenzen . \end{aligned}

Hierbei wurden zur notationellen Vereinfachung die nicht weiter interessierenden Argumente unterdrückt. Stets ist ist $t$ bzw. $x$ gemeint, also z.B. $u (x + h)$ meint $u (t, x + h)$ und so fort. Diese Schreibweise von $u (t, *) = u (*)$ und $u (*, x) = u (x)$ betont die funktionalen Abhängighkeiten. $u$ ist immer eine Funktion zweier Veränderlicher. Vorauszusetzen ist natürlich $u \in C^{4} ([0, t_{e}] \times [0, x_{e}])$ .

Es sei $N := ⌈ t_{e} / τ ⌉$ die Anzahl der Zeitschritte und $M := ⌈ x_{e} / h ⌉$ sei die Anzahl der Ortsschritte. Weiter sei $t_{i} := i τ = 0, τ, 2 τ, \dots$ , für $i = 0, \dots, N$ , und $x_{k} := k h = 0, h, 2 h, \dots$ , für $k = 0, \dots, M$ . Die Näherung für $u (t_{i}, x_{k})$ werde mit $u_{k}^{i}$ bezeichnet. Entsprechend sei $f_{k}^{i}$ die Näherung für $f (t_{i}, x_{k}, u (t_{i}, x_{k}))$ . Offensichtlich ist $t_{N} = t_{e}$ und $x_{M} = x_{e}$ .

1) Mit der Diskretisierung 1) erhält man jetzt das explizite Einschrittverfahren, wenn man $u_{t}$ und $u_{x x}$ entsprechend ersetzt. Durch Zusammenfassung ergibt sich

u_{k}^{i + 1} = (1 - \frac{2 σ τ}{h^{2}}) u_{k}^{i} + \frac{σ τ}{h^{2}} (u_{k + 1}^{i} + u_{k - 1}^{i}) + τ f_{k}^{i}, i = 0, \dots, N - 1.

2) Einsetzen ergibt das explizite Zweischrittverfahren

u_{k}^{i + 1} = u_{k}^{i - 1} + \frac{2 σ τ}{h^{2}} (u_{k + 1}^{i} - 2 u_{k}^{i} + u_{k - 1}^{i}) + 2 τ f_{k}^{i}, i = 1, \dots, N - 1.

Der Wertevektor $u_{k}^{1}$ muß hierbei auf andere Weise erhalten werden, zum Beispiel durch das obige explizite Einschrittverfahren.

3) Einsetzen liefert das implizite Zweischrittverfahren von DuFort/Frankel aus dem Jahre 1953: {DuFort, E.C.}{Frankel, S.P.}%

(1 + 2 α) u_{k}^{i + 1} = 2 α (u_{k + 1}^{i} + u_{k - 1}^{i}) + (1 - 2 α) u_{k}^{i - 1} + 2 τ f_{k}^{i}, i = 1, \dots, N - 1.

4) Einsetzen liefert das implizite Einschrittverfahren von Crank/Nicolson aus dem Jahre 1947:{Crank, J.}{Nicolson, P.}

- \frac{σ τ}{h^{2}} u_{k + 1}^{i + 1} + (\frac{1 + 2 σ τ}{h^{2}}) u_{k}^{i + 1} - \frac{σ τ}{h^{2}} u_{k - 1}^{i + 1} = u_{k}^{i} + τ f_{k}^{i}, i = 0, \dots, N - 1.

5) Einsetzen ergibt das implizite Einschrittverfahren von Crank/Nicolson (II). Dies entspricht also in etwa der Trapezregel:

- \frac{α}{2} u_{k - 1}^{i + 1} + (1 + α) u_{k}^{i + 1} - \frac{α}{2} u_{k + 1}^{i + 1} = \frac{α}{2} u_{k - 1}^{i} + (1 - α) u_{k}^{i} + \frac{α}{2} u_{k + 1}^{i} + τ f_{k}^{i}, i = 0, \dots, N - 1.

Zur Abkürzung wurde oben benutzt $α := σ τ / h^{2}$ , welches auch im folgenden benutzt werden wird. Alle oben angegebenen Verfahren lassen sich in Matrixschreibweise notieren, anhand dessen man dann das Stabilitätsverhalten besser untersuchen kann, als in der Komponentenform. Zur Abkürzung sei daher im weiteren

v^{i} := (\begin{matrix} u_{1}^{i} \\ ⋮ \\ u_{M - 1}^{i} \end{matrix}), v^{0} := (\begin{matrix} η (x_{1}) \\ ⋮ \\ η (x_{M - 1} \end{matrix}), f^{i} := (\begin{matrix} f_{1}^{i} \\ ⋮ \\ f_{M - 1}^{i} \end{matrix})

und seien die Matrizen definiert

\begin{matrix} A_{1} := tridiag (α, 1 - 2 α, α), A_{2} := tridiag (1, - 2, 1), A_{3} := tridiag (1, 0, 1), \\ A_{4} := tridiag (- α, 1 + 2 α, - α), A_{5} := tridiag (- \frac{α}{2}, 1 + α, - \frac{α}{2}), B_{5} := tridiag (\frac{α}{2}, 1 - α, \frac{α}{2}) . \end{matrix}

Mit diesen Vektoren $v^{i}$ , $f^{i}$ und den Matrizen $A_{ν}$ , $B_{5}$ schreiben sich jetzt die alle oben angegebenen Verfahren, wie folgt.

1) $v^{i + 1} = A_{1} v^{i} + τ f^{i}$ , für $i = 0, \dots, N - 1$ .

2) $v^{i + 1} = 2 α A_{2} v^{i} + v^{i - 1} + 2 τ f^{i}$ , für $i = 1, \dots, N - 1$ .

3) $(1 + 2 α) v^{i + 1} = 2 α A_{3} v^{i} + (1 - 2 α) v^{i - 1} + 2 τ f^{i}$ , für $i = 1, \dots, N - 1$ .

4) $A_{4} v^{i + 1} = v^{i} + τ f^{i}$ , für $i = 0, \dots, N - 1$ .

5) $A_{5} v^{i + 1} = B_{5} v^{i} + τ f^{i}$ , für $i = 0, \dots, N - 1$ .

Die Stabilität aller Verfahren ergibt sich damit aus den spektralen Daten der zum Verfahren gehörenden Matrixpolynome. Von tridiagonlen Matrizen, der obigen Gestalt, nämlich Differenzenmatrizen, sind jedoch die Eigenwerte sämtlich angebbar. Es gilt:

1) $A_{1}$ mit den Eigenwerten $λ_{1, i} := 1 - 4 α \sin^{2} \frac{i π}{2 M}$ , für $i = 1, \dots, M - 1$ .

2) $A_{2}$ mit den Eigenwerten $λ_{2, i} := 4 \sin \frac{i π}{2 M} > 0$ , für $i = 1, \dots, M - 1$ .

3) $A_{3}$ mit den Eigenwerten $λ_{3, i} := 2 \cos \frac{i π}{M}$ , für $i = 1, \dots, M - 1$ .

4) $A_{4}$ mit den Eigenwerten $λ_{4, i} := 1 + 4 α \sin^{2} \frac{i π}{2 M}$ , für $i = 1, \dots, M - 1$ .

5) Für $A_{5}$ rechnet man

A_{5} = tridiag (- \frac{α}{2}, 1 + α, - \frac{α}{2}) = \frac{1}{2} [2 I + α tridiag (- 1, 2, - 1)]

und daher $λ_{5 a, i} := \frac{1}{2} (2 + 4 α \sin \frac{i π}{2 M})$ , und

B_{5} = tridiag (\frac{α}{2}, 1 - α, \frac{α}{2}) = \frac{1}{2} [2 I - α tridiag (- 1, 2, - 1)]

mit den Eigenwerten $λ_{5 b, i} := \frac{1}{2} (2 - 4 α \sin \frac{i π}{2 M})$ , für $i = 1, \dots, M - 1$ .

Für die Stabilität ergeben sich nun unter Berücksichtigung der obigen Matrizen, die folgenden Aussagen.

1) Das Matrixpolynom $I μ - A_{1}$ hat die Eigenwerte $μ_{i} = λ_{1, i}$ . Diese sind genau dann betragsmässig kleiner eins, falls $α \leq 1 / 2$ . Wegen $α = σ τ / h^{2}$ , führt dies auf die Begrenzung der Zeitschrittweite $τ$ zu

τ \leq \frac{1}{2 σ} h^{2},

insbesondere ist die Zeitdiskretisierung nicht unabhängig von der Ortsdiskretisierung. Eine sehr feine Ortsdiskretisierung führt somit automatisch auch zu einer sehr restringierten Zeitschrittweite, und dies obwohl vielleicht in der Zeit eine viel größere Schrittweite angemessener wäre. Dies ist ein typisches Phänomen für explizite Verfahren und ein Grund zur Betrachtung impliziter Verfahren.

2) Für das Matrixpolynom $I μ^{2} - 2 α A_{2} μ - I$ ergeben sich nach Durchmultiplikation mit $D$ , wobei $D$ die Transformationsmatrix für $A_{2}$ ist, also $A_{2} = D diag (λ_{2, i}) D^{- 1}$ , die Eigenwerte des Matrixpolynoms zu

det D \cdot det (I μ^{2} - 2 α diag (λ_{2, i}) - I) \cdot det D^{- 1} = 0,

also

μ_{i, 1 / 2} = λ_{2, i} \pm \sqrt{α^{2} λ_{2, i}^{2} + 1}, i = 1, \dots, M - 1.

Der Spektralradius ist also für jedes $α$ größer als 1. Das Verfahren ist demzufolge für alle $τ$ und $h$ instabil und damit nicht global konvergent, insbesondere als Einzelverfahren nicht brauchbar. In der Kombination mit anderen Verfahren, z.B. im Rahmen eines zyklischen Verfahrens, könnte es u.U. konvergent gemacht werden.

3) Das Matrixpolynom lautet $(1 + 2 α) I μ^{2} - 2 α A_{3} μ - (1 - 2 α) I$ . Sei $D$ die Transformationsmatrix auf Diagonalgestalt für die Matrix $A_{3}$ . Damit erhält man die $2 M - 2$ Eigenwerte des Matrixpolynoms als Nullstellen der Gleichung

det D \cdot det (1 + 2 α) I μ^{2} - 2 α \cdot 2 \cos \frac{i π}{M} μ - (1 - 2 α) I) \cdot det D^{- 1} = 0,

also

(1 + 2 α) μ_{i}^{2} - 4 α μ_{i} \cos φ - (1 - 2 α) = 0,

mit $φ := i π / M$ . Mit der Lösungsformel für quadratische Gleichungen der Form $a x^{2} + b x + c = 0$ , $a \neq 0$ , nämlich

x_{1 / 2} = \frac{- b \pm \sqrt{b^{2} - 4 a c}}{2 a},

erhält man

μ_{i, 1 / 2} = \frac{2 α \pm \sqrt{1 - 4 α^{2} \sin^{2} φ}}{1 + 2 α} .

Längere, aber elementare Rechnungen zeigen, daß die beiden Funktionen $μ_{i, ν} : (α, φ) \mapsto μ_{i, ν} (α, φ)$ , $ν = 1, 2$ , auf dem Rechteck $[0, + \infty [\times [0^{\circ}, 180^{\circ}]$ ihre Extrema annehmen für $α = 0$ oder $φ = 0^{\circ}$ bzw. $φ = 180^{\circ}$ . Hierbei sind eine Reihe von Fallunterscheidungen nötig ( $α \to \infty$ , Radikand positiv oder negativ, $\dots$ ). In der Tat also $| μ_{i, 1 / 2} | < 1$ , für alle $α > 0$ und alle $φ \in] 0^{\circ}, 180^{\circ} [$ . Das Verfahren von DuFort/Frankel ist damit unbeschränkt stabil.

4) Das Matrixpolynom $I μ - A_{4}^{- 1}$ hat die Eigenwerte $μ_{i} = λ_{4, i}^{- 1} < 1$ , für alle $α$ , da $λ_{4, i} > 1$ . Das Verfahren von Crank/Nicolson ist damit für alle $τ$ und $h$ stabil. Aufgrund der Konsistenz folgt damit die Konvergenz.

5) Für das entsprechende Matrixpolynom korrespondierend zu $A_{5} μ = B_{5}$ , ergeben sich die Eigenwerte als Quotient der Eigenwerte von $A_{5}$ und $B_{5}$ , also

μ_{i} = \frac{1 - 2 α \sin (i π / 2 M)}{1 + 2 α \sin (i π / 2 M)} .

Damit ist auch dieses Verfahren unbeschränkt stabil, unabhängig also von den beiden Diskretisierungsgrößen $τ$ und $h$ .