7th January 2025, 48 min read

Praktische Gewinnung zyklischer, steif-stabiler Verfahren

Original post is here eklausmeier.goip.de/blog/2025/01-07-praktische-gewinnung-zyklischer-steif-stabiler-verfahren.

• 1. Definition von Stabilitätseigenschaften
• 2. Schreibweise zusammengesetzter Mehrschrittformeln
• 3. Stabilität und Polynome
• 4. Verschärfung von Stabilitätsdefinitionen
• 5. Konstruktive Herleitung steif-stabiler zyklischer Formeln
• 6. Basisformeln der Tendlerschen Zyklen
• 7. Äquilibrierung von Formeln

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_{\mathrm{\infty }}^r[\alpha ]$- und $S_{\mathrm{\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_{\mathrm{\infty }}^r[\alpha ]$- und $S_{\mathrm{\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

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

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,\dots ,\ell$, die durch die folgenden $\ell$ Mehrschrittformeln bestimmt wurden

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=h{f}_i=h{\dot{y}}_i$ benutzt.

Der obige Ausdruck wird zyklisch, falls gilt

Gälte ${\alpha }_{ij}\ne 0$ oder ${\beta }_{ij}\ne 0$ 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, gegenüber 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

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

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

gegen Null gehen ($n\to \mathrm{\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 $\mathcal{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

Nun heißt ein Verfahren $S[\delta ]$-stabil, falls so eine Menge $\mathcal{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 \mathrm{\infty }$, für eine feste positive Schrittweite $h$, falls immer $\lambda h=H\in \mathcal{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 $\mathcal{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 $\mathcal{A}[\alpha ]$ eine zusammenhängende und offene Teilmenge der komplexen Zahlenebene mit Elementen $z$ der Form

Jetzt heißt das Verfahren $A[\alpha ]$-stabil, wenn es w.o. eine derartige Menge $\mathcal{A}[\alpha ]$ gibt derart, daß die berechneten Näherungen dieses Verfahren gegen Null streben, für $n\to \mathrm{\infty }$, wenn man das Verfahren auf die Differentialgleichung $\dot{y}=\lambda y$ anwendet. Diese Abklingeigenschaft soll immer dann gelten, wenn $\lambda h=H\in \mathcal{A}[\alpha ]$ gilt, für feste positive Schrittweite $h$. Für $\alpha \notin [0,{90}^{\circ }]$ 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. $\mathcal{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 $\mathcal{A}[\alpha ]$, bzw. $\mathcal{S}[\delta ]$, “größer”.

Bei den beiden Defintionen beachte man, daß $\mathcal{S}[\delta ]$ die Menge bezeichnet und $S[\delta ]$ die dazugehörige Eigenschaft. Entsprechend gilt dies für $\mathcal{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_{\mathrm{\infty }}[\alpha ]$ und $S_{\mathrm{\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_{\mathrm{\infty }}^0[\alpha ]$- und $S_{\mathrm{\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_{\mathrm{\infty }}^0[\alpha ]$- und $S_{\mathrm{\infty }}^0[\delta ]$-stabil.

(4) Die 10 blockimpliziten Verfahren von Bickart/Picel (1973) sind alle $A[\alpha ]$- und $S[\delta ]$-stabil, sogar $A_{\mathrm{\infty }}^0[\alpha ]$- und $S_{\mathrm{\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

kann man als eine einzige Matrixdifferenzengleichung notieren, nämlich

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 \le 2$, 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,\dots ,\kappa$ lauten:

und

beziehungsweise

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

für $i=0,\dots ,\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

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 $|A\otimes I_n|=\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_{\mathrm{\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

gehörigen Koeffizienten lauten

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

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

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

Die Stabilität definierenden Matrizen sind:

Die charakteristische Gleichung ist

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 $detQ(\mu ,H)=\mathcal{O}(H^{p+1})$ ($\mu =e^H$) scheitert, weil z.B. für DH3 gilt

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}=A{z}_n+h{\phi }_n$, oder $L{z}_{n+1}=U{z}_n+h{\phi }_n$.

3. Stabilität und Polynome

1. Das charakteristische Polynom lautet

Man beweist jetzt, daß das verwendete Verfahren genau dann $A[\alpha ]$-stabil bzw. $S[\delta ]$-stabil ist, wenn die sämtlichen Nullstellen $\mu$ von $detQ(\mu ,H)=0$, alle betragsmässig kleiner eins sind, für alle $H\in \mathcal{A}[\alpha ]$, bzw. alle $H\in \mathcal{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. $\mathcal{A}[\alpha ]$ und $\mathcal{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 $detQ(\mu ,H)=0$ im Komplement von $\mathcal{A}[\alpha ]$, bzw. im Komplement von $\mathcal{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 $detQ(\mu ,H_0)$ für ein festgehaltenes beliebiges $H_0\in \mathcal{A}[\alpha ]$, bzw. $H_0\in \mathcal{S}[\delta ]$, sämtlich betragsmässig kleiner 1 sind und gleichzeitig
2. die Kurve der Punkte $H$, die $detQ(\mu ,H){|}_{|\mu |=1}=0$ erfüllen, im Komplement von $\mathcal{A}[\alpha ]$, bzw. im Komplement von $\mathcal{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 $|H|$ $(\mathrm{Re}H<0)$ immer größer wird, die Nullstellen $\mu$ von $detQ(\mu ,H)=0$ betragsmässig gegen eins streben könnten. Integriert man jedoch $\dot{y}=\lambda y$ ($\mathrm{Re}\lambda <0$) so würde dies heißen, daß es eine Lösung gibt mit $|y_{n+1}/y_n|\to 1$. Diese Eigenschaft spiegelt sicherlich nicht das wahre Verhalten der Lösung wider. Diese Eigenheit hat z.B. die Trapezregel

Aus diesem Grunde sind Verfahren günstiger bei denen die Nullstellen $\mu$ bei $\mathrm{Re}H=-\mathrm{\infty }$ von $detQ(\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 $\mathrm{Re}H\to -\mathrm{\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 $\mathrm{Re}H\to -\mathrm{\infty }$ in der offenen, Kreisscheibe $\{z\in \mathbb{C}:|z|<r\}$ liegen, nennt man $S_{\mathrm{\infty }}^r[\delta ]$-stabil. Für $r=1$ sagt man abkürzend nur $S_{\mathrm{\infty }}[\delta ]$-stabil.

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

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

(4) Der Wert $r$ bei der $S_{\mathrm{\infty }}^r[\delta ]$- bzw. $A_{\mathrm{\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_{\mathrm{\infty }}^0$-Stabilität verallgemeinert offensichtlich den Begriff der L-Stabilität für Einschrittverfahren.

Offensichtlich ist jedes $A_{\mathrm{\infty }}^{r_1}[\alpha ]$-stabiles Polynom auch $A_{\mathrm{\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_{\mathrm{\infty }}^{r_1}[\delta ]$-stabile Polynom ist natürlich auch $S_{\mathrm{\infty }}^{r_2}[\delta ]$-stabil, falls $r_1<r_2$. Insbesondere ist jedes $A_{\mathrm{\infty }}^0[\alpha ]$-stabile Polynom auch $A_{\mathrm{\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_{\mathrm{\infty }}^0[\alpha ]$- und $S_{\mathrm{\infty }}^0[\delta ]$-stabile Verfahren. Aber auch die blockimpliziten Verfahren von Bickart/Picel (1973) sind $A_{\mathrm{\infty }}^0[\alpha ]$- und $S_{\mathrm{\infty }}^0[\delta ]$-stabile Verfahren.

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

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

Für $\xi =-4$ erhält man die BDF3, welche $A_{\mathrm{\infty }}^0[{86.03}^{\circ }]$- und $S_{\mathrm{\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 $detQ(\mu ,H)$ die Gleichung

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

Interessant hierbei ist, daß eine Kombination impliziter Verfahren mit expliziten Verfahren dennoch $A_{\mathrm{\infty }}^0[\alpha ]$-stabil, bzw. $S_{\mathrm{\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

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:

Dies sind also diejenigen Zahlen $H$, sodaß $detQ(\mu ,H)$ verschwindet, dabei $\mu$ aber stets aus der offenen Kreisscheibe $D_r$ ist. Weiter werden jetzt die Mengen $\mathcal{A}[\alpha ]$ und $\mathcal{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

und

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

Mit $\phi$ durchläuft man dann eine diskrete Teilmenge von $[0^{\circ },{180}^{\circ }]$, zum Beispiel $\phi =0^{\circ },{20}^{\circ },\dots ,{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 $|\cdot |$ wähle dann unter den mehreren möglichen “Funktionswerten”, den betragsmässig größten aus, ordnet also letztlich jedem $r{e}^{i\phi }$ eine komplexe Zahl $H$ zu. Das bivariable Polynom $detQ(\mu ,H)$ ist natürlich algebraisch. Lösungsäste algebraischer Funktionen haben höchstens algebraische Singularitäten. $\partial {\mathcal{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 $detQ(\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

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

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

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

in etwas nachlässiger, aber sehr suggestiver Schreibweise. Es ist hier interessant zu vermerken, daß jetzt $A_{\mathrm{\infty }}[\alpha ]$- und $S_{\mathrm{\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_{\mathrm{\infty }}[\alpha ]$-stabil, bzw. $S_{\mathrm{\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_{\mathrm{\infty }}^r[\alpha ]$-stabil, bzw. $S_{\mathrm{\infty }}^r[\delta ]$-stabil genau dann, wenn:

1. die Nullstellen von $detQ(\mu ,\mathrm{\infty })$ liegen in der offenen Einheitskreisscheibe und
2. die Kurve der Punkte $H$, die $detQ(\mu ,H){|}_{|\mu |=r}=0$ erfüllen, im Komplement von $\mathcal{A}[\alpha ]$, bzw. im Komplement von $\mathcal{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 $detQ(\mu ,\mathrm{\infty })=0$ zu Null setzt. Dies ist z.B. der Fall, wenn man fordert:

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\ne 0$, für $i=0,\dots ,\kappa$.

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

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

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

Dann ist

Sei $\mu =1$ einfache Nullstelle von

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

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

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

9. Definition: Ein Matrixpolynom

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

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

Das Matrixpolynom $Q(\mu ,H)$ ist konsistent mit der genauen Ordnung $\underset{i=1}{\overset{\ell }{min}}{p}_i$. Gilt $p_i\ge 1$ ($\mathrm{\forall }i$), so heißt das Matrixpolynom konsistent. Möchte man die besondere Form

besonders betonen, so spricht man von polynomial-konsistent.

Die obige Definition verallgemeinert $\rho (e^h)-h\sigma (e^h)=\mathcal{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 $detQ(\mu ,H_0)=0$ gilt: $|\mu |<1$.

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

heißt Menge der absoluten Stabilität.

Offensichtlich gilt: Enthält ${\mathcal{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 $detQ(\mu ,0)=0$ ist, heißt relativ-stabil für ein $H_0$ genau dann, wenn gilt

${\mu }_1(H)$ ist also die zu ${\mu }_1(0)=1$ gehörende Nullstelle, wobei $Q({\mu }_i(H),H)=0$ ($i\ge 1$). Hat $detQ(\mu ,H)$ in $\mu$ nur den Grad 1, so sei das Matrixpolynom $Q(\mu ,H)$ relativ-stabil.

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

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 $\mathrm{Re}\lambda \ll 0$, 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

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_{\mathrm{\infty }}^r[\delta ]$- und $LR$-$A_{\mathrm{\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 $\mathcal{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. $\dots$
The test for stiff stability of a composite linear multistep method presented herein is entirely graphical. $\dots$
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_{\mathrm{\infty }}^0$-stabile zyklische Verfahren an bis zur Ordnung $p\le 4$. Über der Ordnung $p>4$ sind die weiterhin zweistufigen Zyklen noch $A_{\mathrm{\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

Die BDF erster und zweiter Ordnung

sind beide $A_{\mathrm{\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

Da nur Linearkombinationen von diesen Mehrschrittformeln gebildet werden, ist das Erfülltsein von (LMSV) ebenfalls sofort klar. Dann wurde für jedes $i=1,\dots ,\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:

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

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: