, 26 min read
Die verwendeten zyklischen Formeln im Programm TENDLER
It is entirely possible, and even very likely, in view of the fact that the parameters $C_{ij}$ were “manually” set in an interactive computer search for the maximum $\alpha$, that better and simpler cyclic composite methods exist. An interested reader is invited to try the development of such methods on his own. For those, who are satisfied with the methods obtained thus far, in what follows, we describe how they are implemented in a numerical algorithm.
J.M. Tendler, Th.A. Bickart, Z. Picel (1976)
In diesem Abschnitt werden sämtliche benutzten Formeln vorgestellt, die in dem Programm TENDLER Verwendung gefunden haben. Zuerst kommen die Korrektorformeln und anschliessend die Prädiktorformeln.
Angegeben werden auch mehrere Fehlerfaktoren und Fehlerkonstanten und natürlich auch die Werte, die maßgeblich die Stabilitätseigenschaften quantifizieren, also die Werte $\alpha$ und $\delta$ bei der $A_\infty[\alpha]$-Stabilität und $S_\infty[\delta]$-Stabilität. Hingewiesen wird auch auf Charakteristika und auf den prinzipiellen Aufbau der zyklischen Formeln. Diskutiert wird ebenfalls, welche Bedingungen Stabilität bei unendlich in einfachster Weise garantieren.
Mit angeführt werden auch mehrere Darstellungen für diese Formelpaare; nämlich einmal in der gewöhnlichen Darstellung und das andere Mal in Form von rückwärtsgenommenen Differenzen. Die letztgenannte Form ist auch diejenige Form, die auch in dem Programm TENDLER verwendet wird. Die Darstellung in Form von rückwärtsgenommenen Differenzen bietet Vorteile insbesondere bei der Fehlerkontrolle und bei der Schrittweiten- und Ordnungssteuerung. Dort stellen sie gerade Näherungen für entsprechend benötigte Ableitungswerte dar. Die Korrektoriteration hingegen verhält sich neutral gegenüber der Darstellungsart. Ebenso der Prädiktor. In beiden Darstellungsarten, also sowohl der Darstellung in Form von rückwärtsgenommenen Differenzen als auch in gewöhnlicher Ordinatenform, sind die Prädiktorrechnungen einfach durchführbar. Die Prädiktorrechnungen haben einen häufig unterschätzten Einfluß auf die Gesamtrechenzeit.
- 1. Die sieben zyklischen Korrektorformeln von J.M. Tendler
- 2. Bemerkungen zu den Formeln von J.M. Tendler
- 3. Die Fehlerfaktoren der Korrektorformeln
- 4. Die Stabilitätseigenschaften der Korrektorformeln
- 5. Die sieben Prädiktorformeln
1. Die sieben zyklischen Korrektorformeln von J.M. Tendler
At first it was hoped that high order $A$-stable cyclic methods could be found. This goal was not achieved. However, the new integration methods are applicable to the solution of stiff systems. Moreover they exhibit better stability properties than the backward differentiation formulas in Gear's widely accepted algorithm.
1. Die Verfahren der Ordnung $p=1,\ldots,7$, die Joel Marvin Tendler im Rahmen seiner Dissertation 1973 neu entwickelt hat, werden hier nun angegeben.
Für die an anderer Stelle genauer erläuterte Fehlerkonstante von Henrici
werden auch gleichzeitig die Linkseigenvektoren $v_p$, für $p=1,\ldots,7$ der entsprechenden Matrizen aufgeführt. Der Rechtseigenvektor ist natürlich immer $w_p=(1,\ldots,1)^\top$, für $p=1,\ldots,7$. Zu dieser Konstanten vgl. man auch den Aufsatz von Albrecht (1985).
Die beiden ersten Verfahren, also die der Ordnung 1 und 2, sind nichts anderes als die zyklische Wiederholung ein und der selben Formel, hier also des impliziten Euler-Verfahrens für die Ordnung 1 und BDF2 für das Verfahren der Ordnung 2.
Bibliographisch:
- Joel Marvin Tendler: “A Stiffly Stable Integration Process Using Cyclic Composite Methods”, Ph.D. Diss., Syracuse University, Syracuse, New York, 26-Feb-1973
- Tendler, Joel Marvin und Bickart, Theodore A. und Picel, Zdenek: “STINT: STiff ordinary differential equations INTegrator. Program and User Manual”, Electrical and Computer Engineering Department, Syracuse University, Syracuse, New York 13210, USA, and Department of Electrical Engineering, Delft University of Technology, Delft, The Netherlands, Technical Report TR-76-12, December 1976, {\it ii/}+85~S.
- Theodore A. Bickart (1936—2023), obituary
- Peter Albrecht, 1985: “Numerical Treatment of O.D.E.s: The Theory of $A$-Methods”
- Euler, Leonhard (1707—1783)
- Henrici, Peter Karl Eugen (1923—1987)
- Zdenek Picel
2. Die Formeln der Ordnung $p=1,\ldots,7$ sind stets von der Gestalt
In Matrixschreibweise erhält man z.B. für $\ell=3$ und $k=2$, also für ein dreistufiges Verfahren mit 2 Startwerten, ausgeschrieben:
Aus Übersichtlichkeitsgründen werden die sieben nachfolgenden zyklischen Verfahren in der folgenden transponierten Schreibweise angegeben.
Später wird mit den Tableaux in transponierter Darstellung auch in dieser Form gerechnet.
3. Für die zyklische, dreifache Wiederholung des impliziten Euler-Verfahrens
und für die zyklische, dreifache Wiederholung der BDF2
ergeben sich die beiden transponierten Tableaux für die Verfahren wie untenstehend.
Ordnung 1.
Ordnung 2.
Für die Links- und Rechtseigenvektoren $v_1$, $v_2$, $w_1$ und $w_2$ erhält man
und für die Fehlerkonstanten von Henrici ergeben sich damit
Diese beiden Verfahren werden sowohl in dem Programm STINT als auch in dem neu entwickelten Programm TENDLER in dieser Form benutzt. Da beide Verfahren, also sowohl das implizite Euler-Verfahren und auch BDF2, optimale, also nicht verbesserbare Stabilitätseigenschaften $A_\infty^0[\alpha]$ und $S_\infty^0[\delta]$ besitzen, wurden sie unverändert so übernommen.
4. Für $p=k=3$, $p=k=4$ und $p=k=5$ konnte Tendler die folgenden Werte finden. Bei allen drei Verfahren sind wieder die beiden ersten Stufen der Zyklen BDF3, BDF4, bzw. BDF5. Nur die letzte Stufe sorgt, bzw. die letzten beiden Stufen sorgen für die Vergrößerung des Widlund-Winkels $\alpha$ bei der $A_\infty[\alpha]$-Stabilität und der Verkleinerung der Widlund-Distanz $\delta$ bei der $S_\infty[\delta]$-Stabilität.
Ordnung 3.
Ordnung 4.
Ordnung 5.
Die Vektoren $v_i$ mit $v_i\rho(1)=0$, die Vektoren $w_i$ mit $\rho(1){\mskip 3mu}w_i=0$ und die Fehlerkonstanten von Henrici sind
Wäre beim Verfahren dritter Ordnung auch noch die letzte Stufe gleich BDF3, also das gesamte Verfahren dreimal BDF3 hintereinander, so erhielte man für die Fehlerkonstante natürlich $\tilde C_3=-3/4$. Der Linkseigenvektor von $\rho(1)$, der zugehörige Rechtseigenvektor und die Fehlerkonstante von Henrici zu dem Verfahren 5.ter Ordnung berechnen sich zu
5. Man beachte bei allen diesen obigen Formeln, als auch bei den noch folgenden Formeln, den grundlegenden Aufbau. In gewisser Hinsicht wird das Potential, welches man prinzipiell durch Einführung der zyklischen Verfahren erhält, nicht in voller Gänze ausgeschöpft, jedoch stellt sich die vermeintliche Einschränkung unter anderen Gesichtspunkten (programmiermässige Einfachheit, Rechenbedarf, Speicherplatzbedarf) als gar keine Einschränkung heraus, sondern als echter Vorteil. Stets sind die Formeln von der Gestalt,
D.h. es handelt sich immer um die zyklische Hintereinanderschaltung linearer Mehrschrittverfahren mit fester Anzahl von Startwerten. Insbesondere heißt dies z.B. für die letzte Stufe, daß nicht alle möglichen $k+(\ell-1)$ Werte benutzt werden, sondern tatsächlich nur $k$, oder bei der Ordnung $p=3$ sogar nur $k-1$. Erneut wird deutlich, daß es sich hier um eine echte Verallgemeinerung der BDF auf Mehrstufigkeit handelt. Insbesondere werden also nicht alle möglichen zurückliegenden Ableitungswerte verwendet. Aufgrund des Einbein-Charakters der BDF sind natürlich auch andere Verallgemeinerungen denkbar.
Tischer (1983) in seiner Dissertation und Tischer/Sacks-Davis (1983) verfahren hier etwas anders. Hier werden in der letzten Stufe der Zyklen ein Startwert zusätzlich verwendet. Neben den potentiell besseren Stabilitätseigenschaften von zyklischen Verfahren ist dies ein weiterer, möglicher Vorteil von zyklischen Verfahren gegenüber einstufigen linearen Mehrschrittverfahren.
Beachtenswert ist auch die Benutzung der skalierten Ableitungswerte $hf_{m\ell+i}$. Wie bei den $y_{m\ell+i}$ werden hier nicht alle möglichen Werte wirklich benutzt. Zum maßgeblichen Teil hängt dies an der $A_\infty^0[\alpha]$- bzw. $S_\infty^0[\delta]$-Stabilität. Dennoch muß man hierbei im Auge behalten, daß man diese zusätzlichen Werte auch zu speichern hat. Insbesondere hat man auch für den Fall vorzusorgen, daß ein gesamter Zyklus zurückgewiesen wird und somit noch weitere Werte abgespeichert werden müssen. Darüber hinaus entstehen Rechnungen für die skalierten Ableitungen bei der Interpolation für evtl. benötigte Zwischengitterpunkte bei einem Schrittweitenwechsel.
6. Für die Ordnung $p=k=6$, wo nur die erste Stufe des Zykluses gleich der BDF6 ist und schließlich für $p=k=7$, für die die BDF noch niemals $D$-stabil ist, konnte Tendler die folgenden zyklischen, steif-stabilen Verfahren finden. Bei $p=k=7$ sind die ersten beiden Stufen wieder die BDF7, wobei man dann auch erneut ein Beispiel dafür sieht, daß die ^{zyklische Kombination klassischer (instabiler) Mehrschrittverfahren} das Gesamtverhalten gegenüber dem Einzelverhalten u.U. verbessern kann.
Ordnung 6.
Ordnung 7.
Die entsprechenden Eigenvektoren zu $\rho(1)$ und die Fehlerkonstante von Henrici sind hierbei
Für das letzte siebente Verfahren lauten die Eigenvektoren und die Fehlerkonstante von Henrici
Bibliographisch:
- Tischer, Peter E.: $ldquo;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. + 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. Bemerkungen zu den Formeln von J.M. Tendler
1. Man könnte es als entscheidenden Nachteil ansehen, daß die drei Verfahren der Ordnungen $p=5$, $p=6$ und $p=7$ nun 4 Stufen verlangen. Umfangreiche Testläufe mit dem Programm TENDLER zeigten nun jedoch, daß bei einer geringfügigen Modifikation der Schrittweiten- und Ordnungssteuerung Schrittzurückweisungen innerhalb, oder am Ende, eines Zykluses “fast nie” auftreten. Zudem sind Zurückweisungen in dem Programm TENDLER ohnehin die Ausnahme. Die Testläufe legten sogar nahe noch weiter zu gehen. Zyklen mit 5 Stufen wären völlig unproblematisch in das Programm TENDLER zu integrieren. Diese Erfahrungen stehen im Gegensatz zu den _Entwurfs_zielen bei den zyklischen Verfahren von Tischer (1983) (Dissertation) und Tischer/Sacks-Davis (1983), nicht jedoch zu den praktischen Erfahrungen, die Tischer/Gupta (1983) und Tischer/Gupta (1985) mit dem Programm ODIOUS gemacht haben.
Bibliographisch:
2. Auffällig ist bei den $\alpha_{ij}$-Werten die ^{Vorzeichenalternierung entlang herunter einer Spalte} und bei den $\beta_{ij}$-Werten sticht ins Auge, die “Konzentration” der Minus-Zeichen “in der Mitte”, falls nicht dort Nullen stehen, wiederum entlang herunter einer Spalte.
Stets sind die ersten beiden Stufe nichts anderes als die BDF, außer bei dem Verfahren sechster Ordnung.
Für die Abhängigkeit der Stufenzahl $\ell$ von der Ordnung $p$ erhält man übersichtlich zusammengefasst,
welche man auch mit der Formel $\ell=3+\lfloor p/5\rfloor$ berechnen könnte $(1\le p\le7)$.
3. Für alle Formeln von Tendler (1973) gilt, daß $z_i=hf_i$ für $i\le m\ell$ nicht benötigt werden, man vergl. hier auch das vorher angegebene Sternchen-Schema. Stets ist also $B_0=0$, bei invertierbarer Matrix $B_1$. Hiermit sichert man automatisch, daß die Nullstellen des charakteristischen Polynoms $\det Q(\mu,H)$ für $\mathop{\rm Re}H\to-\infty$ gegen Null streben, wegen
(in nicht ganz korrekter Schreibweise) und somit aufgrund von $B_0=0$ dann
Damit enthält das Spektrum des Matrixpolynoms $Q(\mu,\infty)$ nur die Null, falls $B_1$ invertierbar vorausgesetzt wird.
Dies ist nicht der einzige mögliche Weg sicherzustellen, daß die Eigenwerte bei $\mathop{\rm Re}\nolimits H=-\infty$ verschwinden. Wichtig ist lediglich, daß das Spektrum des Matrixpolynoms
nur aus der Null besteht. Für den Falle eines Polynoms der Form $B_1\mu+B_0$ ergibt sich daher die folgende Überlegung. Man beachte, daß man auf dieses lineare Polynom alle anderen Polynome zurückführen kann, u.U. durch Vergrößern der Dimension, z.B. durch das Begleitpolynom.
Im monischen Falle, also $\det B_1\ne 0$, besteht das Spektrum des besagten linearen Polynomes $B_1\mu+B_0$ genau dann nur aus der Null, wenn $B_1^{-1}B_0$ ähnlich ist zu
also die Matrizen $J_i$ Jordankästen zum Eigenwert Null sind. Der Fall $\det B_\kappa\ne0$ ist stets zu erreichen, wenn alle Stufen implizit, oder sogar blockimplizit sind.
Bei expliziten Stufen oder bei Mischung aus expliziten und impliziten Stufen, ist das reguläre, oder häufig sogar singuläre Matrizenbüschel $B_1\mu+B_0$ entsprechend zu untersuchen, siehe Gantmacher (1986).
Tischer (1983) (Dissertation) und Tischer/Sacks-Davis (1983) z.B. gehen hier den anderen möglichen Weg mit $B_0\ne0$.
Dennoch bedeutet $B_0=0$ im Falle $B_1\mu+B_0$, bzw. $B_i=0$ ($i=0,\ldots,\kappa$) im Falle $\sum_{i=0}^\kappa B_i \mu^i$, einen großen Rechen- und Speichervorteil, wenn der Prädiktor entsprechend gewählt wird. Sehr vereinfacht ausgedrückt heißt dies: ``Je mehr Nullen in den Stufen, desto besser."\space Die Vereinfachung liegt hier maßgeblich u.a. daran, daß der Prädiktor entsprechend gewählt werden muß und, daß Nullen bei den $y_{m\ell+i}$-Werten nicht automatisch auch zu Nullen bei der Implementierung führen, welche ja auf der Darstellung in Form von rückwärtsgenommenen Differenzen beruht. Man beachte, daß bei z.B. 1000 Schritten u.U. 2000–9000 mehr Vektoradditionen und Skalar-Vektor-Multiplikationen ausgeführt werden müssen, wenn man “Nullen verschenkt.”\space Es ist nicht ungewöhnlich, daß eine Funktionsauswertung billiger ist als die Bildung einer längeren Linearkombination von gespeicherten Vektoren.
Adams-Formeln und BDF sind hier typische Vertreter für beide möglichen Extrema, was z.T. ihren hohen Grad an Effizienz plausibel macht.
4. Beispiel: Für das explizite Eulerverfahren $y_{n+1}-y_n=hf_n$ erhält man mit den Vektoren
die vier Matrizen
Die Matrizen $A_1$, $A_0$, $B_1$ und $B_0$ sind natürlich nicht eindeutig. An der Nichtinvertierbarkeit von $B_1$ ändert dies aber nichts. Offensichtlich hat das Matrixpolynom $B_1\mu+B_0$, also hier das verallgemeinerte Eigenwertproblem, einen Eigenwert bei unendlich und damit ist das Verfahren natürlich nicht $A_\infty^0$-stabil. An dieser Stelle werden auch die Überlegung bezüglich des singulären Matrizenbüschels von Wichtigkeit.
5. Es sei nocheinmal darauf hingewiesen, daß für die Ordnung $p=1$ und $p=2$ keine originär neuen Formeln gesucht und damit auch nicht gefunden werden konnte. Es wurde tatsächlich dreimal hintereinander dasselbe Verfahren bei den beiden niedrigsten Ordnungen als zyklisches Verfahren benutzt. Im Hinblick darauf, daß sich die beiden Stabilitätsmerkmale, wie Widlund-$\alpha$-Winkel, also $A_\infty[\alpha]$-Stabilität und Widlund-Distanz, also $S_\infty[\delta]$-Stabilität, nicht weiter verbessern lassen, wird dies verständlich. Zudem bleibt damit auch der Wert $\gamma$={$\beta$}/{$\alpha$}, der bei der Iterationsmatrix $W=I-h\gamma J$ wichtig wird, über den gesamten Zyklus konstant, obwohl dies nicht so wesentlich ist, wie vermutet werden könnte. Kleinere Fehlerkonstanten und Schrittweitenwechsel-Stabilitätseigenschaften blieben bei der Auffindung der neuen zyklischen Formeln unberücksichtigt.
6. Das Tschebyscheffsche-$\alpha$-Äquilibrierungsmaß $\mu$ einer linearen Mehrschrittformel ist definiert durch
Hiermit wird ein erster Eindruck von der unterschiedlichen Gewichtung einzelner Terme in dem Verfahrensausdruck vermittelt.
Die zyklischen Formeln von Tendler weisen annähernd ein Äquilibrierungsmaß von $1:p$ auf, wobei $p$ die Ordnung des Zykluses ist.
Die BDF hingegen weisen annähernd ein Äquilibrierungsmaß von $1:{p\over2}$ auf, welches also lediglich halb so groß ist, wie die Äquilibrierungsmaße der Formeln von J.M. Tendler.
Diesem Punkte wurde beim Aufsuchen der Formeln keinerlei Berücksichtigung zuteil. Anders ist dies etwa bei den Formeln von Filippi/Kraska (1973).
Hier wurden die Formeln durch Lösungen eines linearen Optimierungsproblems bestimmt, wobei die Summe der Beträge der Koeffizienten der Verfahren in der Zielfunktion auftauchten. Die Formeln von Tischer (1983) (Dissertation) und Tischer/Sacks-Davis (1983) weisen besonders ungünstige Äquilibrierungsmaße auf. Die Maße bewegen sich von ca. $1:p^2$ bis zu ca. $1:2p^2$.
7. Zusammenstellend noch die Werte
welche wichtig sind für die Konvergenz der Picard-Iteration. Die Werte sind auf eine Dezimale hinter dem Komma genau.
8. Es seien jetzt wieder nur die Tendlerschen Zyklen betrachtet. Die Matrix $Q$ sei diejenige Matrix, wie sie bei der Umwandlung von Ordinatenwerten in Werte für rückwärtsgenommene Differenzen auftaucht. Durch Vormultiplikation mit der Matrix $VQ^\top V$, passenden Ranges, mit $V=(\delta_{n+1-i,j})$, ergeben sich für die Ordinatenwerte $y_{m\ell+i}$ die Werte in der folgenden Tabelle, also die Werte für die Darstellung der Formeln in Form von rückwärtsgenommenen Differenzen. Dieses sind jedoch noch nicht die Werte, die nachher wirklich gespeichert werden. Die Werte $z_{m\ell+i}$ bleiben hier unberücksichtigt, da sie nicht in der Form rückwärtsgenommener Differenzen verwendet werden. Dies liegt an der günstigen Eigenschaft der Tendlerschen Zyklen, daß die skalierten Ableitungswerte niemals interpoliert zu werden brauchen. Bei einem Schrittweitenwechsel muß ein Ableitungswert lediglich neu skaliert werden. Dies ist eine typische Eigenschaft der Tendlerschen Zyklen. Für anders gebaute Zyklen gilt dies nämlich nicht so ohne weiteres.
Wegen der ersten Zeile der Matrix $Q^\top$ und der Konsistenzbedingung $\rho(1)=0$, wird der Wert für $\nabla^0$ nicht mit vermerkt. Die erste Stufe bei allen zyklischen Formeln ist die BDF und hier gilt
Die Normierung von $\nabla y_{n+1}$ wird hier nicht gesondert durchgeführt.
Bibliographisch:
3. Die Fehlerfaktoren der Korrektorformeln
The new methods presented were obtained solely on the basis of optimizing the Widlund wedge angle associated with the method. It is felt that attention should also be focused on simultaneously minimizing the local discretization error. In addition, higher stiffly stable methods should be derived.
J.M. Tendler (1973)
Nun seien die Beträge der Fehlerfaktoren der sieben zyklischen Verfahren angegeben. Zum Vergleich werden auch die Fehlerfaktoren der BDF angeschrieben. Diese sind gerade die Fehlerfaktoren der ersten Stufe der 7 zyklischen Verfahren. Mit Fehlerfaktor $c_{i,p+1}$ wird hier die erste nichtverschwindende Komponente des Matrix-Vektor-Produktes aus $C_{p+1,k}$ und den Koeffizienten $(\alpha,\beta)^\top$ des Verfahrens bezeichnet, skaliert mit einer entsprechenden Fakultät. Es ist also
1. Die Fehlervektoren $\gamma_p$ für alle Ordnungen $p=1,\ldots,7$ und alle Stufen, lauten exakt ausgerechnet wie in der Tabelle sichtbar.
2. Normiert man stufenweise den führenden Koeffizienten des Verfahrens auf eins, so normieren sich entsprechend die Werte $c_{i,p+1}$, für $i=1,\ldots,\ell$. Für die zyklischen Formeln von Tendler ergeben sich dann nach kurzer Rechnung, die Werte in der entsprechenden Tabelle.
3. Zur Gegenüberstellung nun die gleichen Fehlerfaktoren diesmal jedoch skaliert mit dem Kehrwert der Summe der $\beta$-Koeffizienten, was gerade zur Fehlerkonstante von Henrici führt. Man beachte, daß hier erneut stets stufenweise die Fehlerfaktoren gebildet wurden.
Auffällig ist, daß die Tabellenwerte betragsmässig leicht äquilibrierter erscheinen und im ganzen gesehen, etwas größer geworden sind. Natürlich gilt hier immer $\sigma_{ip}(1)\ne0$, für alle Stufen $i=1,\ldots,\ell$ und alle Ordnungen $p=1,\ldots,7$. Das heißt, 1 ist einfache Wurzel jeder Stufe, nicht notwendig die betragsmässig größte.
4. Übersichtlich zum Vergleich gegenübergestellt, werden jetzt die zyklusweise gebildeten Fehlerkonstanten von Henrici aufgeführt. Man beachte, daß diese Konstante jeweils immer zyklusweise gebildet wird. Die Werte lauten für alle Ordnungen $p=1,\ldots,7$:
Diese Konstanten sind nun betragsmässig am größten, von allen bisher angegebenen Konstanten. Auffallend ist hier die Konstante für die Ordnung $p=3$.
5. Um einen Überblick über die Größenordnungen zu erhalten, seien auch die Fehlerkonstanten der Adams-Moulton-Verfahren in der folgenden Tabelle angeschrieben.
Wegen $\rho'(1)=\sigma(1)=1$ und $\alpha_k=1$ stimmen alle drei verschiedenen Arten an Fehlerkonstanten über-ein. Vergleicht man nun alle bisher angegebenen Konstanten, so muß man beachten vor welcher $h$-Potenz diese Konstanten stehen. Beispielsweise führt bei einem Verfahren siebenter Ordnung eine Halbierung der Schrittweite asymptotisch betrachtet, zu einer Erhöhung der “Genauigkeit” um den Faktor 256. Überdies ist der Unterschied zwischen einem Verfahren $p$-ter und $(p+1)$-ter Ordnung nur gering, falls die Schrittweite variabel sein kann. Ist das Verfahren $p$-ter Ordnung nicht mehr ausreichend genau genug, so ist es gerade die zentrale Aufgabe der Ordnungssteuerung dies zu erkennen und durch Wahl einer passenderen Ordnung, den Genauigkeitsverlust wieder auszugleichen. Weiterhin beachte man die Wachstumseigenschaften der Fehlerkonstante von Henrici bei mehrfacher Wiederholung einer Stufe.
Bibliographisch:
4. Die Stabilitätseigenschaften der Korrektorformeln
While methods suitable for non-stiff problems, $\ldots$, are easy to find and to analyse, methods appropriate for stiff problems seem to be able to hide their presence much more successfully. However, there are such methods $\ldots$
J.C. Butcher (1987)
Zum Vergleich seien die Werte für den Widlund-$\alpha$-Winkel der Formeln von J.M. Tendler mit denen der BDF in der untenstehenden Tabelle gegenübergestellt. Man beachte, daß die BDF7 nicht mehr $D$-stabil ist, erst recht also nicht mehr $A_\infty[\alpha]$-stabil. Rechts daneben ist auch eine Gegenüberstellung der Widlund Distanzen $\delta$, bei der $S_{\infty}[\delta]$-Stabilität.
Widlung-$\alpha$-Winkel. Für die Berechnung der Winkel der BDF vgl. man Lajos Lóczi1 (2019) und Akrivis/Katsoprinakis (2019): Maximum angles of $A(\vartheta)$-stability of backward difference formulae, PDF.
| p | TENDLER | BDF |
|---|---|---|
| 1 | 90.00 | 90.00 |
| 2 | 90.00 | 90.00 |
| 3 | 89.43 | 86.03 |
| 4 | 80.88 | 73.35 |
| 5 | 77.48 | 51.84 |
| 6 | 63.25 | 17.84 |
| 7 | 33.53 | * |
Widlung-Distanz $\delta$.
| p | TENDLER | BDF |
|---|---|---|
| 1 | 0.0000 | 0.0 |
| 2 | 0.0000 | 0.0 |
| 3 | 0.0048 | 0.083 |
| 4 | 0.2400 | 0.67 |
| 5 | 1.4000 | 2.3 |
| 6 | 2.9000 | 6.1 |
| 7 | 10.2000 | * |
Die Definition der Begriffe Widlund-Winkel, Widlund-Distanz, $A_\infty[\alpha]$-Stabilität, $S_\infty[\delta]$-Stabilität findet man bei Tendler (1973) und Tendler/Bickart/Picel (1978), aber sinngemäß auch bei Albrecht (1979), §7.2., und weiteren Autoren.
Bibliographisch:
- Butcher, John Charles (*1933)
- Albrecht, Peter: “Die numerische Behandlung gewöhnlicher Differentialgleichungen — Eine Einführung unter besonderer Berücksichtigung zyklischer Verfahren”, Carl Hanser Verlag, München Wien, 1979, xi+193 S.
5. Die sieben Prädiktorformeln
The difficulty of resolving various other questions concerned with incorporating general linear methods into practical software should not be underestimated.
John C. Butcher (1987)
Es seien jetzt noch die Prädiktorformeln angegeben. Bei ihnen handelt es sich um nichts anderes als die expliziten BDF\null. Diese Formeln sind für Ordnungen $p>2$ nicht mehr $D$-stabil, d.h. sie dürfen nicht alleinstehend ohne einen Korrektor verwendet werden. Das Verfahren erster Ordnung ist natürlich nichts anderes als das explizite Eulerverfahren $y_{n+1}=y_n+hf_n$, und die Prädiktorformel zweiter Ordnung ist die (explizite) Mittelpunktsregel $y_{n+1}=y_{n-1}+2hf_n$ (= explizite ^{Nyström-Formel zweiter Ordnung}).
Die Angabe der Werte geschieht zunächst in der Form
für die Ordnungen $p=1,\ldots,7$. Die Werte $\beta_{n-i}$ für $i>0$ sind sämtlich Null, d.h. heißt es werden vom Prädiktor keine weiteren vergangenen skalierten Ableitungen benutzt.
Zuerst nun die Werte $\alpha_{n-i}$ und $\beta_n$ in der obigen Darstellung mit gewöhnlichen Ordinatenwerten $y_{m\ell+i}$, siehe erste Tabelle.
Die Darstellung der Prädiktorformeln in Form von rückwärtsgenommenen Differenzen, also in der Darstellungsform
für die Ordnungen $p=1,\ldots,7$, entnimmt man der zweiten Tabelle.
In kompakter Form seien noch diejenigen Werte angegeben, die auch tatsächlich im Programm TENDLER benutzt werden, wenn auch nicht in dieser kompakten Form. An späterer Stelle wird nocheinmal kurz auf die Herleitung dieser Koeffizienten eingegangen. Die nachstehende Tabelle enthält in gewisser Hinsicht ein “Gemisch” von Prädiktor- und Korrektorformeln. Dabei werden für $y$ rückwärtsgenommene Differenzen verwendet, für $z$ jedoch nicht.
Hierbei beachte man, daß die Werte in Tendler/Bickart/Picel (1976) für $p=6$, $i=4$ nicht sämtlich richtig sind. In dem Programm STINT und auch in dem Programm DIFJMT werden jedoch die korrekten Werte gespeichert.
Bibliographisch:
- Tendler, J.M. + Bickart, Th. + Picel, Z.: “STINT: STiff ordinary differential equations INTegrator. Program and User Manual”, Electrical and Computer Engineering Department, Syracuse University, Syracuse, New York 13210, USA, and Department of Electrical Engineering, Delft University of Technology, Delft, The Netherlands, Technical Report TR-76-12, December 1976, ii+85 S.