, 58 min read

Konvergenzresultate für feste Schrittweiten

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.

$ \def\diag{\mathop{\rm diag}} \def\col{\mathop{\rm col}} \def\row{\mathop{\rm row}} \def\dcol{\mathop{\rm col\vphantom {dg}}} \def\drow{\mathop{\rm row\vphantom {dg}}} \def\rank{\mathop{\rm rank}} \def\grad{\mathop{\rm grad}} \def\adj#1{#1^*} \def\iadj#1{#1^*} \def\tr{\mathop{\rm tr}} \def\mapright#1{\mathop{\longrightarrow}\limits^{#1}} \def\fracstrut{} $

1. Einführung und grundlegende Begriffe

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

$$ u_{n+1} = Su_n+h\varphi(u_{n-1}),\qquad n=0,1,\ldots,N, \qquad u_i\in\mathfrak{B}, $$

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

$$ N := \left|b-a\over h\right|, \qquad\hbox{also}\qquad \mathopen|h\mathclose| = {\mathopen|b-a\mathclose|\over N}. $$

Subsumiert sind also nicht Verfahren der Vorschrift

$$ A_1u_{n+1} + A_0u_n = h\cdot(B_1F_{n+1} + B_0F_n),\qquad n=0,1,\ldots,N. $$

bzgl. der rein formalen Schreibweise.

Verfahren der leicht allgemeineren Form

$$ u_{n+1} = Su_n+h\varphi(u_{n+1},u_n),\qquad n=0,1,\ldots,N,\tag{*} $$

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 $\varphi$ 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

$$ \eqalign{ \mathopen|\varphi(\hat u_{n+1},\hat u_n) - \varphi(u_{n+1},u_n)\mathclose| &{}\le K_1 \mathopen|\hat u_{n+1}-u_{n+1}\mathclose|,\cr \mathopen|\varphi(\hat u_{n+1},\hat u_n) - \varphi(u_{n+1},u_n)\mathclose| &{}\le K_2 \mathopen|\hat u_n-u_n\mathclose|.\cr } $$

1. Satz: Die Differenzengleichung $(*)$ besitzt, für genügend kleine Schrittweiten $\mathopen|h\mathclose|$, 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 $\varphi$, nach dem Banachschen Fixpunktsatz eine eindeutig bestimmte Lösung und lässt sich durch Fxpunktiteration gewinnen.

$$ \mathopen|h K_2\mathclose| \lt 1, \qquad\hbox{für geeigntes $h$}. $$

    ☐

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

Es seien betrachtet die beiden Verfahren der Form

$$ \eqalign{ \hat u_{n+\ell}+A_{\ell-1}\hat u_{n+\ell-1}+\cdots+A_0\hat u_n &{}= h{\mskip 3mu}\varphi(\hat u_{n+\ell},\hat u_{n+\ell-1},\ldots,\hat u_n)+r_n,\cr u_{n+\ell}+A_{\ell-1}u_{n+\ell-1}+\cdots+A_0u_n &{}= h{\mskip 3mu}\varphi(u_{n+\ell},u_{n+\ell-1},\ldots,u_n).\cr } $$

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 := \pmatrix{I&0&\ldots&0},\qquad R_1 := \pmatrix{0\cr \vdots\cr 0\cr I\cr}, $$

und $\delta_{n+\ell} := \hat u_{n+\ell}-u_{n+\ell}$ und dazu

$$ \hat\delta_{n+\ell} := \varphi(\hat u_{n+\ell},\ldots,\hat u_n) - \varphi(u_{n+\ell},\ldots,u_n). $$

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

Weiter sei

$$ % bold overlined P1 (pee one) bold overlined R1 (err one) \def\bov#1#2{\overline{\bf #1}_{#2}} % boldface and overlined \def\bopo{\bov P1} \def\boro{\bov R1} \def\bfR{{\bf R}} \def\bovy#1{\bov Y{\!#1}} \def\ovbf#1{\overline{\bf #1}} U := \pmatrix{u_0\cr \vdots\cr u_N\cr},\qquad \hat U := \pmatrix{\hat u_0\cr \vdots\cr \hat u_N\cr},\qquad \bf R := \pmatrix{r_0\cr \vdots\cr r_{N+\ell-1}}. $$

Die einzelnen $u_i$ und $\hat u_i$ sind aus dem Vektorraum $\mathbb{R}$, nicht notwendig endlichdimensional. Hierbei sind die $A_i:\mathfrak{B}\to\mathfrak{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. $\mathfrak{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 $\mathfrak{B}=\mathbb{C}^k$ und $\mathfrak{B}=\mathbb{R}^k$, mit $k\ge1$. Gelegentlich gelten die Sätze auch in nicht notwendigerweise kommutativen Ringen $\mathfrak{B}$. Für $\mathfrak{B}$ wird im folgenden stets $\mathbb{C}^k$ gewählt. Die Mengen $\mathbb{C}^{k\ell\times k\ell}$ wären dann entsprechend zu ersetzen durch $\mathbb{R}^{\ell\times\ell}$ 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_\nu$, bei

$$ \sum_{\nu=0}^\ell A_\nu u_{n+\nu} = h{\mskip 3mu}\sum_{\nu=0}^\ell B_\nu{\mskip 3mu}F_{n+\nu}, \qquad n=0,1,\ldots,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_\nu$ 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. $\mathopen|b-a\mathclose|$ wird dann beliebig groß. Der Satz zeigt, daß bei endlichem Integrationsintervall $\mathopen|b-a\mathclose|$, die $\mathopen|\hat U-U\mathclose|$-Norm mit der $\left| \bopo [C_1]^{-1} \boro \bfR\right|$-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)\le w(x)+\int_a^x k(t){\mskip 3mu}h(t){\mskip 3mu}dt,\qquad\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)\le w(x)+\int_a^x \exp\left(\int_t^x k(\tau){\mskip 3mu}d\tau\right)k(t){\mskip 3mu}w(t){\mskip 3mu}dt, \qquad\forall x\in[a,b]. $$

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

$$ H(x) := \int_a^x k(t){\mskip 3mu}h(t){\mskip 3mu}dt. $$

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

$$ H'(x) = k(x){\mskip 3mu}h(x), \qquad H(a)=0, \qquad \forall x\in[a,b]. $$

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

$$ H'(x) = k(x){\mskip 3mu}h(x) \le k(x)\cdot\left(w(x)+H(x)\right), $$

also die lineare Differentialungleichung

$$ H'(x) - k(x){\mskip 3mu}H(x) \le k(x){\mskip 3mu}w(x),\qquad H(a)=0.\tag{*} $$

Multiplikation mit

$$ e^{-K(x)}\gt 0,\qquad K(x) := \int_a^x k(t){\mskip 3mu}dt, $$

führt zu

$$ e^{-K(x)}\left[H'(x) - k(x){\mskip 3mu}H(x)\right] = \left(e^{-K(x)}\cdot H(x)\right)' \buildrel{\displaystyle{(*)\atop\downarrow}}\over\le e^{-K(x)}\cdot k(x){\mskip 3mu}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) \le \int_a^x e^{-K(t)}\cdot k(t){\mskip 3mu}w(t){\mskip 3mu}dt, $$

also wegen $H(a)=0$ somit

$$ H(x)\le\int_a^x e^{K(x)-K(t)}\cdot k(t){\mskip 3mu}w(t){\mskip 3mu}dt $$

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

$$ h(x)\le w(x)+\int_a^x\left(\exp\int_t^x k(\tau){\mskip 3mu}d\tau\right)\cdot k(t){\mskip 3mu}w(t){\mskip 3mu}dt. $$

    ☐

2. Folgerung: Gilt $h(x)\le w+k\int_a^x h(t){\mskip 3mu}dt$, mit festen, nicht-negativen Konstanten $w$ und $k$, so folgt die Abschätzung

$$ h(x)\le w{\mskip 3mu}e^{k\cdot(x-a)},\qquad\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\le\eta_0\le\eta_1\le\ldots\le\eta_m$ vorgegeben. Ferner sei $\delta\ge0$, $h_j\ge0$ und $x_{j+1}=x_j+h_j$. Es gelte die Ungleichung

$$ \varepsilon_0\le\eta_0\qquad\hbox{und}\qquad \varepsilon_{j+1}\le \eta_j + \delta\sum_{\nu=0}^j h_\nu\varepsilon_\nu, \qquad j=0,\ldots,m-1. $$

Dann gilt

$$ \varepsilon_j\le \eta_j{\mskip 3mu}e^{\delta\cdot(x_j-x_0)},\qquad j=0,\ldots,m. $$

Beweis: siehe erneut Helmut Werner und Herbert Arndt in Werner/Arndt (1986). Der Fall $\delta=0$ ist einfach, wegen $e^0=1$. Sei nun $\delta>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 $\delta>0$ vorausgesetzt werden kann. Hier gilt nun

$$ \eqalign{ \varepsilon_{j+1} &{}\le\eta_{j+1}+\delta\sum_{\nu=0}^j h_\nu{\mskip 3mu}\varepsilon_\nu\cr &{}\le\eta_{j+1}+\delta\sum_{\nu=0}^j h_\nu{\mskip 3mu}\eta_\nu{\mskip 3mu}e^{\delta\cdot(x_\nu-x_0)}\cr &{}\le\eta_{j+1}\cdot\left(1+\delta\sum_{\nu=0}^j h_\nu{\mskip 3mu}e^{\delta\cdot(x_\nu-x_0)}\right)\cr &{}\le\eta_{j+1} {\mskip 5mu} e^{\delta\cdot(x_{j+1}-x_0)}.\cr } $$

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

$$ \sum_{\nu=0}^j h_\nu{\mskip 3mu}e^{\delta\cdot(x_\nu-x_0)} \le \int_{x_0}^{x_{j+1}} e^{\delta\cdot(t-x_0)}{\mskip 3mu}dt = {1\over\delta} \left( e^{\delta(x_{j+1}-x_0)}-1 \right). $$

    ☐

3. Notation und Darstellungssatz für Differenzengleichungen

Man vgl. auch Matrixpolynome.

Zur multiplen Lipschitzkonstanten von $\varphi$:

$$ \mathopen|\varphi(u_\ell,\ldots,\hat u_i,\ldots,u_0) - \varphi(u_\ell,\ldots,u_i,\ldots,u_0)\mathclose| \le K_i \cdot \mathopen|\hat u_i-u_i\mathclose|,\qquad i=0,\ldots,\ell. $$

$(\ell+1)$-malige Anwendung der Dreiecksungleichung liefert

$$ \mathopen|\varphi(\hat u_\ell,\ldots,\hat u_0)-\varphi(u_\ell,\ldots,u_0)\mathclose| \le \sum_{i=0}^\ell K_i \mathopen|\hat u_i-u_i\mathclose| = \langle \pmatrix{K_0\cr \vdots\cr K_\ell\cr}, \pmatrix{\mathopen|\hat u_0-u_0\mathclose|\cr \vdots\cr \mathopen|\hat u_\ell-u_\ell\mathclose|\cr} \rangle $$

In der Schreibweise von $\varphi$ seien fortan zahlreiche hier nicht weiter interessierende Argumente der Schreibvereinfachung und der Klarheit wegen weggelassen. Es ist $\varphi(u_\ell,\ldots,u_1) = \varphi(t_\ell,h_\ell,u_\ell,\ldots,u_1)$.

1. Beispiel: Für die Verfahrensvorschrift der Form

$$ A_\ell u_{n+\ell}+\cdots+A_0u_n = h(B_\ell F_{n+\ell}+\cdots+B_0F_n), \qquad F_{n+i} := \pmatrix{f_{N\ell+*}\cr \vdots\cr f_{N\ell+*}\cr}, $$

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

$$ \left|f(t,\hat y)-f(t,y)\right| \le L \mathopen|\hat y-y\mathclose|. $$

Dann gilt für die obigen Lipschitzkonstanten $K_i$ die Verbindung mit der Lipschitzkonstanten der Differentialgleichung zu

$$ K_i = \left|B_i\right|\cdot L,\qquad\hbox{oder ggf.}\qquad K_i = \left|A^{-1}_\ell B_i\right|\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

$$ \left[T\right] := \pmatrix{ I & & & \cr -T & I & & \cr &\ddots&\ddots&\cr & & -T & I\cr}, \qquad \left[T\right]^{-1} = \pmatrix{ I & & & \cr T & I & & \cr \vdots & \vdots & \ddots & \cr T^n & T^{n-1} & \ldots & I\cr}. $$

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 $\mathop{\rm col}$, $\mathop{\rm row}$, $\mathop{\rm diag}$, $[\cdot]$) Es gilt

  1. $\mathop{\rm col} A_\nu B_\nu = \mathop{\rm diag} A_\nu{\mskip 5mu}\mathop{\rm col} B_\nu$.
  2. $\mathop{\rm col} A_\nu B = \left(\mathop{\rm col} A_\nu\right) B$;   Rechtsdistributivität des $\mathop{\rm col}$-Operators.
  3. $\mathop{\rm row} A_\nu B_\nu = \mathop{\rm row} A_\nu{\mskip 5mu}\mathop{\rm diag} B_\nu$.
  4. $\mathop{\rm row} AB_\nu = A{\mskip 3mu} \mathop{\rm row} B_\nu$;   Linksdistributivität des $\mathop{\rm row}$-Operators.
  5. $\mathop{\rm diag} A_\nu B_\nu = \mathop{\rm diag} A_\nu{\mskip 5mu}\mathop{\rm diag} B_\nu$;   multiplikative Distributivität des Bidiagonaloperators.
  6. $\left[S^{-1}TS\right] = \mathop{\rm diag} S^{-1}{\mskip 3mu} \left[T\right]{\mskip 3mu} \mathop{\rm diag} S$.
  7. $\left[S^{-1}TS\right]^{-1} = \mathop{\rm diag} S^{-1}{\mskip 5mu}\left[T\right]^{-1} \mathop{\rm diag} S$.

Beweis: Zu (1):

$$ \mathop{\rm col}_{\nu=0}^n A_\nu B_\nu = \pmatrix{A_0B_0\cr \vdots\cr A_nB_n\cr} = \pmatrix{A_0&&\cr &\ddots&\cr &&A_n\cr}\pmatrix{B_0\cr \vdots\cr B_n\cr}. $$

Zu (3):

$$ \mathop{\rm row}_{\nu=0}^n A_\nu B_\nu = (A_0B_0{\mskip 3mu}\ldots{\mskip 3mu}A_nB_n) = (A_0{\mskip 3mu}\ldots{\mskip 3mu}A_n)\pmatrix{B_0&&\cr &\ddots&\cr &&B_n\cr}. $$

Zu (5)

$$ \mathop{\rm diag}_{\nu=0}^N A_\nu B_\nu = \pmatrix{A_0B_0&&\cr &\ddots&\cr &&A_nB_n\cr} = \pmatrix{A_0&&\cr &\ddots&\cr &&A_n\cr} \pmatrix{B_0&&\cr &\ddots&\cr &&B_n\cr}. $$

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

$$ \pmatrix{ A_{11} & \ldots & A_{1n}\cr \vdots & \ddots & \vdots\cr A_{m1} & \ldots & A_{mn}\cr} \pmatrix{S_1&&\cr &\ddots&\cr &&S_n\cr} = \pmatrix{ A_{11}S_1 & \ldots & A_{1n}S_n\cr \vdots & \ddots & \vdots\cr A_{m1}S_1 & \ldots & A_{mn}S_n\cr}, $$

bzw.

$$ \pmatrix{S_1&&\cr &\ddots&\cr &&S_m\cr} \pmatrix{ A_{11} & \ldots & A_{1n}\cr \vdots & \ddots & \vdots\cr A_{m1} & \ldots & A_{mn}\cr} = \pmatrix{ S_1A_{11} & \ldots & S_1A_{1n}\cr \vdots & \ddots & \vdots\cr S_mA_{m1} & \ldots & S_mA_{mn}\cr}. $$

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

4. Beispiele: Es gilt

$$ \dcol_{i=0}^{\ell-1} (XT^i) = \left(\mathop{\rm diag}_{i=0}^{\ell-1} X\right) \left(\dcol_{i=0}^{\ell-1} T^i\right)\qquad\hbox{und}\qquad \drow_{i=0}^{\ell-1} (T^iY) = \left(\drow_{i=0}^{\ell-1} T^i\right) \left(\mathop{\rm diag}_{i=0}^{\ell-1} Y\right). $$

Im allgemeinen gilt

$$ \mathop{\rm diag}_{i\in U} \mathop{\rm diag}_{k\in V} A_k \ne \mathop{\rm diag}_{k\in V} \mathop{\rm diag}_{i\in U} 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

$$ \left. \eqalign{ \hat u_{n+\ell}+A_{\ell-1}\hat u_{n+\ell-1}+\cdots+A_0\hat u_n &= h{\mskip 3mu}\varphi(\hat u_{n+\ell},\ldots,\hat u_n)+r_{n+\ell}\cr u_{n+\ell}+A_{\ell-1}u_{n+\ell-1}+\cdots+A_0u_n &= h{\mskip 3mu}\varphi(u_{n+\ell},\ldots,u_n)\cr } \right\} \qquad n=0,1,\ldots,N. $$

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

$$ \left.\eqalign{ \delta_{n+\ell} &:= \hat u_{n+\ell} - u_{n+\ell}, \cr \hat\delta_{n+\ell} &:= \varphi(\hat u_{n+\ell},\ldots,\hat u_n) - \varphi(u_{n+\ell},\ldots,u_n) \cr }\right\} \qquad n=0,\ldots,N. $$

Die Differenzengleichung für $\hat u_n$ habe die Startwerte $\hat u_i := u_i + r_i$, für $i=0,\ldots,\ell-1$. Es sei $\delta_i := r_i$, für $i=0,\ldots,\ell-1$, und $r_\nu := \delta_\nu := \hat\delta_\nu := 0$, für $\nu>N$.

Behauptung:

$$ \eqalign{ \delta_n &= P_1 C_1^n \pmatrix{\delta_0\cr \vdots\cr \delta_{\ell-1}\cr} + P_1 \sum_{\nu=0}^{n-\ell} C_1^{n-1-\nu} R_1 \left( r_{\nu+\ell} + h \hat\delta_{\nu+\ell} \right) \cr &= P_1 C_1^n \pmatrix{\delta_0\cr \vdots\cr \delta_{\ell-1}\cr} + P_1 \sum_{\nu=0}^{n-\ell} C_1^{n-1-\nu} R_1 r_{\nu+\ell} + P_1 \sum_{\nu=0}^{n-\ell} C_1^{n-1-\nu} R_1 \hat\delta_{\nu+\ell} \cr } $$

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

$$ x_{n+\ell}+A_{\ell-1}x_{n+\ell-1}+\cdots+A_0x_n = y_n, \qquad n=0,1,\ldots,N $$

lautet

$$ x_n = P_1 C_1^n z_0 + P_1 \sum_{\nu=0}^{n-1} C_1^{n-1-\nu} R_1 y_\nu. $$

Mit den obigen Abkürzungen für $\delta_n$ und $\hat\delta_n$, ergibt sich eine Differenzengleichung für $\delta_n$ zu

$$ \delta_{n+\ell}+A_{\ell-1}\delta_{n+\ell-1}+\cdots+A_0\delta_n = h \hat\delta_n + r_n. $$

Diese Gleichung hat die Lösungsdarstellung

$$ \delta_n = P_1 C_1^n \pmatrix{\delta_0\cr \vdots\cr \delta_{\ell-1}\cr} + P_1 \sum_{\nu=0}^{n-1} C_1^{n-1-\nu} R_1 \left(h\hat\delta_{\nu+\ell} + r_{\nu+\ell}\right) . $$

In der ersten Summe verschwinden die letzten $(\ell-1)$ Terme, wegen $P_1 C_1^i R_1=0$, für $i=0,\ldots,\ell-2$. Daß hier $P_1 C_1^{\ell-1} R_1 = I$, braucht man noch nicht. Für $\nu>n-\ell$ ist $n-1-\nu \le \ell-2$. Daher folgt genau die behauptete Gleichung, wie oben angegeben.     ☐

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

6. Hilfssatz: (Abschätzungssatz) Voraussetzung: $\varphi(\cdot)$ sei in jeder Komponente Lipschitz-stetig mit den Lipschitzkonstanten $K_i$. Die Werte $\delta_{\nu+\ell}$ und $\hat\delta_{\nu+\ell}$ seien wie oben definiert.

Behauptung:

$$ \eqalign{ \sum_{\nu=0}^{n-\ell} \mathopen| \hat\delta_{\nu+\ell} \mathclose| &\le K_\ell\mathopen|\delta_n\mathclose| + \left(\sum_{i=0}^\ell K_i\right) \left(\sum_{\nu=0}^{n-1} \mathopen|\delta_\nu\mathclose|\right)\cr & \le \left(\sum_{i=0}^\ell K_i\right) \left(\sum_{\nu=0}^n \mathopen|\delta_\nu\mathclose|\right)\cr & \le (\ell+1)\cdot\left( \max_{i=0}^\ell K_i\right)\sum_{\nu=0}^n \mathopen|\delta_\nu\mathclose|.\cr } $$

Beweis: Für $\nu=0,\ldots,n-1$ ist

$$ \eqalign{ \mathopen| \hat\delta_{\nu+\ell} \mathclose| &= \left| \varphi(\hat u_{\nu+\ell},\ldots,\hat u_\nu) - \varphi(u_{\nu+\ell},\ldots,u_\nu) \right|\cr &\le K_0 \mathopen|\delta_\nu\mathclose| + K_1 \mathopen|\delta_{\nu+1}\mathclose| + \cdots + K_\ell \mathopen|\delta_{\nu+\ell}\mathclose|.\cr } $$

Sei jetzt, in einer nicht zu Mißverständnissen führenden Doppelbezeichnung, zur Schreibvereinfachung gesetzt $\delta_\nu \gets \mathopen|\delta_\nu\mathclose|$ und $\hat\delta_\nu \gets \mathopen|\hat\delta_\nu\mathclose|$, d.h. die Betragsstriche werden einfach weggelassen. Nun ist hiermit

$$ %\setbox1=\hbox{$\displaystyle{K_1\left(\delta_1+\cdots+\delta_{n-\ell+2} \right)}$} %\dimen1=\wd1 \eqalign{ \sum_{\nu=0}^{n-\ell} \hat\delta_\nu &\le \left(K_0\delta_0+\cdots+K_\ell\delta_\ell\right)+ \left(K_0\delta_1+\cdots+K_\ell\delta_{\ell+1}\right)+\cdots+ \left(K_0\delta_{n-\ell+1}+\cdots+K_\ell\delta_n\right)\cr &= K_0\left(\delta_0+\cdots+\delta_{n-\ell+1}\right)\cr & \qquad+K_1\left(\delta_1+\cdots+\delta_{n-\ell+2} \right)\cr & \qquad\qquad+\qquad\cdots\cr & \qquad\qquad\qquad+K_\ell\left(\delta_\ell+\cdots+\delta_n\right)\cr } $$

Summation und Abschätzung bzgl. der Spalten im obigen Schema zeigt sofort die erste Abschätzung, wenn man die allerletzte Spalte mit $K_\ell$ und $\delta_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 := \pmatrix{I&0&\ldots&0\cr} \in\mathbb{C}^{k\times k\ell},\qquad R_1 := \pmatrix{0\cr \vdots\cr 0\cr I\cr} \in\mathbb{C}^{k\ell\times k}. $$

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

$$ C_1 := \pmatrix{ 0 & I & 0 & \ldots & 0\cr 0 & 0 & I & \ldots & 0\cr \vdots & \vdots & \vdots & \ddots & \vdots\cr & & & \ldots & I\cr -A_0 & -A_1 & & \ldots & -A_{\ell-1}\cr} \in\mathbb{C}^{k\ell\times k\ell}. $$

Desweiteren sei

$$ \bopo := \mathop{\rm diag}_{\nu=0}^N P_1 = \pmatrix{ I & 0 & \ldots & 0 &&&&&&&&&\cr & & & & I & 0 & \ldots & 0 &&&&&\cr & & & & & & & & \ddots &&&&\cr & & & & & & & & & I & 0 & \ldots & 0\cr} \in\mathbb{C}^{(N+1)k\times(N+1)k\ell}, $$

und

$$ \boro := \mathop{\rm diag}\left(I_{k\ell\times k\ell},{\mskip 3mu} \mathop{\rm diag}_{\nu=1}^N R_1\right) = \pmatrix{ I_{k\ell\times k\ell} &&&\cr & 0 &&\cr & \vdots &&\cr & 0 &&\cr & I &&\cr & & \ddots &\cr & && 0\cr & && \vdots\cr & && 0\cr & && I\cr} \in\mathbb{C}^{(N+1)k\ell\times(N+\ell)k}, $$

wobei

$$ \bfR := \mathop{\rm col}_{\nu=0}^{N+\ell-1} r_\nu = \pmatrix{r_0\cr \vdots\cr r_{N+\ell-1}\cr} \in\mathbb{C}^{(N+\ell)k}. $$

Für das Produkt gilt: $[C_1]^{-1} \boro \in \mathbb{C}^{(N+1)k\ell \times (N+\ell)k}$. Es sei $(X,T,Y)$ ein beliebiges Standard-Tripel. Weiter sei

$$ \ovbf X := \mathop{\rm diag}_{\nu=0}^N X = \pmatrix{ X&&&\cr &X&&\cr &&\ddots&\cr &&&X\cr} $$

und

$$ \ovbf Y := \mathop{\rm diag}\left[\left(\dcol_{i=0}^{\ell-1} XT^i\right)^{-1}, \mathop{\rm diag}_{\nu=1}^N Y\right] = \pmatrix{ \left(\mathop{\rm col}_{i=0}^{\ell-1} XT^i\right)^{-1} &&&\cr &Y&&\cr &&\ddots&\cr &&&Y\cr}. $$

Es ist, aufgrund der Biorthogonalitätsbeziehung,

$$ \left( \mathop{\rm col}_{i=0}^{\ell-1} XT^i \right)^{-1} = \left( \mathop{\rm row}_{i=0}^{\ell-1} T^iY \right) B, $$

mit der Block-Hankel-Matrix $B$ zu

$$ B = \pmatrix{ A_1 & \ldots & A_\ell\cr \vdots & \unicode{x22F0} & \cr A_\ell & & \cr }, \qquad A_\ell = I. $$

Die Sonderbehandlung der Blockmatrix bei $\boro$ und $\ovbf Y$ in dem ersten “Diagonalelement” hat seinen Ursprung in der Lösungsdarstellung einer Differenzengleichung für Matrixpolynome der Form

$$ x_n = XJ^n\left(\mathop{\rm col}_{i=0}^{\ell-1} XJ^i\right)\pmatrix{y_0\cr \vdots\cr y_{\ell-1}\cr} + X \sum_{\nu=0}^{n-1} J^{n-1-\nu} Y y_{\nu+\ell}. $$

Für den Fall $\ell=1$, also $\rho(\mu)=I\mu-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 := {\bf0} \in \mathbb{C}^{k\ell\times k\ell}$, 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

$$ \left| [C_1]^{-1} \boro \bfR \right| \sim \left| \bopo [C_1]^{-1} \boro \bfR \right|. $$

Beweis: Der Beweis wird in zwei Teile aufgespalten. Man schätzt beide Stabilitätsfunktionale gegeneinander ab. Die Abschätzung $\left|\bopo [C_1]^{-1} \boro \bfR\right| \le \left|\bopo\right| \left|[C_1]^{-1} \boro \bfR\right|$ ist klar, wobei die Zeilensummennorm von $\bopo$ 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

$$ \left| C_1^n z_0 \right| \le \left| C_1^{\ell-1} \right| {\mskip 3mu} \left| C_1^{n-\ell+1} z_0 \right| = \left| C_1^{\ell-1} \right| {\mskip 3mu} \max_{i=0}^{\ell-1} \left| P_1 C_1^{n+i-\ell+1} z_0 \right|, $$

wegen

$$ \left| C_1^n z_0 \right| = \max_{i=0}^{\ell-1} \left| P_1 C_1^{n+i} z_0 \right|. $$

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

$$ \left| C_1^{n-1-\nu} R_1 r_{\nu+\ell} \right| \le \left| C_1^{\ell-1} \right| {\mskip 3mu} \left| C_1^{n-\ell-\nu} R_1 r_{\nu+\ell} \right| = \left| C_1^{\ell-1} \right| {\mskip 3mu} \max_{i=0}^{\ell-1} \left| P_1 C_1^{n-\ell-\nu+i} R_1 r_{\nu+\ell} \right|, $$

wegen $r_\nu := 0$, für $\nu>N$.     ☐

2. Beispiel: Sei $\ell=2$ und sei $N:=n:=3$. Es ist $\rho(\mu)=I\mu^2+A_1\mu+A_0\in\mathbb{C}^{k\times k}$ und die Potenzen der ersten Begleitmatrix $C_1$ lauten $C_1^\nu$, für $\nu=1,\ldots,N$:

$$ C_1 = \pmatrix{ 0 & I\cr -A_0 & -A_1\cr},\quad C_1^2 = \pmatrix{ -A_0 & -A_1\cr A_1A_0 & -A_0+A_1^2\cr},\quad C_1^3 = \pmatrix{ A_1A_0 & -A_0+A_1^2\cr A_0^2-A_1^2A_0 & A_0A_1+A_1A_0-A_1^3\cr}. $$

Es war

$$ P_1 = \pmatrix{I&0\cr} \in\mathbb{C}^{k\times 2k},\qquad R_1 = \pmatrix{0\cr I\cr}\in\mathbb{C}^{2k\times k}. $$

Die Matrizen $\bopo$ und $\boro$ haben das Aussehen

$$ \bopo = \pmatrix{ I&0 && &&\cr && I&0 &&\cr && && I&0\cr}\in\mathbb{C}^{3k\times6k},\qquad \boro = \pmatrix{ I&&&\cr &I&&\cr &&0&\cr &&I&\cr &&&0\cr &&&I\cr}\in\mathbb{C}^{6k\times4k}. $$

Man berechnet

$$ [C_1]^{-1} \boro = \pmatrix{ I & \cr C_1 & R_1 &\cr C_1^2 & C_1R_1 & R_1 &\cr C_1^3 & C_1^2R_1 & C_1R_1 & R_1\cr } \in \mathbb{C}^{8k\times5k} $$

zu

$$ \begin{pmatrix} \matrix{I&0\cr 0&I\cr} &\\[1em] %\noalign{\vskip 9pt} \matrix{0&I\cr -A_0&-A_1\cr} & \matrix{0\cr I\cr} &\\[1em] %\noalign{\vskip 9pt} \matrix{-A_0&-A_1\cr A_1A_0&-A_0+A_1\cr} & \matrix{I\cr -A_1\cr} & \matrix{0\cr I\cr} &\\[1em] %\noalign{\vskip 9pt} \matrix{A_1A_0&-A_0+A_1^2\cr A_0^2-A_1^2A_0&A_0A_1+A_1A_0-A_1^3\cr}& \matrix{-A_1\cr -A_0+A_1^2\cr} & \matrix{I\cr -A_1\cr} & \matrix{0\cr I\cr} \cr \end{pmatrix} $$

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

3. Beispiel: Es sei nun $\ell=3$ und $N=3$. Es ist $\rho(\mu)=I\mu^3+A_2\mu^2+A_1\mu+A_0 \in \mathbb{C}^{k\times k}$. Nun berechnet man $C_1^\nu$ für $\nu=1,\ldots,N$:

$$ C_1 = \pmatrix{ 0 & I & 0\cr 0 & 0 & I\cr -A_0 & -A_1 & -A_2\cr},\qquad C_1^2 = \pmatrix{ 0 & 0 & I\cr -A_0 & -A_1 & -A_2\cr A_2A_0 & -A_0+A_2A_1 & -A_1+A_2^2\cr} $$

und

$$ C_1^3 = \pmatrix{ -A_0 & -A_1 & -A_3\cr A_2A_0 & -A_0+A_2A_1 & -A_1+A_2^2\cr -A_2A_0^2 & A_2A_0+A_1^2-A_2A-1 & -A_0+A_1A_2+A_2A_1-A_2^3\cr}. $$

Weiter ist $\bopo\in\mathbb{C}^{4k\times4k\ell}$ und $\boro\in\mathbb{C}^{4k\ell\times6k}$ mit

$$ \bopo = \pmatrix{ I&0&0 &&& &&& &&&\cr &&& I&0&0 &&& &&&\cr &&& &&& I&0&0 &&&\cr &&& &&& &&& I&0&0\cr},\qquad \boro = \pmatrix{ I && &&&\cr & I & &&&\cr && I &&&\cr &&& 0 &&\cr &&& 0 &&\cr &&& I &&\cr &&& & 0 &\cr &&& & 0 &\cr &&& & I &\cr &&& && 0\cr &&& && 0\cr &&& && I\cr}. $$

Nun ist $[C_1]^{-1} \boro \in \mathbb{C}^{12k\times6k}$ mit

$$ [C_1]^{-1} \boro = \pmatrix{ I &&&\cr C_1 & R_1 &&\cr C_1^2 & C_1R_1 & R_1 &\cr C_1^3 & C_1^2R_1 & C_1R_1 & R_1 } \in \mathbb{C}^{4k\ell\times6k}, $$

also

$$ \begin{pmatrix} \matrix{I&&\cr &I&\cr &&I\cr} &&&\\[1em] %\noalign{\vskip 9pt} \matrix{0 & I & 0\cr 0 & 0 & I\cr -A_0 & -A_1 & -A_2\cr} & \matrix{0\cr 0\cr I\cr}&&\\[1em] %\noalign{\vskip 9pt} \matrix{ 0 & 0 & I\cr -A_0 & -A_1 & -A_2\cr A_2A_0 & -A_0+A_2A_1 & -A_1+A_2^2\cr} & \matrix{0\cr I\cr -A_2\cr} & \matrix{0\cr 0\cr I\cr} &\\[1em] %\noalign{\vskip 9pt} \matrix{ -A_0 & -A_1 & -A_3\cr A_2A_0 & -A_0+A_2A_1 & -A_1+A_2^2\cr -A_2A_0^2 & A_2A_0+A_1^2-A_2A-1 & -A_0+A_1A_2+A_2A_1-A_2^3\cr} & \matrix{I\cr -A_2\cr -A_1+A_2^2\cr} & \matrix{0\cr I\cr -A_2\cr} & \matrix{0\cr 0\cr I\cr}\cr \end{pmatrix} $$

Das zugrunde liegende Schema ist hier

$$ \begin{matrix} \matrix{1&1&1\cr 2&2&2\cr 3&3&3\cr}\\[1em] %\noalign{\vskip 9pt} \matrix{2&2&2\cr 3&3&3\cr 4&4&4\cr}\\[1em] %\noalign{\vskip 9pt} \matrix{\vdots & \vdots & \vdots\cr}\cr \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

$$ \rho(\mu) := I\mu^\ell+A_{\ell-1}\mu^{\ell-1}+\cdots+A_1\mu+A_0, $$

vom Grade $\ell\ge1$. Die Funktion $\varphi$ sei Lipschitz-stetig in jeder Komponente mit den Lipschitzkonstanten $K_i$, also

$$ \left|\varphi(u_\ell,\ldots,\hat u_i,\ldots,u_0)- \varphi(u_\ell,\ldots,u_i,\ldots,u_0)\right| \le K_i\cdot\left|\hat u_i-u_i\right|,\quad\hbox{für}\quad i=0,\ldots,\ell. $$

Die Potenzen der Matrix $C_1$ seien beschränkt durch die obere Schranke $D$, also $\left|C_1^\nu\right|\le D$, $\forall\nu\in\mathbb{N}$. Seien $\xi$ und $\hat\xi$ definiert durch

$$ \xi := \left|P_1\right| D \left|R_1\right| K_\ell,\qquad \hat\xi := \left|P_1\right| D \left|R_1\right| \left(\sum_{i=0}^\ell K_i\right). $$

Die Größe $\hat\xi$ ist eine Funktion von mehreren Veränderlichen, es ist $\hat\xi=\hat\xi(P_1,D,R_1,K_0,\ldots,K_\ell)$. Es sei $\mathopen|b-a\mathclose|\ne0$. Die Schrittweite $h$ sei so gewählt, daß erstens $\mathopen|b-a\mathclose| / \mathopen|h\mathclose|$ natürlich ist und zweitens gleichzeitig gilt

$$ \mathopen|h\mathclose| \lt \cases {1/\xi,&falls $\xi\gt 0$;\cr \infty,&falls $\xi=0$.\cr} $$

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

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

$$ \eqalign{ \hat u_{n+\ell}+A_{\ell-1}\hat u_{n+\ell-1}+\cdots+A_0\hat u_n &= h{\mskip 3mu}\varphi(\hat u_{n+\ell},\ldots,\hat u_n)+r_{n+\ell}\cr u_{n+\ell}+A_{\ell-1}u_{n+\ell-1}+\cdots+A_0u_n &= h{\mskip 3mu}\varphi(u_{n+\ell},\ldots,u_n)\cr } $$

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

(2) Für die maximale normmässige Abweichung $\left|\hat u_n-u_n\right|$ gilt die beidseitige Abschätzung bzgl. der additiven Störglieder $r_n$, wie folgt

$$ c_1 \left| \bopo [C_1]^{-1} \boro \bfR \right| \le \left| \hat U-U\right| \le c_2 \left| \bopo [C_1]^{-1} \boro \bfR \right| \le c_3 N \left| \bfR \right| . $$

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

$$ c_1 = {1\over1+\hat\xi\mathopen|b-a\mathclose|},\qquad c_2 = {1\over1-\mathopen|h\mathclose|\xi} \exp{\hat\xi\mathopen|b-a\mathclose| \over 1-\mathopen|h\mathclose|\xi},\qquad c_3 = c_2 \left|P_1\right| D \left|R_1\right|. $$

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

$$ \bov X1 [T_1]^{-1} \bovy1 \bfR = \bov X2 [T_2]^{-1} \bovy2 \bfR, $$

für zwei beliebige Standard-Tripel $(X_1,T_1,Y_1)$ und $(X_2,T_2,Y_2)$ zum Matrixpolynom $\rho$.

(5) Das verkürzte Funktional $\left|[C_1]^{-1} \boro \bfR \right|$ ist ebenfalls Stabilitätsfunktional und zum unverkürzten Funktional äquivalent, unabhängig von $N$, d.h. es gilt

$$ \left| \bopo [C_1]^{-1} \boro \bfR \right| \sim \left| [C_1]^{-1} \boro \bfR \right|. $$

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

$$ \left| [T_1]^{-1} \bovy1 \bfR \right| \sim \left| [T_2]^{-1} \bovy2 \bfR \right|. $$

Beweis: Zur Abkürzung werde wieder benutzt

$$ \delta_{n+\ell} := \hat u_{n+\ell} - u_{n+\ell},\qquad \hat\delta_{n+\ell} := \varphi(\hat u_{n+\ell},\ldots,\hat u_n) - \varphi(u_{n+\ell},\ldots,u_n). $$

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

$$ \hat u_{n+\ell} = \hat F(\hat u_{n+\ell} := h\varphi(\hat u_{n+\ell},\ldots{\mskip 5mu})+\hat\psi, \qquad\hbox{bzw.}\qquad u_{n+\ell} = F(u_{n+\ell}) := h\varphi(u_{n+\ell},\ldots{\mskip 5mu})+\psi, $$

sind kontrahierend, falls $\mathopen|h\mathclose| K_\ell < 1$. Durch die oben vorausgesetzte Einschränkung an $h$, nämlich $\mathopen|h\mathclose|\xi<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

$$ \eqalignno{ \left|P_1 C_1^n \pmatrix{r_0\cr \vdots\cr r_{\ell-1}\cr} + P_1 \sum_{\nu=0}^{n-1} C_1^{n-1-\nu} R_1 r_{\nu+\ell}\right| &\le \mathopen|\delta_n\mathclose| + \mathopen|h\mathclose| \left|P_1\right| D \left|R_1\right| \sum_{\nu=0}^{n-\ell} \left|\hat\delta_{\nu+\ell}\right| \cr &\le \mathopen|\delta_n\mathclose| + \mathopen|h\mathclose| \left|P_1\right| D \left|R_1\right| \left(\sum_{i=0}^\ell K_i\right) \sum_{\nu=0}^{n-1} \mathopen|\delta_\nu\mathclose| \cr &\le \mathopen|\delta_n\mathclose| + \mathopen|b-a\mathclose| \left|P_1\right| D \left|R_1\right| \left(\sum_{i=0}^\ell K_i\right) \sup_{\nu=0}^{n-1} \mathopen|\delta_\nu\mathclose| \cr &\le \sup_{\nu=0}^n \mathopen|\delta_\nu\mathclose| + \mathopen|b-a\mathclose| \left|P_1\right| D \left|R_1\right| \left(\sum_{i=0}^\ell K_i\right) \sup_{\nu=0}^n \mathopen|\delta_\nu\mathclose| \cr &= \left( 1+\hat\xi\mathopen|b-a\mathclose| \right) \sup_{\nu=0}^n \mathopen|\delta_\nu\mathclose| . \cr } $$

Hierbei wurde die Abschätzung

$$ \sum_{\nu=0}^{n-\ell} \left|\hat\delta_{\nu+\ell}\right| \le K_\ell \mathopen|\delta_n\mathclose| + \left(\sum_{i=0}^{\ell-1} K_i\right) \left(\sum_{\nu=0}^{n-1} \mathopen|\delta_\nu\mathclose|\right) $$

des Abschätzungssatzes benutzt, desweiteren die Gleichung $N\mathopen|h\mathclose| = \mathopen|b-a\mathclose|$ und schließlich die Abschätzung $\sum_{\nu=0}^n \mathopen|\delta_\nu\mathclose| \le N\sup_{\nu=0}^{n-1} \mathopen|\delta_\nu\mathclose|$. Durchmultiplikation mit

$$ { 1 \over 1 + \hat\xi \mathopen|b-a\mathclose| } $$

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

$$ \eqalign{ \mathopen|\delta_n\mathclose| &\le \mathopen|h\mathclose|{\mskip 3mu}\left|P_1\right| D \left|R_1\right| \sum_{\nu=0}^{n-\ell} \left|\hat\delta_{\nu+\ell}\right| + \left| P_1 C_1^n \pmatrix{r_0\cr \vdots\cr r_{\ell-1}\cr} + P_1 \sum_{\nu=0}^{n-1} C_1^{n-1-\nu} R_1 r_{\nu+\ell} \right| \cr &\le \mathopen|h\mathclose| {\mskip 3mu} \underbrace{ \left|P_1\right| D \left|R_1\right| \left( \sum_{i=0}^\ell K_i \right) }_{\displaystyle{{}=\hat\xi}} \sum_{\nu=0}^{n-1} \mathopen|\delta_\nu\mathclose| + \mathopen|h\mathclose| {\mskip 3mu} \underbrace{ \left|P_1\right| D \left|R_1\right| K_\ell } _{\displaystyle{{}=\xi}} \mathopen|\delta_n\mathclose| + \left| \bopo [C_1]^{-1} \boro \bfR \right|, \cr } $$

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

$$ \left(1-\mathopen|h\mathclose|\xi\right) \mathopen|\delta_n\mathclose| \le \mathopen|h\mathclose|\hat\xi \sum_{\nu=0}^{n-1} \mathopen|\delta_\nu\mathclose| + \left| \bopo [C_1]^{-1} \boro \bfR \right|. $$

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

$$ \varepsilon_{j+1} \gets \mathopen|\delta_n\mathclose|,\qquad \eta_j \gets {\left|\bopo[C_1]^{-1}\boro\bfR\right| \over1-\mathopen|h\mathclose|\xi}, \qquad \delta \gets {\hat\xi \over 1-\mathopen|h\mathclose|\xi},\qquad h_\nu \gets \mathopen|h\mathclose|, $$

erhält man jetzt die Abschätzung

$$ \mathopen|\delta_n\mathclose| \le {\left|\bopo[C_1]^{-1}\boro\bfR\right| \over 1-\mathopen|h\mathclose|\xi} \exp {\hat\xi\mathopen|b-a\mathclose| \over 1-\mathopen|h\mathclose|\xi}. $$

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_1S,\qquad T_2=S^{-1}T_1S,\qquad Y_2=S^{-1}Y_1, $$

oder

$$ X_1 = X_2 S^{-1}, \qquad T_1 = S T_2 S^{-1}, \qquad Y_1 = S Y_2. $$

Nun ist

$$ \eqalignno{ \bov X1 [T_1]^{-1} \bovy1 \bfR &= \left(\mathop{\rm diag}_{\nu=0}^N X_1\right) [T_1]^{-1} \mathop{\rm diag}\left[ \left(\drow_{i=0}^{\ell-1} T_1^iY\right) B, {\mskip 3mu} \mathop{\rm diag}_{\nu=1}^N Y_1 \right] \bfR \cr &= \left(\mathop{\rm diag}_{\nu=0}^N(X_2S^{-1})\right) [ST_2S^{-1}]^{-1} \mathop{\rm diag}\left\{ \drow_{i=0}^{\ell-1}\left[ \left(ST_2S^{-1}\right)^i SY_2\right] B, {\mskip 3mu} \mathop{\rm diag}_{\nu=1}^N (SY_2) \right\} \bfR \cr &= \left(\mathop{\rm diag}_{\nu=0}^N X_2\right) [T_2]^{-1} \mathop{\rm diag}\left[ \left(\drow_{i=0}^{\ell-1} T_2^iY\right) B, {\mskip 3mu} \mathop{\rm diag}_{\nu=1}^N Y_2 \right] \bfR \cr &= \bov X2 [T_2]^{-1} \bovy2 \bfR. \cr } $$

Hierbei wurden die Recheneigenschaften der Operatoren $\mathop{\rm diag}$, $\mathop{\rm 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

$$ \eqalign{ [T_1]^{-1} \bovy1 \bfR &= [ST_2S^{-1}]^{-1} \mathop{\rm diag}\left\{ \drow_{i=0}^{\ell-1}\left[ \left(ST_2S^{-1}\right)^i SY_2\right] B, {\mskip 3mu} \mathop{\rm diag}_{\nu=1}^N SY_2 \right\} \bfR \cr &= \left(\mathop{\rm diag}_{\nu=0}^N S\right) [T_2]^{-1} \mathop{\rm diag}\left[ \drow_{i=0}^{\ell-1} \left(T_2^iY_2\right) B,{\mskip 3mu} \mathop{\rm diag}_{\nu=1}^N Y_2 \right] \bovy2 \bfR \cr &= \left(\mathop{\rm diag}_{\nu=0}^N S\right) [T_2]^{-1} \bovy2 \bfR . \cr } $$

Durch Multiplikation von links mit $\mathop{\rm diag}_{\nu=0}^N S^{-1}$ folgt sofort

$$ [T_2]^{-1} \bovy2 \bfR = \left(\mathop{\rm diag}_{\nu=0}^N S^{-1}\right) [T_1]^{-1} \bovy1 \bfR. $$

Damit sind beide verkürzten Funktionale äquivalent.     ☐

5. Bemerkungen: Zur Voraussetzung: Aufgrund der Einschränkung der Schrittweite $h$ in Abhängigkeit der Konstanten $\xi$, 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 $\xi$, $\hat\xi$, $c_1$, $c_2$, $c_3$ schnell unangemessen groß. Die Konstanten $\xi$ und $\hat\xi$ enthalten direkt als multiplikativen Faktor die obere Schranke $D$ für die Matrixpotenzen. $\xi$ seinerseits geht exponentiell in die Abschätzung ein. Die Aufspaltung in zwei sehr ähnliche Konstanten $\xi$ und $\hat\xi$ geschieht nur, weil $\xi$ i.a. kleiner ist als $\hat\xi$ und damit schärfere Schranken liefert. Man könnte mit $\hat\xi$ alleine auskommen. Dabei würde man $\xi$ vollständig durch $\hat\xi$ ersetzen. $\hat\xi$ lässt sich wiederum ersetzen durch $\left|P_1\right| D \left|R_1\right| \left(\ell+1\right) \left(\max_{i=0}^\ell K_i\right)$, man vergl. hier den entsprechenden Hilfssatz mit den diesbezüglichen Abschätzungen, siehe Abschätzungssatz. Erfüllt die erste Begleitmatrix $C_1$ die Bedingung $\left|C_1^\nu\right|\le D$, $\forall\nu\in\mathbb{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 $\mathopen|h\mathclose| \xi < 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 $\left|C_1^i\right|$ hängt u.U. ab von der Dimension der zugrundeliegenden Differentialgleichung $\dot y=f(t,y)$, aufgrund der Beziehung $\left|C_1\otimes I\right|=\left|C_1\right|$, 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 $\mathopen|h\mathclose|$ 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 $\left|\bopo [C_1]^{-1} \bopo \bfR\right|$, $\left|[C_1]^{-1} \boro \bfR\right|$ und $\left|\bfR\right|$ sind unabhängig von $\varphi(\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 $\mathopen|b-a\mathclose|$. 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+hf_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+hf_{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 $\mu$ definiert zu

$$ \mu(A) := \lim_{\varepsilon\to0,{\mskip 3mu}\varepsilon\gt 0} {\left|I+\varepsilon A\right| - 1 \over \varepsilon}, A \in \mathbb{C}^{k\times k}. $$

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

$$ \eqalignno{ \mu(A) &= \lambda_{\rm max}, \qquad \lambda_{\rm max} \hbox{ größter Eigenwert von } {1\over2}(A^\top+A),\cr \mu(A) &= \max_{k=1}^n \left(\mathop{\rm Re}\nolimits a_{ii} + \sum_{i\ne k} \left|a_{ik}\right|\right),\cr \mu(A) &= \max_{i=1}^n \left(\mathop{\rm Re}\nolimits a_ii + \sum_{k\ne i} \left|a_{ki}\right|\right). } $$

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_1B^{-1},C_2,BR_1)$, mit der Block-Hankel-Matrix zu

$$ B := \pmatrix{ A_1 & A_2 & \ldots & A_{\ell-1} & I\cr A_2 & \unicode{x22F0} & \unicode{x22F0} & \unicode{x22F0} & \cr \vdots & \unicode{x22F0} & \unicode{x22F0} & & \cr A_{\ell-1} & I & & & \cr I & & & & \cr}. $$

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

$$ \left(\mathop{\rm row}\_{i=0}^{\ell-1} T^iY\right) {\mskip 3mu} B {\mskip 3mu} \left(\mathop{\rm col}_{i=0}^{\ell-1} XT^i\right) = I_{k\ell\times k\ell}, $$

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

6. Corollar: Voraussetzung: $\xi$, $c_1$, $c_2$ wie beim Hauptsatz.

Behauptung: $\xi\to0$, $c_1\to1$, $c_2\to1$, falls $\displaystyle\max_{i=0}^\ell K_i \to 0$.

D.h. die beidseitige Ungleichungskette entartet zu einer Gleichungskette, falls alle Lipschitzkonstanten $K_\rho$ 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) dt$ als Spezialfall.

In Komponentenschreibweise liest sich der Hauptsatz wie folgt.

7. Hauptsatz: Voraussetzung: Es sei

$$ \mathopen|h\mathclose| \lt {1\over \left|P_1\right| D \left|R_1\right| K_\ell} $$

Behauptung: (2) Für die maximale normmässige Abweichung $\left|\hat u_n-u_n\right|$ gilt die beidseitige Abschätzung bzgl. der additiven Störglieder $r_n$, wie folgt

$$ \eqalign{ c_1 \sup_{n=0}^N \left| P_1 C_1^n \pmatrix{r_0\cr \vdots\cr r_{\ell-1}\cr} + P_1 \sum_{\nu=0}^{n-1} C_1^{n-1-\nu} R_1 r_{\nu+\ell} \right| &\le \sup_{n=0}^N \left|\hat u_n-u_n\right| \cr &\le c_2 \sup_{n=0}^N \left| P_1 C_1^n \pmatrix{r_0\cr \vdots\cr r_{\ell-1}\cr} + P_1 \sum_{\nu=0}^{n-1} C_1^{n-1-\nu} R_1 r_{\nu+\ell} \right| \cr &\le c_3 N \sup_{n=0}^N \left|r_n\right| . \cr } $$

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

$$ \displaylines{ \sup_{n=0}^N \left| \hat X \hat T^n \left(\mathop{\rm row}\_{i=0}^{\ell-1} \hat T^i \hat Y\right) B \pmatrix{r_0\cr \vdots\cr r_{\ell-1}\cr} + \hat X \sum_{\nu=0}^{n-1} \hat T^{n-1-\nu}\hat Y r_{\nu+\ell} \right| \cr {}= \sup_{n=0}^N \left| X T^n \left(\mathop{\rm row}\_{i=0}^{\ell-1} T^i Y\right) B \pmatrix{r_0\cr \vdots\cr r_{\ell-1}\cr} + X \sum_{\nu=0}^{n-1} T^{n-1-\nu} Y r_{\nu+\ell} \right| , \cr } $$

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; $\nu=n-\ell$, für $\nu>n-\ell$, also $n-1-\nu\le\ell-2$, man denke an $P_1C_1^{n-1-\nu}R_1=0$. M.a.W. für $\nu=n-\ell$ kann $r_n$ nicht in den Kernbereich von $P_1C_1^{n-1-\nu}R_1$ gelangen. Im Falle von Diskretisierungsverfahren, wo die $r_\nu$ die lokalen Diskretisierungsfehlervektoren darstellen, hat dies zur Konsequenz: Die Ordnung in $h$ der $r_\nu$ kann nie überschritten werden. Durch Summation kann selbstverständlich eine Reduktion der Ordnung anfallen. Ist beispielsweise $r_\nu={\cal O}(h^{p+1})$, so ist $\left|\bopo [C_1]^{-1} \boro \bfR\right| = {\cal O}(h^{p+1+\varepsilon})$ mit $\varepsilon>0$ unmöglich. Für ein Diskretisierungsverfahren ist dennoch ein Ordnungssprung größer 1 möglich, falls gewisse Komponenten von $r_\nu\in\mathbb{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\varphi}$ [$\varphi\ne0 \pmod{2\pi}$] 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) $\mu=1$ trägt. $T$ sei die Transformationsmatrix von $C_1$ auf Jordannormalform, also

$$ C_1=TJT^{-1},\qquad J=\mathop{\rm diag}(1,\ldots,1,\hbox{weitere Jordanblöcke}). $$

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

$$ E := T\hat JT^{-1}, \qquad \hat J := \mathop{\rm diag}(1,\ldots,1,0,\ldots,0). $$

Es zeigt sich, daß $E$ nicht speziell von der Matrix $C_1$ abhängt.

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

Die Matrix $E$ hat nun eine Reihe von Recheneigenschaften, die in nachstehendem Satz zusammengefaßt sind. $\mathbb{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_1H$. Dann gelten

  1. $E$ ist idempotent, d.h. $E^2=E$, also allgemein $E^i=E$, für alle $i\in\mathbb{N}$. $E$ ist damit ein Projektor.
  2. $E$ ist unabhängig von $C_1$. $E$ hängt nur ab von $V$, beim Standard-Tripel $(X,V,Y)$ zu $C_1$.
  3. $SE=E=ES$. $S^\nu E=E=ES^\nu$, $\forall\nu\in\mathbb{N}$.
  4. $S^n=E+(S-E)^n$, $\forall n\in\mathbb{N}$.
  5. $[E]^{-1}[S]=[S-E]$, also $\left|[E]^{-1}[S]\right| = 1+\left|S-E\right|$.
  6. $[S]^{-1}[E]=[S-E]^{-1}$, also $\left|[S]^{-1}[E]\right| = 1 +\left|S-E\right| + \cdots + \left|(S-E)^N\right|.$
  7. Es gilt
$$ [S-E]^{-1} \mathop{\rm diag}_{\nu=0}^N E = \mathop{\rm diag}_{\nu=0}^N E = \left(\mathop{\rm diag}_{\nu=0}^N E\right) [S-E]^{-1} $$

und

$$ [S-E] \mathop{\rm diag}_{\nu=0}^N E = \mathop{\rm diag}_{\nu=0}^N E = \left(\mathop{\rm diag}_{\nu=0}^N E\right) [S-E], \qquad [S]^{-1} = [S-E]^{-1} + [E]^{-1} - I. $$

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

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

Zu (3): $SE=(TJT^{-1})(T\hat JT^{-1})=TJ\hat JT^{-1}=T\hat JT^{-1}=E$. Für $E=ES$ 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_{\nu=0}^n {n\choose\nu} (-1)^\nu S^{n-\nu} E^\nu = S^n + \sum_{\nu=1}^n {n\choose\nu} (-1)^\nu E = S^n - E, $$

wegen

$$ 0 = (1-1)^n = \sum_{\nu=0}^n {n\choose\nu} (-1)^\nu \qquad\hbox{und}\qquad {n\choose 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-ES,$ 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

$$ \pmatrix{ I & 0 & 0 & 0\cr E & I & 0 & 0\cr E^2 & E & I & 0\cr E^3 & E^2 & E & I\cr} \pmatrix{ I & 0 & 0 & 0\cr -S & I & 0 & 0\cr 0 & -S & I & 0\cr 0 & 0 & -S & I\cr} = \pmatrix{ I & 0 & 0 & 0\cr E-S & I & 0 & 0\cr E-ES & E-S & I & 0\cr E-ES & E-ES & E-S & I\cr} $$

Es ist $E-ES=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 $\limsup(\cdot)$ gelangt.

3. Satz: Voraussetzungen: $\left|S^\nu\right| \le D$, $\forall\nu\in\mathbb{N}$, $S=XJY$, 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 $SX=XJ$ und $YS=JY$. Weiter gilt $S=XJx^{-1}=Y^{-1}JY=XJY$.

Behauptungen: (1) $\displaystyle \hat c_1 \left| [S]^{-1} \bfR \right| \le \left| \hat U-U \right| \le \hat c_2 \left| [S]^{-1} \bfR \right| \le \hat c_3 N \left|\bfR\right|$.

(2) $\displaystyle \left|[S]^{-1}\bfR\right| \sim \left|[E]^{-1}\bfR\right| \sim \left|[J]^{-1} \ovbf Y \bfR\right|$.

(3) Die Konstanten $\hat c_1$, $\hat c_2$ und $\hat c_3$ sind gegeben durch

$$ \hat c_1 = {c_1\over1+\left|S-E\right|},\qquad \hat c_2 = c_2 \sum_{\nu=0}^\infty \left|(S-E)^\nu\right|,\qquad \hat c_3 = \hat c_2 N \max\left(1,\left|E\right|\right). $$

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 $\left|[E]^{-1}\bfR\right|$ und $\left|[S]^{-1}\bfR\right|$ gegeneinander ab. Es ist

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

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

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

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

$$ [S]^{-1} \bfR = \ovbf X [J]^{-1} \ovbf Y \bfR. $$

    ☐

In Komponentenschreibweise liest sich der obige Satz wie folgend.

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

Behauptungen: (1) Die Stabilitätsnormen

$$ \sup_{n=0}^N\left|\sum_{\nu=0}^n S^{n-\nu} r_\nu\right| \qquad\hbox{und}\qquad \sup_{n=0}^N\left|r_n + E \sum_{\nu=0}^{n-1} r_\nu\right| $$

sind zueinander äquivalent.

(2) Es gilt die zweiseitige Fehlerabschätzung

$$ \hat c_1\sup_{n=0}^N\left|r_n + E \sum_{\nu=0}^{n-1} r_\nu\right| \le \sup_{n=0}^N \left|\hat u_n-u_n\right| \le \hat c_2 \sup_{n=0}^N\left|r_n + E \sum_{\nu=0}^{n-1} r_\nu\right| \le \hat c_3 \sup_{n=0}^N \left|r_n\right|. $$

5. Satz: Voraussetzungen: $(P_1,C_1,R_1)$ sei das erste Begleitertripel, $(X,J,Y)$ sei das Jordan-Tripel zum Matrixpolynom

$$ \rho(\mu) = A_\ell\mu^\ell + A_{\ell-1}\mu^{\ell-1} + \cdots + A_0 \in \mathbb{C}^{k\times k}, $$

$C_1^\nu\le D$, $\forall\nu\in\mathbb{N}$. Die Matrix $\hat J$ filtere aus $J$ nur die Jordanblöcke zum Eigenwert $\mu=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 \mathop{\rm col}_{i=0}^{\ell-1} XJ^i = \left(\mathop{\rm col}_{i=0}^{\ell-1} XJ^i\right) 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ß

$$ \left| [\hat C_1]^{-1} \boro \bfR \right| \sim \left| \bopo [\hat C_1]^{-1} \boro \bfR \right| \sim \left| [C_1]^{-1} \boro \bfR \right| \sim \left| \bopo [C_1]^{-1} \boro \bfR \right|. $$

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

$$ \left| [\hat J]^{-1} \ovbf Y \bfR \right| \sim \left| \ovbf X [\hat J]^{-1} \ovbf Y \bfR \right| \sim \left| [J]^{-1} \ovbf Y \bfR \right| \sim \left| \ovbf X [J]^{-1} \ovbf Y \bfR \right|. $$

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

$$ \left| [\hat C_1]^{-1} \boro \bfR \right| = \left| [\hat C_1]^{-1} [C_1] [C_1]^{-1} \boro \bfR \right| = \left| [C_1-\hat C_1] [C_1]^{-1} \boro \bfR \right| $$

Für die Rückabschätzung rechnet man

$$ \left| \bopo [C_1]^{-1} \boro \bfR \right| = \left| \bopo [C_1]^{-1} [\hat C_1] [\hat C_1]^{-1} \boro \bfR \right| = \left| \bopo [C_1-\hat C_1]^{-1} [\hat C_1]^{-1} \boro \bfR \right|. $$

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 $\left|C_1^\nu\right|\le D$, $\forall\nu\in\mathbb{N}$. Die Beschränktheit von $\left|\bopo\right|$, 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 $\rho(\mu)=I\mu-S$ beschränkt, also überall lediglich den Fall $\ell=1$ betrachtet. Die untersuchten Verfahrenstypen bleiben dabei die gleichen, man verliert also letztlich nichts an Allgemeinheit. Die Werte $\left|P_1\right|$, $\left|R_1\right|$ 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 $\rho(1)=0$ für jede Stufe eines zusammengesetzten Verfahrens, gilt $Sw=w$ ($S\in\mathbb{C}^{k\times k}$), mit $w=(1,\ldots,1)^\top \in\mathbb{C}^k$, wenn man das Verfahren in der Form $u_{n+1}=Su_n+h\varphi(u_n,u_{n+1})$ notiert. Bei zyklischer oder auch nicht-zyklischer Kombination mehrerer Verfahren gelangt man zu $u_{n+1}=S_\nu S_{\nu-1} \ldots S_2 S_1 u_n + \ldots{\mskip 3mu}$. Als notwendige und hinreichende Bedingung für Stabilität erhält man

$$ \left\|S_i S_{i+1} \ldots S_{j-1} S_j\right\| \hbox{ beschränkt,}\qquad \forall i\gt j. $$

Ein typischer Fall ist die Benutzung eines Grundverfahrens $u_{n+1}=Su_n+h\varphi$ 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:

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

Den ersten Fall kann man so erreichen, daß man die $\alpha_{ij}$-Koeffizienten eines Verfahrens vorgibt und anschließend die $\beta_{ij}$-Koeffizienten berechnet. Das Verfahren ist offensichtlich stabil (die zu $\alpha_{ij}$ gehörende Matrix $S$ war ja stabil) und es konvergiert mit der gleichen Konvergenzordnung wie $S$, wenn man die $\beta_{ij}$ 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 = (\varepsilon_1 w, \varepsilon_2 w, \ldots, \varepsilon_k w), \qquad \varepsilon_1,\ldots,\varepsilon_k \in \mathbb{C} $$

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

$$ \varepsilon w = 1, \qquad \varepsilon=(\varepsilon_1,\ldots,\varepsilon_k), $$

im Falle von $w=(1,\ldots,1)$ also $\varepsilon_1+\cdots+\varepsilon_k=1$.

Egal wie $w$ aussieht, $T$ ist ein Projektor, also $T^2=T$, oder was dasselbe ist: $\ker T$ und $\mathop{\rm Im} T$ sind zueinander komplementäre % Unterräume des $\mathbb{C}^k$ ($A_1\cap A_2=\emptyset$, $A_1+A_2=\mathbb{C}^k$). Offensichtlich ist aber nicht jeder konsistente Projektor von der Gestalt $T=(\varepsilon_1 w,\ldots,\varepsilon_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_\nu \in \mathbb{C}^{k\times k}$ ($\nu=1,\ldots,k-1$) mit

$$ T_1 = \pmatrix{ 0&&&&\cr &0&&&\cr &&\ddots&&\cr &&&0&\cr &&&&1\cr }, \quad T_2 = \pmatrix{ 0&1&&&\cr &0&&&\cr &&\ddots&&\cr &&&0&\cr &&&&1\cr }, \quad \ldots, \quad T_{k-1} = \pmatrix{ 0&1&&&\cr &0&1&&\cr &&0&&\cr &&&\ddots&\cr &&&&1\cr }. % \in \mathbb{C}^k. $$

Behauptung: Es gilt

$$ \left. \eqalign{ &T=(\varepsilon_1 w,\ldots,\varepsilon_k w)\cr &\varepsilon w=1\cr} \right\} \iff T\sim T_1 = T_\nu^\nu % (\nu=1,\ldots,k) $$

Beweis: “$\Rightarrow$”: Die Matrix $T=(\varepsilon_1 w,\ldots,\varepsilon_k w)$ hat den Rang genau 1: Jeder Minor mit 2 oder mehr Zeilen verschwindet (Spalten Vielfaches voneinander oder Nullspalte); aus $w\ne0$ folgt $T\ne0$. Damit ist $T$ ähnlich zu einer der Matrizen $T_\nu$, $\nu=1,\ldots,k-1$. Da eine Projektoreigenschaft invariant unter einem Basiswechsel ist % ($P^2=P$ $\Rightarrow$ $S^{-1}PSS^{-1}PS=S^{-1}PS$) 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=XT_1Y$, mit

$$ X=(*,\ldots,*,w), \qquad Y=\pmatrix{*\cr\vdots\cr*\cr v\cr}. $$

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

$$ XT_1=(0,\ldots,0,w), \qquad T_1Y=\pmatrix{0\cr\vdots\cr0\cr v\cr}. $$

Offensichtlich hat $T=(XT_1)(T_1Y)$ dann die verlangte Gestalt $T=(v_1w,\ldots,v_kw)$ (dyadisches Produkt), mit $vw=1$ aufgrund der Biorthogonalitätsbeziehung $XY=YX=I$.     ☐

Der Beweis zeigt gleichzeitig, daß $\varepsilon T=\varepsilon$, also $\varepsilon$ 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_\nu$ nur sooft wiederholen, wie der Nilpotenzgrad von 0 angibt, also $\nu$-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\mathbb{C}^{n\times n}$ habe die Eigenwerte $\lambda_1=1$ und $\left|\lambda_i\right|<1$ ($i=2,\ldots,n$). Es gelte $Tw=w$, $v_1T=v_1T$, $v_1w=1$, $v_1=\mathop{\rm row}(v_{1i})$ und es sei $T^\infty := (v_{11}w,\ldots,v_{1n}w)$. Dann gilt: $T^\nu\to T^\infty$ ($\nu\to\infty$).

Beweis: Für $S:=T-T^\infty$ gilt wegen $TT^\infty=T^\infty=(T^\infty)^\nu=T^\infty T$ offensichtlich $S^\nu=T^\nu-T^\infty$ ($\nu\in\mathbb{N}$). Mit den Linkseigenvektoren $v_2,\ldots,v_n$ zu $\lambda_2,\ldots,\lambda_n$ ergibt sich $v_iTT^\infty=v_iT^\infty=\lambda_iv_iT^\infty$, also $v_iT^\infty=0$, somit $v_i(T-T^\infty)=v_iS=\lambda_iv_i$, für $i=2,\ldots,n$. Weiter ist $v_1T^\infty=v_1=v_1T$, daher $v_1(T-T^\infty)=v_1S=0$, folglich hat $S$ die Eigenwerte $0,\lambda_2,\ldots,\lambda_n$, ergo $\rho(S)<1$.     ☐

Erneut muß $w$ nicht gleich $(1,\ldots,1)^\top$ sein. $v_1w=1$ lässt sich immer erreichen. Bei stark stabilen konsistenten Matrizen ist entweder $T\sim T_1\nu$, was äquivalent ist mit $T^\nu=(\varepsilon_1w,\ldots,\varepsilon_nw)$, 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 = \mathop{\rm tridiag}(c,a,b) := \pmatrix{ a & b & & \cr c & \ddots & \ddots & \cr & \ddots & \ddots & b\cr & & c & a\cr } \in \mathbb{R}^{n\times n} , $$

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

$$ \lambda_i = a + 2\sqrt{bc}{\mskip 3mu}\cos{i\pi\over n+1}, \qquad i=1,\ldots n. $$

Sei

$$ E = \mathop{\rm tridiag}(1,0,1) = \pmatrix{ 0 & 1 & & \cr 1 & \ddots & \ddots & \cr & \ddots & \ddots & 1\cr & & 1 & 0\cr } \in \mathbb{R}^{n\times n}. $$

$E$ ist diagonalisierbar mit den Eigenwerten

$$ \lambda_i = 2\cos{i\pi\over n+1},\qquad i=1,\ldots n. $$

Weiter gelten

$$ B = \mathop{\rm tridiag}(-1,a,-1) = \pmatrix{ a & -1 & & \cr -1 & \ddots & \ddots & \cr & \ddots & \ddots & -1\cr & & -1 & a\cr } \in \mathbb{C}^{n\times n},\quad \lambda_i=a-2\cos{i\pi\over n+1},\quad i=1,\ldots n. $$

und

$$ T = \mathop{\rm tridiag}(-1,2,-1) = \pmatrix{ 2 & -1 & & \cr -1 & \ddots & \ddots & \cr & \ddots & \ddots & -1\cr & & -1 & 2\cr } \in \mathbb{R}^{n\times n},\quad \lambda_i=4\sin{i\pi\over 2(n+1)},\quad i=1,\ldots n. $$

Mit $T$ gilt dann für $C:=1-\alpha T=\mathop{\rm tridiag}(\alpha,1-2\alpha,\alpha)$, also

$$ C = \pmatrix{ 1-2\alpha & \alpha & & \cr \alpha & \ddots & \ddots & \cr & \ddots & \ddots & \alpha\cr & & \alpha & 1-2\alpha\cr } \in \mathbb{R}^{n\times n},\quad \lambda_i=1-4\alpha\sin{i\pi\over 2(n+1)},\quad i=1,\ldots 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 = \sigma u_{xx}+f(t,x,u), \qquad \sigma\gt 0 \hbox{ (Materialkonstante)} $$

mit den Rand- und Anfangsdaten

$$ \eqalign{ u(0,x) &= \eta(x),\cr u(t,0) &= u(t,x_e),\cr }\quad \eqalign{ &\hbox{für alle}\quad x\in[0,x_e],\cr &\hbox{für alle}\quad t\in[0,t_e].\cr } $$

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

$$ \eqalign{ u_t &= {u(t+\tau)-u(t)\over\tau}+{\cal O}(\tau), \qquad\hbox{Vorwärtsdifferenz}\cr u_{xx} &= {u(x+h)-2u(x)+u(x-h)\over h^2}+{\cal O}(h),\qquad \hbox{zentrale Differenz}\cr } $$

2)

$$ \eqalign{ u_t &= {u(t+\tau)-u(t-\tau)\over2\tau}+{\cal O}(\tau^2),\qquad \hbox{zentrale Differenz}\cr u_{xx} &= {u(x+h)-2u(x)+u(x-h)\over h^2}+{\cal O}(h),\qquad \hbox{zentrale Differenz}\cr } $$

3)

$$ \eqalign{ u_t &= {u(t+\tau)-u(t)\over\tau}+{\cal O}(\tau), \qquad\hbox{Vorwärtsdifferenz}\cr u_{xx} &= {u(x+h)-\bigl(u(t+\tau,x)+u(t-\tau,x)\bigr)+u(x-h)\over h^2} +{\cal O}(h^2).\cr } $$

4)

$$ \eqalign{ u_t &= {u(t+\tau)-u(t)\over\tau}+{\cal O}(\tau), \qquad\hbox{Vorwärtsdifferenz}\cr u_{xx} &= {u(t+\tau,x+h)-2u(t+\tau,x)+u(t+\tau,x-h)\over h^2} +{\cal O}(h), \qquad\hbox{zentrale Differenz bei $(t+\tau,x)$}.\cr } $$

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

$$ \eqalign{ u_t &= {u(t+\tau)-u(t)\over\tau}+{\cal O}(\tau), \qquad\hbox{Vorwärtsdifferenz}\cr u_{xx} &= \left({u(x+h)-2u+u(x-h)\over h^2} +{u(t+\tau,x+h)-2u(t+\tau,x)+u(t+\tau,x-h)\over h^2}\right)\bigg/2 +{\cal O}(h^2),\fracstrut\cr &\qquad\hbox{Mittelwert zweier zentraler Differenzen}.\cr } $$

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:=\lceil t_e/\tau\rceil$ die Anzahl der Zeitschritte und $M:=\lceil x_e/h\rceil$ sei die Anzahl der Ortsschritte. Weiter sei $t_i:=i\tau=0,\tau,2\tau,\ldots$, für $i=0,\ldots,N$, und $x_k:=kh=0,h,2h,\ldots$, für $k=0,\ldots,M$. Die Näherung für $u(t_i,x_k)$ werde mit $u^i_k$ bezeichnet. Entsprechend sei $f^i_k$ 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_{xx}$ entsprechend ersetzt. Durch Zusammenfassung ergibt sich

$$ u^{i+1}_k=\left(1-{2\sigma\tau\over h^2}\right)u^i_k +{\sigma\tau\over h^2}\left(u^i_{k+1}+u^i_{k-1}\right)+\tau f^i_k, \qquad i=0,\ldots,N-1. $$

2) Einsetzen ergibt das explizite Zweischrittverfahren

$$ u^{i+1}_k=u^{i-1}_k+{2\sigma\tau\over h^2}\left(u^i_{k+1}-2u^i_k+u^i_{k-1}\right) +2\tau f^i_k,\qquad i=1,\ldots,N-1. $$

Der Wertevektor $u^1_k$ 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\alpha)u^{i+1}_k=2\alpha\left(u^i_{k+1}+u^i_{k-1}\right) +(1-2\alpha)u^{i-1}_k+2\tau f^i_k,\qquad i=1,\ldots,N-1. $$

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

$$ -{\sigma\tau\over h^2}u^{i+1}_{k+1}+\left(1+2\sigma\tau\over h^2\right)u^{i+1}_k -{\sigma\tau\over h^2}u^{i+1}_{k-1}=u^i_k+\tau f^i_k, \qquad i=0,\ldots,N-1. $$

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

$$ -{\alpha\over2}u^{i+1}_{k-1}+(1+\alpha)u^{i+1}_k-{\alpha\over2}u^{i+1}_{k+1} ={\alpha\over2}u^i_{k-1}+(1-\alpha)u^i_k+{\alpha\over2}u^i_{k+1}+\tau f^i_k, \qquad i=0,\ldots,N-1. $$

Zur Abkürzung wurde oben benutzt $\alpha:=\sigma\tau/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 := \pmatrix{u^i_1\cr \vdots\cr u^i_{M-1}\cr},\qquad v^0 := \pmatrix{\eta(x_1)\cr \vdots\cr \eta(x_{M-1}\cr},\qquad f^i := \pmatrix{f^i_1\cr \vdots\cr f^i_{M-1}\cr} $$

und seien die Matrizen definiert

$$ \displaylines{ A_1 := \mathop{\rm tridiag}(\alpha,1-2\alpha,\alpha),\quad A_2 := \mathop{\rm tridiag}(1,-2,1),\quad A_3 := \mathop{\rm tridiag}(1,0,1),\cr A_4 := \mathop{\rm tridiag}(-\alpha,1+2\alpha,-\alpha),\quad A_5 := \mathop{\rm tridiag}\left(-{\alpha\over2},1+\alpha,-{\alpha\over2}\right),\quad B_5 := \mathop{\rm tridiag}\left({\alpha\over2},1-\alpha,{\alpha\over2}\right). \fracstrut\cr } $$

Mit diesen Vektoren $v^i$, $f^i$ und den Matrizen $A_\nu$, $B_5$ schreiben sich jetzt die alle oben angegebenen Verfahren, wie folgt.

1) $v^{i+1}=A_1v^i+\tau f^i$, für $i=0,\ldots,N-1$.

2) $v^{i+1}=2\alpha A_2v^i+v^{i-1}+2\tau f^i$, für $i=1,\ldots,N-1$.

3) $(1+2\alpha)v^{i+1}=2\alpha A_3v^i+(1-2\alpha)v^{i-1}+2\tau f^i$, für $i=1,\ldots,N-1$.

4) $A_4v^{i+1}=v^i+\tau f^i$, für $i=0,\ldots,N-1$.

5) $A_5v^{i+1}=B_5v^i+\tau f^i$, für $i=0,\ldots,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 $\lambda_{1,i} := 1-4\alpha\sin^2{i\pi\over2M}$, für $i=1,\ldots,M-1$.

2) $A_2$ mit den Eigenwerten $\lambda_{2,i} := 4\sin{i\pi\over2M} > 0$, für $i=1,\ldots,M-1$.

3) $A_3$ mit den Eigenwerten $\lambda_{3,i} := 2\cos{i\pi\over M}$, für $i=1,\ldots,M-1$.

4) $A_4$ mit den Eigenwerten $\lambda_{4,i} := 1+4\alpha\sin^2{i\pi\over2M}$, für $i=1,\ldots,M-1$.

5) Für $A_5$ rechnet man

$$ A_5=\mathop{\rm tridiag}\left(-{\alpha\over2},1+\alpha,-{\alpha\over2}\right) ={1\over2}\left[2I+\alpha\mathop{\rm tridiag}(-1,2,-1)\right] $$

und daher $\lambda_{5a,i} := {1\over2}\left(2+4\alpha\sin{i\pi\over2M}\right)$, und

$$ B_5=\mathop{\rm tridiag}\left({\alpha\over2},1-\alpha,{\alpha\over2}\right) ={1\over2}\left[2I-\alpha\mathop{\rm tridiag}(-1,2,-1)\right] $$

mit den Eigenwerten $\lambda_{5b,i} := {1\over2}\left(2-4\alpha\sin{i\pi\over2M}\right)$, für $i=1,\ldots,M-1$.

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

1) Das Matrixpolynom $I\mu-A_1$ hat die Eigenwerte $\mu_i=\lambda_{1,i}$. Diese sind genau dann betragsmässig kleiner eins, falls $\alpha\le1/2$. Wegen $\alpha=\sigma\tau/h^2$, führt dies auf die Begrenzung der Zeitschrittweite $\tau$ zu

$$ \tau\le{1\over2\sigma}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\mu^2-2\alpha A_2\mu-I$ ergeben sich nach Durchmultiplikation mit $D$, wobei $D$ die Transformationsmatrix für $A_2$ ist, also $A_2=D\mathop{\rm diag}(\lambda_{2,i})D^{-1}$, die Eigenwerte des Matrixpolynoms zu

$$ \det D\cdot\det\left(I\mu^2-2\alpha\mathop{\rm diag}(\lambda_{2,i})-I\right) \cdot\det D^{-1}=0, $$

also

$$ \mu_{i,1/2}=\lambda_{2,i}\pm\sqrt{\alpha^2\lambda_{2,i}^2+1}, \qquad i=1,\ldots,M-1. $$

Der Spektralradius ist also für jedes $\alpha$ größer als 1. Das Verfahren ist demzufolge für alle $\tau$ 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\alpha)I\mu^2-2\alpha A_3\mu-(1-2\alpha)I$. Sei $D$ die Transformationsmatrix auf Diagonalgestalt für die Matrix $A_3$. Damit erhält man die $2M-2$ Eigenwerte des Matrixpolynoms als Nullstellen der Gleichung

$$ \det D\cdot\det\left(1+2\alpha)I\mu^2 -2\alpha\cdot2\cos{i\pi\over M}\mu-(1-2\alpha)I\right)\cdot\det D^{-1}=0, $$

also

$$ (1+2\alpha)\mu_i^2-4\alpha\mu_i\cos\varphi-(1-2\alpha)=0, $$

mit $\varphi := i\pi/M$. Mit der Lösungsformel für quadratische Gleichungen der Form $ax^2+bx+c=0$, $a\ne0$, nämlich

$$ x_{1/2}={-b\pm\sqrt{b^2-4ac}\over2a}, $$

erhält man

$$ \mu_{i,1/2}={2\alpha\pm\sqrt{1-4\alpha^2\sin^2\varphi} \over 1+2\alpha}. $$

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

4) Das Matrixpolynom $I\mu-A_4^{-1}$ hat die Eigenwerte $\mu_i=\lambda_{4,i}^{-1} < 1$, für alle $\alpha$, da $\lambda_{4,i} > 1$. Das Verfahren von Crank/Nicolson ist damit für alle $\tau$ und $h$ stabil. Aufgrund der Konsistenz folgt damit die Konvergenz.

5) Für das entsprechende Matrixpolynom korrespondierend zu $A_5\mu=B_5$, ergeben sich die Eigenwerte als Quotient der Eigenwerte von $A_5$ und $B_5$, also

$$ \mu_i = {1-2\alpha\sin(i\pi/2M) \over 1+2\alpha\sin(i\pi/2M) }. $$

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