, 57 min read

Das Fehlerverhalten zusammengesetzer linearer Mehrschrittformeln

Bei zusammengesetzten Verfahren, also Verfahren mit mehr als einer Stufe, besitzt ersteinmal jede Stufe für sich eine eigene Fehlerkonstante im herkömmlichen Sinne. Dennoch zeigt z.B. die zyklische Hintereinanderausführung des impliziten und expliziten Euler-Verfahrens, daß das Einzelverhalten der Stufen nicht unbedingt auch das Gesamtverhalten des Zykluses wiedergibt. Das implizite Euler-Verfahren für sich alleine betrachtet hat die Konvergenzordnung 1, ebenso hat das explizite Euler-Verfahren für sich alleine betrachtet die Konvergenzordnung 1. Das zusammengesetzte Verfahren hat allerdings schon die Konvergenzordnung 2. Es ist nun naheliegend zu fragen, ob noch höhere Konvergenzordnungssprünge möglich sind. Desweiteren wird man für diesen Sprung der Konvergenzordnung eine Erklärung wünschen.

Allerdings wird man nicht in so unstetigen Übergängen denken wollen. Bei den klassischen Verfahren, wie linearen Mehrschrittformeln oder Runge-Kutta-Verfahren, ist bekannt, daß ein sehr genaues Verfahren der Ordnung $p$, sich ähnlich verhält wie ein sehr ungenaues Verfahren der eins höheren Ordnung $p+1$. Man erwartet also vielmehr einen gleitenden Übergang zwischen Verfahren. Paradebeispiel ist hierfür übrigens das $\vartheta$-Verfahren

$$ y_{n+1}-y_n = h(\vartheta f_{n+1}+(1-\vartheta)f_n),\qquad n=0,1,\ldots $$

Für $\vartheta=1/2$ erhält man die implizite Trapezregel mit der Ordnung 2 und in allen anderen Fällen nur Verfahren der Ordnung 1, insbesondere für $\vartheta=0$ ein explizites Verfahren. Es zeigt sich nun, daß der primäre dominante lokale Fehler eine erste Auskunft gibt über das Fehlerverhalten. Verschwindet der primäre dominante lokale Fehler, liegt also der Fall der annullierten Dominanz vor, so gibt der sekundäre dominante lokale Fehler ein weiteres Bild über das Fehlerverhalten.

Zuerst erscheint zur Klarstellung der Bezeichnungen, die Festlegung des Konsistenzbegriffes. Anschliessend werden eine Reihe von zueinander äquivalenten Beschreibungsmöglichkeiten für hohe Konsistenzordnung gegeben. Diese Beschreibungen sind direkt anwendbar zur Berechnung neuer Verfahren. Eine Reihe von Fehlerkonstanten werden miteinander verglichen und Gemeinsamkeiten deutlich gemacht.

1. Konsistenz, Konsistenzordnung und Fehlerkonstanten

Es sei

$$ \alpha_0y_n+\alpha_1y_{n+1}+\cdots+\alpha_ky_{n+k} = h\left(\beta_0f_n+\beta_1f_{n+1}+\cdots+\beta_kf_{n+k}\right), \qquad\alpha_k\ne0, $$

ein lineares $k$-Schrittverfahren. An die Koeffizienten des linearen $k$-Schrittverfahrens sind gewisse Einschränkungen zu stellen, damit die Lösungen, die durch das $k$-Schrittverfahren berechnet werden, auch etwas mit der Lösung der Differentialgleichung zu tun haben. Man unterscheidet zweierlei Bedingungsarten: einmal Konsistenzbedingungen und zum anderen Stabilitätsbedingungen. Es zeigt sich nachher, daß die Stabilitätsbedingung die einschränkendere Bedingung ist. Die Konsistenzbedingungen führen auf lineare Gleichungssyteme. Die Stabilitätsbedingungen führen auf nicht-lineare und Gleichungs- und Ungleichungssysteme.

Die Konsistenzbedingungen sind:

$$ C_{p,k}\cdot\pmatrix{\alpha\cr \beta\cr} = \pmatrix{A{\mskip 3mu}|{\mskip 3mu}B\cr} \pmatrix{\alpha\cr \beta\cr} = 0. $$

Die Konsistenzmatrix $C_{p,k}$ für lineare $k$-Schrittverfahren der Konsistenzordnung $p$ lautet im Verein mit dem Koeffizientenvektor des Verfahrens und der entsprechenden Bedingung

$$ \left( \begin{array}{ccccc|ccccc} 1& 1& 1& 1& \ldots& 1&& 0& 0& 0& 0& \ldots& 0\cr 0& 1& 2& 3& \ldots& k&&-1& -1& -1& -1& \ldots& -1\cr 0& 1^2& 2^2& 3^2& \ldots& k^2&& 0& -2\cdot1& -2\cdot2& -2\cdot3& \ldots& -2\cdot k\cr 0& 1^3& 2^3& 3^3& \ldots& k^3&& 0& -3\cdot1^2& -3\cdot2^2& -3\cdot3^2& \ldots& -3\cdot k^2\cr \vdots & \vdots & \vdots & \vdots & \ddots&\vdots && \vdots & \vdots & \vdots & \vdots & \ddots&\vdots\cr 0& 1^p& 2^p& 3^p& \ldots& k^p&& 0& -p1^{p-1}& -p2^{p-1} & -p3^{p-1}& \ldots& -pk^{p-1}\cr \end{array} \right) \cdot \pmatrix{\alpha_0\cr \alpha_1\cr \vdots\cr \alpha_{k-1}\cr \alpha_k\cr \beta_0\cr \beta_1\cr \vdots\cr \beta_{k-1}\cr \beta_k\cr} = \pmatrix{0\cr 0\cr 0\cr 0\cr \vdots\cr 0\cr}. $$

Beispielsweise lautet die Konsistenzmatrix $C_{p,k}$ bei spezieller Wahl von $p$ und $k$, wie folgt:

1. Beispiel: Konsistenzmatrix $C_{p+1,k}$ für $p=3$ und $k=3$:

$$ C_{4,3} = \pmatrix{ 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0\cr 0 & 1 & 2 & 3 & -1 & -1 & -1 & -1\cr 0 & 1 & 4 & 9 & 0 & -2 & -4 & -6\cr 0 & 1 & 8 & 27 & 0 & -3 & -12 & -27\cr 0 & 1 & 16 & 81 & 0 & -4 & -32 & -108\cr } $$

und für $p=5$, $k=5$ lautet $C_{p+1,k}$ mithin

$$ C_{6,5} = \pmatrix{ 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\cr 0 & 1 & 2 & 3 & 4 & 5 & -1 & -1 & -1 & -1 & -1 & -1\cr 0 & 1 & 4 & 9 & 16 & 25 & 0 & -2 & -4 & -6 & -8 & -10\cr 0 & 1 & 8 & 27 & 64 & 125 & 0 & -3 & -12 & -27 & -48 & -75\cr 0 & 1 & 16 & 81 & 256 & 625 & 0 & -4 & -32 & -108 & -256 & -500\cr 0 & 1 & 32 & 243 & 1024 & 3125 & 0 & -5 & -80 & -405 & -773 & -3125\cr 0 & 1 & 64 & 729 & 4096 & 15625 & 0 & -6 & -192 & -1458 & -6144 & -18750\cr } $$

2. Ein lineares $k$-Schrittverfahren mit mindestens der Konsistenzordnung $p$ muß also im Kern der Matrix $C_{p,k}\in\mathbb{Z}^{(p+1)\times(2k+2)}$ liegen, also

$$ \pmatrix{\alpha\cr \beta\cr}\in\ker C_{p,k}. $$

Erweitert man die Matrix $C_{p,k}$ in offensichtlicher Weise unten um eine weitere Zeile und damit zur Matrix $C_{p+1,k}\in\mathbb{Z}^{(p+2)\times(2k+2)}$, und ist dann das Matrix-Vektor Produkt nicht mehr der Nullvektor, so hat das Verfahren die genaue Konsistenzordnung $p$, und die von Null verschiedene Komponente des Ergebnises ist ein unskalierter Fehlerfaktor $c_{p+1}$. Wenn auf die Unskalierung besonders hingewiesen werden soll auch $\lambda c_{p+1}$, mit $\lambda\in\{\alpha_k,\beta_k,\sigma(1),\ldots\}$. Man hat also

$$ C_{p+1,k}\pmatrix{\alpha\cr\beta\cr}=\pmatrix{0\cr\vdots\cr0\cr c_{p+1}\cr}. $$

Den Wert $c_{p+1}$ teilt man jetzt noch durch $(p+1)!$, aus später ersichtlichen Gründen, die mit einer Taylorentwicklung zu tun haben.

Übliche Skalierungsgrößen sind nun $\alpha_k$ oder $\sum_{i=0}^k \beta_i$. Skaliert man mit der letztgenannten Summe, so heiße der resultierende Faktor auch Henrici's Fehlerkonstante. Sie tritt in natürlichster Art und Weise auf bei der Behandlung von differential-algebraischen Gleichungen, die man mit Verfahren löst, wie man sie bei gewöhnlichen Differentialgleichungen einsetzt.

Bibliographisch: Henrici, Peter Karl Eugen (1923--1987).

3. Die Minus-Zeichen in der Konsistenzmatrix $C_{p,k}\in\mathbb{Z}^{(p+1)\times(k+1)(m+1)}$ (bisher $m=1$) tauchen nicht mehr auf, wenn man statt

$$ \sum \alpha_i y_{n+i} = h \sum \beta_i \dot y_{n+i} $$

einfach

$$ \sum \alpha_i y_{n+i} + h\beta_i \dot y_{n+i} + h^2\gamma_i \ddot y_{n+i} = 0 $$

schreibt (oben für $m=2$).

4. Bei Diskretisierungen zur Lösung von Gleichungen der From $F(t,y,\dot y)=0$, hier insbesondere linearen Mehrschrittverfahren, wird man unmittelbar dazu geführt, die Gleichung

$$ {1\over h}\sum_{i=0}^k\alpha_iy_i = \sum_{i=0}^k\beta_i\dot y_i $$

wie folgt zu interpretieren. Die linke Summe stellt eine Näherung für die Ableitung $\dot y$ an einer hier nicht weiter interessierenden Zwischenstelle dar. Die rechte Summe ist genau dann ein gewichtetes Mittel der Werte $\dot y_i$, wenn sich die $\beta_i$-Summanden genau zu eins aufsummieren. Das letzte kann man aber stets erreichen durch geeignete Vormultiplikation der obigen Gleichung mit dem Kehrwert der Summe $\sum_{i=0}^k\beta_i$. Aufgrund der sofort offensichtlichen Linearität der Konsistenzbedingungen, erscheint die entsprechende Summe dann auch in dem Fehlerfaktor $c_{p+1}$.

Genauso ist aber auch

$$ \sum_{i=0}^{k-1} \beta_i \dot y_{n+i} + \sum_{i=0}^k \alpha_i y_{n+i} = \dot y_{n+k} $$

interpretierbar als Näherung eben an der Stelle $t_{n+k}$.

5. Es gibt mehrere weitere Möglichkeiten die Fehlerkonstante von Henrici abzuleiten.

Insbesondere durch Überlegungen bzgl. des Einflußes der lokalen Fehler auf den globalen Fehler. Dies sind die Überlegungen, wie man sie in den Aufsätzen von Skeel (1976) und Albrecht (1985) findet. In dem letztgenannten Aufsatz werden diese Überlegungen in allgemeinster Form unter Berücksichtigung von mehrfachen Eigenwerten $\mu=1$ durchgeführt. Ähnliche Erwägungen werden in der Dissertation von Tischer (1983) durchgeführt. U.a. findet man eine Darstellung und Ableitung der Fehlerkonstanten von Henrici in dem Buche von Hairer/Wanner/Nørsett (1987) und natürlich in dem Buche von Henrici (1962) selber. Die Fehlerkonstante von Henrici ist nicht originär von Henrici erfunden worden. Sie taucht ebenfalls bei zahlreichen anderen Autoren auf, wie z.B. bei Hull/Newberry (1961).

Bibliographisch: Henrici, Peter Karl Eugen (1923--1987), Hairer, Ernst (*1949), Wanner, Gerhard (*1942), Nørsett, Syvert Paul, Thomas E. Hull, A.C.R. Newberry, Robert David Skeel, Peter Albrecht.

Peter E. Tischer: "The Cyclic Use of Linear Multistep Formulas for the Solution of Stiff Differential Equations", Ph.D. Thesis, Department of Computer Science, Monash University, Clayton, Victoria, Australia, August 1983, x+180 pages.

Bei linearen Verfahren, die der starken Wurzelbedingung genügen, ergibt sich die Fehlerkonstante $C_A$ zu

$$ C_A = {v{\mskip 3mu}\gamma\over vA_1w}, \qquad \left\{\eqalign{v{\mskip 3mu}(A_0+A_1)&=0,\quad v\ne0\cr (A_0+A_1)w&=0.\quad w\ne 0\cr}\right. $$

Hierbei ist also $v$ der Linkseigenvektor zum Matrixpolynom $A_1\mu+A_0$ zum Eigenwert $\mu=1$ und $w$ ist der Rechtseigenvektor zum gleichen Matrixpolynom und gleichem Eigenwert. Um eine eindeutige Konstante zu erhalten, wird der letzte Vektor $w$ noch normiert z.B. durch die Norm

$$ \left\|w\right\| = {1\over s}\sqrt{w_1^2+\cdots+w_s^2}. $$

Das starke Wurzelkriterium garantiert, daß $vA_1w$ nicht verschwindet. Genauer: die Tatsache, daß $\mu=1$ einziger Eigenwert inklusive Multiplizität ist, garantiert das Nichtverschwinden.

2. Die Anwendung linearer Mehrschrittverfahren bei DAE

1. Um mit Einbein-Verfahren differential-algebraische Gleichungen der Form

$$ \eqalign{ F(t,\dot x, x, y) &= 0,\cr G(t,x,y) &= 0,\cr } $$

zu lösen, rechnet man

$$ \eqalign{ F(\tau_n,{1\over h_n}\rho_nx_n, \sigma_nx_n,\sigma_ny_n) &= 0,\cr G(t_n,x_n,y_n) &= 0.\cr } $$

Die differential-algebraische Gleichung enthalte vermittels $G$ eindeutig als rein algebraische Restriktionen identifizierbare Gleichungen. Das Einbein-Verfahren für gewöhnliche Differentialgleichungen der Form $\dot x=f(t,x)$, ist gegeben durch

$$ {1\over h_n}\sum_{\nu=0}^k \alpha_{\nu,n}x_{n-k+\nu} - f\left(\tau_n,{\mskip 3mu}\sum_{\nu=0}^k \beta_{\nu,n}x_{n-k+\nu}\right)=0, $$

und

$$ \tau_n=\sum_{\nu=0}^k \beta_{\nu,n} t_{n-k+\nu} \qquad\hbox{wobei}\quad \sum_{\nu=0}^k\beta_{\nu,n}=1,\enspace\forall n $$

Die Operatoren $\rho_n$ und $\sigma_n$ sind wie üblich

$$ \rho_n z_n = \sum_{\nu=0}^k \alpha_{\nu,n} z_{n-k+\nu},\qquad \sigma_n z_n = \sum_{\nu=0}^k \beta_{\nu,n} z_{n-k+\nu}. $$

Ebenfalls mögliche Diskretisierungen des Gleichungspaares $F(t,\dot x,x,y)=G(t,x,y)=0$ sind

$$ \eqalign{ F(\tau_n,{1\over h}\rho_nx_n, \sigma_nx_n,\sigma_ny_n) &= 0,\cr G(\tau_n,\sigma_nx_n,\sigma_ny_n) &= 0,\cr } $$

was sich besonders dann anbietet, falls bei einer differential-algebraischen Gleichung, die Trennung zwischen “reiner” Differentialgleichung und rein algebraischen Restriktionen, nicht klar zutage tritt. Dies ist beispielsweise bei dem Problem P9 der Fall.

2. Für lineare Mehrschrittverfahren der Form

$$ {1\over h_n}\left(\sum_{\nu=0}^k \alpha_{\nu,n}x_{n-k+\nu}\right) - \sum_{\nu=0}^k \beta_{\nu,n}\dot x_{n-k+\nu} = 0, \qquad n=k,k+1,\ldots, $$

stets mit der Normierung $\sum_{\nu=0}^k\beta_{\nu,n}=1$, gewinnt man eine Diskretisierung für die Gleichung $F(t,\dot x,x,y)=G(t,x,y)=0$, wie folgt. Man fügt temporär eine zusätzliche Variable $w:=\dot x$ dem System $F=G=0$ hinzu und erhält dann nach Diskretisierung, unter Beachtung von $w_n:=\dot x_n$, das vergrößerte System

$$ \eqalign{ F(t_n,w_n,x_n,y_n) &= 0,\cr G(t_n,x_n,y_n) &= 0,\cr {1\over h_n}\rho_nx_n - \sigma_nw_n &= 0,\cr } $$

mit den Unbekannten $x_n$, $y_n$ und $w_n$. Die letzte Gleichung ist allerdings sehr leicht nach $w_n$ aufzulösen. Man erhält

$$ w_n = {1\over\beta_{k,n}}\left({1\over h_n}\rho_nx_n-\sigma'_nw_n\right), $$

mit

$$ \sigma'_n = \sum_{\nu=0}^{k-1} \beta_{\nu,n}w_{n-k+\nu}. $$

Der Operator $\sigma'_n$ wirkt also lediglich auf schon zurückliegende berechnete Werte. Einsetzen des so aufgelösten $w_n$, führt auf

$$ \eqalign{ F(t_n,{1\over\beta_{k,n}}\left({1\over h_n}\rho_nx_n-\sigma'_nw_n\right), x_n,y_n) &= 0,\cr G(t_n,x_n,y_n) &= 0,\cr } $$

zur Bestimmung von $(x_n,{\mskip 3mu}y_n)$.

3. Prinzipiell sind auch explizite Operatoren $(\rho_n,\sigma_n$) denkbar. Entsprechend erhält man dann

$$ \eqalign{ F (t_{n-1}, {1\over\beta_{k-1,n}}\left({1\over h_n}\rho_nx_n- \sigma''_nw_n\right),x_{n-1}, y_{n-1}) &= 0,\cr G (t_{n-1}, x_{n-1}, y_{n-1}) &= 0,\cr } $$

wobei $\sigma''_n$ noch einen Term weniger enthält. Der Vorteil der leichteren Auflösbarkeit oder gar der Vermeidung von Nichtlinearitäten solcher Diskretisierungen, ist jedoch im Falle von differential-algebraischen Gleichungen nicht mehr vorhanden. Es muß in jedem Falle in jedem Zeitschritt ein i.d.R. nicht-lineares Gleichungssystem gelöst werden.

Zu dieser Diskretisierung schreibt Liniger (1979):

although to the author's knowledge, thus far such methods have not been used practically …

Der Einsatz von linearen Mehrschrittverfahren und Einbein-Verfahren, und die Ableitung von Diskretisierungen bei differential-algebraischen Gleichungen, wird bei Liniger (1979) untersucht. Insbesondere Fragen der Konsistenzordnung der Diskretisierung findet man bei Liniger (1979). Stabilitätsuntersuchungen sind für differential-algebraische Gleichungen schwieriger und werden daher auch dort nicht behandelt. Einen Konvergenzbeweis für die BDF findet man bei Petzold/Lötstedt (1986a). Eine gute übergreifende Darstellung bietet Griepentrog/März (1986).

Bibliographisch: Werner Liniger (1927--2017), Linda Ruth Petzold (*1954), Eberhard Griepentrog (1933--2023), Roswitha März (*1940), Per Lötstedt.

3. Mehrere Charakterisierungen der Konsistenzordnung

An dieser Stelle sollen eine Reihe von gleichwertigen Charakterisierungen angegeben werden, die garantieren, daß ein lineares Mehrschrittverfahren der Form

$$ \sum_{i=0}^k \alpha_iy_{n+i} = h\sum_{i=0}^k \beta_if_{n+i}, $$

mindestens die Konsistenzordnung $p$ hat. Unter Umständen hat das Verfahren sogar eine noch höhere Ordnung.

1. Es sei

$$ L(y,t_0,h) = \sum_{i=0}^k \left(\alpha_i y(t_0+ih) - h\beta_i\dot y(t_0+ih)\right). $$

und

$$ \rho(\mu) := \sum_{i=0}^k \alpha_i\mu^i,\qquad\hbox{ferner}\qquad \sigma(\mu) := \sum_{i=0}^k \beta_i\mu^i. $$

2. Satz: Nun sind die folgenden 9 Aussagen paarweise zueinander äquivalent.

(1) $C_{p,k}{\alpha\choose\beta}= {\bf0}$, also $(\alpha,\beta)^\top\in\ker C_{p,k}$.

(2) $\sum_{i=0}^k \alpha_i i^q = q\sum_{i=0}^k \beta_i i^{q-1}$, für $q=0,1,\ldots,p$.

(3) $\rho(e^h)-h\sigma(e^h)={\cal O}(h^{p+1})$, für $h\to 0$.

(4) $\zeta=1$ ist mindestens $p$-fache Nullstelle der Funktion

$$ h\mapsto{\rho(\zeta)\over\ln\zeta}-\sigma(\zeta), $$

also $\rho(\zeta)/\ln\zeta=\sigma(\zeta)+{\cal O}\bigl((\zeta-1)^p\bigr)$, für $\zeta\to 1$.

(5) $L(f,t,h)={\cal O}(h^{p+1}),\quad\forall f\in C^{p+2}(G,\mathbb{R}{})$.

(6) Die Monome bis zum Grade $p$ liegen im Kern des $h=1$ Schnittes von $L$, also $L(t^i,t_0,1)=0$, für $i=0,1,\ldots,p$.

(7) $L(f,t_0,h)={\cal O}(h^{p+1})$, für die spezielle Funktion $t\mapsto f(t)=e^t$.

(8) $y(t_0+kh)-y_k={\cal O}(h^{p+1})$, falls man die Startwerte als exakte Werte wählt.

(9) $L(y,t_0,h)=c_{p+1}h^{p+1}y^{(p+1)}(t_0)+{\cal O}(h^{p+2})$ mit $c_{p+1}=\sum_{i=0}^k\bigl(\alpha_ii^{p+1}-(p+1)\beta_ii^p\bigr)/(p+1)!$.

Diese Liste liesse sich fortsetzen. Einige der Nummern sind lediglich Umformulierungen anderer Nummern. Dennoch ist es gelegentlich nützlich eine der möglichen Formeln in der oben zitierten Form parat zu haben. Die Bedingung der Konsistenz und die Konsistenzordnung ist unabhängig von der Entwicklungsstelle der Taylorentwicklung. Zur praktischen Berechnung von Mehrschrittformeln, die zuerst nicht unbedingt stabil sein müssen, wählt man häufig (1), für den Beweis der Dahlquistbarriere werden (3) und (4) herangezogen. Einzelne Eigenschaften der obigen Angaben tragen in der Literatur häufig auch gesonderte Namen. Für die Beweise sei beispielsweise verwiesen auf Hairer/Wanner/Nørsett (1987) oder auch {Werner, Helmut}{Arndt, Herbert}Werner/Arndt (1986).

Wählt man, wie bei (8) die Startwerte exakt, also $y_0=y(t_0)$, $y_1=y(t_0+h)$ und so fort bis $y_{k-1}=y(t_0+(k-1)h)$, so gilt für den dann entstehenden Fehler zwischen dem nun durch die Formel bestimmten Wert $y_k$ und dem exakten Wert $y(t_0+kh)$, die Beziehung

$$ y(t_0+kh)-y_k = {1\over\alpha_k}W^{-1}L(y,t_0,h). $$

Hierbei ist $W=I-h\gamma J$, $\gamma=\beta_k/\alpha_k$ und $J$ ist die Jacobimatrix von $f$, deren Zeilen ausgewertet wurden an Stellen auf der Strecke zwischen $y(t_k)$ und $y_k$. Insbesondere gilt nicht notwendig $J=f_y(y_k)$, jedoch gilt dies näherungsweise. Im Falle steifer Differentialgleichungen und differential-algebraischer Gleichungen wird auch tatsächlich $W^{-1}L(y,t_0,h)$, bzw. eine Näherung hiervon, zum Gebrauch als Fehlerkonstante empfohlen, man vgl. hier Petzold (1982) und die dort angeführte Literatur. Im Falle von expliziten Formeln, also $\beta_k=0$, ist natürlich $W=I$.

4. Die erste Dahlquist-Barriere

You know, I am a multistep-man … and don't tell anybody, but the first program I wrote for the first Swedish computer was a Runge-Kutta code …
(G. Dahlquist, 1982, after some glasses of wine; printed with permission), Hairer/Wanner/Nørsett (1987)

Das Bemerkenswerte am Hopfschen Problem ist, daß eine einfach zu formulierende algebraische Frage eine einfache algebraische Antwort findet, daß indessen zur Lösung nichttriviale Methoden der Topologie erforderlich sind; hier erscheint zum ersten Mal der “topologische Stachel im Fleisch der Algebra”, der bis auf den heutigen Tag von vielen Algebraikern als so schmerzhaft empfunden wird.
R. Remmert, M. Koecher (1983)

Bibliographisch: Hairer, Ernst (*1949), Wanner, Gerhard (*1942), Nørsett, Syvert Paul, Remmert, Reinhold (1930--2016), Koecher, Max (1924--1990),

Der Kern der Konsistenzmatrix $C_{p,k}$ hat eine Fülle von sehr speziellen Eigenschaften, die hier zusammengestellt werden sollen. Anhand des sehr strukturierten Aufbaus dieser Matrix sind einige dieser Eigenschaften nicht verwunderlich. Insbesondere einige Symmetrieeigenschaften der Matrix übertragen sich auf den Kern. Nimmt man als Nebenbedingung noch das Erfülltsein der Wurzelbedingung hinzu, so gelten weitreichende Konsequenzen, hier die beiden Dahlquist-Barrieren. Diese beiden Barrieren sind auch eine der Gründe für das verstärkte Interesse an zusammengesetzten Verfahren.

1. Definition: Ein lineares $k$-Schrittverfahren heißt symmetrisch, falls

$$ \alpha_i=-\alpha_{k-i}\qquad\hbox{und}\qquad\beta_i=\beta_{k-i} \qquad\hbox{für alle}\quad i=0,\ldots,k. $$

2. Beispiel: Das Milne-Simpson-Verfahren $y_{n+1}-y_{n-1}=h\cdot(f_{n+1}+4f_n+f_{n-1})/3$ ist symmetrisch, während hingegen alle Adams-Verfahren, implizit oder explizit, nicht symmetrisch sind. Die Nullstellen des charakteristischen Polynoms $\rho$ des Milne-Simpson-Verfahren liegen bei $(+1)$ und bei $(-1)$.

Für symmetrische lineare Mehrschrittverfahren gilt $\rho(\mu)=-\mu^k\rho(1/\mu)$, aufgrund der obigen Definition. Mit $\mu$ ist auch gleichzeitig $1/\mu$ Nullstelle von $\rho$.

Für ein stabiles lineares Mehrschrittverfahren liegen somit alle Wurzeln auf dem Einheitskreis und sind einfach. An solchen Verfahren ist man eher weniger interessiert, da bei Schrittweiten $h\ne0$, diese Wurzeln vom Betrage $1$ aus dem Einheitskreis nach aussen wandern. Allerdings hängt bei Schrittweiten $h\ne0$ dieses Verhalten stark vom Prädiktor ab. In der Kombination als Picard-Prädiktor-Korrektor-Verfahren oder als Stufen zyklischer Verfahren haben jedoch symmetrische lineare Mehrschrittverfahren gegenüber anderen Verfahren keine Nachteile.

Ein Verfahren der Konsistenzordnung $p$ liefert für eine Differentialgleichung mit Polynomen des Grades $p$ als Lösung, unter völliger Vernachlässigung von Rundungsfehlern, die exakte Lösung. Ist das Verfahren zusätzlich noch stabil, so konvergiert das Verfahren.

3. Satz: Es gilt

  1. Es existiert kein lineares $k$-Schrittverfahren der Konsistenzordnung $2k+1$.
  2. Es gibt genau ein explizites $k$-Schrittverfahren der Konsistenzordnung $2k-1$. Aus $\beta_k=0$ folgt also automatisch $p\le2k-1$.
  3. Es gibt genau ein implizites lineares $k$-Schrittverfahren der Konsistenzordnung $2k$.
  4. Dieses eindeutig bestimmte implizite lineare $k$-Schrittverfahren der maximalen Konsistenzordnung $2k$, ist symmetrisch.
  5. Symmetrische lineare Mehrschrittverfahren haben immer eine gerade Konsistenzordnung. Dies heißt, hat man von einem symmetrischen Mehrschrittverfahren die Konsistenzordnung $2\nu-1$ nachgewiesen, so hat das Verfahren auch schon automatisch die eins höhere Konsistenzordnung $2\nu$.
  6. Zu vorgegebenen Polynomgraden $k_\alpha$ und $k_\beta$, mit $k_\beta\le k_\alpha$, kann man Polynome $\rho$ und $\sigma$ finden, mit $\mathop{\rm grad}\rho=k_\alpha$ und $\mathop{\rm grad}\sigma=k_\beta$, die ein lineares $k=k_\alpha$-Schrittverfahren festlegen und zwar mit der Konsistenzordnung $1+k_\alpha+k_\beta$.
  7. Der Rang der Konsistenzmatrix $C_{p,k}$ ist für alle $p$ und $k$ maximal. Es gilt also
$$ \mathop{\rm rank} C_{p,k}=\min(p+1,k+1),\qquad\forall p,k\in\mathbb{N}. $$

Es war $C_{p,k}\in\mathbb{Z}^{(p+1)\times(2k+2)}$ und mit dem ^{Dimensionssatz für lineare Gleichungssysteme} ergibt sich daher $\dim\ker C_{p,k}=(2k+2)-\min(p+1,k+1)$.
Den führenden Koeffizienten $\alpha_k$ kann man als Normierungsfaktor auffassen, und damit gibt es genau eine $2k-p$ parameterabhängige Schar von linearen $k$-Schrittverfahren, falls $p\ge k$. Von dieser Schar ist aber nach der ersten Dahlquist-Barriere nur ein kleiner Teil als einstufiges Verfahren konvergent.

Es zeigt sich, daß die maximal erreichbare Konsistenzordnung eines linearen Mehrschrittverfahrens der Form $\sum_{i=0}^k (\alpha_i y_{n+i} - h\beta_i f_{n+i})$ durch die Wurzelbedingung an $\rho$, begrenzt ist. Ist das charakteristische Polynom $\rho$ stabil, so ist dies gleichzeitig ein Hemmschuh für die maximal erreichbare Konsistenzordnung, d.h. Konsistenzordnung und Stabilität beissen sich gegeneinander, oder in noch anderer Formulierung, mehr geometrischer Sprechweise:

Die Menge aller derjenigen Vektoren

$$ \left(\alpha_0,\ldots,\alpha_\kappa,\beta_0,\ldots,\beta_\kappa\right) \in\mathbb{R}^{2\kappa+2}, $$

welche zu einem stabilen Polynom $\rho(\mu)=\sum_{i=0}^\kappa \alpha_i\mu^i$ führen, beschreibbar durch Ungleichungen vom Routh-Hurwitz-Typ, {Routh, E.J.}{Hurwitz, Adolf (1859--1919} stellt einen nicht-linearen Kegel im $\mathbb{R}^{2\kappa+2}$ dar. Die lineare Untermannigfaltigkeit beschrieben durch $C_{p,\kappa}{\alpha\choose\beta}=0$ schneidet den nicht-linearen Kegel überhaupt nicht für $p>\kappa+2$ und berührt ihn für $p=\kappa+1+{1\over2}\left[1+(-1)^\kappa\right]$. Es gilt der

4. Satz: (^{Erste Dahlquist-Barriere}) Das lineare $k$-Schrittverfahren sei wenigstens konsistent der Ordnung $1$ und das charakteristische Polynom $\rho$ erfülle die Wurzelbedingung.

(a) Dann unterliegt die Konsistenzordnung und damit gleichzeitig einhergehend die Konvergenzordnung $p$, einer Schranke nach oben wie folgt

$$ p\le\cases{ k+2, & falls $k$ gerade ist;\cr k+1, & falls $k$ ungerade ist;\cr k, & falls $\beta_k/\alpha_k\le0$, insbesondere falls das Verfahren explizit ist.\cr } $$

(b) Weiterhin gilt: Stabile lineare Mehrschrittverfahren der maximalen Konsistenzordnung $k+2$ ($k$ also gerade) sind symmetrisch. Bei einem solchen linearen $k$-Schrittverfahren mit geradem $k$, besitzt $\rho$, also wie oben bemerkt, die Wurzeln $(+1)$, $(-1)$ und $(k-2)/2$ Paare verschiedener unimodularer Wurzeln. Zu jedem $\rho$ mit dieser Eigenschaft der Wurzeln, existiert genau ein einziges $\sigma$, sodaß $(\rho,\sigma)$ ein lineares $k$-Schrittverfahren der maximalen Konsistenzordnung $k+2$ ist.

5. Wegen dem großen Interesse, dem man diesem Satz beimißt, findet man den sehr schönen, funktionentheoretischen Beweisgang in mehreren ($\ge8$) Büchern. Verwiesen sei hier lediglich auf die schon mehrfach zitierten beiden Bücher von Werner/Arndt (1986) und Hairer/Wanner/Nørsett (1987) sowie Werner/Schaback (1972). Man vgl. auch die Bücher von Gear (1971c), Henrici (1962) und Stetter (1973) (unvollständiger Beweis).

Darüberhinaus ist der Satz erweitert worden, um für eine noch größere Klasse von Verfahren seine entsprechende Gültigkeit zu behalten, man vgl. hier Jeltsch/Nevanlinna (1986). Eine Kurzübersicht über Ordnungsbeschränkungen findet man in dem Tagungsaufsatz von Wanner (1987). G. Dahlquist bewies diesen Satz 1956.

Bibliographisch: Jeltsch, Rolf, Nevanlinna, Olavi, Hairer, Ernst (*1949), Wanner, Gerhard (*1942), Nørsett, Syvert Paul, Dahlquist, Germund, Stetter, Hans Jörg (*1930), Schaback, Robert (*1945), Werner, Helmut (1931--1985), Arndt, Herbert, Gear, Charles William (1935--2022), Henrici, Peter Karl Eugen (1923--1987).

6. Konsequenzen dieser ersten Dahlquist-Barriere sind:

6.1. Es gibt kein stabiles lineares $3$-Schrittverfahren mit der Konsistenzordnung 6. Es gibt allerdings mehrere lineare zyklische Mehrschrittformeln der Konvergenzordnung 6 mit nur 3 Startwerten, nämlich die Verfahren DH2 und DH3.{Donelson III, John}{Hansen, Eldon} Allerdings erhält man die Lösungswerte stets im “3er-Pack”. Die letzte Stufe dieser beiden dreistufigen Zyklen benutzt für sich alleine nur 3 Startwerte, davon stammen zwei aus dem aktuellen Zyklus, insgesamt jedoch hat man dann 6 Lösungswerte mit äquidistanten Gitterabstand vorliegen.

6.2 Die Adams-Moulton-Verfahren

$$ y_{n+1}-y_n=h\sum_{i=0}^\kappa \beta_i f_{n+i-\kappa} $$

haben die Konsistenzordnung $\kappa$ und die Konvergenzordnung $\kappa+1$, welches für ungerade $\kappa$ von keinem anderem konvergenten linearen einstufigen Mehrschrittverfahren überboten werden kann.

Als ein Beispiel für die Verallgemeinerungen, sei der Satz von Reimer aus dem Jahre 1968 angegeben, ohne Beweis. Während die erste Dahlquist-Barriere lediglich für lineare Mehrschrittverfahren galt, verallgemeinerte der Satz von Reimer, 12 Jahre später nach der ersten Dahlquist-Barriere, den Sachverhalt auf lineare Verfahren mit beliebig vielen Ableitungen.

Bibliographisch: Reimer, Manfred.

7. Satz: (Satz von Reimer über den maximalen Grad stabiler Differenzenformeln.) Die lineare Differenzenform

$$ L(h,y,\ldots,y^{(m)}) := \sum_{\nu=0}^\kappa \sum_{\mu=0}^m h^\mu \alpha_\nu^{(\mu)} y_\nu^{(\mu)}, \qquad\hbox{mit}\quad\kappa\ge1, m\ge1 $$

und der Normierung $\alpha_\kappa^{(0)}=1$ habe den Grad $p$ und sei stabil. Dann gilt

$$ p \le N+{1\over2}\left[1+(-1)^{N+1}\right], \qquad\hbox{mit}\quad N=(\kappa+1)m. $$

Die Schranke wird für jedes Paar $(\kappa,m)\in\mathbb{N}^2$ angenommen. Weiterhin wird der maximal mögliche Grad $p=N+1$ genau dann erreicht, wenn $N=(\kappa+1)m$ ungerade ist und das Verfahren symmetrisch ist.

Beweis: siehe Reimer (1982).     ☐

8. Die hier häufig auftauchenden BDF$i$ haben das folgende Stabilitätsverhalten. Die BDF$i$ ist ein lineares $i$-Schrittverfahren der Konsistenzordnung $i$ und definiert durch die Formeln

$$ \sum_{\nu=1}^i {1\over\nu}\nabla^\nu y_{n+1} = hf_{n+1},\qquad n=0,\ldots . $$

9. Satz: Die BDF$i$ sind stabil für $i\in\{1,\ldots,6\}$ und instabil für alle $i\ge7$.

Den funktionentheoretischen Beweis findet man in dem Buche von Hairer/Wanner/Nørsett (1987). Die Tatsache, daß die BDF$i$ für alle $i\ge7$ instabil sind, wurde zuerst 1972 von C.W. Cryer bewiesen. Drei Jahre später erschien ein alternativer Beweis von D.M. Creedon und J.J.H. Miller (alle BIT). Der Satz gilt nicht für zusammengesetzte Verfahren, wo die BDF$i$ als Stufen auftauchen, wie z.B. das zyklische Verfahren siebenter Ordnung von Tendler (1973) deutlich macht. Man siehe auch Tendler/Bickart/Picel (1978).

Bibliographisch: Colin Walker Cryer, Theodore A. Bickart (1936--2023), obituary.

Joel Marvin Tendler: "A Stiffly Stable Integration Process Using Cyclic Composite Methods", Ph.D. Diss., Syracuse University, Syracuse, New York, 26.Feb.1973, viii+iv+172 pages.

5. Die zweite Dahlquist-Barriere

1. Die erste Dahlquist-Barriere schränkt die höchstmögliche Konvergenzordnung ein, aber allerdings nicht allzu drastisch. Sehr drakonisch wird jedoch die Vielfalt $A$-stabiler linearer Mehrschrittverfahren durch die zweite Dahlquist-Barriere eingeengt. Es stellt sich heraus, daß man über die Ordnung 2 grundsätzlich nicht hinaus kommt. Hier gilt dann also

2. Satz: Zweite Dahlquist-Barriere.

(1) Ein lineares $A$-stabiles Mehrschrittverfahren ist stets implizit und besitzt höchstens die Konsistenzordnung und damit die Konvergenzordnung 2.

(2) Unter allen linearen $A$-stabilen Mehrschrittverfahren der Konsistenzordnung 2, besitzt die Trapezregel $y_{n+1}=y_n+h\cdot(f_{n+1}+f_n)/2$ ($\vartheta=1/2$-Einschrittverfahren) die kleinstmögliche Fehlerkonstante.

Die BDF2 ist das einzige lineare $2$-Schrittverfahren, welches $A_\infty^0$-stabil ist.

Für den funktionentheoretischen Beweis sei auf die Originalarbeit von Dahlquist (1963) verwiesen, Germund G. Dahlquist (1925--2005). Der Beweis beruht maßgeblich auf den folgenden beiden Sachverhalten und dem Vorzeichenverhalten gewisser Terme.

3. Lemma: Ein $k$-Schrittverfahren ist $A$-stabil genau dann, wenn $\rho(\mu)/\sigma(\mu)$ das Äußere des Einheitskreises auf die komplexe linke Halbebene abbildet, also $\mathop{\rm Re}\nolimits \left[\rho(\mu)/\sigma(\mu)\right]<0$, für $\left|\mu\right|>1$.

4. Satz: Satz von Riesz-Herglotz.

Voraussetzungen: Die Funktion $\varphi(z)$ sei holomorph für $\mathop{\rm Re}\nolimits z>0$. Desweiteren gelte für alle Argumente in der rechten Halbebene $\mathop{\rm Re}\nolimits \varphi(z)\ge0$, also

$$ \mathop{\rm Re}\nolimits z\gt 0 {\mskip 5mu}\Longrightarrow{\mskip 5mu} \mathop{\rm Re}\nolimits \varphi(z)\ge0. $$

Ferner genüge $\varphi$ auf der positiven reellen Achse der Beschränktheitsbedingung $\sup\left\{\left|x\varphi(x)\right|: 0<x<\infty\right\}$.

Behauptung: $\varphi$ hat für alle Argumente aus der rechten Halbebene die Darstellung

$$ \varphi(z) = \int_{-\infty}^\infty {d\omega(t)\over z-it}, \qquad\hbox{für}\quad\mathop{\rm Re}\nolimits z\gt 0, $$

wobei $\omega(t)$ beschränkt und nicht fallend ist.

Bibliographisch: Gustav Ferdinand Maria Herglotz (1881--1953), Friedrich Riesz (1880--1956).

Aus dieser einschränkenden Begrenzung der Konvergenzordnung, beziehen Begriffe wie $A[\alpha]$-, $S[\delta]$-Stabilität, etc., überhaupt ihre Berechtigung. Gäbe es $A$-stabile lineare $k$-Schrittverfahren beliebig hoher Ordnung, so wären dies die idealen Verfahren zur Lösung steifer Differentialgleichungen, unter der Voraussetzung daß nicht andere Eigenschaften, wie z.B. Schrittzahl $k$, Fehlerkonstanten, u.s.w., erheblich verschlechtert würden. Beispielsweise wäre ein $A_\infty^0$-stabiles, lineares 4-Schrittverfahren, bei Fehlerkonstanten im Bereich von ca. $1\over10$, ein ideales Verfahren zur Lösung steifer Differentialgleichungen. Die zweite Dahlquist-Barriere besagt, daß es solch ein Verfahren prinzipiell nicht geben kann.

Alternative Beweisgänge werden in Wanner (1987), s.u., angedeutet, allerdings nicht streng bewiesen. Die Trapezregel ist nicht $A_\infty^0$-stabil. Die Dissertation Tischer (1983) und Tischer/Sacks-Davis (1983) geben $A_\infty^0$- und $S_\infty^0$-stabile zyklische zweistufige Verfahren an, mit den Konvergenzordnungen $p=2,3,4$ und benötigten Startwerten von $k=2,3,4$. Allerdings sind bei allen Stufen der zyklischen Formeln von Tischer (1983) (Dissertation) und Tischer/Sacks-Davis (1983), die Äquilibrierungsmaße vergleichsweise hoch.

Bibliographisch: Peter E. Tischer, Ron Sacks-Davis.

Peter E. Tischer: "The Cyclic Use of Linear Multistep Formulas for the Solution of Stiff Differential Equations", Ph.D. Thesis, Department of Computer Science, Monash University, Clayton, Victoria, Australia, August 1983, x+180 pages.

Wanner, Gerhard (*1942): "Order Stars and Stability", in ``The State of the Art in Numerical Analysis", Proceedings of the joint IMA/SIAM conference held at the University of Birmingham, 14--18 April 1986, Edited by A. Iserles and M.J.R. Powell, Clarendon Press, Oxford, 1987, 451--471

Die Tatsache, daß ein explizites lineares Mehrschrittverfahren nicht $A$-stabil sein kann, ist sehr leicht einzusehen. Auch eine beliebige Kombination expliziter linearer Mehrschrittverfahren, mit äquidistanter Gitterweite $h$, kann nicht $A$-stabil sein. Sehr wohl kann jedoch eine Kombination von impliziten und expliziten Verfahren sogar $A_\infty^0$-stabil sein.

Zusätzliches Licht auf die Beziehung zwischen Konsistenzordnung und Stabilitätseigenschaften wirft der Satz von Jeltsch/Nevanlinna (1982).

Bibliographisch: Rolf Jeltsch, Olavi Nevanlinna.

5. Satz: siehe Jeltsch/Nevanlinna (1982). Es gilt

$$ \forall k\in\mathbb{N}:\forall\alpha\in[0^\circ,90^\circ):\quad \hbox{existiert mindestens ein $A[\alpha]$-stabiles Verfahren} $$

und

$$ \forall k\in\mathbb{N}:\forall\delta\gt 0:\quad\hbox{existiert mindestens ein $S[\delta]$-stabiles Verfahren}. $$

Weiter gibt es funktionale Zusammenhänge zwischen Fehlerkonstanten und erreichbarer Höchstordnung.

6. Satz: siehe Jeltsch/Nevanlinna (1986).

Voraussetzung: Es sei $c_{p+1}$ der Fehlerfaktor von

$$ (L_hy)(t) := \sum_{i=0}^k \sum_{j=0}^m \alpha_{ij} h^j y^{(j)}(t+ih), \qquad c_{p+1}={1\over y^{(p+1)}(0)} (L_1y)(0). $$

Die Wurzeln von $\rho_0(\zeta)$ seien $\zeta_1=1,\zeta_2,\ldots,\zeta_k$, welche in der Einheitskreisscheibe liegen:

$$ \left|\zeta_i\right|\le R,\qquad 0\le R\le1. $$

Behauptung: (1) Ist die Formel explizit und von der Ordnung $p=mk$, so gilt

$$ c_{p+1} \ge {1\over(m+k)!} 2^{1-k} \sum_{j=0}^{k-1} {k-1\choose j} \left(1-R\over1+R\right)^j f_{jkm}. $$

(2) Im impliziten Falle und der Ordnung $p=(k+1)m$ gilt

$$ (-1)^m c_{p+1} \ge {(-1)^m\over\left[(k+1)m\right]!} 2^{1-k} \sum_{j=0}^{k-1} {k-1\choose j} \left(1-R\over1+R\right)^j e_{jkm}. $$

(3) In beiden Fällen, also in Fall (1) und (2) ist Gleichheit möglich bei der Formel maximaler Ordnung mit dem charakteristischen Polynom

$$ \rho_0(\zeta) = (\zeta-1)(\zeta+R)^{k-1}, $$

welches stabil ist, für $R<1$ oder $R=1$ und $k\le2$.

Die Größen $f_{jkm}$ und $e_{jkm}$ sind gegeben durch

$$ f_{jkm} := \sum_{i=0}^{k-1} T_{ij} W_{i,k-1,m}, \qquad W_{jkm} := \int_j^{j+1} \prod_{\nu=0}^k (t-\nu) dt, $$

und

$$ e_{jkm} := \sum_{i=0}^{k-1} T_{ij} W_{i,k+1,m}, \qquad T_{ij} := \sum_{\mu=\max(0,i-j)}^{\min(n-j,i)} {j\choose i-\mu} {n-j\choose\mu} (-1)^{j-i+\mu}, $$

für $j=0,\ldots,k-1$.

Beweis: Jeltsch/Nevanlinna (1986).     ☐

Nach dem Satz von Jeltsch/Nevanlinna (1982) gibt es also “fast $A$-stabile” lineare Mehrschrittverfahren mit beliebiger Anzahl von Startwerten. Dennoch sind diese Verfahren nicht $A$-stabil und schon gar nicht $A_\infty^0$- oder $S_\infty^0$-stabil, nach der zweiten Dahlquist-Barriere. Vielmehr rücken die Wurzeln betragsmässig immer mehr der eins näher, für $H\to\infty$, und die Fehlerkonstanten werden mit größer werdenden $k$ immer größer. Beispielsweise verwendet Gupta (1985) $k$-Schrittverfahren mit den in der Tabelle angegebenen Eigenschaften.

Bibliographisch: Gopal K. Gupta.

p $\alpha$ $\delta$ $c_{p+1}$ $\mu_\infty$ k
1 90.00 0.0 0.5 0.0 1
2 90.00 0.0 0.083 1.0 1
3 86.46 0.075 0.242 0.32 3
4 80.13 0.282 0.374 0.43 5
5 73.58 0.606 0.529 0.567 7
6 67.77 1.218 0.724 0.878 9
7 65.53 1.376 1.886 0.898 12
8 64.96 1.149 7.686 0.790 16
9 62.78 2.086 16.737 0.989 19
10 63.74 1.223 133.955 0.878 26

Hierbei bedeutet $p$ die Konvergenzordnung, $k$ die Anzahl der Startwerte, $c_{p+1}$ der Fehlerfaktor, $\alpha$ der Widlund-Winkel, und $\delta$ ist der entsprechende Wert bei der $S[\delta]$-Stabilität. $\mu_\infty$ ist der Betrag der betragsmässig größten Wurzel bei $\infty$.

Es zeigt sich, daß das Programm DSTIFF, welches diese Formeln benutzt, häufig doppelt so viele Funktionsauswertungen und doppelt so lange Rechenzeiten beansprucht, wie das Programm LSODE, welches auf den BDF$i$, mit $i\in\{1,\ldots,5\}$ basiert. Man beachte allerdings, daß hier Implementierungen von Formeln und Heuristik miteinander verglichen wurden. Dieser Vergleich, den Gupta (1985) anstellte, ist also keine endgültige Wertung von Formeln, sondern eine Wertung von Programmen. Eine geschickte Programmierung und eine durchdachte Heuristik sind von nicht zu unterschätzender Wichtigkeit.

Aufgrund einer modifizierten Strategie für die Korrektoriteration in dem Programm DSTIFF, sind die Anzahlen für $LU$-Zerlegungen (= Jacobimatrixauswertungen) leicht geringer, als für das Programm LSODE. Lediglich erwartungsgemäss für das Problem B5 benötigt das Programm DSTIFF bedeutend weniger Schritte, als das Programm LSODE. Bei diesem Problem werden die Stabilitätseigenschaften besonders gefordert. Wüßte man im voraus um die Lage der Eigenwerte der konstanten Jacobimatrix, so könnte man bei dem Programm LSODE gleich von vorne herein eine passende Höchstordnung wählen und damit würden sich beide Programme wieder angleichen. Das Verhältnis der Schritte der beiden Programme zueinander beträgt bei B5 grob $1:4$ (DSTIFF:LSODE). Dieses Verhältnis übersetzt sich allerdings nicht in gleichem Maßstab auf die Rechenzeit. Die Rechenzeit ist lediglich um einen wesentlich geringeren Betrag angestiegen.

Eine gewisse Sonderstellung nehmen die BDF$i$ ein, wegen $\sigma(\mu)=\mu^i$.

7. Bemerkung: Die BDF$i$ sind die einzigen linearen $i$-Schrittverfahren der Konsistenzordnung $i$, für $i=1,\ldots,6$, die $A_\infty^0[\alpha]$- bzw. $S_\infty^0[\delta]$-stabil sind.

Es gibt weitere lineare Mehrschrittverfahren ($\ne$BDF$i$), die $A_\infty^0[\alpha]$- bzw. $S_\infty^0[\delta]$-stabil sind, jedoch ist dann die Konsistenzordnung $i$ nicht mehr mit $i$ Startwerten zu erreichen. Für mehrstufige Verfahren gilt die Bemerkung nicht mehr, wie z.B. die zyklischen Formeln von Tendler (1973) deutlich machen, siehe auch Tendler/Bickart/Picel (1978).

6. Annullierte Dominanz und Totalannullation

Wie üblich bedeute “$\cong$” Gleichheit bis auf ${\cal O}(h^{p+2})$. Jede Stufe eines zusammengesetzten Verfahrens wird nun in einen Taylorabschnitt (Taylor, Brook (1685--1731)) zerlegt und man erhält hierfür

$$ \gamma \cong \pmatrix{ c_{11}\dot yh+\cdots+c_{1,p+1}y^{(p+1)}h^{p+1}\cr c_{21}\dot yh+\cdots+c_{2,p+1}y^{(p+1)}h^{p+1}\cr \vdots \qquad \cdots \qquad \vdots\cr c_{s1}\dot yh+\cdots+c_{s,p+1}y^{(p+1)}h^{p+1}\cr } = \underbrace{\pmatrix{ c_{11} & \ldots & c_{1p}\cr c_{21} & \ldots & c_{2p}\cr \vdots & \ddots & \vdots\cr c_{s1} & \ldots & c_{sp}\cr}}_{\in\mathbb{C}^{s\times (p+1)}}{\mskip 3mu} \underbrace{\pmatrix{\dot yh\cr \ddot yh^2\cr \vdots\cr y^{(p)}h^p\cr}} _{\in\mathbb{C}^p} $$

$\gamma$ kann aufgespalten werden in eine Summe

$$ \gamma \cong \pmatrix{c_{11}\cr\vdots\cr c_{s1}\cr}\dot yh +\pmatrix{c_{21}\cr\vdots\cr c_{s,2}\cr}\ddot yh +\cdots+\pmatrix{c_{1,p+1}\cr\vdots\cr c_{s,p+1}\cr}y^{(p+1)}h^{p+1} $$

und jeder Summand wird einzeln auf annullierte Dominanz, oder Totalannullation geprüft. Bei Verfahren, bei denen alle Stufen die gleiche Konsistenzordnung haben, sind $c_{ij}=0$, für $j\le p$, und $i=1,\ldots,s$, und die obige Gleichung reduziert sich dann auf

$$ \gamma \cong \pmatrix{ c_{1,p+1}y^{(p+1)}h^{p+1}\cr \vdots\cr c_{s,p+1}y^{(p+1)}h^{p+1}\cr} = \pmatrix{c_{1,p+1}\cr \vdots\cr c_{s,p+1}\cr} y^{(p+1)}h^{p+1}. $$

1. Beispiel: Es werde das explizite Euler-Verfahren als Prädiktor verwendet und mit der BDF2 werde zweimal anschliessend iteriert.

Schritt Formel Verfahren
Prädiktor $y^0_{n+1}=y_n+z_n$ explizites Euler-Verfahren
Korrektor $y_{n-1}-4y_n+3y^1_{n+1}=2z^0_{n+1}$ BDF2
Korrektor $y_{n-1}-4y_n+3y_{n+1}=2z^0_{n+1}$ BDF2

Mit $u_n=(y^0_n,{\mskip 3mu}y^1_n,{\mskip 3mu}y_n)$ ergibt sich für die sechs Matrizen

$$ A_0=\pmatrix{0&0&0\cr 0&0&1\cr 0&0&1\cr},\quad A_1=\pmatrix{0&0&-1\cr 0&0&-4\cr 0&0&-4\cr},\quad A_2=\pmatrix{1&0&0\cr 0&3&0\cr 0&0&3\cr},\quad B_0=B_1={\bf 0},\quad B_2=\pmatrix{0&0&0\cr 2&0&0\cr 0&2&0\cr}. $$

Für das Matrixpolynom $\rho(\mu)=A_0+A_1\mu+A_2\mu^2$ ergibt sich

$$ \rho(\mu) = \pmatrix{ \mu^2 & 0 & 0\cr 0 & 3\mu^2 & 1-4\mu\cr 0 & 0 & 3\mu^2-4\mu+1\cr} $$

Auffällig ist die obere Dreiecksgestalt und die Verteilung der $\mu^\kappa$-Terme auf der Diagonalen. Weiterhin erscheint in der rechten unteren Ecke das charakteristische Polynom des Korrektors. Dieser Sachverhalt gilt ganz allgemein. Der Nullensatz für Prädiktor-Korrektor-Verfahren lautet:

2. Satz: Sei $\rho_c$ das charakteristische Polynom des Korrektors. Für das Matrixpolynom eines $P(EC)^i\{E\}$-Verfahrens, mit

$$ \rho(\mu) = \sum_{\nu=0}^\kappa A_\nu\mu^\nu,\qquad \deg\rho = \kappa,\qquad A_\nu\in\mathbb{C}^{(i+1)\times(i+1)},\quad\nu=0,\ldots,\kappa, $$

hat man für alle $\mu\in\mathbb{C}$ die Darstellung

$$ \rho(\mu) = \pmatrix{ \alpha_{11}\mu^\kappa & 0 & \ldots & *\cr 0 & \alpha_{22}\mu^\kappa & \ldots & *\cr \vdots & \vdots & \ddots & \vdots\cr 0 & 0 & \ldots & \rho_c(\mu)\cr} \in \mathbb{C}^{(i+1)\times(i+1)}. $$

Beweis: (Nullensatz) Es ist

$$ u_{n+\nu} = \pmatrix{y^0_{n+\nu}\cr \vdots\cr y^{i-1}_{n+\nu}\cr y_{n+\nu}\cr}, \qquad \nu=-\kappa+1,\ldots,1. $$

Die ersten $i$ Komponenten der Vektoren $u_{n+\nu}$ kommen in der letzten Korrektorstufe nicht vor. Die Matrizen $A_0,\ldots,A_{\kappa-1}$ tragen Elemente lediglich auf der letzten Spalte, sind also alle von der Form

$$ A_0\sim A_1\sim\ldots A_{\kappa-1}\sim\pmatrix{ 0 & \ldots & 0 & *\cr \vdots & \ddots & \vdots & \vdots\cr 0 & \ldots & 0 & *\cr} $$

während hingegen $A_\kappa$ Diagonalgestalt hat, also

$$ A_\kappa = \pmatrix{ \alpha_{11} & & \llap{0}\cr & \ddots & \cr \rlap{0} & \ldots & \alpha_\kappa\cr}, $$

wobei $\alpha_{11}=\alpha_\kappa^P$, der führende Koeffizient des Prädiktor-Verfahrens ist und $\alpha_{22}=\cdots=\alpha_{i+1,i+1}=\alpha_\kappa$ gleich dem führenden Koeffizient des Korrektors ist. Summation der $A_\nu\mu^\nu$ ergibt die Behauptung.     ☐

Nur $y^*_{n+1}$, also die Lösung der impliziten Korrektorgleichung für die Zeit $t_{n+1}$, ist gesucht. Die anderen vergangenen Werte sind schon gefunden. Die Zwischenwerte der vergangenen Iterationen werden nicht mehr verwendet, anders als bei semiiterativen Verfahren. Würde man diese dennoch verwenden, so könnte sich natürlich auch das Spektrum ändern, weil sich sich dann auch die obere Dreiecksgestalt ändert. Vorausgesetzt wird hier ebenfalls, daß immer mit dem gleichen Korrektor iteriert wird, und daß nur ein einziger Prädiktor genommen wird. Z.B. verwenden off-step-point Verfahren u.U. mehrere Prädiktoren, so auch Filippi/Kraska (1973).

Bibliographisch: Siegfried Filippi (1929--2022), Ernst Kraska (1932--2021), Todesanzeige.

3. Folgerungen: (1) Das charakteristische Polynom $Q(\mu,0)=\det\rho(\mu)$ sowohl des $P(EC)^iE$-, als auch des $P(EC)^i$-Verfahrens hat mindestens $\kappa i$ Nullen im Spektrum und die restlichen $\kappa$ Eigenwerte stimmen mit den Wurzeln des charakteristischen Polynomes $\rho_c$ des Korrektors überein.

(2) Insbesondere ist ein $P(EC)^iE$- bzw. $P(EC)^i$-Verfahren genau dann stabil, wenn der Korrektor stabil ist. Die Stabilität der Prädiktorformel ist völlig unerheblich für die $D$-Stabilität des $P(EC)^iE$- bzw. $P(EC)^i$-Verfahrens. Sehr wohl hat natürlich der Prädiktor Einfluß auf das Aussehen des Stabilitätsgebietes.

(3) Die Linkseigenvektoren $v_m$ mit $m=1,\ldots$, zu den nicht zu Null gehörenden Eigenwerten, also somit $m\le\kappa$, sind alle von der Form

$$ v_m = (0,{\mskip 3mu}\ldots,{\mskip 3mu}0,{\mskip 3mu}1)\in\mathbb{C}^{1\times(i+1)},\qquad m\le\kappa. $$

Folgerung (2) kann man auch auf andere Art und Weise einsehen. Seien die beiden Funktionen $\varphi$ und $\psi$ Lipschitz-stetig, dann ist auch die verkettete Funktion $\varphi\circ\psi$ Lipschitz-stetig. Ist also

$$ \left|\varphi(x)-\varphi(y)\right|\le K\left|x-y\right|\qquad\hbox{und}\qquad \left|\psi(u)-\psi(v)\right|\le L\left|u-v\right|, $$

so gilt

$$ \left|\varphi(\psi(u))-\varphi(\psi(v))\right| \le K\left|\psi(u)-\psi(v)\right|\le KL\left|u-v\right|. $$

Denkt man sich nun ein Picard-Prädiktor-Korrektor-Verfahren nicht als einziges mehrstufiges, i.d.R. lineares Verfahren, sondern denkt man es sich als ein ineinander verschachteltes, i.d.R. nicht-lineares Verfahren, so sieht man ebenfalls sofort, daß für die Stabilitätseigenschaften bzgl. $h\to0$, nur der Korrektor maßgeblich ist und der Prädiktor unmaßgeblich ist.

4. Sei der Prädiktor $\hat z_{n+\kappa}$ gegeben durch die Matrix-Differenzengleichung

$$ \hat A_0z_n+\cdots+\hat A_{\kappa-1}z_{n+\kappa-1}+\hat A_\kappa\hat z_{n+\kappa} = h\varphi(z_n,\ldots,z_{n+\kappa-1}),\qquad \hat A_\kappa=I $$

und der Korrektorwert ergebe sich als Lösung der Matrix-Differenzengleichung

$$ A_0z_n + \cdots + A_\kappa z_{n+\kappa-1} + A_\kappa = h \psi(z_n,\ldots,z_{n+\kappa-1},z_{n+\kappa}), \qquad A_\kappa = I. $$

Picard-Iteration besteht nun darin, daß man den Wert $z_{n+\kappa}$ in der Funktion $\psi$ ersetzt durch die Iterierte der vorherigen Iteration, also wird $z_{n+\kappa}$ in $\psi$ ersetzt durch $\hat z_{n+\kappa}$. Direkte Substitution der Bestimmungsgleichung für den Prädiktorwert $\hat z_{n+\kappa}$ in die Funktion $\psi$, ergibt dann

$$ \sum_{i=0}^\kappa A_i z_{n+i} = h\vartheta(z_n,\ldots,z_{n+\kappa-1}), $$

wobei die Funktion $\vartheta$ Lipschitz-stetig ist und damit erhält man mit den üblichen Sätzen bei gegebener Konsistenz und Stabilität dann die Konvergenz. Über die Entstehungsgeschichte von $\vartheta$ braucht man nichts zu wissen, außer eben, daß $\vartheta$ Lipschitz-stetig bzgl. seiner Argumente ist. Insbesondere die Nullstabilität hängt jetzt offensichtlich nur noch von den Matrizen $A_i$ ab. Iteriert man häufiger als einmal, also Picard-$P(EC)^i\{E\}$ (i>1), so wird die Verschachtelungtiefe nur höher, am Prinzip ändert sich nichts.

Aus der Folgerung (3) ergibt sich jetzt sofort, daß die Fehlervektoren des Prädiktors und die Fehlervektoren der Zwischeniterierten von den Linkseigenvektoren $v_m$ vollständig weggefiltert werden. Dies heißt, die Eigenwerte nicht gleich Null bekommen nur den Fehlerfaktor des Korrektors zu Gesicht und die restlichen Fehlerfaktoren werden von den Nullen im Spektrum total annulliert. Überhaupt gibt es Parallelen zwischen Runge-Kutta-Verfahren und Prädiktor-Korrekor-Verfahren. Prädiktor-Korrektor-Verfahren steht hier sowohl für $P(EC)^iE$- also auch für $P(EC)^i$-Verfahren, kurz $P(EC)^i\{E\}$-Verfahren. Diesen Effekt der Totalannullation kann man anhand eines Beispiels besonders deutlich nachvollziehen.

5. Beispiel: Als Prädiktor-Korrektor-Verfahren werde verwendet

Schritt Formel Verfahren
Prädiktor $y^0_{n+1}=y_n+z_n$ explizites Euler-Verfahren
Korrektor $y_{n-1}-4y_n+3y_{n+1}=2z^0_{n+1}$ BDF2

Mit $u_n=(y^0_n,{\mskip 3mu}y_n)$ erhält man die sechs Matrizen

$$ A_0=\pmatrix{0&0\cr 0&1\cr},\quad A_1=\pmatrix{0&-1\cr 0&-4\cr},\quad A_2=\pmatrix{1&0\cr 0&3\cr},\qquad B_0={\bf 0},\quad B_1=\pmatrix{0&1\cr 0&0\cr},\quad B_2=\pmatrix{0&0\cr 2&0\cr} . $$

Als Lösung für die Differenzengleichung

$$ A_2u_{n+2}+A_1u_{n+1}+A_0u_n=\Gamma\hat Z,\qquad n=0,1,\ldots $$

erhält man nach Durchmultiplikation mit $A_2^{-1}$ für den dominanten Term

$$ \pmatrix{1&*&*&1\cr 0&*&*&1\cr} \pmatrix{0&&&\cr &0&&\cr &&1/3&\cr &&&1\cr} \pmatrix{*&*\cr *&*\cr 0&1\cr 0&1\cr} \pmatrix{1/2&*\cr 0&-2/3\cr} \hat Z. $$

Die ersten Fehlerterme des Prädiktors sind $h^2/2\ddot y+{\cal O}(h^3)$ und für den Korrektor lauten sie $-2h^3/3y^{III}+{\cal O}(h^4)$.

Die Fehlervektoren des Prädiktors liegen nun gerade so, daß sie genau auf diese Nullen heraufpassen, das heißt die Fehlervektoren stehen senkrecht auf den nicht zu Null gehörenden Jordanvektoren. Die vom Prädiktor gelieferten niedrigen Konsistenzordnungen, werden deswegen total annulliert (Totalannullation). Dieses Verhalten ist völlig analog dem Verhalten bei Runge-Kutta-Verfahren, wo die Stufen mit niedrigen Konsistenzordnungen von den Nullen in der Jordanmatrix

$$ J = \pmatrix{ 0 & & & \llap0\cr & \ddots & & \cr & & 0 & \cr \rlap0 & & & 1\cr} $$

vollständig bedämpft werden und somit keinerlei Wirkung zeigen. Dies gilt zumindestens im asymptotischen Falle $h\to 0$, wo allein einzig die $A_\nu$ entscheidend wirken und die Matrizen $B_\nu$ keine Rolle spielen.

6. Beispiel: Runge-Kutta-Verfahren mit insgesamt 4 Stufen.

$$ A = \pmatrix{ 0 & 0 & 0 & 1\cr 0 & 0 & 0 & 1\cr 0 & 0 & 0 & 1\cr 0 & 0 & 0 & 1\cr}, \quad X = \pmatrix{ 1 & 0 & 0 & 1\cr 0 & 1 & 0 & 1\cr 0 & 0 & 1 & 1\cr 0 & 0 & 0 & 1\cr}, \quad J = \pmatrix{ 0 & & & \llap0\cr & 0 & & \cr & & 0 & \cr \rlap0 & & & 1\cr}, \quad Y = \pmatrix{ 1 & -1 & 0 & 0\cr 0 & 1 & -1 & 0\cr 0 & 0 & 1 & -1\cr 0 & 0 & 0 & 1\cr} $$

Die Konsistenzordnung kann pro Stufe um eine Einheit steigen. Die Matrix $C$ hat somit die Form

$$ C = \pmatrix{ c_1 & * & * & *\cr 0 & c_2 & * & *\cr 0 & 0 & c_3 & *\cr 0 & 0 & 0 & c_4\cr} $$

7. Beispiel: Die hier zutage tretende Ähnlichkeit zwischen Runge-Kutta-Verfahren und $P(EC)^i\{E\}$-Verfahren gilt sogar soweit, daß manche $P(EC)^i\{E\}$-Verfahren mit bestimmten expliziten Runge-Kutta-Ver-fahren völlig gleichwertig sind. Zum Beispiel gilt für das verbesserte Euler-Verfahren mit dem Parametertableau

$$ \begin{array}{c|cc} 1 && 1 & \cr \hline && {1\over2} & {1\over2}\cr \end{array} \qquad\qquad \eqalign{k_0 &= f(t_0,y_0)\cr k_1 &= f(t_0+h,y_0+hk_0)\cr \hline y_1 &= y_0+{h\over2}\left(k_0+k_1\right)\cr } $$

daß es völlig identisch ist mit dem impliziten Trapezverfahren, wobei das explizite Euler-Verfahren als Prädiktor verwendet wird:

$$ % sich ausdehnender Rechtspfeil, s.Knuth, S.325 % ev. noch ein \smash einfügen, sodaß \matrix "nichts" merkt %\def\mapright#1{ % \setbox0=\hbox{$#1$} % temporär ablegen und dann ausmessen % \dimen0=\wd0 % weil man \wd0 nicht "advancen" kann % \advance\dimen0 by 0.7cm % bißchen mehr Platz % \mathop{ % \limits will dies halt so % \hbox to \dimen0{\rightarrowfill} % } \limits^{#1} % darüber die Information %} y_{n+1} = y_n + {h\over2}\left\{f_n+f_{n+1}\right\} %\mapright{PECE} \overset{PECE}\longrightarrow y_n + {h\over2}\left\{f_n+f(t_{n+1},y_n+hf_n)\right\}. $$

Die Analyse beider Verfahren geschieht häufig völlig getrennt. Die Konsistenzordnung 2 des verbesserten Euler-Verfahrens weist man häufig durch Taylorentwicklung direkt nach. Beim Prädiktor-Korrektor Verfahren wendet man die Konsistenzsätze an, weist Stabilität des Korrektors nach und zeigt schließlich mit Hilfe des Satzes von Liniger (1971), daß die Konsistenzordnung des Korrektors erhalten bleibt, wenn man ausreichend lange iteriert.

8. Wichtig für das Konvergenzverhalten ist die spektrale Struktur der drei Matrizen $X$, $J$ und $Y$:

$$ \eqalign{ X: & \qquad XJ^\nu Y\gamma=\ldots,\cr J: & \qquad J^\nu Y\gamma=\ldots,\cr Y: & \qquad Y\gamma=(\ldots 0\ldots)^\top.\cr } $$

Seien $v_1,\ldots,v_r$ die Linkseigenvektoren zum Matrixpolynom $\rho$ zum Eigenwert $\lambda=1$, so gilt als Bedingung der annullierten Dominanz

$$ v_i{\mskip 3mu}\rho(1)=0,\quad v_i{\mskip 3mu}\gamma=0,\quad v_i\ne{\bf0}^\top,\qquad i=1,\ldots,r. $$

Für den Fall $r=1$ erhält man also die geometrische Bedingung, daß die Spalten der Matrix $(\rho(1),{\mskip 3mu}\gamma)$ aus dem orthogonalen Komplement von $v$ sein müssen, somit

$$ (\rho(1),{\mskip 3mu}\gamma)\in v^\bot. $$

Ist der Linkskern von $\rho(1)$ nicht mehr 1-dimensional, sondern $r$-dimensional, so hat die Bedingung zu gelten

$$ \def\spn{\mathop{\rm span}} \left(\rho(1),\gamma\right)\in\left(\spn(v_1,\ldots,v_r)\right)^\bot $$

Für Eigenwerte $\left|\mu\right|=1$ gilt ganz entsprechend

$$ \left(\rho(\mu),\gamma\right)\in\left(\spn(v_1,\ldots,v_r)\right)^\bot, $$

wobei $r$ jetzt die Vielfachheit des Eigenwertes zu $\left|\mu\right|=1$ ist und entsprechend $v_1,\ldots,v_r$ die Linkseigenvektoren zu diesem Eigenwert sind. Algebraische Vielfachheit (=Multiplizität der Nullstelle des charakteristischen Polynoms) und geometrische Vielfachheit (=Dimension des invarianten Unterraumes) müssen bei dominanten Eigenwerten natürlich gleich sein. Anhand der oben schon angegebenen Darstellung für die Lösung der Matrixdifferenzengleichung, und zwar in der Form

$$ u_{m+1} = XT^{m+1}c + X\sum_{i=0}^m T^{m-i}Yf_i, \qquad m=0,1,\ldots, $$

ist das Auftauchen der Linkseigenvektoren sofort offenkundig. Die Bedingung der annullierten Dominanz ist eine stetige Invariante, da bei einfachen Eigenwerten die Eigenvektoren stetig von kleinen Änderungen abhängen. Bei mehrfachen Eigenwerten muß dies nicht unbedingt gelten.

7. Das $n$-dimensionale äußere Produkt für $n-1$ Vektoren

Man vgl. auch On Differential Forms.

1. Da die Determinante in jeder Spalte linear ist, stellt

$$ x\mapsto\det\left(a_1,\ldots,a_{n-1},x\right) $$

für fest gegebene $a_1,\ldots,a_{n-1}\in\mathbb{R}^n$, eine Linearform des $\mathbb{R}^n$ dar. Nach dem Darstellungssatz von Riesz:

Sei $H$ ein Hilbertraum, und sei $f:H\to\mathbb{C}$ ein stetiges lineares Funktional, dann

$$ \dot\exists b\in H:\enspace\forall x\in H:\quad f(x)=\langle b,x\rangle$$

und weiter ist $|b|=|f|$.

Daher gibt es genau einen Vektor $b\in\mathbb{R}^n$, sodaß die Linearform als Skalarprodukt geschrieben werden kann:

$$ \dot\exists b:\enspace\forall x:\quad \det\left(a_1,\ldots,a_{n-1},x\right)=\langle b,x\rangle \qquad(a_i\hbox{ fest}). \tag{*} $$

Bibliographisch: Riesz, Friedrich (1880--1956).

2. Diesen, implizit durch das Skalarprodukt, eindeutig bestimmten Vektor $b$ nennt man das Vektorprodukt (auch Kreuzprodukt oder äußeres Produkt) und schreibt hierfür

$$ b =: a_1\wedge\cdots\wedge a_{n-1} =: \bigwedge_{i=1}^{n-1} a_i, $$

oder auch

$$ a_1\times\cdots\times a_{n-1}=\mathop{\times}_{i=1}^{n-1} a_i. $$

Es gilt also

$$ \det\left(a_1,\ldots,a_{n-1},x\right) = \left\langle\bigwedge_{i=1}^{n-1}a_i,x\right\rangle.\tag{**} $$

3. Hieraus liest man ab

$$ \eqalign{ \bigwedge_{i=1}^{n-1}a_i=0\quad &\Longleftrightarrow\quad a_i\hbox{ linear abhängig},\cr a_1\wedge\cdots\wedge a_i\wedge\cdots\wedge a_k\wedge\cdots\wedge a_{n-1}&= -\left(a_1\wedge\cdots\wedge a_k\wedge\cdots\wedge a_i\wedge\cdots\wedge a_{n-1}\right),\cr a_1\wedge\cdots\wedge a_i+\hat a_i\wedge\cdots\wedge a_{n-1} &= \left(a_1\wedge\cdots\wedge a_i\wedge\cdots\wedge a_{n-1}\right) + \left(a_1\wedge\cdots\wedge \hat a_i\wedge\cdots\wedge a_{n-1}\right),\cr a_1\wedge\cdots\wedge\lambda a_i\wedge\cdots\wedge a_{n-1} &= \lambda\left(a_1\wedge\cdots\wedge a_i\wedge\cdots\wedge a_{n-1}\right),\cr \left\langle\bigwedge_{i=1}^{n-1}a_i,a_k\right\rangle &= 0, \quad k=1,\ldots,n-1.\cr } $$

Die letzte Gleichung sagt, daß das äußere Produkt senkrecht auf jedem “Einzelfaktor” steht. Weiter kann man jetzt noch die Jacobische und die Grassmannsche Identität leicht nachrechnen. Die obigen Gleichungen gelten auch für $n=2$, wobei dann $\bigwedge_{i=1}^1a_i=a_1$ ist.

Bibliographisch: Grassmann, Hermann (1809--1877), Jacobi, Carl Gustav (1804--1851).

Das oben eingeführte äußere Produkt ist ein spezielles äußeres Produkt. Es gibt weitere äußere Produkte. Bei diesen ist der Bildbereich nicht mehr notwendig gleich $\mathbb{C}^n$, sondern $\mathbb{C}[{n\choose m}]$, bei einem $m$-fachen Produkt. Für $m=n-1$ ergibt sich natürlich genau das oben angegebene Produkt, bis auf Proportionalität.

4. Die Komponenten des Vektors $b$ bei $(*)$, ergeben sich durch sukzessives Einsetzen der $n$ Einheitsvektoren $e_i$ zu

$$ b_i = \left\langle b,e_i\right\rangle = \det\left(a_1,\ldots,a_{n-1},e_i\right), \qquad i=1,\ldots,n. $$

Für den Betrag des äußeren Produktes gilt

$$ \left|\bigwedge_{i=1}^{n-1}a_i\right|^2 = \left|\matrix{ \langle a_1,a_1\rangle & \ldots & \langle a_1,a_{n-1}\rangle\cr \vdots & \ddots & \vdots\cr \langle a_1,a_{n-1}\rangle & \ldots & \langle a_{n-1},a_{n-1}\rangle\cr }\right|, $$

weges des Satzes über die Gramsche Determinante (Gram, Jorgen Pedersen (1850--1916))

$$ \det(a_1,\ldots,a_n){\mskip 3mu}\det(b_1,\ldots,b_n) = \left|\matrix{ \langle a_1,b_1\rangle & \ldots & \langle a_1,b_n\rangle\cr \vdots & \ddots & \vdots\cr \langle a_n,b_1\rangle & \ldots & \langle a_n,b_n\rangle\cr }\right|. $$

5. Die Definition des Vektorproduktes kann auch in der folgenden Form geschehen:

$$ x\mapsto\det(a_1,\ldots,a_{i-1},x,a_{i+1},\ldots,a_n) $$

ist eine Linearform für feste $a_i$, u.s.w. Die $a_1,\ldots,a_{i-1},x,,a_{i+1},\ldots,a_n$ bilden ein Rechtssystem.

6. Beispiel: $n=3$. Gesucht sind die Komponenten des Vektorproduktes $a\times b$, mit $a=(\alpha_1,{\mskip 3mu}\alpha_2,{\mskip 3mu}\alpha_3)$ und $b=(\beta_1,{\mskip 3mu}\beta_2,{\mskip 3mu}\beta_3)$. Zu berechnen sind drei Determinanten,

$$ a\times b = \pmatrix{ \left|\matrix{\alpha_1&\beta_1&1\cr \alpha_2&\beta_2&0\cr \alpha_3&\beta_3&0\cr}\right|\\[0.5em] \left|\matrix{\alpha_1&\beta_1&0\cr \alpha_2&\beta_2&1\cr \alpha_3&\beta_3&0\cr}\right|\\[0.5em] \left|\matrix{\alpha_1&\beta_1&0\cr \alpha_2&\beta_2&0\cr \alpha_3&\beta_3&1\cr}\right|\cr} = \pmatrix{ \alpha_2\beta_3-\alpha_3\beta_2\cr \alpha_3\beta_1-\alpha_1\beta_3\cr \alpha_1\beta_2-\alpha_2\beta_1\cr}. $$

7. Beispiel: $n=2$. Gesucht sind die Komponenten des Vektors $a^\bot$, welcher senkrecht steht auf $a$, mit $a={\alpha_1\choose\alpha_2}$. Das äußere Produkt liefert gerade solch einen Vektor. In diesem Fall hat das Produkt nur einen Faktor. Zu berechnen sind hier $n=2$ Determinanten und zwar

$$ a^\bot = \pmatrix{ \left|\matrix{\alpha_1 & 1\cr \alpha_2 & 0\cr}\right| \\[0.5em] \left|\matrix{\alpha_1 & 0\cr \alpha_2 & 1\cr}\right| \cr} = \pmatrix{-\alpha_2\cr \alpha_1\cr}. $$

Eine Einführung in das Vektorprodukt findet man beispielsweise in den Büchern von Walter (1986) oder Koecher (1985). Besonders hervorzuheben ist hierbei die ausführliche Darstellung von Gröbner (1966).

Bibliographisch: Rolf Walter (1937--2022), Max Koecher (1924--1990), Wolfgang Gröbner (1899--1980).

8. Äußeres Produkt und Fehlerkonstanten

Für die Fehlerkonstante von Henrici

$$ C := {v{\mskip 3mu}\gamma\over v{\mskip 3mu}\rho'(1){\mskip 3mu}w},\qquad\left\{\eqalign{ v{\mskip 3mu}\rho(1)&=0, \quad v\ne0,\cr \rho(1){\mskip 3mu}w&=0, \quad w\ne0,\cr }\right. $$

erhält man nun das folgende Resultat. Da das äußere Produkt $\bigwedge_{i=1}^{n-1}a_i$ senkrecht steht auf $a_i$, für $i=1,\ldots,n-1$, ist dieses Produkt also Linkseigenvektor von $\rho(1)$, wenn man die Spalten der Matrix $\rho(1)$ mit $a_i$ bezeichnet und einen Spaltenvektor, sagen wir $a_n$, herausstreicht. Wenn man einmal von Umnumerierungen absieht, so hat man damit alle Fälle abgedeckt. Die Restmatrix sei

$$ \widehat{\rho(1)}\in\mathbb{R}^{n\times(n-1)}. $$

Ist $\bigwedge_{i=1}^{n-1}a_i=0$, so hat $\rho(1)$ einen mehrfachen Eigenwert $\mu=1$; man erhält hier also zugleich ein leichtes Kriterium, unter der Voraussetzung starker Stabilität. Dies liegt daran, daß das Vektorprodukt genau dann verschwindet, falls die Faktoren linear abhängig sind, siehe Gröbner (1966). Da die Berechnung von $\bigwedge_{i=1}^{n-1}a_i$ allerdings häufig über Determinanten geschieht, ist dieses Kriterium von der praktischen Rechnung nicht immer günstig. Wie starke Stabilität nachgewiesen wurde, sei hier dazu noch nicht einmal berücksichtigt. Für die Fehlerkonstante ergibt sich wegen $(**)$

$$ C = {\det(\widehat{\rho(1)},\gamma)\over \det(\widehat{\rho(1)},\rho'(1){\mskip 3mu}w)}. $$

Dies heißt, das Volumen der durch die linear unabhängigen Spaltenvektoren von $\rho(1)$ und dem Vektor $\gamma$ aufgespannten Körpers, ist der Zähler der Henricischen Fehlerkonstante.

Verschwindet der Zähler der Henricischen Fehlerkonstante, so liegt annullierte Dominanz vor. Albrecht (1979) nennt dies die Ordnungsbedingung. Durch Berechnen von $\det\left(\widehat{\rho(1)},\gamma\right)$ kann man dies also überprüfen. Diese Prüfung auf annullierte Dominanz kann man natürlich auch ohne den Umweg über das Kreuzprodukt, wie folgt herleiten. Aus

$$ v\ne0,\qquad v{\mskip 3mu}\rho(1)=0,\qquad v{\mskip 3mu}\gamma=0 $$

folgt, daß die zusammengesetzte Matrix $\left(\rho(1),\gamma\right)$ nicht maximalen Rang haben kann, also

$$ \mathop{\rm rank}\left(\rho(1),\gamma\right) = \kappa-n_1,\qquad \hbox{somit}\qquad \det\left(\widehat{\rho(1)},\gamma\right)=0. $$

Hierbei war $n_1$ die Vielfachheit des Eigenwertes $\mu=1$. Für Eigenwerte $\left|\mu\right|=1$ gilt allgemein $\mathop{\rm rank}\left(\rho(\mu),\gamma\right)=\kappa-n_\mu$, mit $n_\mu$ Vielfachheit des Eigenwertes $\left|\mu\right|=1$. Falls $n_\mu\ge1$ dann $\det\left(\widehat{\rho(1)},\gamma\right)=0$.

Damit sind alle denkbaren Fälle erschöpft, wenn man von Umnumerierungen absieht.

1. Beispiel: Zweistufiges, zyklisches lineares Mehrschrittverfahren mit zwei Startwerten. Die erste Stufe, mit Fehlerfaktor $\gamma_1$, sei

$$ \alpha_0^1y_{2n}+\alpha_1^1y_{2n+1}+\alpha_2^1y_{2n+2}=\ldots $$

und die zweite Stufe, mit Fehlerfaktor $\gamma_2$, sei

$$ \alpha_0^2y_{2n}+\alpha_1^2y_{2n+1}+\alpha_2^2y_{2n+2}+\alpha_3^2y_{2n+3}=\ldots{\mskip 3mu}. $$

Dann ist

$$ \rho(1) = \pmatrix{ \alpha_0^1+\alpha_2^1 & \alpha_1^1\cr \alpha_0^2+\alpha_2^2 & \alpha_1^2+\alpha_3^2\cr}. $$

2. Bedingung der annullierten Dominanz für zweistufige Verfahren. Für zweistufige Verfahren hat man

$$ \def\sumalf#1#2#3{\displaystyle\sum_{i\mathbin\%#3={#1}}\alpha_i^{#2}} \rho(1) = \pmatrix{ \sumalf012 & \sumalf112\cr \sumalf022 & \sumalf122\cr} =: \pmatrix{ \alpha_g^1 & \alpha_u^1\cr \alpha_g^2 & \alpha_u^2\cr} $$

und der Fehlervektor sei $\gamma=(\gamma_1,{\mskip 3mu}\gamma_2)$. Aufgrund der Konsistenz $\rho(1){\mskip 3mu}{1\choose 1}={0\choose0}$, ist

$$ \alpha_u^1+\alpha_g^1=0, \qquad \alpha_g^2+\alpha_u^2=0. $$

Ist nun die Matrix $\rho(1)$ nicht die Nullmatrix, so erhält man als Linkseigenvektor zu $\rho(1)$ natürlich

$$ \pmatrix{\alpha_g^1\cr \alpha_g^2\cr}^\bot = \pmatrix{-\alpha_g^2\cr \alpha_g^1\cr} $$

und damit als Bedingung für annullierte Dominanz

$$ \alpha_g^1\gamma_2 = \alpha_g^2\gamma_1. $$

Wäre jetzt $(\alpha_g^1,{\mskip 3mu}\alpha_g^2)=(0,{\mskip 3mu}0)$, so wäre die Matrix $\rho(1)$ gleich der Nullmatrix und die Bedingung der annullierten Dominanz führte zu der Bedingung, daß der Fehlervektor $\gamma$ sowohl auf $1\choose 0$, als auch auf $0\choose 1$ senkrecht stehen müßte. Damit wäre $\gamma_1=\gamma_2=0$, die Bedingung also leer. Die Konsistenzordnung im modifizierten Sinne wäre schon eine Ordnung höher als die wirklich erreichte Konvergenzordnung.

3. Beispiel: Verwendet man als erste Stufe das ^{Verfahren von Milne-Simpson} der Ordnung 4, mit

$$ 3y_{n+1} - 3y_{n-1} = z_{n+1} + 4z_n + z_{n-1} $$

und als zweite Stufe ein beliebiges Verfahren dritter Ordnung, so ist die Bedingung der annullierten Dominanz automatisch erfüllt, wegen

$$ \alpha_g^1 = 0, \qquad \gamma_1 = 0. $$

Das so gebildete zweistufige Verfahren konvergiert dann insgesamt mit der Ordnung 4.

4. Bedingung der annullierten Dominanz für dreistufige Verfahren. Sei der Fehlervektor des Verfahrens bezeichnet mit $\gamma=(\gamma_1,{\mskip 3mu}\gamma_2,{\mskip 3mu}\gamma_3)$ und der Matrix $\rho(1)$ sei gegeben durch

$$ \rho(1) = \pmatrix{ \sumalf013 & \sumalf113 & \sumalf213\cr \sumalf023 & \sumalf123 & \sumalf223\cr \sumalf033 & \sumalf133 & \sumalf233\cr} =: \pmatrix{ * & m_1 & n_1\cr * & m_2 & n_2\cr * & m_3 & n_3\cr} $$

wobei

$$ \eqalign{ m_i: &\quad\hbox{Summe der $\alpha_i$-Koeffizienten mit $i=3k+1$,}\cr n_i: &\quad\hbox{Summe der $\alpha_i$-Koeffizienten mit $i=3k+2$ .}\cr } $$

Als Bedingung für annullierte Dominanz ergibt sich nun

$$ \gamma_1(m_2n_3-m_3n_2)+\gamma_2(m_3n_1-m_1n_3)+\gamma_3(m_1n_2-m_2n_1)=0, $$

unter der Voraussetzung, daß $\rho(1)$ den Rang 2 hat.

Die Verallgemeinerung auf den $r$-stufigen Fall ergibt unmittelbar

$$ \rho(1) = \pmatrix{ \sumalf01r & \ldots & \sumalf{r-1}1r\cr \vdots & \ddots & \vdots\cr \sumalf0rr & \ldots & \sumalf{r-1}rr\cr } $$

9. Rechenregeln für Fehlerkonstanten

1. Henrici, Peter Karl Eugen (1923--1987). Hier wird nun allgemeiner eine Klasse von Fehlerkonstanten vorgestellt und die Beziehung zueinander werden aufgezeigt. Liegt das Verfahren $z_{n+1}=Az_n+h\varphi_n$ zugrunde mit Fehlervektor $\gamma$, so wird eine Fehlerkonstante definiert durch

$$ C_A = {\tilde v{\mskip 3mu}\tilde\gamma\over\tilde vw}, \qquad\left\{ \eqalign{\tilde vA&=\tilde v,\quad\tilde v\ne0,\cr Aw&=w,\quad w\ne0.\cr}\right. $$

Sei jetzt leicht allgemeiner $Lz_{n+1} = Uz_n + h\varphi_n$. Dann gilt

$$ C_A = {v{\mskip 3mu}\gamma\over vLw}, \qquad\left\{ \eqalign{v(L+U)&=0,\quad v\ne0\cr (L+U)w&=0,\quad w\ne0.\cr}\right. $$

Dies gilt wegen $A=-L^{-1}U$, daher

$$ A-I=-(L^{-1}U+I)=-L^{-1}(L+U). $$

Aus $\tilde v(A-I)=0$ folgt

$$ \tilde vL^{-1}(L+U) = (vL) L^{-1}(L+U)=0, $$

und damit ist $v$ Linkseigenvektor von $L\mu+U$ zum Eigenwert $\mu=1$, während $\tilde v=vL^{-1}$ und $\tilde\gamma=L^{-1}\gamma$ war. Vorausgesetzt ist natürlich, daß das Matrixpolynom $L\mu+U$ monisch ist, also $L$ invertierbar ist. Wegen der Null-Stabilität des Verfahrens ist das natürlich der Fall. Erkennbar ist auch, daß der Nenner nicht verschwinden kann, da die Matrix $A$ zur Klasse M gehört, siehe Ortega (1972).

Bibliographisch: Ortega, James McDonough.

Die obige Fehlerkonstante verallgemeinert sich sinngemäß bei mehrfachen Eigenwerten $\mu=1$.

Da die zyklischen Verfahren in der Dissertation von Tendler (1973) oder Tendler/Bickart/Picel (1978), in der Dissertation von Tischer (1983) und Tischer/Sacks-Davis (1983) und schließlich auch alle zyklischen Verfahren von Donelson/Hansen (1971) jedoch nur einen einfachen Eigenwert bei $\mu=1$ besitzen, wird dieser Fall hier nicht weiter verfolgt.

Bibliographisch: Peter E. Tischer, Ron Sacks-Davis, Donelson III, John (1941--2010), biography, Hansen, Eldon Robert (*1927), wiki.

Joel Marvin Tendler: "A Stiffly Stable Integration Process Using Cyclic Composite Methods", Ph.D. Diss., Syracuse University, Syracuse, New York, 26.Feb.1973, viii+iv+172 pages.

Peter E. Tischer: "The Cyclic Use of Linear Multistep Formulas for the Solution of Stiff Differential Equations", Ph.D. Thesis, Department of Computer Science, Monash University, Clayton, Victoria, Australia, August 1983, x+180 pages.

Im folgenden werden zwei Eigenschaften der Fehlerkonstanten von Henrici gezeigt. Zum einen ist die Fehlerkonstante von Henrici unabhängig von einer Skalierung des Verfahrens und zum anderen kann man sich auf den Fall einer Linearisierung des Matrixpolynomes beschränken. Für die praktische Rechnung von Linkseigenvektoren und weiteren Größen ist es natürlich günstiger das Verfahren in Form eines Matrixpolynomes mit möglichst geringer Dimension darzustellen. An anderer Stelle wiederum ist es angebrachter die Linearisierung zu betrachten, um nur mit einer einzigen Matrix zu hantieren. Daher ist es günstig für beide Darstellungen äquivalente Beschreibungen zur Verfügung zu haben.

Als nächstes wird nun also gezeigt, daß die Fehlerkonstante von Henrici unanhängig von einer Skalierung der Stufen ist. Jede Stufe darf beliebig mit einem Faktor $(\ne\!0)$ multipliziert werden, der gesamte Zyklus darf sogar einer nichtsingulären Skalierung unterzogen werden. Weiterhin ersieht man hieraus, daß die Fehlerkonstante von Henrici unabhängig von der Reihenfolge der Stufen ist. Für die umgedrehte Reihenfolge der Stufen wählt man beispielsweise einfach die Hankelmatrix $D=(\delta_{i,s+1-j})_{i,j=1}^s\in\mathbb{C}^{s\times s}$. Eine Vertauschung von Stufen innerhalb eines Zykluses kann natürlich sehr wohl die Anzahl der Startwerte ändern, u.U. kann sich also auch sogar die Anzahl der Matrizen im Matrixpolynom ändern. Dieser Fall ist dennoch mitberücksichtigt, da man ja das “alte” Verfahren mit Nullmatrizen ergänzen kann.

2. Satz: Sei $D\in\hbox{GL}(\mathbb{C},s)$ und sei $\hat\rho(\mu)=D\rho(\mu)$ das skalierte Polynom. Das zu dem Matrixpolynom $\rho$ gehörige Verfahren habe die Fehlerkonstante

$$ C={v{\mskip 3mu}\gamma\over v{\mskip 3mu}\rho'(1){\mskip 3mu}w}\qquad\left\{\eqalign{ v{\mskip 3mu}\rho(1)&=0,\quad v\ne0,\cr \rho(1){\mskip 3mu}w&=0,\quad w\ne0.\cr}\right. $$

Der Vektor $\gamma$ enthält hierbei die Fehlerfaktoren der Stufen. $\hat C$ sei die Fehlerkonstante von Henrici des skalierten Verfahrens. Dann gilt

$$ C = \hat C. $$

Beweis: Es ist $\hat\rho(1)=D\rho(1)$, $\hat\gamma=D\gamma$, und $\hat v:=vD^{-1}$ ist Linkseigenvektor von $\hat\rho(1)$. Die Ableitung des Matrixpolynomes $\rho$ an der Stelle 1 skaliert sich ebenfalls entsprechend, also $\hat\rho'(1)=D\rho'(1)$. Der Rechtseigenvektor $w$ ist auch gleichzeitig Rechtseigenvektor von $\hat\rho(1)$. Nun ist

$$ \hat C = {\hat v\hat\gamma\over\hat v{\mskip 3mu}\hat\rho'(1){\mskip 3mu}\hat w} = {vD^{-1}D\gamma\over vD^{-1}D\rho'(1){\mskip 3mu}w} = C. $$

    ☐

3. Sei jetzt allgemeiner statt des Matrixpolynoms $\rho(\mu)=L\mu+U$, betrachtet der Fall des Matrixpolynoms

$$ \rho(\mu) := A_0+A_1\mu+\cdots+A_\kappa\mu^\kappa. $$

Dann ist

$$ C_H = {v{\mskip 3mu}\gamma\over v{\mskip 3mu}\rho'(1){\mskip 3mu}w}, $$

eine Fehlerkonstante. Hierbei sind $v$ und $w$ entsprechend die Links- und Rechtseigenvektoren des Matrixpolynomes $\rho(\mu)$ zum dominanten Eigenwert $\mu=1$, es ist also

$$ v{\mskip 3mu}\rho(1)=0,\quad v\ne0\qquad\quad\hbox{und}\qquad\quad \rho(1){\mskip 3mu}w=0,\quad w\ne0. $$

Durch Wahl einer speziellen Norm und entsprechende Normierung des Vektors $w$ kann man dann auch den bestimmten Artikel benutzen.

In natürlicher und offensichtlicher Weise wird hiermit die klassische Fehlerkonstante von Henrici verallgemeinert. Auch hier kann der Nenner nicht verschwinden, da, wie unten gezeigt wird, diese Konstante mit der oben angegebenen Konstante äquivalent ist. Sind die Koeffizienten des Polynomes $A_i$ nicht von der Form $\alpha_i\otimes I$, so gilt nicht notwendig $\rho'(1)=\sigma(1)$, wie man anhand des folgenden Beispiels einsieht.

4. Beispiel: Die zyklische, zweimalige Hintereinanderausführung der BDF2 führt auf die Matrix mit den Einträgen

$$ \pmatrix{A_0, &A_1 &| &B_0, &B_1 \cr} $$

und den Elementen

$$ \left( \begin{array}{cccc|cccc} 1 & -4 & 3 & 0 & \tt & 0 & 0 & 2 & 0\cr 0 & 1 & -4 & 3 & \tt & 0 & 0 & 0 & 2\cr \end{array} \right). $$

Hier ist $\rho(\mu)=A_0+A_1\mu$ und $\sigma(\mu)=B_1\mu$. Offensichtlich gilt jetzt nicht $\rho'(1)=\sigma(1)$, da $\rho'(1)\equiv A_1\ne\sigma(1)\equiv B_1$ und dies obwohl jede Stufe die gleiche Konsistenzordnung hat, ja sogar alle Stufen gleich sind.

Nun wird gezeigt, daß alle angegebenen Fehlerkonstanten äquivalent sind. Um die nachstehenden Überlegungen durchsichtiger zu gestalten, soll anhand eines einfachen Beispieles unter anderem einige Eigenschaften der Begleitmatrix gezeigt werden.

5. Beispiel: Es sei das Polynom

$$ \def\aa{\alpha_0} \def\ab{\alpha_1} \def\ac{\alpha_2} \rho(\mu)=\aa+\ab\mu+\ac\mu^2+\mu^3 $$

vorgelegt, und es sei $\rho(1)=0$. Die Koeffizienten dieses Polynoms seien aus einem beliebigen Ring, nicht notwendig kommutativ, wobei 1 das neutrale Element bezeichne. Die Begleitmatrix zu $\rho$ sei

$$ C_1 = \pmatrix{ 0 & 1 & 0\cr 0 & 0 & 1\cr -\aa & -\ab & -\ac\cr}, \qquad\hbox{also}\qquad I-C_1 = \pmatrix{ 1 & -1 & 0\cr 0 & 1 & -1\cr \aa & \ab & 1+\ac\cr}. $$

Jetzt ist

$$ v = (\ab+\ac+1, \ac+1, 1) $$

Linkseigenvektor des Matrixpolynoms $I\mu-C_1$ zu $\mu=1$, wegen

$$ v(I-C_1) = \bigl( (\ab+\ac+1)+\aa,{\mskip 3mu}-(\ab+\ac+1)+(\ac+1)+\ab,{\mskip 3mu}-(\ac+1)+(\ac+1)\bigr) = (0,{\mskip 3mu}0,{\mskip 3mu}0), $$

da ja $\aa+\ab+\ac+1=0$, aufgrund $\rho(1)=0$. Wichtig ist noch zu vermerken, daß die Summe der Komponenten des Vektors $v$, gerade die Ableitung des Polynoms an der Stelle 1 ist, also es gilt

$$ v\pmatrix{1\cr 1\cr 1\cr} = (\ab+\ac+1)+(\ac+1)+1 = 3+2\ac+1\ab = \rho'(1). $$

Wegen $\rho(1)=0$ ist selbstverständlich $w^\top=(1,1,1)$ Rechtseigenvektor der Matrix $C_1$. Den Linkseigenvektor zu $C_1\in\mathbb{C}^{s\times s}$ kann man natürlich auch über das äußere Produkt erhalten. Aus der Matrix $(I-C_1)$ streicht man eine Spalte und ersetzt diese Spalte sukzessive $s$-mal durch den $i$-ten Einheitsvektor, für $i=1,\ldots,s$ und berchnet die $s$ Determinanten, also die Komponenten des äußeren Produktes.

Interessant ist in diesem Zusammenhang der nachstehende Zusammenhang zwischen Rechtseigenvektoren und Begleitmatrix, siehe Schäfke/Schmidt (1973), S.94.

Bibliographisch: Schäfke, Friedrich Wilhelm (1922--2010), Schmidt, Dieter (*1941).

6. Satz: Sei $C_1$ die Begleitmatrix des monischen Polynoms $\rho$ des Grades $n$. Ist $0\ne\mu\in\mathbb{C}$ Nullstelle von $\rho(\mu)$ der genauen Ordnung $r$, so liefert die vektorwertige Funktion $w_0\colon\mathbb{C}\to\mathbb{C}^n$ definiert durch

$$ w_0(\mu) := \pmatrix{1\cr \mu\cr \vdots\cr \mu^{n-1}\cr} $$

mit

$$ w_i(\mu) := {1\over i!}w^{(i)}(\mu), \qquad i=0,1,\ldots,r-1 $$

ein System von Rechts-Jordanvektoren zum Eigenwert $\mu$ von $C_1$, für welches also gilt

$$ (C_1-\mu I)w_i = \cases{ 0, &für $i=1$,\cr w_{i-1}, &für $i=2,3,\ldots,r$.\cr } $$

Beweis: Man geht aus von der Identität $(\lambda I-C_1)=\rho(\lambda)e_n$, wobei $e_n=(0,\ldots,0,1)\in\mathbb{C}^n$. $i$-malige Differentiation liefert

$$ \left(\lambda I-C_1\right) w_0^{(i)}(\lambda) = \rho^{(i)}(\lambda) e_n - i w_0^{(i-1)}(\lambda) $$

Einsetzen von $\lambda=\mu$, für $i=0,1,\ldots,r-1$ ergibt mit $\rho^{(r)}(\mu)=0$ (algebraische Vielfachheit von $\mu$), daß die $r$ linear unabhängigen Vektoren $w_i$ Hauptvektoren zu $\mu$ von $C_1$ sind.     ☐

Um nun eine gewisse Äquivalenz der Fehlerkonstanten zu zeigen, verfährt man wie nachstehend. Bei gewissen Einschränkungen an die Links- und Rechtseigenvektoren, kann man tatsächlich Gleichheit erzielen. Zumindestens Proportionalität ist stets gewährleistet.

7. Satz: Voraussetzungen: Es sei

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

($\ell{\buildrel\land\over=}$Stufenzahl), ferner sei $C_1\in\mathbb{C}^{\ell\kappa\times\ell\kappa}$ die erste Begleitmatrix zu $\rho(\mu)$, also

$$ 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_{\kappa-1}\cr } . $$

$v$ und $w$ seien beliebige aber feste Links- und Rechtseigenvektoren von $\rho(\mu)$ zu $\mu=1$ und

$$ \eqalignno{ v_c &:= v{\mskip 3mu}(A_1+A_2+\cdots+I, A_2+\cdots+I, \ldots, I) \in \mathbb{C}^{\kappa\ell}, \cr w_c &:= \pmatrix{I\cr \vdots\cr I\cr} w \in \mathbb{C}^{\kappa\ell}, \qquad I=I_\ell\in\mathbb{C}^{\ell\times\ell} . \cr } $$

$\gamma\in\mathbb{C}^\ell$ sei gänzlich beliebig und $\gamma_c = (0,\ldots,0,\gamma)^\top \in \mathbb{C}^{\kappa\ell}$, ($\gamma,\gamma_c{\buildrel\land\over=}$Fehlervektoren).

Behauptungen: (1) $v_c$ und $w_c$ sind Links- und Rechtseigenvektoren von $\rho_c(\mu) := I\mu - C_1$ zu $\mu=1$.

(2) Es gilt

$$ {v{\mskip 3mu}\gamma\over v{\mskip 3mu}\rho'(1){\mskip 3mu}w} = {v_c{\mskip 3mu}\gamma_c\over v_c{\mskip 3mu}w_c} = {v_c{\mskip 3mu}\gamma_c\over v_c{\mskip 3mu}\rho_c'(1){\mskip 3mu}w_c} . $$

Beweis: zu (1): Man sieht schnell, daß tatsächlich $v_c\rho_c(1)=0$ und $\rho_c(1)w_c=0$, mit

$$ \rho_c(1) = I - C_1 = \pmatrix{ I & -I & 0 & \ldots & 0\cr 0 & I & -I & \ldots & 0\cr \vdots & \vdots & \vdots & \ddots & \vdots\cr & & & \ldots & -I\cr A_0 & A_1 & & \ldots & I+A_{\kappa-1}\cr } . $$

zu (2): Gezeigt wird, daß Zähler und Nenner jeweils gleich sind. Für die Zähler ist dies unmittelbar klar. Für die Nenner rechnet man leicht nach, daß

$$ v_c{\mskip 3mu}\rho_c'(1){\mskip 3mu}w_c = v_c{\mskip 3mu}I_{\kappa\ell\times\kappa\ell}{\mskip 3mu}w_c = v_c{\mskip 3mu}w_c = v{\mskip 3mu}\rho'(1){\mskip 3mu}w . $$

    ☐

Die hier durchgeführten Überlegungen gelten sinngemäß in beliebigen, nicht notwendigerweise kommutativen Ringen. Hierzu ersetzt man $\mathbb{C}^\ell$ durch $\mathbb{R}$. Der obige Satz rechtfertigt in gewisser Hinsicht

8. Definition: Die Linearform $\gamma\mapsto v\gamma/v\rho'(1)w$ heißt Henrici-Linearform, mit $v$, $w$ wie oben. Insbesondere für einen Fehlervektor (spezieller Vektor des $\mathbb{C}^\ell$) heißt der Wert dann Henrici-Fehlerkonstante.

Selbstverständlich wird nicht behauptet, daß $v\gamma/v\rho'(1)w = v_c\gamma/v_c\rho_c'(1)w_c$ für beliebige Links- und Rechtseigenvektoren $v$, $w$, bzw. $v_c$, $w_c$. Dies erkennt man unmittelbar, falls man einen der Vektoren beliebig streckt oder staucht. Für den Zähler waren gewisse Unterraumeigenschaften von $\gamma_c$, nämlich $\gamma_c=(0,\ldots,0,*)^\top$ bedeutsam.

Die weiteren Verallgemeinerungen führen dann direkt zu den Begriffen der annullierten Dominanz und der Totalannullation. Um nun die Verbindung mit der klassischen Fehlerkonstante von Henrici weiter aufzuzeigen, sei auf den folgenden Sachverhalt hingewiesen.

Bei $m$-facher Wiederholung ein und desselben Verfahrens, multipliziert sich die oben angegebene Fehlerkonstante mit $m$. Dieses Ergebnis ergibt sich sofort, wenn man erkennt, daß

$$ (1,\ldots,1) \in \mathbb{C}^{1\times m} $$

Linkseigenvektor von $\rho(1)$ ist. Dann steht im Zähler $(1,\ldots,1){\mskip 3mu}\gamma$ und da jede Komponente von $\gamma$ natürlich gleich ist, erhält man sofort das verlangte Resultat, wenn man noch weiß, daß der Nenner natürlich bei welcher Dimension auch immer, gleich bleibt. Auch hier gelten wieder, w.o. schon bemerkt, diese Ergebnisse in beliebigen Ringen, nicht notwendig kommutativ.