, 48 min read

Praktische Gewinnung zyklischer, steif-stabiler Verfahren

In diesem Text werden die theoretischen Grundlagen, die u.a. zu den zyklischen Verfahren von Tendler (1973) führten, angeschnitten. Es werden die Begriffe $A[\alpha]$- und $S[\delta]$-Stabilität eingeführt. Die naheliegenden weiteren Definitionen von $A_\infty^r[\alpha]$- und $S_\infty^r[\delta]$-Stabilität werden ebenso diskutiert. Diese Definitionen weichen gelegentlich von den Definitionen ab, wie sie in der Literatur zu finden sind. Die Begriffe $A_\infty^r[\alpha]$- und $S_\infty^0[\delta]$-Stabilität werden hier nur so benutzt, wie sie an dieser Stelle erklärt werden. Diese Stabilitätsmerkmale ermöglichen einen Vergleich verschiedener Verfahren gleicher Ordnung und damit eine Abschätzung über ihre Brauchbarkeit bei einem Einsatz in einem Programmpaket. Dennoch setzen diese Stabilitätsdefinitionen gewisse idealisierte Bedingungen voraus. Aus diesem Grunde ist es ratsam, ihnen nicht unangemessen großes Gewicht beizumessen. Erfüllt also ein Verfahren diese Eigenschaft besonders “gut”, so heißt dies nicht automatisch, daß es sich in einem Rechenprogramm auch immer “gut” verhalten würde.

Schließlich wird dargelegt, wie die neuen zyklischen Formeln generiert wurden. Das letzte ist von Interesse, wenn versucht werden sollte, die von J.M. Tendler gefundenen Formeln weiter zu verbessern. Von seinen prinzipiellen Grundbestandteilen her sind die Tendlerschen Zyklen vergleichsweise einfach aufgebaut, nämlich als lineares Komposituum linear unabhängiger Mehrschrittformeln. Dennoch ist die tatsächliche Erzeugung, insbesondere die systematische Generierung nicht ganz einfach. Es handelt sich hier nicht um einen einfach erlernbaren Kalkül, sondern um einen iterativen Suchvorgang.

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

1. Definition von Stabilitätseigenschaften

$\ldots$ wird sich diese Störung allerdings erst später einstellen, aber auf alle Fälle wird die Integration nach ausreichend vielen Schritten gründlich verpfuscht. Man nennt diese Erscheinung numerische Instabilität. Sie ist nicht etwa eine Folge der Akkumulation von Rundungsfehlern, denn die Formel beschreibt den Integrationsprozeß unter der Annahme, daß vollständig exakt gerechnet wird.

Eduard Stiefel (1968)

Ausschließlich Differentialgleichungen vom Typ

$$ \dot y = f(t,y), \quad\hbox{mit dem Anfangswert}\quad y(t_0) = y_0 $$

werden hier betrachtet. Es sei für diesen und alle weiteren Abschnitte angenommen, daß $f$ und damit $y$ so oft differenzierbar sind, wie es für die jeweilige Betrachtung nötig ist. Es wird nun $t_n=t_0+nh$ gesetzt, und $y_n$ bezeichnet die berechneten Näherung für $y(t_n)$, für $n=m\ell+j$ und $j=1,\ldots,\ell$, die durch die folgenden $\ell$ Mehrschrittformeln bestimmt wurden

$$ \sum_{j=-k+1}^{\ell} \left( \alpha_{ij}y_{m\ell+j} - h\beta_{ij}\dot y_{m\ell+j} \right) = 0, \qquad i=1,\ldots,\ell, $$

wobei $\dot y_i = f(t_i,y_i)$ ist.

Ausführlich aufgeschrieben lautet der obige Ausdruck z.B. für ein dreistufiges Verfahren mit drei Startwerten, also $\ell=3$, und $k=3$, wie untenstehend. Hier und im weiteren wird häufig die Abkürzung $z_i=hf_i=h\dot y_i$ benutzt.

$$ \eqalignno{ & \alpha_{1,-2}y_{3m-2}+\alpha_{1,-1}y_{3m-1}+\alpha_{10}y_{3m}+ \alpha_{11}y_{3m+1}+\alpha_{12}y_{3m+2}+\alpha_{13}y_{3m+3} \cr = & \beta_{1,-2}z_{3m-2}+\beta_{1,-1}z_{3m-1}+\beta_{10}z_{3m}+ \beta_{11}z_{3m+1}+\beta_{12}z_{3m+2}+\beta_{13}z_{3m+3} \cr & \alpha_{2,-2}y_{3m-2}+\alpha_{2,-1}y_{3m-1}+\alpha_{20}y_{3m}+ \alpha_{21}y_{3m+1}+\alpha_{22}y_{3m+2}+\alpha_{23}y_{3m+3} \cr = & \beta_{2,-2}z_{3m-2}+\beta_{2,-1}z_{3m-1}+\beta_{20}z_{3m}+ \beta_{21}z_{3m+1}+\beta_{22}z_{3m+2}+\beta_{23}z_{3m+3} \cr & \alpha_{3,-2}y_{3m-2}+\alpha_{3,-1}y_{3m-1}+\alpha_{30}y_{3m}+ \alpha_{31}y_{3m+1}+\alpha_{32}y_{3m+2}+\alpha_{33}y_{3m+3} \cr = & \beta_{3,-2}z_{3m-2}+\beta_{3,-1}z_{3m-1}+\beta_{30}y_{3m}+ \beta_{31}z_{3m+1}+\beta_{32}z_{3m+2}+\beta_{33}z_{3m+3}.\cr } % Ende eqalign $$

Der obige Ausdruck wird zyklisch, falls gilt

$$ \alpha_{ij} = 0 \quad\land\quad \beta_{ij} = 0 \qquad (1\le i \lt j \le \ell). $$

Gälte $\alpha_{ij}\ne0$ oder $\beta_{ij}\ne0$ für ein $i<j$, so wäre das zusammengesetzte Verfahren blockimplizit. Der Aufwand zur “Beseitigung” der Implizitheit der Formeln wäre dann größer als bei (einstufigen) linearen Mehrschrittformeln.

Bei zyklischen Verfahren ist der Aufwand zur Berechnung der Näherungen nicht größer als wie bei Mehrschrittverfahren. Andererseits hat man jedoch nun die Möglichkeit, da nicht alle Stufen gleich gewählt werden müssen, Stabilitätseigenschaften zu erreichen die u.U. den Stufen überlegen sind. Die Stufen müssen untereinander nicht verschieden sein. Einige oder auch alle dürfen übereinstimmen. Praktische Ergebnisse deuten auf den merkwürdigen Umstand hin, daß die dreifache Wiederholung der BDF2 als zyklisches Verfahren entsprechend programmiert, ggü. der dreifachen Anwendung von 3 Schritten mit der BDF2 mit einem Programm konzipiert für einstufige Verfahren, Vorteile bieten kann.

Während bis jetzt noch jede Stufe des Zykluses $2(k+i)$ Parameter, d.h. noch nicht bestimmte Koeffizienten, enthält, sei noch zusätzlich angenommen, um die die Anzahl der der Bedingungen zu reduzieren, daß gelte

$$ \alpha_{ij} = 0 \quad\land\quad \beta_{ij} = 0 \qquad (-k+i\gt j,\quad i=1,\ldots,\ell). \tag{LMSV} $$

D.h., daß jede Stufe ein lineares $k$-Schrittverfahren ist. Es verbleiben jetzt noch $2(k+1)$ Parameter. Aus Vereinfachungsgründen wird nun verlangt, daß jede Stufe die gleiche Konsistenzordnung $p=k$ hat. Berücksichtigt man die Konsistenzbedingungen, so verbleiben noch $k+1$ Parameter. Diese werden dazu benutzt die gewünschten Stabilitätseigenschaften des Verfahrens zu erreichen.

Zur Untersuchung des Lösungsverhaltens der Differenzengleichung, welche man zur numerischen Lösung von Differentialgleichungen einsetzt, betrachtet man die sehr spezielle Testgleichung

$$ \dot y = \lambda y, \qquad y(t_0) = y_0. $$

Da man durch Linearisierung viele Differentialgleichung zumindestens in einer kleinen Umgebung auf diese obige Form zurückführen kann, liefert die Untersuchung dieser speziellen Gleichung einen ersten Eindruck was man wohl im allgemeineren Falle zu erwarten hat. Desweiteren kann man damit gewisse Ineffizienzen von Programmen verstehen und erklären. Zur Wiederholung seien in sehr salopper Form die beiden bekanntesten Stabilitätsdefinitionen angegeben:

Eine Formel heißt stabil, wenn die sämtlichen berechneten Lösungen für $H=0$ nicht unbeschränkt wachsen können; $H=h\lambda$.

Die Formel heißt $A$-stabil, wenn sämtliche berechneten Lösungen eine Nullfolge bilden, für alle $H\in\mathbb{C}^-$. Ausführlicher

1. Definition: nach Germund Dahlquist, Dahlquist (1963). Ein lineares $k$-Schrittverfahren heißt $A$-stabil, falls alle Lösungen von

$$ \sum_{i=0}^k \alpha_i y_{n+i} = \sum_{i=0}^k \beta_i f_{n+i} $$

gegen Null gehen ($n\to\infty$), wenn das Verfahren mit fester Schrittweite $h$ auf Differentialgleichungen der Form $\dot y=\lambda y$ angewendet wird, mit einer festen komplexen Konstante $\lambda$, mit negativen Realteil.

Die Beschränkung auf lineare $k$-Schrittverfahren ist unwesentlich. Die Beschäftigung mit zyklischen linearen Verfahren hat zu einem nicht unmaßgeblichen Teil ihren Ursprung in der zweiten Dahlquist-Barriere. Diese besagt, daß $A$-stabile, lineare Mehrschrittformeln keine beliebig hohe Ordnung erreichen können: maximal die Ordnung zwei. Durch die drastische Einschränkung der Vielfalt $A$-stabiler, linearer Mehrschrittverfahren werden zwei Erweiterungen nahegelegt: Abschwächung des Begriffes der $A$-Stabilität und Hinzufügung größerer Klassen von Verfahren. Zur Quantifizierung von Stabilitätseigenschaften einerseits und zur Abschwächung der $A$-Stabilität andererseits dienen die nachstehenden Definitionen.

2. Definition: Bezeichne ${\cal S}[\delta]$ eine zusammenhängende und offene Teilmenge der komplexen Ebene mit dem Ursprung als Randpunkt, die mindestens einen unendlich großen Teil der linken Halbebene enthält mit Elementen $z$ der Form

$$ \mathop{\rm Re}\nolimits z \lt -\delta \le 0, \qquad \delta\ge0 $$

Nun heißt ein Verfahren $S[\delta]$-stabil, falls so eine Menge ${\cal S}[\delta]$ existiert derart, daß wenn man das Verfahren auf die lineare und zeitinvariante Differentialgleichung $\dot y=\lambda y$ anwendet, die mit dem Verfahren berechneten Lösungen $y_n$ gegen Null streben, wenn $n\to\infty$, für eine feste positive Schrittweite $h$, falls immer $\lambda h=H \in {\cal S}[\delta]$. Existiert kein unendlich großer Bereich der linken Halbebene, so werde von nicht $S[\delta]$-stabil gesprochen. Im Falle der Existenz heißt das betragsmässig kleinste $\delta$ die Widlund-Distanz, nach Olof B. Widlund.

Diese Definition hat eine sehr anschauliche Interpretation, wenn man das Gebiet ${\cal S}[\delta]$ aufzeichnet. Von besonderem Interesse ist hier das betragsmässig kleinstmögliche $\delta$. Nur von diesem kleinstmöglichen $\delta$ wird im weiteren überhaupt noch die Rede sein. Gleiche Anschaulichkeit hat auch die folgende Definition.

3. Definition: Sei ${\cal A}[\alpha]$ eine zusammenhängende und offene Teilmenge der komplexen Zahlenebene mit Elementen $z$ der Form

$$ \mathopen|\arg \{-z\}\mathclose| \lt \alpha \le 90^\circ, \qquad \alpha\gt 0. $$

Jetzt heißt das Verfahren $A[\alpha]$-stabil, wenn es w.o. eine derartige Menge ${\cal A}[\alpha]$ gibt derart, daß die berechneten Näherungen dieses Verfahren gegen Null streben, für $n\to\infty$, wenn man das Verfahren auf die Differentialgleichung $\dot y = \lambda y$ anwendet. Diese Abklingeigenschaft soll immer dann gelten, wenn $\lambda h=H \in {\cal A}[\alpha]$ gilt, für feste positive Schrittweite $h$. Für $\alpha\notin\left[0,90^\circ\right]$ werde von nicht $A[\alpha]$-stabil gesprochen. Im Falle der Existenz nennt man den größtmöglichen Winkel $\alpha$ auch den Widlund-$\alpha$-Winkel, ebenfalls nach Olof B. Widlund.

In sinnfälliger Verallgemeinerung nennt man ein Verfahren $A[0]$-stabil, wenn es ein $\alpha$ gibt, sodaß das Verfahren $A[\alpha]$-stabil ist. ${\cal A}[0]$ ist in $\mathbb{C}$ nicht offen. Wie die obige Definition kann man diese Eigenschaft wiederum sehr anschaulich deuten. Hier gilt je größer der Winkel $\alpha$, desto besser. Bei der vorherigen Definition galt je kleiner $\delta$, desto günstiger. Durch ein größeres $\alpha$, bzw. kleineres $\delta$, werden die entsprechenden Mengen ${\cal A}[\alpha]$, bzw. ${\cal S}[\delta]$, “größer”.

Bei den beiden Defintionen beachte man, daß ${\cal S}[\delta]$ die Menge bezeichnet und $S[\delta]$ die dazugehörige Eigenschaft. Entsprechend gilt dies für ${\cal A}[\alpha]$ und $A[\alpha]$. Man stellt schnell fest, daß ein Verfahren stets $A[\alpha]$-stabil ist, für ein gewisses $\alpha(\delta)$, wenn es $S[\delta]$-stabil ist. Die Umkehrung ist nicht notwendig richtig. Bei den noch folgenden Eigenschaften $A_{\infty}[\alpha]$ und $S_{\infty}[\delta]$ wird dies jedoch der Fall sein. Ist $\delta = 0$ oder $\alpha = 90^\circ$, so ist dies das gleiche wie $A$-Stabilität. Beide Werte führen zur selben Menge, nämlich der gesamten linken komplexen Halbebene ohne die imaginäre Achse.

4. Beispiele: (1) Die impliziten BDF der Ordnung $p=1$ bis $p=6$ erfüllen mit gewissen $\alpha$ und $\delta$ die obigen Stabilitätseigenschaften.

(2) Sämtliche zyklischen Verfahren von Tendler (1973), oder Tendler/Bickart/Picel (1978) sind ebenfalls Beispiele für Verfahren, die $A[\alpha]$- bzw. $S[\delta]$-stabil sind. Zusätzlich sind die Formeln von Tendler auch $A_\infty^0[\alpha]$- und $S_\infty^0[\delta]$-stabil. (Definition folgt.)

(3) Genauso sind auch alle zyklischen Formeln von Peter E. Tischer (1983) und Tischer/Sacks-Davis (1983) alle $A[\alpha]$- bzw. $S[\delta]$-stabil, sogar $A_\infty^0[\alpha]$- und $S_\infty^0[\delta]$-stabil.

(4) Die 10 blockimpliziten Verfahren von Bickart/Picel (1973) sind alle $A[\alpha]$- und $S[\delta]$-stabil, sogar $A_\infty^0[\alpha]$- und $S_\infty^0[\delta]$-stabil.

(5) Explizite Verfahren mit konstanter Schrittweite erfüllen die obigen Eigenschaften nicht. Das explizite Euler-Verfahren mit konstanter Schrittweite ist also weder $A[\alpha]$- noch $S[\delta]$-stabil und zwar für kein $\alpha$ und kein $\delta$.

Bibliographisch: Tischer, Peter E.: "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 S.

Tischer, Peter E. und Sacks-Davis, Ron: "A New Class of Cyclic Multistep Formulae for Stiff Systems", SIAM Journal on Scientific and Statistical Computing, Vol 4, No 4, December 1983, pp.733—747.

Bickart, Theodore A. und Picel, Zdenek: "High Order Stiffly Stable Composite Multistep Methods for Numerical Integration of Stiff Differential Equations", BIT, Vol 13, 1973, pp.272—286

2. Schreibweise zusammengesetzter Mehrschrittformeln

1. Wann nun ein Verfahren diese Eigenschaften hat, wird durch sein charakteristisches Polynom bestimmt. Dieses Polynom in zwei komplexen Veränderlichen erhält man nun wie folgt.

Den gesamten $\ell$-stufigen Zyklus

$$ \sum_{j=-k+1}^{\ell} \left[ \alpha_{ij}y_{m\ell+j} - h\beta_{ij}\dot y_{m\ell+j} \right] = 0,\qquad\hbox{für}\quad i=1,\ldots,\ell, $$

kann man als eine einzige Matrixdifferenzengleichung notieren, nämlich

$$ \sum_{i=0}^\kappa \left[ A_i Y_{m+i} - h B_i \dot Y_{m+i} \right] = 0, % \eqno \hbox{(Matrix)} $$

mit $\kappa=\lceil k/\ell \rceil$. $\kappa$ ist hierbei die kleinste ganze Zahl, die größer oder gleich $k/\ell$ ist. Die insgesamt $\kappa+1$ Matrizen $A_i$, $B_i$ haben die Größe $\ell\times\ell$. Es genügen genau dann zwei Matrizen $A_1$ und $A_0$, bzw. $B_1$ und $B_0$, wenn $\kappa=1$ ist, also $k\le\ell$.

Vergrößert man die Dimensionen der Matrizen $A_i$ und $B_i$ mit Hilfe von sogenannten Trivialstufen ($1=1$ Stufen), so läßt sich natürlich die Anzahl der Matrizen verringern, z.B. auf genau zwei, also $A_0$ und $A_1$, bzw. $B_0$ und $B_1$. Bei den zyklischen Verfahren von J.M. Tendler ist stets $\kappa\le2$, d.h. man kommt mit jeweils höchstens drei Matrizen zur Beschreibung der Matrixdifferenzengleichung aus.

2. Die Matrizen $A_i$ und $B_i$, für $i=0,\ldots,\kappa$ lauten:

$$ A_\kappa = \pmatrix{ \alpha_{11} & \ldots & \alpha_{1\ell}\cr \vdots & \ddots & \vdots\cr \alpha_{\ell1} & \ldots & \alpha_{\ell\ell}\cr }, \quad B_\kappa = \pmatrix{ \beta_{11} & \ldots & \beta_{1\ell}\cr \vdots & \ddots & \vdots\cr \beta_{\ell1} & \ldots & \beta_{\ell\ell}\cr }, \quad Y_{m+\kappa} = \pmatrix{y_{m\ell+1}\cr \vdots\cr y_{m\ell+\ell}}, $$

und

$$ A_{\kappa-1} = \pmatrix{ \alpha_{1,-\ell+1} & \ldots & \alpha_{10}\cr \vdots & \ddots & \vdots\cr \alpha_{\ell,-\ell+1} & \ldots & \alpha_{\ell0}\cr }, \quad B_{\kappa-1} = \pmatrix{ \beta_{1,-\ell+1} & \ldots & \beta_{10}\cr \vdots & \ddots & \vdots\cr \beta_{\ell,-\ell+1} & \ldots & \beta_{\ell0}\cr }, \quad Y_{m+\kappa-1} = \pmatrix{y_{m\ell-\ell+1}\cr \vdots\cr y_{m\ell}}, $$

beziehungsweise

$$ A_{\kappa-2} = \pmatrix{ \alpha_{1,-2\ell+1} & \ldots & \alpha_{1,-\ell}\cr \vdots & \ddots & \vdots\cr \alpha_{\ell,-2\ell+1} & \ldots & \alpha_{\ell,-\ell}\cr }, \quad B_{\kappa-2} = \pmatrix{ \beta_{1,-2\ell+1} & \ldots & \beta_{1,-\ell}\cr \vdots & \ddots & \vdots\cr \beta_{\ell,-2\ell+1} & \ldots & \beta_{\ell,-\ell}\cr }, \quad Y_{m+\kappa-2} = \pmatrix{y_{m\ell-2\ell+1}\cr \vdots\cr y_{m\ell-\ell}\cr} $$

und so weiter bis zu den Matrizen $A_0$, $B_0$ und Vektoren $Y_m$. Allgemein also

$$ A_i = \pmatrix{ \alpha_{1,(i-\kappa)\ell+1} & \ldots & \alpha_{1,(i-\kappa)\ell+\ell}\cr \vdots & \ddots & \vdots\cr \alpha_{\ell,(i-\kappa)\ell+1} & \ldots & \alpha_{\ell,(i-\kappa)\ell+\ell}\cr }, \qquad Y_{m+i} = \pmatrix{y_{(m+i-\kappa)\ell+1}\cr \vdots\cr y_{(m+i-\kappa)\ell+\ell}\cr}, \qquad \dot Y_{m+i} = \pmatrix{f_{(m+i-\kappa)\ell+1}\cr \vdots\cr f_{(m+i-\kappa)\ell+\ell}\cr} $$

für $i=0,\ldots,\kappa$. Die $B_i$ analog wie die $A_i$.

3. Bei einer vektorwertigen Funktion mit $n$ Komponenten hat man die oben angegebenen Matrizen zu ersetzen durch

$$ A_i\gets A_i\otimes I_n \in \mathbb{R}^{n\ell\times n\ell}, \qquad B_i\gets B_i\otimes I_n \in \mathbb{R}^{n\ell\times n\ell}. $$

Hierbei ist $I_n$ die $n\times n$-Einheitsmatrix und $\otimes$ bezeichnet das Kroneckerprodukt, nach Leopold Kronecker (1823—1891). In den meisten Fällen kann man auf diese gesonderte Notierung verzichten, u.a. weil gilt $\left|A\otimes I_n\right| = \left|A\right|$.

4. Beispiele: (1) Im Falle einstufiger Zyklen, also $\ell=1$, reduziert sich die obige Matrixdifferenzengleichung zu der gewöhnlichen skalaren Differenzengleichung für lineare Mehrschrittverfahren, mit $A_i=\alpha_i$, $B_i=\beta_i$ und $Y_i=y_i$.

(2) Das lineare zweistufige zyklische $A_\infty^0$-stabile Verfahren der Ordnung 4, also $k=4$, $p=4$ und $\ell=2$, von Tischer (1983) und Tischer/Sacks-Davis (1983), soll in der obigen Schreibweise angegeben werden. Die zu den Werten

$$ \left( \begin{array}{cccccc|cccccc} y_{2n-3} & y_{2n-2} & y_{2n-1} & y_{2n} & y_{2n+1} & \hbox{---} & z_{2n-3} & z_{2n-2} & z_{2n-1} & z_{2n} & z_{2n+1} & \hbox{---}\cr y_{2n-3} & y_{2n-2} & y_{2n-1} & y_{2n} & y_{2n+1} & y_{2n+2} & z_{2n-3} & z_{2n-2} & z_{2n-1} & z_{2n} & z_{2n+1} & z_{2n+2}\cr \end{array} \right) $$

gehörigen Koeffizienten lauten

$$ \left( \begin{array}{cccccc|cccccc} 0.6 & -2.7 & 4.7 & -3.5 & 1 & 0 & 0.9 & -4.1 & 7.1 & -5.6 & 1.7 & 0\cr 0 & 0 & -0.3 & 1.2 & -2.0 & 1 & -2.6 & 11.9 & -20.9 & 17.9 & -7.6 & 1.7\cr \end{array} \right) $$

wobei alle Koeffizienten auf eine Stelle hinter dem Komma gerundet worden sind, also der “Konsistenzschnelltest” $\rho(1)=0$ hier immer $\pm0.1$ ergibt.

In der Schreibweise als Matrixdifferenzengleichung w.o. ergibt sich nun unmittelbar

$$ \begin{array}{cccccc} & \pmatrix{0.6 & -2.7\cr 0 & 0\cr}\pmatrix{y_{2n-3}\cr y_{2n-2}\cr} &+&\pmatrix{4.7 & -3.5\cr -0.3 & 1.2\cr}\pmatrix{y_{2n-1}\cr y_{2n}\cr} &+&\pmatrix{1 & 0\cr -2.0 & 1\cr}\pmatrix{y_{2n+1}\cr y_{2n+2}\cr}\cr \cr %\noalign{\phantom{:}} =& \pmatrix{0.9 & -4.1\cr -2.6 & 11.9\cr}\pmatrix{z_{2n-3}\cr z_{2n-2}\cr} &+&\pmatrix{7.1 & -5.6\cr -20.9 & 17.9\cr}\pmatrix{z_{2n-1}\cr z_{2n}\cr} &+&\pmatrix{1.7 & 0\cr -7.6 & 1.7\cr}\pmatrix{z_{2n+1}\cr z_{2n+2}\cr}\cr \end{array} $$

Eine typische Eigenschaft der Formeln von Tischer (1983) und Tischer/Sacks-Davis (1983) ist die relativ geringe Äquilibrierung der $\beta_{ij}$-Werte. Bei diesem Verfahren beträgt sie über 1:20. Für die Ordnung $p=8$ hat man sogar Verhältnisse von 1:130. Maßgeblich betrifft dies stets die letzte Stufe der immer zweistufigen Zyklen.

Die Formeln von Tendler (1973) besitzen diese Eigenschaften nicht. Bei gewöhnlichen Quadratur- und Kubaturformeln für $y=\int_a^b f(x) dx$ ist bekannt, daß diese Eigenschaft unerwünscht ist.

(3) Das zweistufige, zyklische Verfahren DH3 von Donelson/Hansen (1971)

Donelson III, John und Hansen, Eldon: "Cyclic Composite Multistep Predictor-Corrector Methods", SIAM Journal on Numerical Analysis, Vol 8, 1971, pp.137—157

$$ \begin{array}{r|rr} p=5 & 1 & 2\cr \hline -2 & 0 & 0\cr -1 & -57 & 31\cr 0 & 24 & -12\cr 1 & 33 & -39\cr 2 & 0 & 20\cr \hline -2 & -1 & 0\cr -1 & 24 & -10\cr 0 & 57 & -39\cr 1 & 10 & 12\cr 2 & 0 & 7\cr \end{array} $$

Die Stabilität definierenden Matrizen sind:

$$ \eqalign{ A_0 &= \pmatrix{0&0\cr 0&0\cr},\cr &\cr A_1 &= \pmatrix{-57&24\cr 31&-12\cr},\cr &\cr A_2 &= \pmatrix{33&0\cr -39&20\cr},\cr } \qquad \eqalign{ B_0 &= \pmatrix{0&-1\cr 0&0\cr},\cr &\cr B_1 &= \pmatrix{24&57\cr -10&-39,\cr}\cr &\cr B_2 &= \pmatrix{10&0\cr 12&7\cr},\cr } $$

Die charakteristische Gleichung ist

$$ \det\left(A_0+A_1\mu+A_2\mu^2\right) = -660(\mu-1)(\mu+1/11)\mu^2. $$

Keinesfalls müssen zyklische Verfahren nur zweistufig sein. Die Stufenzahl eines zyklischen Verfahrens unterliegt keiner Schranke, außer, daß diese Anzahl nicht unendlich ist und daß diese Anzahl dabei fest ist. Gleichzeitig erkennt man, daß eine formale Übertragung $\det Q(\mu,H)={\cal O}(H^{p+1})$ ($\mu=e^H$) scheitert, weil z.B. für DH3 gilt

$$ \det Q(\mu,H) = -720H - 3060H^2 - 5838H^3 - 7090H^4 - 6303H^5 - {52985\over12}H^6 + {\cal O}(H^7) \ne {\cal O}(H^6), $$

nach Substitution von $\mu=e^H$.

5. Es ist natürlich ebenso möglich die Linearisierung des obigen Matrixpolynoms zu betrachten, jedoch werden dann die Dimensionen der Matrizen größer, falls $\kappa>1$, also $k>\ell$ ist. Das Verfahren schreibt sich dann in der Form $z_{n+1}=Az_n+h\varphi_n$, oder $Lz_{n+1}=Uz_n+h\varphi_n$.

3. Stabilität und Polynome

1. Das charakteristische Polynom lautet

$$ \det Q(\mu,H) = \det\Biggl\{\sum_{i=0}^\kappa (A_i - H B_i) \mu^i \Biggr\}. % \eqno \hbox{(char.Pol.)} $$

Man beweist jetzt, daß das verwendete Verfahren genau dann $A[\alpha]$-stabil bzw. $S[\delta]$-stabil ist, wenn die sämtlichen Nullstellen $\mu$ von $\det Q(\mu,H)=0$, alle betragsmässig kleiner eins sind, für alle $H \in {\cal A}[\alpha]$, bzw. alle $H \in {\cal S}[\delta]$. Die Eigenschaft eines Polynomes $Q(\mu,H)$ für ein spezielles $H_0$ sämtliche Wurzeln in der offenen Einheitskreisscheibe zu haben, nennt man absolut-stabil. ${\cal A}[\alpha]$ und ${\cal S}[\delta]$ sind also Mengen der absoluten Stabilität.

Äquivalent hierzu ist, daß das Verfahren genau dann $A[\alpha]$-stabil ist, bzw. $S[\delta]$-stabil ist, falls alle Nullstellen $H$ von $\det Q(\mu,H)=0$ im Komplement von ${\cal A}[\alpha]$, bzw. im Komplement von ${\cal S}[\delta]$, liegen, für alle $\mu$ aus dem Komplement der offenen Einheitskreisscheibe. Man erhält somit den folgenden Test.

2. Test A: Ein Verfahren ist genau dann $A[\alpha]$-stabil, bzw. $S[\delta]$-stabil, falls

  1. die Nullstellen $\mu$ von $\det Q(\mu,H_0)$ für ein festgehaltenes beliebiges $H_0 \in {\cal A}[\alpha]$, bzw. $H_0 \in {\cal S}[\delta]$, sämtlich betragsmässig kleiner 1 sind und gleichzeitig
  2. die Kurve der Punkte $H$, die $\det Q(\mu,H)|_{\mathopen|\mu\mathclose|=1} = 0$ erfüllen, im Komplement von ${\cal A}[\alpha]$, bzw. im Komplement von ${\cal S}[\delta]$ liegen.

Leider sind Verfahren, die die obigen beiden Eigenschaften erfüllen, noch nicht wirklich gut zur Integration von steifen Differentialgleichungen geeignet. Dies liegt u.a. daran, daß wenn $\mathopen|H\mathclose|$ $(\mathop{\rm Re}\nolimits H<0)$ immer größer wird, die Nullstellen $\mu$ von $\det Q(\mu,H)=0$ betragsmässig gegen eins streben könnten. Integriert man jedoch $\dot y = \lambda y$ ($\mathop{\rm Re}\nolimits \lambda<0$) so würde dies heißen, daß es eine Lösung gibt mit $\mathopen|y_{n+1}/y_n\mathclose| \to 1$. Diese Eigenschaft spiegelt sicherlich nicht das wahre Verhalten der Lösung wider. Diese Eigenheit hat z.B. die Trapezregel

$$ y_{m+1} = y_m + (h/2)(f_{m+1} + f_m). $$

Aus diesem Grunde sind Verfahren günstiger bei denen die Nullstellen $\mu$ bei $\mathop{\rm Re}\nolimits H=-\infty$ von $\det Q(\mu,H)=0$ in der offenen Einheitskreisscheibe liegen, also keinerlei Wurzeln des Betrages eins auftauchen. Damit hätte man die Trapezregel als zwar $A$-stabiles Verfahren, letztlich ausgesondert. Um diesen Wunsch zu präzisieren und um handliche Begriffe bereitzustellen, seien diese Anforderungen nun nomenklaturmässig genauer festgelegt.

4. Verschärfung von Stabilitätsdefinitionen

Suppose, for example, that it is known that the planetary system is stable “in the past”. If it captures a new body, say, a speck of dust arriving from “infinity”, then the newly formed system of bodies looses its stability: with probability one, either a collision occurs, or one of the bodies escapes again at infinity. Moreover, it is by no means necessary that the speck of dust is the one that leaves the Solar system: Jupiter, or even the Sun may equally well leave.

V.I. Arnold, V.V. Kozlov, A.I. Neishtadt (1988), Mathematical Aspects of Classical and Celestial Mechanics

Bibliographisch: Vladimir Arnold (1937—2010), Valery Vislevich Kozlov (*1950), Anatoly Neishtadt.

Es liegt nun nahe, die Begriffe $A[\alpha]$-Stabilität und $S[\delta]$-Stabilität zu verfeinern. Hierzu wird das Verhalten der Differenzengleichung für $\mathop{\rm Re}\nolimits H\to-\infty$ mit in die Nomenklatur integriert. Dies führt jetzt zur

1. Definition: (1) Ein $S[\delta]$-stabiles Matrixpolynom $Q(\mu,H)$ mit der Eigenschaft, daß die sämtlichen Nullstellen $\mu$ für $\mathop{\rm Re}\nolimits H\to-\infty$ in der offenen, Kreisscheibe $\{z\in\mathbb{C}: \mathopen|z\mathclose|<r\}$ liegen, nennt man $S_\infty^r[\delta]$-stabil. Für $r=1$ sagt man abkürzend nur $S_\infty[\delta]$-stabil.

(2) Analog definiert man $A_\infty^r[\alpha]$-Stabilität als Erweiterung der $A[\alpha]$-Stabilität. Erneut bezeichnet $A_\infty[\alpha]$-Stabilität die Abkürzung für $A_\infty^1[\alpha]$-Stabilität.

(3) Für $r=0$ ist $r\to0$ gemeint, also: Die sämtlichen Nullstellen $\mu$ für $\mathop{\rm Re}\nolimits H\to-\infty$ verschwinden. Hierfür schreibt man $A_\infty^0[\alpha]$- bzw. $S_\infty^0[\delta]$-stabil.

(4) Der Wert $r$ bei der $S_\infty^r[\delta]$- bzw. $A_\infty^r[\alpha]$-Stabilität heißt Widlund-Endradius. I.a. wird ein kleinstmöglicher Wert angegeben und von dem Widlund-Endradius gesprochen.

Hat ein Verfahren das Matrixpolynom $Q(\mu,H)$, so vererbt sich der entsprechende Name auf das Verfahren. Genauso hätte man auch die $A[\alpha]$- bzw. $S[\delta]$-Stabilität einführen können. Die Werte $\alpha$ und $\delta$ bei den Stabilitätsdefinitionen werden in der Schreibweise unterdrückt, wenn die Wert $\alpha=90^\circ$, bzw. $\delta=0$ betragen, genau w.o. auch.

Die Definition von $A_\infty^0$-Stabilität verallgemeinert offensichtlich den Begriff der L-Stabilität für Einschrittverfahren.

Offensichtlich ist jedes $A_\infty^{r_1}[\alpha]$-stabiles Polynom auch $A_\infty^{r_2}[\alpha]$-stabil, falls $r_1<r_2$, nicht jedoch immer umgekehrt. Gleiches gilt für die entsprechend andere Stabilitätseigenschaft, also jedes $S_\infty^{r_1}[\delta]$-stabile Polynom ist natürlich auch $S_\infty^{r_2}[\delta]$-stabil, falls $r_1<r_2$. Insbesondere ist jedes $A_\infty^0[\alpha]$-stabile Polynom auch $A_\infty[\alpha]$- und $A[\alpha]$-stabil.

I.d.R. wird man ein kleinstmögliches $r$ bei größtmöglichen $\alpha$, bzw. ein kleinstmögliches $r$ bei kleinstmöglichem $\delta$, anstreben.

Welche der Eigenschaften insgesamt am günstigsten ist, kann allgemein nicht gesagt werden. Man hat hier einen Zielkonflikt. Dies hängt auch von der zu lösenden Differentialgleichung ab. Weiß man im voraus, daß beispielsweise die zu lösende lineare Differentialgleichung $\dot y=Ay$ eine Jacobimatrix $A$ mit nur negativen reellen Eigenwerten hat, so ist ein großer Widlund-$\alpha$-Winkel unerheblich. Aber selbst einer beliebigen festen Matrix anzusehen, ob sie nur reelle Eigenwerte hat, ist keine triviale Aufgabe.

2. Beispiele: (1) Die BDF der Ordnung $p=1$ bis $p=6$, die Formeln von Peter Tischer aller Ordnungen und auch die Formeln von Tendler aller Ordnungen sind $A_\infty^0[\alpha]$- und $S_\infty^0[\delta]$-stabile Verfahren. Aber auch die blockimpliziten Verfahren von Bickart/Picel (1973) sind $A_\infty^0[\alpha]$- und $S_\infty^0[\delta]$-stabile Verfahren.

(2) Liniger/Gagnebin (1974) zeigten u.a., daß

$$ (40-\xi)y_{n+3}-(60-\xi)y_{n+2}+(24-3\xi)y_{n+1}-(4-\xi)y_n=24h\dot y_{n+3} $$

$A_\infty^0=S_\infty^0$-stabil ist für $0\le\xi\le16$. Der lokale Fehler ist hierbei

$$ -(4+\xi)h^3y^{(3)}(t) + {\cal O}(h^4). $$

Für $\xi=-4$ erhält man die BDF3, welche $A_\infty^0[86.03^\circ]$- und $S_\infty^0[0.083]$-stabil ist.

Werner Liniger (1927—2017) und Gagnebin, Thierry: "Construction of a Family of Second Order, $A$-Stable $k$-Step Formulas Depending on the Maximum Number, $2k-2$, of Parameters", in "Stiff Differential Systems", Editor Ralph A. Willoughby, Plenum Press, New York and London, 1974, pp.217—227

(3) Für ein nicht-triviales Beispiel, wo diese Eigenschaften erfüllt, aber auch möglicherweise nicht erfüllt sind, kann man wie folgt vorgehen. Kombiniert man in zyklischer Reihenfolge $i$-mal das explizite und $j$-mal das implizite Eulerverfahren, so erhält man sofort für die Nullstelle $\mu$ von $\det Q(\mu,H)$ die Gleichung

$$ \mu = {(1+H)^i\over (1-H)^j}, $$

der man augenblicklich ansieht, daß für betragsmässig große $H$ mit negativem Realteil, die Nullstelle $\mu$ sich wie untenstehend verhält:

$$ \mathopen|\mu\mathclose| \stackrel{\mathop{\rm Re}\nolimits H\to-\infty}{\longrightarrow} %\mapright^{\mathop{\rm Re}\nolimits H\to-\infty} \cases{\infty,& falls $i\gt j$;\cr 1, & falls $i=j$;\cr 0,& falls $i\lt j$.\cr} $$

Interessant hierbei ist, daß eine Kombination impliziter Verfahren mit expliziten Verfahren dennoch $A_\infty^0[\alpha]$-stabil, bzw. $S_\infty^0[\delta]$-stabil sein kann. Bei partiellen Differentialgleichungen, wo die Implizitheit von Formeln häufig noch schwerlastiger auf die gesamte Rechenzeit einwirken, werden sogar tatsächlich Kombinationen von expliziten und impliziten Verfahren verwendet, beispielsweise beim leapfrog-DuFort/Frankel-Verfahren, oder beim odd-even-hopscotch-Verfahren.

Schwächt man noch weiter ab, daß also nicht mehr konstante Schrittweiten zu nehmen sind, sondern auch variable Schrittweiten, so kann man sogar alleine mit expliziten Verfahren prinzipiell auskommen. Der Hinweis auf das Prinzipielle deswegen, weil die resultierenden Verfahren nicht notwendig auch im praktischen Gebrauch gut einsetzbar sind.

Beim expliziten Eulerverfahren liegen die Verhältnisse besonders durchsichtig vor. Hier ist

$$ \mu = (1+H_1)(1+H_2)\cdots(1+H_\nu) $$

und offensichtlicherweise kann man durch ein geeignet gewähltes $H_i$, mit $1\le i\le\nu$ das obige Produkt zum Verschwinden bringen. Bei vektorwertigen Differentialgleichungen liegen die Verhältnisse nicht mehr so einfach. John (Jack) D. Lambert in Lambert (1986) äußert sich hier in größerem Umfange zu dieser Problematik. Hans Jörg Stetter reißt dieses Thema kurz in seinem Buch "Analysis of Discretization Methods for Ordinary Differential Equations" von 1973 an, in dem Abschnitt über die Stabilitätsbereiche zyklischer Verfahren.

3. Weiteren Aufschluß über die Stabiltät eines linearen Verfahrens erhält man, wenn man versucht zu quantifizieren, wie schnell die vom Verfahren gelieferten Näherungswerte abklingen. Für genügend rasch abfallende Lösungen ist die folgende Stabilitätsmenge von Interesse:

$$ {\cal Y}_r := \left\{{\mskip 3mu} H\in\mathbb{C}: \det Q(\mu,H)=0 \: \land\: \mathopen|\mu\mathclose|\lt r {\mskip 3mu}\right\}. $$

Dies sind also diejenigen Zahlen $H$, sodaß $\det Q(\mu,H)$ verschwindet, dabei $\mu$ aber stets aus der offenen Kreisscheibe $D_r$ ist. Weiter werden jetzt die Mengen ${\cal A}[\alpha]$ und ${\cal S}[\delta]$ verallgemeinert zu Mengen, die den Radius $r$ mit berücksichtigen. Hier ist insbesondere natürlich $0<r<1$ von Interesse. Es seien also die Mengen betrachtet

$$ {\cal A}_r[\alpha] := \left\{{\mskip 3mu} H\in\mathbb{C}^-: \left|\arg(-H)\right| \lt \alpha \le90^\circ \hbox{ und } H\in{\cal Y}_r {\mskip 3mu}\right\} $$

und

$$ {\cal S}_r[\delta] := \left\{{\mskip 3mu} H\in\mathbb{C}^-:\mathop{\rm Re}\nolimits H\lt \delta\le0 \hbox{ und } H\in{\cal Y}_r {\mskip 3mu}\right\}. $$

Zur Zeichnung des Randes $\partial{\cal Y}_r$ löst man genügend oft das verallgemeinerte Eigenwertproblem in $H$ zu

$$ (A_1\mu+A_0){\bf x} = H(B_1\mu+B_0){\bf x}, \qquad \mu = r e^{i\varphi}. $$

Mit $\varphi$ durchläuft man dann eine diskrete Teilmenge von $[0^\circ,180^\circ]$, zum Beispiel $\varphi = 0^\circ, 20^\circ, \ldots, 160^\circ, 180^\circ$. Mit $r$ wandert man von $r=1$ beispielsweise in $(-1/10)$ Schritten nach $r=0.1$. Damit erhält man also die Höhenlinien der “Funktion” $H(r)$, bzw. von $H(\mu)$. Die “Funktion” $H(r)$ ist implizit definiert. Die Betragsfunktion $\left|\cdot\right|$ wähle dann unter den mehreren möglichen “Funktionswerten”, den betragsmässig größten aus, ordnet also letztlich jedem $re^{i\varphi}$ eine komplexe Zahl $H$ zu. Das bivariable Polynom $\det Q(\mu,H)$ ist natürlich algebraisch. Lösungsäste algebraischer Funktionen haben höchstens algebraische Singularitäten. $\partial{\cal Y}_r$ kann aus mehreren Zusammenhangskomponenten bestehen, insbesondere für kleiner werdende $r$ ziehen sich die “kreisförmigen” Komponenten immer mehr zu den Nullstellenpunkten von $\det Q(\mu,0)=0$ zusammen.

Ein mehr geometrisch-anschaulich motivierter Ansatz geht von der Betrachtung von $\mu(H)$ aus. Dann löst man das verallgemeinerte Eigenwertproblem in $\mu$ zu

$$ \left[ (A_1-HB_1)\mu + (A_0-HB_0) \right] {\bf x} = {\bf 0}, $$

auf einem Gitter von $\mathbb{C}$ für genügend viele $H$ und zeichnet entweder nur die betragsmässig größten $\mu$, oder zeichnet mehrere Schichten ein. Eine räumliche Veranschaulichung dieses Sachverhaltes findet man in Wanner (1987).

Wanner, Gerhard: "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

Für das verallgemeinerte Eigenwertproblem empfehlen sich die Unterprogramme cqzhes(), cqzval() und cqzvec() von Burton S. Garbow (1978), und Garbow (1984), welche eine Implementierung des QZ-Algorithmuses von Moler/Stewart (1973) darstellen. Aufgrund der Reellität der Matrixkoeffizienten, also

$$ \overline{\sum A_i}=\sum A_i, \hbox{ und wegen } \overline{e^{i\varphi}} = e^{-i\varphi}, $$

ist das Stabilitätsgebiet stets achsensymmetrisch zur reellen Achse $\mathbb{R}\subset\mathbb{C}$. Die Eigenvektoren mittels cqzvec() werden für obige Stabilitätsüberlegungen nicht benötigt.

Bibliographisch: Cleve Moler und G.W. "Pete" Stewart. Verfahrensbeschreibung: An Algorithm for the Generalized Matrix Eigenvalue Problem Ax = Lambda Bx. Ferner Burton S. Garbow.

4. Bei dem, wie schon oben angegebenen charakteristischen Polynom

$$ \det Q(\mu,H) = \det\left\{\sum_{i=0}^\kappa (A_i - H B_i) \mu^i\right\}, $$

erhält man für immer mehr dominierendes $H$ schließlich im Grenzfalle das Polynom $\det Q(\mu,\infty)$ zu

$$ \det Q(\mu,\infty) = \det\left\{ \sum_{i=0}^\kappa B_i\mu^i \right\}, $$

in etwas nachlässiger, aber sehr suggestiver Schreibweise. Es ist hier interessant zu vermerken, daß jetzt $A_{\infty}[\alpha]$- und $S_{\infty}[\delta]$-Stabilität zueinander äquivalent sind, für geignete $\alpha$ und $\delta$. Dies ist ein weitergehendes Ergebnis, als man es von der $A[\alpha]$- und $S[\delta]$-Stabilität her kannte.

5. Um nun zu überprüfen ob ein Verfahren sogar $A_{\infty}[\alpha]$-stabil, bzw. $S_{\infty}[\delta]$-stabil ist, erweitert man die 1.te Bedingung im Test A zu einer neuen Bedingung (1) und erhält damit den

6. Test B: Ein Verfahren ist $A_\infty^r[\alpha]$-stabil, bzw. $S_\infty^r[\delta]$-stabil genau dann, wenn:

  1. die Nullstellen von $\det Q(\mu,\infty)$ liegen in der offenen Einheitskreisscheibe und
  2. die Kurve der Punkte $H$, die $\det Q(\mu,H)|_{\mathopen|\mu\mathclose|=r} = 0$ erfüllen, im Komplement von ${\cal A}[\alpha]$, bzw. im Komplement von ${\cal S}[\delta]$ liegen.

7. Eine Möglichkeit die Bedingung (1) bei Test B zu erfüllen ist, wie es u.a. Gear (1968) vorgeschlagen hat, daß man alle Nullstellen $\mu$ von $\det Q(\mu,\infty)=0$ zu Null setzt. Dies ist z.B. der Fall, wenn man fordert:

$$ \beta_{ij} = 0\quad (j\le0), \qquad \beta_{ii}\ne0. \tag {NULL} $$

Eine einfache, notwendige und hinreichende Bedingung wird später angegeben. Von der Sicht des Speicherplatzes, des Rechenbedarfs und der programmiermässigen Einfachheit, sind diese Nullen besonders vorteilhaft einsetzbar. Diese Vorteile werden in dem Programm TENDLER genutzt. Das Programm ODIOUS kann diese Vorteile natürlich nicht nutzen, aufgrund $B_i\ne0$, für $i=0,\ldots,\kappa$.

Zusammen mit der Bedingung (LMSV) ergibt nun die Bedingung (NULL), daß

$$ \det Q(\mu,\infty) = \pm\mu^{\kappa\ell} \prod_{i=1}^{\ell} \beta_{ii}. $$

Nun enthält jede Stufe des Zykluses der Ordnung $p=k$, die der Bedingung (LMSV) und (Null) genügt, noch $\min\{i,k+1\}$ Parameter. Jede Stufe kann aber noch skaliert werden, ohne die Stabilitätseigenschaften zu verändern, d.h. von jeder Stufe muß man noch einen freien Parameter abziehen. Es bleiben dann $\min\{i-1,k\}$ völlig frei bestimmbare Parameter für die $i$-te Stufe, also ist

$$ \eqalignno{ {\ell(\ell-1)\over2}, &\qquad\hbox{falls}\quad\ell\le k+1\cr {k\left[ 2\ell - (k+1)\right]\over2}, &\qquad\hbox{falls}\quad\ell\ge k+1.\cr} $$

Wählt man die $\beta_{ij}$ für $0<j<i$, alle zu Null, so erhält man wieder die BDF. In diesem Sinne kann man die zyklischen Formeln von J.M. Tendler als echte Verallgemeinerung der einstufigen BDF auffassen. In der Tat unterscheiden sich die Formeln auch nicht allzu sehr in ihren Eigenschaften und in ihren Ergebnissen, die sie liefern, insbesondere ist bei den zyklischen Verfahren von Tendler die erste Stufe auch immer gleich der entsprechenden BDF.

8. Die sämtlichen bisherigen Stabilitätsdefinitionen gehören zur “Begriffswelt” der absoluten Stabilität, nicht jedoch zur relativen Stabilität. Für den Begriff der relativen Stabilität benötigt man die folgenden Überlegungen. Sei

$$ L(t,y,h) = \sum_{\nu=0}^m \sum_{i=0}^k \alpha_{i\nu} h^\nu y^{(\nu)}(t+ih). $$

Dann ist

$$ \det Q(\mu,H) = \sum_{\scriptstyle{0\le\nu\le m}\atop\scriptstyle{0\le i\le k}} \alpha_{i\nu} H^\nu \mu^i = \sum_{0\le i\le k} a_i(H)\mu^i = a_k(H) (\mu-\xi_1(H))\cdot\ldots\cdot(\mu-\xi_k(H)). $$

Sei $\mu=1$ einfache Nullstelle von

$$ \rho_0 = \sum_{0\le i\le k} \alpha_{i0} \mu^i. $$

Da dann $\mu=1$ separiert ist, folgt aufgrund der stetigen Abhängigkeit der Wurzeln für genügend kleine $h$

$$ \mathopen|e^H - \xi_\nu(H)\mathclose| \ge C, \qquad \nu=2,\ldots,k, $$

das ist die Fortpflanzung der Separiertheit für genügend kleine $h$. Jetzt muß aber

$$ \xi_1(H) = e^H + {\cal O}(H^{p+1}) = e^H + {\cal O}(h^{p+1}) $$

gelten, damit $Q(\mu,H)={\cal O}(h^{p+1})$ für $\mu=e^H$, aufgrund der Konsistenz von $L$.

9. Definition: Ein Matrixpolynom

$$ Q(\mu,H) = \sum_{\scriptstyle{0\le\nu\le m}\atop\scriptstyle{1\le i\le k}} A_{i\nu} H^\nu \mu^i $$

heißt konsistent mit mindestens der Konsistenzordnung $p$ genau dann, wenn die formale Ersetzung der Skalare $\mu^i$ ($i=0,\ldots,k$) durch die Vektoren

$$ \pmatrix{1\cr \vdots\cr \mu^{\ell-1}\cr},\quad\ldots\quad \pmatrix{\mu^{i\ell}\cr \vdots\cr \mu^{i\ell+\ell-1}\cr},\quad\ldots\quad \pmatrix{\mu^{k\ell}\cr \vdots\cr \mu^{k\ell+\ell-1}\cr} $$

($\mu=e^H$) dann liefert

$$ Q(\mu,H) = \pmatrix{{\cal O}(H^{p_1+1})\cr \vdots\cr {\cal O}(H^{p_\ell+1})\cr}, \qquad (p\ge p_i). $$

Das Matrixpolynom $Q(\mu,H)$ ist konsistent mit der genauen Ordnung $\min_{i=1}^\ell p_i$. Gilt $p_i\ge1$ ($\forall i$), so heißt das Matrixpolynom konsistent. Möchte man die besondere Form

$$ [\mu^{i\ell},\ldots,\mu^{i\ell+\ell-1}], \qquad (\mu=e^H) $$

besonders betonen, so spricht man von polynomial-konsistent.

Die obige Definition verallgemeinert $\rho(e^h)-h\sigma(e^h)={\cal O}(h^{p+1})$ auf $\ell$-Mehrstufigkeit mit $m$-facher Ableitung für äquidistante Gitter. Die Konsistenz nach obiger Festlegung ist eine Eigenschaft des Matrixpolynomes. Sie vererbt sich sinngemäß auf das Verfahren, welches dem Matrixpolynom i.d.R. zugrunde liegt.

Zur besseren Unterscheidung zwischen absoluter und relativer Stabilität seien beide einander gegenübergestellt.

10. Definition: (10.1) Ein Matrixpolynom $Q\mu,H)$ heißt absolut-stabil für ein $H_0$ genau dann, wenn für die sämtlichen Nullstellen $\mu$ von $\det Q(\mu,H_0)=0$ gilt: $\mathopen|\mu\mathclose|<1$.

(10.2) Die Menge ${\cal Y}_1$ derjenigen $H$, für die das Matrixpolynom $Q(\mu,H)$ absolut-stabil ist, also

$$ {\cal Y}_1 = \{H\in\mathbb{C}: \det Q(\mu,H)=0\: \land\: \mathopen|\mu\mathclose|\lt 1\} $$

heißt Menge der absoluten Stabilität.

Offensichtlich gilt: Enthält ${\cal Y}_1$ die offene, linke komplexe Halbebene $\mathbb{C}^-$, so ist das Verfahren $A$-stabil.

11. Definition: (11.1) Ein konsistentes Matrixpolynom $Q(\mu,H)$, für welches $\mu=1$ einfache Nullstelle von $\det Q(\mu,0)=0$ ist, heißt relativ-stabil für ein $H_0$ genau dann, wenn gilt

$$ \mathopen|\mu_i(H)\mathclose| \lt \mathopen|\mu_1(H)\mathclose| \quad(i\ge2), \qquad \mu_1(0)=1. $$

$\mu_1(H)$ ist also die zu $\mu_1(0)=1$ gehörende Nullstelle, wobei $Q(\mu_i(H),H)=0$ ($i\ge1$). Hat $\det Q(\mu,H)$ in $\mu$ nur den Grad 1, so sei das Matrixpolynom $Q(\mu,H)$ relativ-stabil.

(11.2) Die Menge $\cal Z$ derjenigen $H$, für die das Matrixpolynom $Q(\mu,H)$ relativ-stabil ist, also

$$ {\cal Z} = \{H\in\mathbb{C}: \det Q(\mu,H)=0\: \land\: \mathopen|\mu_i\mathclose|\lt \mathopen|\mu_1\mathclose|\: (i\gt 2)\} $$

heißt Menge der relativen Stabilität.

Im Falle eines linearen Differentialgleichungssystems $\dot y=Ay\in\mathbb{R}^n$ mit diagonalähnlicher Matrix $A$ entstehen bei Transformation $n$ skalare Differentialgleichungen. Eigenwerte von $A$ mit stark negativen Realteil, also $\mathop{\rm Re}\nolimits \lambda\ll0$, führen zu schnell abklingenden Lösungen. In diesem Falle zu forden, daß die Nebenwurzeln $\mu_i(H)$ die Hauptwurzel $\mu_1(H)$ nicht betragsmässig übersteigen ist nur wichtig, wenn man an sehr hoher relativer Genauigkeit interessiert ist, sonst allerdings nicht. Wesentlicher ist in diesem Falle die absolute Stabilität. Man ist an Exponentialtermen nur so lange interessiert, wie sie einen im Rahmen der Rechengenauigkeit signifikanten Einfluß auf das Gesamtergebnis besitzen. Spielen stark abklingende Terme keine Rolle mehr, so genügt absolute Stabilität. Dies motiviert die

12. Definition: (12.1) Ein konsistentes Matrixpolynom $Q\mu,H)$ heißt $LR$-$S[\delta]$-stabil genau dann, wenn es $S[\delta]$-stabil ist und das um den Ursprung angeordnete Rechteck

$$ {\cal R} = \{H\in\mathbb{C}: a_1\le\mathop{\rm Re}\nolimits H\le a_2\: \land\: b_1\le\mathop{\rm Im} H\le b_2\} $$

mit $a_1<0<a_2$, $b_1<0<b_2$ Teilmenge der Menge der relativen Stabilität ist.

(12.2) Analog definiert man $LR$-$S_\infty^r[\delta]$- und $LR$-$A_\infty^r[\alpha]$-Stabilität.

$LR$ soll an lokal-relativ erinnern. Die Festlegung auf ein Rechteck ist nicht wesentlich. Genauso gut hätte man eine Kreisscheibe, eine elliptische Scheibe, ein Hexagon, etc. nehmen können, solange es um den Ursprung angeordnet ist. Je umfassender $\cal R$ ist, desto günstiger. Die Begrenzung in $b_1$ und $b_2$ rührt daher, daß Exponentialterme mit stark imaginären Anteil, aber geringen negativen Realteil als Exponent, noch mehr Eigenschaften benötigen als relative Stabilität. Für Differentialgleichungen mit periodischen Lösungen sind gesonderte Formeln besser geeignet.

Abweichend von der Definition in Gear (1968) und auch abweichend von der Definition in Werner/Arndt (1986), aber völlig sinnentsprechend wie in Tendler/Bickart/Picel (1978), in Bickart/Picel (1973) und in Albrecht (1979), sei wie folgt definiert.

13. Definition: Ein konsistentes Matrixpolynom $Q(\mu,H)$, welches einer oder mehrerer der Bedingungen der $A[\alpha]$-, $S[\delta]$-, oder $LR$-$A[\alpha]$-, $LR$-$S[\delta]$-Stabilität genügt, heißt steif-stabil.

Bibliographisch: C. William Gear (1935—2022), Theodore A. Bickart (1936—2023), obituary, Werner, Helmut (1931--1985), Arndt, Herbert, Gear, Charles William (1935--2022), Peter Albrecht.

Bickart, Theodore A. und Picel, Zdenek: "High Order Stiffly Stable Composite Multistep Methods for Numerical Integration of Stiff Differential Equations", BIT, Vol 13, 1973, pp.272—286

5. Konstruktive Herleitung steif-stabiler zyklischer Formeln

A computer program for plotting both the Lambda and Zeta Loci was presented and described. Using an interactive version of this program new and efficient cyclic composite linear multistep methods of orders three through seven were obtained. $\ldots$
The test for stiff stability of a composite linear multistep method presented herein is entirely graphical. $\ldots$
Interactive communication was via a teletypewriter, the only device available for this purpose at the computer facility used.

J.M. Tendler (1973)

1. In diesem Abschnitt wird darauf eingegangen, wie die in dem Programm TENDLER verwendeten Formeln, von Tendler (1973) entwickelt wurden. Neben der schon mehrfach erwähnten Dissertation von J.M. Tendler, findet man eine Darstellung der Parametrisierung von Verfahren in dem Aufsatz von Donelson/Hansen (1971). U.a. weitere Hinweise zur Konstruktion generell zyklischer Verfahren und auch steif-stabiler zyklischer Verfahren, einschließlich entsprechender Parametrisierungen, findet man in dem Buch von Albrecht (1979).

Die beiden Aufsätze von Mihelčić (1977), Mihelčić (1978), die Veröffentlichung von Mihelčić/Wingerath (1981), und die Disserattion von Tischer (1983), die Arbeit Tischer/Sacks-Davis (1983)](https://doi.org/10.1137/0904051) bzw. Tischer/Gupta (1985), liefern ebenfalls Aufschluß über die Gewinnung zyklischer, steif-stabiler Verfahren. Tischer/Gupta (1985) äußern sich dort mehr zu den Erfahrungen mit den von Tischer (1983)2 und Tischer/Sacks-Davis (1983)3 gewonnen Formeln, wobei auch weitere zyklische Formeln in die Bewertung mit eingeschlossen werden.

Die Dissertation von Tischer (1983) und die Arbeit von Tischer/Sacks-Davis (1983) geben sogar zweistufige $A_\infty^0$-stabile zyklische Verfahren an bis zur Ordnung $p\le4$. Über der Ordnung $p>4$ sind die weiterhin zweistufigen Zyklen noch $A_\infty^0[\alpha]$-stabil, mit $\alpha$=$86^\circ$, $76^\circ$, $57^\circ$, und $22^\circ$. Die Formeln erfüllen sogar noch zahlreiche weitere wünschenswerte Eigenschaften, allerdings haben sie auch gewisse Nachteile.

Bibliographisch. Matija (Miško) Mihelčić (1935—2011), Peter E. Tischer.

Matija Mihelčić: "Fast $A$-stabile Donelson-Hansensche zyklische Verfahren zur numerischen Integration von 'stiff'-Differentialgleichungssystemen", Angewandte Informatik, 1977, Vol 19, Heft 7, pp.299—305

Matija Mihelčić: "$A(\alpha)$-Stable Composite Multistep Methods of Order 5", Computing, Vol 20, 1978, pp.267—272

Matija Mihelčić und K. Wingerath: "$A(\alpha)$-stable Cyclic Composite Multistep Methods of Orders 6 and 7 for Numerical Integration of Stiff Ordinary Differential Equations", ZAMM, Band 61, 1981, pp.261—264

Tischer, Peter E.: "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 S.

Tischer, Peter E. und Sacks-Davis, Ron: "A New Class of Cyclic Multistep Formulae for Stiff Systems", SIAM Journal on Scientific and Statistical Computing, Vol 4, No 4, December 1983, pp.733—747.

2. Für die gestellte Aufgabe zur Gewinnung neuerer und besserer Verfahren für die Integration steifer Differentialgleichungen ging Tendler (1973) wie folgt vor: Man entwerfe zyklische Formeln der Konvergenzordnung 1 bis 6, die Ordnung für Ordnung bessere Stabilitätseigenschaften besitzen als die BDF. Falls möglich sollten sogar Verfahren mit einer Ordnung höher als 6 entwickelt werden. Als Vergleichsmaßstab wurde der Widlund-$\alpha$-Winkel gewählt. Andere Kriterien spielten keine Rolle.

Die Bedingungen (Null) und (LMSV) an die Koeffizienten $\alpha_{ij}$ und $\beta_{ij}$, sorgen dafür, daß der Zyklus stets die folgende Form hat

$$ \underbrace{ \begin{array}{cccccc} *&*&*&*& & \cr &*&*&*&*& \cr & &*&*&*&* \end{array} }_{ \displaystyle\alpha_{ij} } \qquad \quad \underbrace{ \begin{array}{cccccc} 0&0&0&*& & \cr &0&0&*&*& \cr & &0&*&*&* \end{array} }_{ \displaystyle\beta_{ij} } $$

Die BDF erster und zweiter Ordnung

$$ \eqalign{ y_{n+1} &= y_n + h\dot y_{n+1}, \cr 3y_{n+1} &= 4y_n - y_{n-1} + 2h\dot y_{n+1} } $$

sind beide $A_{\infty}^0[\alpha]$-stabil mit optimalen Widlund-Winkel von $\alpha=90^\circ$ und optimaler Widlund-Distanz $\delta=0$, sodaß damit eine weitere Verbesserung nicht mehr möglich ist, wenn man wie gesagt nur den Widlund-$\alpha$-Winkel als Vergleichsmaßstab nimmt. Diese beiden Verfahren werden daher auch für $p=1$ und entsprechend für $p=2$ in dem Programm TENDLER verwendet.

Die $i$ Parameter in der $i$-ten Stufe deuten darauf hin, daß $i$ linear unabhängige Formeln in geeigneter Kombination ausreichen, um ein neues Verfahren zu konstruieren. Für die Ordnungen $p=3$ und höher, wurde daher wie folgt die Suche nach besseren Formeln durchgeführt:

  1. Im ersten Schritt wurden zyklische Formeln gewonnen, die die Bedingung (LMSV) und (Null) erfüllen und
  2. im zweiten Schritt wurde der Widlund-$\alpha$-Winkel numerisch berechnet und mit dem entsprechenden Winkel für die BDF verglichen.
  3. War dann der Winkel größer und damit besser, wurde die Suche abgebrochen, andernfalls wurden die Parameter weiter variiert und der erste Schritt, ggf. mit einer anderen Stufenzahl $\ell$, erneut durchgeführt.

Beim ersten Schritt ging man wie folgt vor. Zuerst wurde eine Folge von $k+1$ linear unabhängigen Mehrschrittformeln aufgestellt, sodaß die $\beta$-Koeffizienten eine Diagonalmatrix bilden. Wegen der linearen Unabhängigkeit ist dies immer zu erreichen. Diese Diagonalform erleichtert das Generieren von Zyklen mit der Eigenschaft

$$ \beta_{ij}=0 \quad (j\le0), \qquad \beta_{ii}\ne0. $$

Da nur Linearkombinationen von diesen Mehrschrittformeln gebildet werden, ist das Erfülltsein von (LMSV) ebenfalls sofort klar. Dann wurde für jedes $i=1,\ldots,\ell$ eine Linearkombination dieser Formeln berechnet um daraus die $i$-te Stufe des Zykluses zu bilden. Die Stufenzahl $\ell$ wurde versucht so klein wie möglich zu halten.

Um das Vorgehen zu illustrieren ein längeres Beispiel.

3. Beispiel: Es seien die $k+1=4$ linear unabhängigen Mehrschrittformeln $F1$ bis $F4$ der Konsistenzordnung 3 betrachtet:

$$ \matrix{F1:\cr F2:\cr F3:\cr F4:\cr}\qquad \pmatrix{ -2 & 9 & -18 & 11 & 0 & 0 & 0 & -6\cr 1 & -6 & 3 & 2 & 0 & 0 & -6 & 0\cr -2 & -3 & 6 & -1 & 0 & -6 & 0 & 0\cr -11 & 18 & -9 & 2 & -6 & 0 & 0 & 0 } \pmatrix{ y_{n+i-3}\cr y_{n+i-2}\cr y_{n+i-1}\cr y_{n+i}\cr h\dot y_{n+i-3}\cr h\dot y_{n+i-2}\cr h\dot y_{n+i-1}\cr h\dot y_{n+i} } = \pmatrix{ 0\cr 0\cr 0\cr 0 }. $$

Jede Zeile der Matrix ist dabei das entsprechende Mehrschrittverfahren, also die erste Zeile ist F1, die zweite ist F2 und so fort. Das Verfahren $F1$ ist implizit, während hingegen $F2$, $F3$ und $F4$ explizite Verfahren sind. Wie so häufig gilt auch hier, daß die einzelnen Stufen instabil sind, mit Ausnahme der ersten, welche gerade die BDF3 ist, mit den Nullstellen 1 und $0.3\overline{18}\pm 0.284i$. Die Nullstellen der zweiten “Zeile” liegen bei 1, $-2.686$ und $0.186$. Die dritte “Zeile” hat die Nullstellen 1, $5.372$ und $-0.372$ und schließlich die letzte vierte Zeile hat die Nullstellen 1 und $1.75\pm 1.561i$. Diese Instabilität der Stufen bzw. von Komponenten der Stufen, überträgt sich nicht automatisch auf den Gesamtzyklus. Die Stabilität des Gesamtzyklus wird getrennt überprüft.

Von diesem Arsenal an Formeln wurde nun der Zyklus für das Verfahren dritter Ordnung, wie folgt aufgebaut: Jede Stufe wurde als Linearkombination der $F1$ bis $F4$ erhalten, d.h. die $i$-te Stufe ist gleich

$$ c_{i1}F1 + c_{i2}F2 + c_{i3}F3 + c_{i4}F4. $$

Um nun die Bedingung (NULL) zu erfüllen, muß offensichtlich $c_{ij}=0$ für $j>i$ gelten. Der Widlund-$\alpha$-Winkel wurde jetzt für die Werte maximiert:

$$ %\newcommand{\md}[2]{\mathchoice{\textstyle\nicefrac{#1}{#2}}{\nicefrac{#1}{#2}}{\nicefrac{#1}{#2}}{\nicefrac{#1}{#2}}} \def\md#1{#1} \begin{array}{lllll} \ell=3: & i=1: & c_{11}=1, & & \cr & i=2: & c_{21}=1, & c_{22}=0, & \cr & i=3: & c_{31}=\md{1/3}, & c_{32}=-\md{2/3}, & c_{33}=-\md{2/3}.\cr \end{array} $$

Der so gebildete Zyklus der Ordnung $p=3$ kann daher wie folgt geschrieben werden:

$$ \begin{array}{ccccccccccccccccc} -2y_{m\ell-2}&+&9y_{m\ell-1}&-&18y_{m\ell}&+&11y_{m\ell+1}&&&&% &=& 6h\dot y_{m\ell+1},\cr &-&2y_{m\ell-1}&+&9y_{m\ell}&-&18y_{m\ell+1}&+&11y_{m\ell+2}&&% &=&&& 6h\dot y_{m\ell+2},\cr &&&&&&9y_{m\ell+1}&-&12y_{m\ell+2}&+&3y_{m\ell+3}% &=& h\bigl(-4\dot y_{m\ell+1}&-&4\dot y_{m\ell+2}+2\dot y_{m\ell+3}\bigr).\cr \end{array} $$

Weil der Widlund-$\alpha$-Winkel grösser und damit besser als der entsprechende Winkel für BDF3 ausgefallen war, wurde der Zyklus akzeptiert. Die Formeln der Konvergenzordnung $p=4$ bis $p=7$ wurden auf genau die oben beschriebene Art und Weise erhalten. Daß bei diesem Beispiel die ersten beiden Stufen gleich sind, nämlich hier gleich BDF3, ist mehr oder minder Zufall. Mit der Bedingung $c_{ij}=0$ für $i>j$ hat dies nichts zu tun. Beim zyklischen Verfahren sechster Ordnung sind auch tatsächlich alle Stufen untereinander verschieden. Bei allen anderen Verfahren hingegen, außer eben bei $p=6$, sind die ersten beiden Stufen gleich und zwar gleich der entsprechenden BDF.

4. Es waren 3 Stufen nötig für die Verfahren 3.ter und 4.ter Ordnung und 4 Stufen für die Formeln 5.ter, 6.ter und 7.ter Ordnung. Die genauen Zahlenwerte für die Koeffizienten aller Zyklen werden später in einem Schema angegeben. Es sei angemerkt, daß immer die ersten beiden Stufen des $i$-ten Zykluses mit BDF$i$ beginnen. Dies gilt bei allen Verfahren, außer bei der Ordnung $p=6$. Dort ist es nur die 1.te Stufe.

5. Zum Vergleich seien einmal die Werte für den Widlund-$\alpha$-Winkel für die neuen zyklischen Mehrschrittverfahren, die in dem Programm TENDLER verwendet werden, mit denen der BDF in der untenstehenden Tabelle gegenübergestellt. Auch findet man die entsprechenden Werte für die Formeln, die in der Dissertation Tischer (1983), bzw. Tischer/Sacks-Davis (1983) entwickelt haben, unter der Spalte ODIOUS. BP73 bezeichnet hier die blockimpliziten Verfahren von Bickart/Picel (1973). Die BDF$i$, für alle $i\ge7$, sind nicht $D$-stabil, also natürlich erst recht nicht $A[\alpha]$- bzw. $S[\delta]$-stabil. Einen Beweis und alternative Beweisansätze findet man in dem Buch von Hairer/Wanner/Nørsett (1987) oder hier Hairer/Wanner (1983): On the Instability of the BDF Formulas.

Zyklische Verfahren über der Ordnung $p=7$ werden in dem Programm TENDLER nicht verwendet. Die Programme LSODE, LSODA und LSODAR benutzen die BDF lediglich bis zur Ordnung $p=5$. Rechts daneben ist auch eine Gegenüberstellung der $\delta$-Werte bei der $S_{\infty}[\delta]$-Stabilität.

Widlund-$\alpha$-Winkel

p TENDLER BDF ODIOUS BP73
1 90.00 90.00 90.00 90.0
2 90.00 90.00 90.00 90.0
3 89.43 86.03 90.00 88.9
4 80.88 73.35 90.00 87.7
5 77.48 51.84 86.64 85.5
6 63.25 17.84 76.32 82.7
7 33.53 * 57.66 79.5
8 * * 22.15 76.0
9 * * * 72.5
10 * * * 69.6

Widlund-Distanz $\delta$

p TENDLER BDF
1 0 0
2 0 0
3 0.0048 0.083
4 0.24 0.67
5 1.4 2.3
6 2.9 6.1
7 10.2 *

Die Widlund-Winkel der BDF können exakt bestimmt werden: Gander/Wanner (2020): Exact BDF stability angles with maple und Akrivis/Katsoprinakis (2019): Maximum angles of $A(\vartheta)$-stability of backward difference formulae, PDF.

6. Die besseren Stabilitätseigenschaften, die sich durch den größeren Widlund-$\alpha$-Winkel bemerkbar machen, erlauben nun einem Programm, welches die hergeleiteten Formeln benutzt, früher die Schrittweite zu erhöhen, da sie nicht durch Stabilitätsbeschränkungen verkleinert werden muß. Insbesondere für Differentialgleichungen, bei denen die Eigenwerte der Jacobimatrix nahe der imaginären Achse sind, sollten besser integriert werden können.

Das Suchverfahren nach steif-stabilen zyklischen Formeln wird beschrieben in der Dissertation von Tendler (1973). Dort wird auch ein Programm zur Berechnung des Widlund-$\alpha$-Winkels angegeben, basierend auf einem Polynomlöser. Das hier beschriebene Suchverfahren hängt von einem Programm zur Berechnung des Widlund-Winkels ab.

Man beachte, daß alle Verfahren interaktiv gefunden wurden: “Probieren-gucken”, “probieren-gucken”, $\ldots$
Wenn also hier von “maximiert” gesprochen wird, ist damit gemeint, daß durch Suchen, nicht jedoch kalkülmässig, in geschlossener Form, eine Maximierung stattgefunden hat. Sieht man einmal von grundsätzlichen Einstellungen ab, so kann man diesem pragmatischen Ansatz seine Bedeutung nicht absprechen.

Es wurde auch versucht ein Verfahren der Konvergenzordnung $p=8$ zu erhalten, welches steif-stabil ist, doch gelang dies nicht. Dies führt Tendler (1973) darauf zurück, daß eben nur interaktiv gesucht wurde. Eine andere Möglichkeit wird darin gesehen, daß die Begrenzung der Anzahl der Stufen auf $\ell=4$ zu restriktiv gewesen sei.

Testläufe haben ergeben, daß die Auffindung von zyklischen Verfahren mit mehr als 4 Stufen durchaus sinnvoll sein kann. Oder m.a.W., die Erfahrungen mit dem Programm TENDLER deuten darauf hin, daß Zyklen mit mehr als 4 Stufen nicht krass den praktischen Bedürfnissen widerstreben. Ob diese Erfahrung grundsätzlich für geeignete vielstufige Zyklen sich bestätigen könnte, kann selbstverständlich nicht pauschal postuliert werden. Ein Einbau solcher Formeln in das Programm TENDLER zumindestens wäre völlig unproblematisch. Da jede Stufe Element einer ringzyklischen Liste ist, macht es keinen Aufwand hier eine weiteres Ringelement einzufügen. Zudem ist die Konstante 4 als Programmkonstante MAXYRING abgelegt.

Außer Konstantenänderungen sind keine Änderungen nötig, insbesondere keine algorithmischen, vorausgesetzt es handelt sich um echte Verallgemeinerungen der BDF, wie es beispielsweise die Tendlerschen Zyklen darstellen. Man beachte hierzu die Ausführungen bzgl. der Prädiktorformeln. Zu berücksichtigen ist jedoch, daß der Speicherplatzbedarf dann größer wird, allerdings nur in linearen Termen der Dimension $n$ der Differentialgleichung. Dies nur, falls die Ordnung 7 übersteigt. Dies liegt an der Notwendigkeit zur Bildung rückwärtsgenommener Differenzen bis zur Ordnung $p+1$ in dem Programm TENDLER.

6. Basisformeln der Tendlerschen Zyklen

1. Die Bestandteile der zyklischen Formeln von Tendler (1973) lauten:

$$ \begin{array}{c|c} p & i=1,\ldots,\ell\cr \hline j & \alpha_{ij}\cr \hline j & \beta_{ij}\cr \end{array} $$

Die in jeder Tabelle ganz rechts stehende Formel, ist immer die BDF. Die restlichen Formeln sind explizite Verfahren. Alle unten aufgeführten Formeln benötigen $k$ Startwerte und besitzen die Konsistenzordnung $p=k$. Die im folgenden angegebenen Basisformeln sind also ersteinmal nur für die Konstruktion von Zyklen der Konvergenzordnung $p=k$ geeignet. Wünscht man andere Konvergenzordnungen, bei anderer Schrittzahl, so gelten die Überlegungen jedoch entsprechend. Berücksichtigt man den Effekt der annullierten Dominanz, so sind die angegebenen Basisformeln auch für zyklische Verfahren der Konvergenzordnung $p=k+1$ geeignet. Unter den Verfahren befindet sich stets der entsprechende, gänzlich unskalierte Fehlerfaktor, der berechnet wurde nach der Vorschrift

$$ c_{p+1} := {1\over(p+1)!} C_{p+2,k} \pmatrix{\alpha\cr \beta\cr}. $$

Die Basisformeln im einzelnen.

$$ \def\x#1/#2.{\displaystyle{#1\over #2}} \def\y#1/#2.{\displaystyle{-{#1\over #2}}} \begin{array}{r|rrrr} p=3 & 1 & 2 & 3 & 4\cr \hline -2 & -11 & -2 & 1 & -2\cr -1 & 18 & -3 & -6 & 9\cr 0 & -9 & 6 & 3 & -18\cr 1 & 2 & -1 & 2 & 11\cr \hline -2 & 6 & 0 & 0 & 0\cr -1 & 0 & 6 & 0 & 0\cr 0 & 0 & 0 & 6 & 0\cr 1 & 0 & 0 & 0 & 6\cr \hline c_{i,4} & \x3/2. & \y1/2. & \x1/2. & \y3/2.\cr \end{array} $$
$$ \begin{array}{r|rrrrr} p=4 & 1 & 2 & 3 & 4 & 5\cr \hline -3 & -25 & -3 & 1 & -1 & 3\cr -2 & 48 & -10 & -8 & 6 & -16\cr -1 & -36 & 18 & 0 & -18 & 36\cr 0 & 16 & -6 & 8 & 10 & -48\cr 1 & -3 & 1 & -1 & 3 & 25\cr \hline -3 & 12 & 0 & 0 & 0 & 0\cr -2 & 0 & 12 & 0 & 0 & 0\cr 1 & 0 & 0 & 12 & 0 & 0\cr 0 & 0 & 0 & 0 & 12 & 0\cr 1 & 0 & 0 & 0 & 0 & 12\cr \hline c_{i,5} & \y12/5. & \x3/5. & \y2/5. & \x3/5. & \y12/5.\cr \end{array} $$
$$ \begin{array}{r|rrrrrr} p=5 & 1 & 2 & 3 & 4 & 5 & 6\cr \hline -4 & -137 & -12 & 3 & -2 & 3 & -12\cr -3 & 300 & -65 & -30 & 15 & -20 & 75\cr -2 & -300 & 120 & -20 & -60 & 60 & -200\cr -1 & 200 & -60 & 60 & 20 & -120 & 300\cr 0 & -75 & 20 & -15 & 30 & 65 & -300\cr 1 & 12 & -3 & 2 & -3 & 12 & 137\cr \hline -4 & 60 & 0 & 0 & 0 & 0 & 0\cr -3 & 0 & 60 & 0 & 0 & 0 & 0\cr -2 & 0 & 0 & 60 & 0 & 0 & 0\cr 1 & 0 & 0 & 0 & 60 & 0 & 0\cr 0 & 0 & 0 & 0 & 0 & 60 & 0\cr 1 & 0 & 0 & 0 & 0 & 0 & 60\cr \hline c_{i,6} & 10 & -2 & 1 & -1 & 2 & -10\cr \end{array} $$
$$ \begin{array}{r|rrrrrrr} p=6 & 1 & 2 & 3 & 4 & 5 & 6 & 7\cr \hline -5 & -147 & -10 & 2 & -1 & 1 & -2 & 10\cr -4 & 360 & -77 & -24 & 9 & -8 & 15 & -72\cr -3 & -450 & 150 & -35 & -45 & 30 & -50 & 225\cr -2 & 400 & -100 & 80 & 0 & -80 & 100 & -400\cr -1 & -225 & 50 & -30 & 45 & 35 & -150 & 450\cr 0 & 72 & -15 & 8 & -9 & 24 & 77 & -360\cr 1 & -10 & 2 & -1 & 1 & -2 & 10 & 147\cr \hline -5 & 60 & 0 & 0 & 0 & 0 & 0 & 0\cr -4 & 0 & 60 & 0 & 0 & 0 & 0 & 0\cr -3 & 0 & 0 & 60 & 0 & 0 & 0 & 0\cr -2 & 0 & 0 & 0 & 60 & 0 & 0 & 0\cr -1 & 0 & 0 & 0 & 0 & 60 & 0 & 0\cr 0 & 0 & 0 & 0 & 0 & 0 & 60 & 0\cr 1 & 0 & 0 & 0 & 0 & 0 & 0 & 60\cr \hline c_{i,7} & \y60/7. & \x10/7. & \y4/7. & \x3/7. & \y4/7. & \x10/7. & \y60/7.\cr \end{array} $$
$$ \begin{array}{r|rrrrrrrr} p=7 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8\cr \hline -6 & -1089 & -60 & 10 & -4 & 3 & -4 & 10 & -60\cr -5 & 2940 & -609 & -140 & 42 & -28 & 35 & -84 & 490\cr -4 & -4410 & 1260 & -329 & -252 & 126 & -140 & 315 & -1764\cr -3 & 4900 & -1050 & 700 & -105 & -420 & 350 & -700 & 3675\cr -2 & -3675 & 700 & -350 & 420 & 105 & -700 & 1050 & -4900\cr -1 & 1764 & -315 & 140 & -126 & 252 & 329 & -1260 & 4410\cr 0 & -490 & 84 & -35 & 28 & -42 & 140 & 609 & -2940\cr 1 & 60 & -10 & 4 & -3 & 4 & -10 & 60 & 1089\cr \hline -6 & 420 & 0 & 0 & 0 & 0 & 0 & 0 & 0\cr -5 & 0 & 420 & 0 & 0 & 0 & 0 & 0 & 0\cr -4 & 0 & 0 & 420 & 0 & 0 & 0 & 0 & 0\cr -3 & 0 & 0 & 0 & 420 & 0 & 0 & 0 & 0\cr -2 & 0 & 0 & 0 & 0 & 420 & 0 & 0 & 0\cr -1 & 0 & 0 & 0 & 0 & 0 & 420 & 0 & 0\cr 0 & 0 & 0 & 0 & 0 & 0 & 0 & 420 & 0\cr 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 420\cr \hline c_{i,8} & \x105/2. & \y15/2. & \x5/2. & \y3/2. & \x3/2. & \y5/2. & \x15/2. & \y105/2.\cr \end{array} $$
$$ \begin{array}{r|rrrrrrrrr} p=8 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9\cr \hline -7 & -2283 & -105 & 15 & -5 & 3 & -3 & 5 & -15 & 105\cr -6 & 6720 & -1338 & -240 & 60 & -32 & 30 & -48 & 140 & 960\cr -5 & -11760 & 2940 & -798 & -420 & 168 & -140 & 210 & -588 & 3920\cr -4 & 15680 & -2940 & 1680 & -378 & -672 & 420 & -560 & 1470 & -9408\cr -3 & -14700 & 2450 & -1050 & 1050 & 0 & -1050 & 1050 & -2450 & 14700\cr -2 & 9408 & -1470 & 560 & -420 & 672 & 378 & -1680 & 2940 & -15680\cr -1 & -3920 & 588 & -210 & 140 & -168 & 420 & 798 & -2940 & 11760\cr 0 & 960 & -140 & 48 & -30 & 32 & -60 & 240 & 1338 & -6720\cr 1 & -105 & 15 & -5 & 3 & -3 & 5 & -15 & 105 & 2283\cr \hline -7 & 840 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\cr -6 & 0 & 840 & 0 & 0 & 0 & 0 & 0 & 0 & 0\cr -5 & 0 & 0 & 840 & 0 & 0 & 0 & 0 & 0 & 0\cr -4 & 0 & 0 & 0 & 840 & 0 & 0 & 0 & 0 & 0\cr -3 & 0 & 0 & 0 & 0 & 840 & 0 & 0 & 0 & 0\cr -2 & 0 & 0 & 0 & 0 & 0 & 840 & 0 & 0 & 0\cr -1 & 0 & 0 & 0 & 0 & 0 & 0 & 840 & 0 & 0\cr 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 840 & 0\cr 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 840\cr \hline c_{i,9} & \y280/3. & \x35/3. & \y10/3. & \x5/3. & \y4/3. & \x4/3. & \y10/3. & \x35/3. & \y280/3.\cr \end{array} $$
$$ \begin{array}{r|rrrrrrrrrr} p=9 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10\cr \hline -8 & -7129 & -280 & 35 & -10 & 5 & -4 & 5 & -10 & 35 & -280\cr -7 & 22680 & -4329 & -630 & 135 & -60 & 45 & -54 & 105 & -360 & 2835\cr -6 & -45360 & 10080 & -2754 & -1080 & 360 & -240 & 270 & -504 & 1680 & -12960\cr -5 & 70560 & -11760 & 5880 & -1554 & -1680 & 840 & -840 & 1470 & -4704 & 35280\cr -4 & -79380 & 11760 & -4410 & 3780 & -504 & -2520 & 1890 & -2940 & 8820 & -63504\cr -3 & 63504 & -8820 & 2940 & -1890 & 2520 & 504 & -3780 & 4410 & -11760 & 79380\cr -2 & -35280 & 4704 & -1470 & 840 & -840 & 1680 & 1554 & -5880 & 11760 & -70560\cr -1 & 12960 & -1680 & 504 & -270 & 240 & -360 & 1080 & 2754 & -10080 & 45360\cr 0 & -2835 & 360 & -105 & 54 & -45 & 60 & -135 & 630 & 4329 & -22680\cr 1 & 280 & -35 & 10 & -5 & 4 & -5 & 10 & -35 & 280 & 7129\cr \hline -8 & 2520 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\cr -7 & 0 & 2520 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\cr -6 & 0 & 0 & 2520 & 0 & 0 & 0 & 0 & 0 & 0 & 0\cr -5 & 0 & 0 & 0 & 2520 & 0 & 0 & 0 & 0 & 0 & 0\cr -4 & 0 & 0 & 0 & 0 & 2520 & 0 & 0 & 0 & 0 & 0\cr -3 & 0 & 0 & 0 & 0 & 0 & 2520 & 0 & 0 & 0 & 0\cr -2 & 0 & 0 & 0 & 0 & 0 & 0 & 2520 & 0 & 0 & 0\cr -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2520 & 0 & 0\cr 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2520 & 0\cr 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2520\cr \hline c_{i,10} & 252 & -28 & 7 & -3 & 2 & -2 & 3 & -7 & 28 & -252\cr \end{array} $$

2. Die Werte zur Bildung der Linearkombination der zyklischen Formeln von Tendler lauten in der Reihenfolge $p=3,4,5$ und $p=6,7$ wie folgend:

$$ \displaylines{ \pmatrix{1&&\cr 1&0&\cr 1/3&-2/3&-2/3\cr},\qquad \pmatrix{1&&\cr 1&0&\cr 4&-4&5\cr},\qquad \pmatrix{1&&&\cr 1&0&&\cr 10&0&-1&\cr 10&-51/2&-31&1/2\cr}, \cr \pmatrix{1&&&\cr 20&-1&&\cr 20&-1&7&\cr 20&-78&-68&3\cr},\qquad \pmatrix{1&&&\cr 1&0&&\cr 20/7&-31/7&-10/5&\cr 10&-20&-5&2\cr}. \cr } $$

7. Äquilibrierung von Formeln

Es folgt die Definition von Äquilibrierungsmaßen für eine konsistente Formel $F$ der Form

$$ F = (\alpha,{\mskip 3mu}\beta) = (\alpha_0,\ldots,\alpha_\kappa, {\mskip 5mu}\beta_0,\ldots,\beta_\kappa). $$

1. Definition: Eine Formel $F$ hat

(1) Gaußsches-$\alpha$-Äquilibrierungsmaß

$$ \mu_\alpha^G := {1\over\alpha_\kappa^2}\sum_{i=0}^\kappa\left(\alpha_i^2+ \beta_i^2\right)={1\over\alpha_\kappa^2}\lVert F\rVert_2^2, $$

(2) Gaußsches-$\beta$-Äquilibrierungsmaß $\mu_\beta^G := \lVert F\rVert_2^2/\beta_\kappa^2$,

(3) Gauß-Henrici-Äquilibrierungsmaß $\mu_H^G := \lVert F\rVert_2^2 / (\sigma(1))^2$,

(4) Tschebyscheff-$\alpha$-Äquilibrierungsmaß

$$ \mu_\alpha^T := {1\over\lVert alpha_\kappa\rVert}\max_{i=0}^\kappa \left(|\alpha_i|,|\beta_i|\right) = {1\over\lVert \alpha_\kappa\rVert}\lVert F\rVert_\infty, $$

(5) Tschebyscheff-$\beta$-Äquilibrierungsmaß $\mu_\beta^T := \lVert F\rVert_\infty / \lVert\beta_\kappa\rVert$,

(6) Tschebyscheff-Henrici-Äquilibrierungsmaß $\mu_H^T := \lVert F\rVert_\infty / \lVert\sigma(1)\rVert$.

Sollten Zyklen mit geringem Äquilibrierungsmaß in den führenden Koeffizienten stark flukturieren, so kann man daran denken, bei der Suche nach günstigen Formeln die folgende Deviation der führenden Koeffizienten klein zu halten. Sei $\gamma:=\beta_{ii}/\alpha_{ii}$. Es gibt $({\ell\atop2})$ zyklusinterne Deviationen. Man begrenzt dann beispielsweise die Gaußsche Deviation der Koeffizienten, nämlich $\sum_{i<j}(\gamma_i-\gamma_j)^2$.

2. Beispiel: Sei $\ell=3$. Es gibt ${\ell\choose2} = {3\choose2}=3$ Deviationen. Man fordert Beschränktheit, oder in besonderen Situationen auch Minimalität, durch

$$ (\gamma_3-\gamma_2)^2+(\gamma_3-\gamma_1)^2+(\gamma_2-\gamma_1)^2 \cases{\to\min&\cr \le\rm const&\cr} $$

Es sollen nun einige Eigenschaften der Äquilibrierungsmaße hergeleitet werden. Es handelt sich hier um direkt aus der Definition herleitbare Aussagen.

3. Satz: Es gelten

(1) $\mu_\beta^T\ge1$, $\mu_\alpha^T\ge1$. Die Ungleichungen sind scharf, beispielsweise für die beiden Euler-Verfahren

$$ F=(-1,1,{\mskip 5mu}1,0) \qquad\hbox{und}\qquad F=(-1,1,{\mskip 5mu}0,1). $$

(2) $\mu_\alpha^G>1$, $\mu_\beta^G>1$.

(3) $\mu_H^T$ und $\mu_H^G$ sind stets definiert für $D$-stabile Verfahren, wegen $\rho'(1)=\sigma(1)\ne0$, da 1 keine mehrfache Nullstelle sein darf.

Die $\alpha$-Äquilibrierungsmaße treten vornehmlich im Zusammenhang mit expliziten gewöhnlichen Differentialgleichungen der Form $\dot y=f(t,y)$ auf. Die $\beta$-Äquilibrierungsmaße treten bei DAE auf, die man mit linearen Mehrschrittverfahren löst, und die Henrici-Äquilibrierungsmaße treten gehäuft im Zusammenhang mit Einbein-Verfahren auf. Man ist grundsätzlich, ersteinmal ohne Einschränkungen, an möglichst kleinen Äquilibrierungsmaßen interessiert. Bei Quadraturformeln ist aus Erfahrung bekannt, daß große Äquilibrierungsmaße ungünstig sind, da diese Quadraturformeln zu Oszillationen neigen. Jedoch bedeuten kleine Äquilibrierungsmaße nicht auch automatisch kleine Fehlerkonstanten. Gegenbeispiele anhand numerisch gewonnener Formeln sind bekannt.