4th April 2026, 41 min read

TENDLER: 6. Die Schrittweiten- und Ordnungssteuerung

Fortsetzung der TENDLER Programmbeschreibung.

1. TENDLER: 1. Grobaufbau und prinzipielle Überlegungen
2. TENDLER: 2. Benutzung des Programmes
3. TENDLER: 3. Der Prädiktor
4. TENDLER: 4. Die Korrektoriteration
5. TENDLER: 5. Die Fehlerkontrolle

Inhalt.

• 1. Prinzipien der Schrittweiten- und Ordnungssteuerung
• 2. Die Schrittweitensteuerung im Programm TENDLER
• 3. Die Idee des Shampine-Testes
• 4. Der Shampine-Test im Programm TENDLER
• 5. Die Auswahl der Ordnung im Programm TENDLER
• 6. Die beiden Sperren wait[] und problem im Programm TENDLER
• 7. Die Schrittweiten- und Ordnungssteuerung im Programm DE/STEP

For general use, linear multistep formulas are of little value without a means of selecting values of the step size $h$ and method order $q$ for which the method is reasonably accurate and efficient. The algorithm for doing this constitute a major distinction between modern ODE software and its obsolete counterparts. These algorithms are the result of considerable research and development efforts, which are still continuing, and their impact on the accuracy and economy achievable with modern ODE codes is often quite dramatic.

G.D. Byrne, A.C. Hindmarsh (1987)

Hier wird sowohl die Schrittweiten- und Ordnungssteuerung in dem Programm TENDLER als auch in dem Programm DE/STEP beschrieben. Die Schrittweiten- und Ordnungssteuerung in dem Programm TENDLER wird nach dem Durchlaufen aller Stufen ausgeführt. Die Schrittweiten- und Ordnungssteuerung wählt eine passende Schrittweite $h$ und Ordnung $p$ für den nächsten Zyklus aus. Hierduch wird versucht, sowohl den geforderten Genauigkeitswünschen des Benutzers zu entsprechen, als auch die Rechenzeit klein zu halten. Es wird auch der Shampine-Test beschrieben, der dazu dient, zu erkennen, ob eine Differentialgleichung nun im Ganzen gesehen als steif oder als nicht-steif zu klassifizieren ist. Besondere Bedeutung kommt zwei Variablen zu, die nicht nur alleine in der Schrittweiten- und Ordnungssteuerung gesetzt werden, sondern auch von der Korrektoriteration mit gesteuert werden. Ihren Einfluß entfalten sie jedoch erst in dem Schrittweiten- und Ordnungssteuerungssegment. Auf alternative Zugänge zu den Betrachtungen in dem Schrittweiten- und Ordnungssteuerungssegment wird hingewiesen. Grundlegende Prinzipien, wie sie in gängigen Programmen basierend auf linearen Mehrschrittverfahren immer wieder zur Anwendung gelangen, werden den Überlegungen vorangestellt.

Biographisch:

1. Byrne, Georg Dennis
2. Hindmarsh, Alan C. (*1942)

1. Prinzipien der Schrittweiten- und Ordnungssteuerung

For a general problem it is most unreasonable to expect the user to have any idea what is a reasonable order or step size.

C.W. Gear (1971)

1. Hier soll die grundlegende Idee der Schrittweiten- und Ordnungssteuerung erklärt werden. Für ein Programm basierend auf zyklischen Verfahren und Benutzung von rückwärtsgenommenen Differenzen, wird man i.d.R. in einer der unten beschriebenen Art und Weise vorgehen. Durch die rückwärtsgenommenen Differenzen ist nahe gelegt, wie man die Schrittweite und Ordnung für den nächsten Zyklus schätzt und durch die Zyklushaftigkeit werden bestimmte Sicherheitsfaktoren nahegelegt. Zum einen könnte man die Differenz von Prädiktornäherung und Korrektorwert als ein erstes Indiz für die Genauigkeit benutzen. Dieser Weg wird jedoch in dem Programm TENDLER nicht beschritten, da man bei steifen Differentialgleichungen dem Prädiktor noch weniger Vertrauen schenkt, als man dies bei nicht-steifen Gleichungen tut. Hingegen das Programm DE/STEP tut dies; es ist allerdings auch nur für nicht-steife Differentialgleichungen konzipiert. Für nicht-steife Differentialgleichungen ist es sehr erfolgreich.

Biographisch:

1. Gear, Charles William (1935—2022)

2. Es sei $\mathbb{P}$ die Menge der möglichen Ordnungen für den nächsten Zyklus. Für $j\in\mathbb{P} := \{p-1,p,p+1\} \cap \{1,\ldots,7\}$ berechnet man die Größen

wobei für die Konstante $\sigma_{pj}$ gilt, daß $0<\sigma_{pj}\le1$. Die Werte $\overline\eta_{pj}$, die in der folgenden Tabelle aufgelistet erscheinen, sind Schrittweitenbegrenzer, welche von der Stufenanzahl $\ell$ und der Ordnung $p$ abhängen. Diese Begrenzer dienen dazu, die Schrittweite in der Art zu begrenzen, daß man ausschließlich mit schon berechneten Werten für die Interpolation auskommt. Die Schrittweite kann also nur so verändert werden, wenn genügend viele alte Werte vorhanden sind, um die benötigten Startwerte für den neuen Zyklus mit der neuen Schrittweite, durch Interpolation zu erhalten. Andernfalls, wenn dies nicht der Fall ist, so wird die Schrittweite nicht geändert oder zumindestens nicht maximal möglich geändert. Der theoretische Hintergrund für die obige Rechenvorschrift für $\eta_j$ ist die Formel

Hierbei ist wieder $d_{p+1}=A_\kappa^{-1}c_{p+1}$, mit dem Diskretisierungsfehlervektor $c_{p+1}$. Die entsprechend hohe Ableitung multipliziert mit der entsprechenden $h$-Potenz wird wie üblich durch die rückwärtsgenommene Differenz ersetzt.

Man kann natürlich keine optimale Schrittweitenfolge erwarten; interessiert ist man an einer günstigen und brauchbaren Schrittweitenfolge. Die optimale Schrittweitenfolge für eine feste Differentialgleichung und feste Toleranz wäre diejenige Folge von Schrittweiten, die minimalen Rechenaufwand erfordert, ohne daß der Fehler der gefordert wird, überschritten wird. Diese optimale Schrittweitenfolge ist natürlich u.a. von der Rechenanlage abhängig, denn das Verhältnis von benötigter Rechenzeit von Gleitkommaaddition zu Gleitkommamultiplikation ist in jedem Falle von der zugrundeliegenden Rechenmaschine abhängig.

3. Durch die Beschränkung auf die Menge $\mathbb{P}$ wird also nur die um eins niedriegere, die aktuelle und die um eins höhere Ordnung betrachtet. Ein Wechsel z.B. von der Ordnung 2 direkt auf die Ordnung 5 ist damit nicht möglich. In dem Programm TENDLER ist es daher prinzipiell nicht möglich einen solchen sprungartigen Anstieg der Ordnung durchzuführen. Allerdings ist es möglich in langsamerer Progression von der Ordnung 2 auf die Ordnung 5 überzuwechseln. Es gibt mindestens ein Programm basierend auf linearen Mehrschrittverfahren, indem es tatsächlich möglich ist von einer Ordnung $p_1$ auf eine Ordnung $p_2$ überzuwechseln, mit $|p_1-p_2|>1$, insbesondere ist ein Ordnungssprung nach oben um beispielsweise 4 Ordnungen möglich. Üblich ist dies jedoch nicht. Aufgrund der Differenzentabelle, die im Speicher aufgebaut wird, ist eine langsamere Erhöhung, insbesondere in Einerschritten, einfacher zu programmieren und erfordert schließlich auch weniger Rechenaufwand. Folgende Ordnungsfolge kann jedoch ebenfalls nicht vorkommen: 1, 2, 3, 4, 5, 6, 7. Ein direktes Aufsteigen zu der Ordnung $p=7$ in dieser linearen Form ist somit nicht möglich.

Zur Vermeidung dieses Sachverhalts sorgt die Sperrvariable wait[]. Wohl aber möglich ist 1, 2, 3, 4, 4, 5, 5, 5, 6, 6, 6, 6, 7, 7, $\ldots$ . Das letzte ist nur dann wahrscheinlich bei sehr strenger Genauigkeitsanforderung $\varepsilon$. Denkbar ist auch 7, 7, 2, 1; dies muß man jedoch als Ausnahme werten. Dies kann auch nur dann auftreten, falls sich große Probleme während der Integration einstellen. Ein abrupter Wechsel von der höchsten Ordnung $p=7$ zu der niedrigsten Ordnung $p=1$ ist Folge einer besonders schwerwiegenden Problematik, die das Programm TENDLER dadurch versucht zu begegnen, indem versucht wird diesen Schwierigkeiten durch rasche Veränderung der Ordnung beizukommen. Die letzte Folge von Ordnungen 7, 7, 2, 1, wird nicht von der Schrittweiten- und Ordnungs-Steuerung hervorgerufen. Im steifen Modus werden im Programm TENDLER die niedrigeren Ordnungen ggü. den höheren Ordnungen bevorzugt. Im nicht-steifen Modus ist dieser Trend nicht so ausgeprägt.

Niedrigere Ordnungen haben neben ihren i.a. besseren Stabilitätseigenschaften auch den sehr bedeutsamen Vorteil des sehr viel geringeren Rechenaufwandes bzgl. Interpolation und Bildung des Prädiktors. Zudem sind die Veränderungen der Schrittweite, die durchgeführt werden können, größer als diejenigen, die bei den höheren Ordnungen durchführbar sind. Man vgl. hier die nachfolgende Tabelle mit den Schrittweitenbegrenzern $\overline\eta_{pj}$.

4. Es wird nun prinzipiell wie folgt verfahren:

1. Ist $\max \{\eta_j:j\in\mathbb{P}\}{\mskip 3mu}<{\mskip 3mu}\texttt{UPTHRESH}$, so wird die Schrittweite und die Ordnung nicht gewechselt. Dies liegt daran, daß der Aufwand in diesem Falle als zu hoch angesehen wird.
2. Ist jedoch $\max \{\eta_j:j\in\mathbb{P}\}{\mskip 3mu}\ge{\mskip 3mu}\texttt{UPTHRESH}$, so wird $\eta$ auf den Wert $\max \{\eta_j:j\in\mathbb{P}\}$ gesetzt und $\hat p=p$, falls $\eta_p=\eta$ und sonst wird $\hat p$ auf den kleinsten Index $j$ derart gesetzt, sodaß $\eta_j=\eta$ ist.
3. Zum Schluß wird die neue Schrittweite als $\hat h = \eta \cdot h$ gewählt.

Die obigen drei Regeln sind das grundlegende Prinzip der Schrittweiten- und Ordnungssteuerung. Die wirklichen Regeln sind verschachtelter. Es sei für den Moment vorausgesetzt, daß der Schrittweitenbegrenzer $\overline\eta_{pj}$ keine Rolle spielt, also $\overline\eta_{pj}=+\infty$ ist. Wird nun der Abschnitt für die Schrittweiten- und Ordnungssteuerung betreten, so sei vereinfachend angenommen, daß $\left\Vert\nabla^{\hat p+1}y_{m\ell+\ell}\right\Vert_w$ der Wert des entsprechenden Ausdrucks für die zweite Stufe im nächsten Zyklus sei. Die Vereinfachung liegt hier also darin, daß man eine Konstanz der entsprechend hohen Ableitung der Lösungsfunktion samt Schrittweite voraussetzt und damit eine entsprechende Konstanz der rückwärtsgenommenen Differenz. Diese Annahme birgt keine allzu großen Gefahren in sich, da ja bei einer sehr raschen Änderung der hohen Ableitung, die Fehlerkontrolle versagen würde und damit die in diesem Falle ungerechfertigte Annahme abfangen würde.

Nun würde (Phase 2) der Fehlerkontrolle soeben mit Gleichheit passieren, unter der Annahme, daß $\sigma_{pj}$ gleich 1 sei. In diesem speziellen Falle wäre also $\hat h$ die größtmögliche Schrittweite. Nun ist aber i.d.R. $\left\Vert\nabla^{\hat p+1}y_{m\ell+\ell}\right\Vert_w$ nicht gleich $\left\Vert\nabla^{\hat p+1}y_{(m+1)\ell+2}\right\Vert_w$ und den weiteren entsprechenden Werten. Deswegen wurde der Faktor $\sigma_{pj}$ eingeführt, um eine gewisse “Sicherheitszone” zu erreichen.

Durch ausgiebige Probeläufe wurde beobachtet, daß Sicherheitsfaktoren in dem Bereich 0.6 bis 0.9 akzeptable Resultate liefern. Diese Grenzen sind allerdings nicht kritisch. Faktoren im Bereich 0.3 bis 0.6 sind ebenfalls denkbar, aber nicht sehr empfehlenswert. Hierduch würde nämlich dann die Schrittweitensteuerung immer mehr zum Hemmschuh für eine große Schrittweite ausarten. Die genauen Werte sind in den Konstanten nssaffac[][] (non-stiff safety-factors) und stsaffac[][] (stiff safety-factors) abgelegt.

Für die Fehlerkontrolle, die ja ebenfalls eine zukünftige Schrittweite voraussagen muß, ergaben sich Sicherheitsfaktoren in einer ähnlichen Größenordnung. Hierbei muß jedoch dem Rechnen an der Grenzgenauigkeit der Maschine gut Rechnung getragen werden. Es sei betont, daß dieser Wertebereich längeren Probeläufen entnommen wurde und nicht etwa einem Beweis entspringt. Das Programm STINT verwendet vergleichsweise große Sicherheitsfaktoren und schenkt den Besonderheiten des Rechnens am Rande der Maschinengenauigkeit keinerlei besondere Beachtung.

5. Man beachte, daß die Schrittweite und damit die Ordnung nur dann gewechselt werden, wenn sich die neue Schrittweite von der alten um mindestens $100(\texttt{UPTHRESH}-1)$ Prozent unterscheidet. Der Wert der Konstanten $\texttt{UPTHRESH}$ beträgt z.Z. 1.16. Der Grund für die Einführung eines solchen Mindestanstieg der Schrittweite liegt darin, daß der Aufwand für die Interpolation der vergangenen Werte für das neue Gitter, sich als lohnend erweisen sollte. In den Programmen DIFSUB, GEAR und LSODE wird der Faktor 11/10 verwendet.

Gear (1971b) gibt an, daß dadurch auch der Einfluß von Rundungsfehlern gedämpft wird. Das Programm DSTIFF hingegen verwendet einen Faktor von $1.025$ und dies, um schneller zu den höheren Ordnungen zu gelangen.

Dennoch zeigen die Testläufe von Gupta (1985), daß das Programm DSTIFF bei hohen Genauigkeitsanforderungen $\varepsilon$ besonders wenig erfolgreich arbeitet. Dies ist ein Hinweis auf die eingangs gemachte Erwähnung, daß höhere Ordnungen mehr Rechenbedarf erfordern und es sich damit anbietet die hohen Ordnungen so weit wie möglich zu meiden, zumindestens im steifen Falle. Das Programm DSTIFF ist maßgeblich für steife Differentialgleichungen gedacht und nicht schaltfähig.

Bibliographisch:

1. Gupta, Gopal K.: “Description and Evaluation of a Stiff ODE Code DSTIFF”, SIAM Journal on Scientific and Statistical Computing, Vol 6, No 4, October 1985, pp.939–950

6. Wie oben schon erwähnt, werden die vergangenen Werte für das neue Gitter durch Interpolation bestimmt. In der Absicht den Interpolationsfehler möglichst gering zu halten, wurde die folgende Einschränkung eingeführt: Die vergangenen Werte für das neue Gitter müssen im Intervall $[t_{m\ell-p+1}{\mskip 3mu} ,{\mskip 5mu} t_{m\ell+\ell}]$ liegen. Dieses Intervall enthält nun genau die vergangenen Werte des aktuellen Zykluses und die gerade berechneten Näherungen. Diese Einschränkung ergibt die Ungleichung

Die größte Schrittweite $\hat h$, die nun dies erfüllt, ist

Die Werte $\overline \eta_{p\hat p}$ sind nun in der nachstehenden Tabelle, wie schon weiter oben angekündigt, aufgelistet. Sie ergeben sich alle durch die oben angegebene Formel, wenn man dabei noch beachtet, daß die Stufenanzahl $\ell$ von der Ordnung $p$ abhängt. Die Abänderung der Unendlichkeitswerte auf 6 und 4$1\over2$ geschah, weil es sonst zu einem völligen Überschiessen kommen könnte. Dieses Verhalten konnte auch tatsächlich beobachtet werden. Dies gilt insbesondere für den Wert 6.

Während der Startphase sind die Schrittweiten u.U. für das gegebene Problem noch so “klein”, daß die Genauigkeit der Schätzung der Schrittweite nicht überbetont werden darf. Den beiden abgeänderten Unendlichkeitsstellen haftet aber dennoch eine gewisse Willkürlichkeit an. Werte entsprechender Größe, also kleiner 10, sind natürlich ebenfalls einsetzbar. Die Werte sind im Programm TENDLER nicht kritisch.

Es ist $\ell(p) = 3 + \lfloor p/5\rfloor$ ($1\le p\le 7$) und dies sogar für die Ordnungen $p=1$ und $p=2$. Das Verfahren für die Ordnung $p=1$ ist, wie vormals schon erwähnt, die dreimalige Wiederholung des impliziten Euler-Verfahrens. Für die Ordnung $p=2$ wird dreimal hintereinander die BDF2 verwendet. Bei diesen beiden Ordnungen handelt es sich, wenn man so will, nicht um “echte” Zyklen. Weil die Schrittweiten- und Ordnungs-Steuerung mindestens $p+3=k+3$ Werte für die Berechnung von $\left\Vert\nabla^{p+2}y_{m\ell+\ell}\right\Vert_w$ benötigt, muß die Zykluslänge für das Verfahren erster Ordnung mindestens 3 betragen, wenn man nur die Lösungswerte aus dem augenblicklichen Zyklus benützen will.

Diese Bemerkung gilt für alle Ordnungen, insbesondere auch für die beiden niedrigsten Ordnungen, $p=1$ und $p=2$. Diese Bedingung kann nun leicht dadurch erfüllt werden, in dem auch die Ordnungen 1 und 2 die Zykluslänge 3 haben. Bei diesen Ordnungen ist dann also natürlich jede Stufe gleich. Die Programme DIFJMT, STINT und TENDLER verwenden alle die völlig gleichen Formeln, insbesondere werden bei den zwei niedrigsten Ordnungen $p=1$ und $p=2$, dreimal hintereinander jeweils die BDF1 und die BDF2 benutzt.

7. Der Aufbau des Differenzenschemas ist

Die eingeklammerten Werte werden nur von der Schrittweiten- und Ordnungssteuerung benutzt, stehen dem Benutzer aber ebenfalls zur Verfügung, anders als in den Programmen STINT und DIFJMT. Bei diesen beiden letzten Programmen werden diese Werte nur indirekt bei einer Normberechnung bestimmt. In dem Programm TENDLER stehen diese Differenzen direkt als Vektoren dem Benutzer zur Verfügung. Würde man auf diese Vektoren verzichten, so würden $5n$ weniger Speicherzellen benötigt. Das Programm TENDLER würde damit dann echt weniger Speicherplatz beanspruchen als die beiden Programme STINT und DIFJMT und dies trotz erweiterter Funktionalität. Die Ursache hierfür liegt in einem anderen verwendeten Datenschema.

Bei Problemen mit impliziter Zeitvorgabe könnten diese hohen Differenzen u.U. von Vorteil sein.

Zur Erinnerung: Die Berechnung von $\nabla^iy_k$ benötigt die Werte $y_k, y_{k-1}, \ldots, y_{k-i}$, also $i+1$ Ordinatenwerte. Da die zyklischen Verfahren ab der Ordnung $p\ge3$ sowieso mindestens dreistufig sind, entfällt eine Sonderbehandlung dieser hohen Ordnungen. Lässt man die Bedingung, daß bei dem Verfahren erster Ordnung nur zykluseigene Werte benutzt werden dürfen, fallen, so kann man durchaus bei den beiden niedrigsten Ordnungen zweistufige Verfahren wählen. Das Programm TENDLER lässt sich hier durch Austausch nur einer einzigen Konstanten leicht abändern.

Jedoch zeigten Testläufe von Gaffney (1984), daß das Programm STINT bei der Beschränkung auf die beiden niedrigsten Ordnungen, also BDF1 und BDF2, sogar effizienter arbeitet als das Programm LSODE, welches die völlig gleichen Formeln bei diesen beiden niedrigsten Ordnungen benutzt, allerdings in Nordsieck-Darstellung. In dem Programm TENDLER wurde daher vorerst davon abgesehen diese Zyklenlängenverringerung durchzuführen. Auch speicherplatzmässig liesse sich durch eine solche Veränderung keinerlei Gewinn erzielen, es sei denn, man begrenzt die Höchstordnung auf 2 oder auf 4. Bei diesen beiden Ordnungen vollzieht sich ja der Sprung der Stufenanzahl $\ell$. Daß Programme, basierend auf gleichen Formeln, unterschiedliche Wirkungsgrade aufweisen, ist ja schon mehrfach erwähnt worden. Daß allerdings eine, wenn man so will, künstliche Zyklenverlängerung, hier von 1 auf 3 sich positiv auswirkt, ist bemerkenswert. Es wurde schon anderer Stelle darauf hingewiesen, daß diese Erscheinung nicht notwendig ihre kausale Begründung in einer Zyklenverlängerung findet. Jedoch wenn es einen Nachteil für längere Zyklen gäbe, so könnte dieser mögliche Nachteil dennoch in dem allgemeinen Wechselspiel der verschiedenen Segmente keine maßgebliche Rolle spielen. Auf das letzte kommt es aber gerade an.

Bibliographisch:

1. Gaffney, Parick W.: “A Performance Evaluation of Some FORTRAN Subroutines for the Solution of Stiff Oscillatory Ordinary Differential Equations”, ACM Transactions on Mathematical Software, Vo. 10, No. 1, March 1984, pp. 58-72

8. J.M. Tendler gewann seine zyklischen Formeln durch Optimieren des Widlund-Winkels $\alpha$ und der Widlund-Distanz $\delta$ bei der $A_\infty[\alpha]$- und $S_\infty[\delta]$-Stabilität, unter der Prämisse konstanter Schrittweite.

Die in füheren Abschnitten erwähnten Stabilitätseigenschaften, wie $A_\infty^0[\alpha]$- und $S_\infty^0[\delta]$-Stabilität, wurden unter der Voraussetzung gewonnen, daß die Schrittweite sich nicht erhöht oder erniedrigt. Wenn sich nun trotzdem die Schrittweite oder die Ordnung im Laufe der Integration verändert, so muß man u.U. damit rechnen, daß die gewünschten Stabilitätseigenschaften nicht mehr anzutreffen sind.

Gear/Tu/Watanabe (1974) beobachteten, daß bei starken Schwankungen der Schrittweite, die $D$-Stabilität, insbesondere bei höheren Ordnungen, verloren ging. Gear/Tu/Watanabe (1974) berichten über Schrittweitenänderungen jeweils immer um den Faktor 10 bei der BDF2. Wird also bei der BDF2 erst die Schrittweiteite verzehnfacht und dann sofort anschliessend wieder gezehntelt, also

so kann bewiesen werden, daß dieses so gebildete Verfahren nicht D-stabil ist. Genau das gleiche Resultat erhält man auch für die Adams-Formeln. Die $(k-1)$-fache Wiederholung des $k$-Schritt Adams-Verfahrens (mit konstanter Schrittweite) wirkt hingegen stabilisierend.

Mit Hinblick auf die obige Tabelle der Schrittweitenbegrenzer kann eine solche sehr extreme Änderung der Schrittweite in dem Programm TENDLER grundsätzlich nicht vorkommen. Auch in Programmen basierend auf Runge-Kutta-Verfahren werden solche starken Schwankungen von der dortigen Schrittweitensteuerung abgedämpft, man vgl. hier z.B. das Programm DOPRI5 von Hairer/Wanner/Nørsett (1987). Dennoch wird eine grundsätzliche Problematik krass beleuchtet.

Liniger (1979) schreibt hierzu:

The one-leg implementation of $A2$ is $A_0$-contractive with respect to $\left\Vert\cdot\right\Vert_\infty$ (and thus $A_0$-stable), relative to the variable coefficient test equation $\dot x=\lambda(t)x$, for any variable step sequence whatsoever. By contrast, any multistep method (which is not at the same time a one-leg method) can be destabilized by the use of certain variable step sequences when applied to the variable test equation $\dot x=\lambda(t)x$ with real valued $\lambda(t)$.

Es ist nun anzunehmen, daß diese Instabilitätserscheinungen auch für die zyklischen Formeln von Tendler (1973) gelten. Um diesen Schwierigkeiten zumindestens ansatzweise aus dem Wege zu gehen, wird verlangt, daß ein Zyklus immer mit konstanter Schrittweite durchlaufen wird (was in gewisser Hinsicht selbstverständlich sein muß) und ab den Ordnungen 4 und höher, mindestens zwei oder drei Zyklen mit fester Schrittweite benutzt werden.

Die genaue Anzahl der “Zwangs-Wartezyklen” für die Ordnung $p$ wird durch die Sperrvariable $\texttt{wait[}p\texttt{]}$ festgelegt.

Ein weiterer Punkt, dem Aufmerksamkeit gewidmet werden sollte bei der Schrittweiten- und Ordnungssteuerung in dem Programm TENDLER, ist, daß sowohl für den steifen als auch nicht-steifen Teil, verschiedene Schrittweiten- und Ordnungssteuerungsunterprogramme aufgerufen werden. Im steifen Teil werden niedrige Ordnungen sehr stark favorisiert, während im nicht-steifen Teil diese Tendenz abgeschwächt wurde, da bei nicht-steifen Differentialgleichungen hohe Ordnungen in dem Programm TENDLER wesentlich günstiger sind. Ein Austausch des Schrittweiten- und Ordnungssteuerungsunterprogramm und insbesondere die Hinzufügung weiterer Unterprogramme ist einfach (es handelt sich erneut um einen Vektor von Unterprogrammen). Das Programm LSODAR verwendet nur ein Schrittweiten- und Ordnungssteuerungssegment. Das Programm DEASY enthält natürlich 2 solcher Segmente: eines für das Runge-Kutta-Verfahren und eines für die BDF. Dies hat natürlich seine Ursachen darin, daß zwei ansonsten völlig selbstständige Programme zusammengefaßt wurden zu einem einzigen neuen.

Biographisch:

1. Werner Liniger (1927–2017)
2. Hairer, Ernst (*1949)
3. Wanner, Gerhard (*1942)
4. Nørsett, Syvert Paul
5. Gear, Charles William (1935—2022)
6. Tu, Kai-Wen
7. Watanable, Daniel S.

2. Die Schrittweitensteuerung im Programm TENDLER

The need to increase step size for efficiency and the simultaneous need to maintain strict error control run somewhat counter to each other. Yet, this is precisely the type of performance we look for in high quality ODE software.

G.D. Byrne, A.C. Hindmarsh (1987)

Vor einer allzu häufigen und allzu starken Änderung der Schrittweite sei jedoch gewarnt; hier sind noch nicht alle Fragen der Stabilität bei Schrittweitenänderungen befriedigend beantwortet worden.

H. Werner, H. Arndt (1986)

1. Unter der Annahme konstanter Schrittweiten werden die theoretischen Begründungen für dieses Segment abgeleitet. Hier wird beschrieben wie die Schrittweite bestimmt wird. Die Auswahl der Ordnung folgt später. Sie ist abhängig von der Schrittweitenwahl.

Berechnet werden nur die wirklich benötigten Werte, also diejenigen Werte, die auch wirklich an der Entscheidungsfindung für eine neue Ordnung und Schrittweite partizipieren.

Biographisch:

1. Helmut Werner (1931—1985)
2. Herbert Arndt

2. Mit der Variablen same berechnet man denjenigen Faktor, mit dem man die Schrittweite $h$ für den nächsten Rechenschritt multiplizieren könnte. Der Wert same wird immer berechnet und zwar nach der Rechenvorschrift:

3. Die Berechnung des Wertes down folgt der Formel:

4. Der Wert up wird nur dann berechnet, falls eine Ordnungserhöhung überhaupt möglich ist. Dies ist gegeben, falls

d.h., keine der Sperren ist gesetzt, die Maximalordnung $p_{\max}$ wird auch bei einer Ordnungserhöhung nicht überschritten, und es wäre möglich, auch bei der aktuellen Ordnung einen Schrittweitengewinn, bzw. zumindestens keinen Verlust zu erleiden. Nun wird up berechnet wie folgend

Für den nicht-steifen Teil werden weniger restriktive Anforderungen gestellt. Zudem ist sf[UP] für den nicht-steifen Teil größer. Die Matrizen stsaffac[][] und nssaffac[][] entsprechen jeweils $\sigma_{p,j}$. Im nicht-steifen Teil wird unter down, same und up schließlich der größte Wert herausgesucht, während hingegen im steifen Teil eine Tendenz über mehrere Ordnungen hinweg notwendig ist, damit die Ordnung erhöht werden kann.

5. Für zwei reelle Zahl $a$, $b$ ist $a\land b=\inf(a,b)$, also hier natürlich $a\land b=\min(a,b)$. Die Konstanten sf[SAME], sf[DOWN] und sf[UP] sind Sicherheitsfaktoren (sf: safety factor), wobei sf entweder stsaffac[][] oder nssaffac[][] ist.

Der Vektor errconst[] enthält die Konstanten für denjenigen Fehlertest, der zuerst der Schrittweiten- und Ordnungssteuerung nachfolgt.

6. Die hierbei mehrfach aufgetretene Funktion rsmooth() (ratio smooth) ist definiert durch

Diese Funktion ist nicht global-stetig. Man kann sie auffassen als nicht-normalisierte unscharfe Funktion (fuzzy function).

Die Konstanten DOWNTHRESH und UPTHRESH sind gegenwärtig auf die Werte $0.885$ und $1.16$ gesetzt. Das Programm TENDLER verhält sich nicht sensitiv ggü. diesen beiden Faktoren, jedoch führt eine geschickte Wahl dieser beiden Werte zu einer entsprechenden Verringerung der Rechenarbeit. Die Funktion rsmooth()

1. unterstützt die Vermeidung eines Schrittweitenwechsels im Bereich $\xi\in[\texttt{DOWNTHRESH},\texttt{UPTHRESH}]$,
2. reicht das Element $\xi$ unverändert weiter für Werte größer als UPTHRESH und begrenzt ansonsten das Absinken auf $\frac{1}{4}$.

Die Endbewertung des Sicherheitsfaktors ratio geschieht hauptsächlich über die Variablen stsaffac[][] und nssaffac[][].

7. Man kann den Prozeß der Schrittweiten- und Ordnungssteuerung auffassen als ein stochastisches, multikriterielles Optimierungsproblem. Es soll also eine optimale Strategie gefunden werden, jedoch sind die Auswirkungen dieser Entscheidung, dem Entscheidungsträger im voraus nicht bekannt. Die Optimierungsaufgabe ist multikriteriell, da mehrere, konfliktiöse Zielfunktionen gegeneinander konkurrieren. U.a. unter den Nebenbedingungen

sind die drei folgenden Zielfunktionen simultan zu optimieren:

Die augenblickliche Lösungsstrategie ist stark heuristischer Natur. Insbesondere garantiert sie nicht, das globale Optimum zu finden. Sie verspricht auch nicht per se schrittweitenwechsel-stabil zu sein und dies obwohl Überlegungen bzgl. der Schrittweitenwechselstabilität letztlich gewisse Entscheidungen beeinflussten. Unter veränderten Gesichtspunkten, nämlich kontrolltheoretischen Überlegungen, kommen Gustafsson/Lundh/Söderling (1988) zu einer alternativen Lösungsstrategie.

Von den genannten Autoren werden diese Versuche anhand des Programmes DOPRI5, siehe Hairer/Wanner/Nørsett (1987), durchgeführt.

Bibliographisch:

1. Gustafsson, K., Lundh, M. & Söderlind, G. API stepsize control for the numerical solution of ordinary differential equations. BIT 28, 270–287 (1988).
2. Kjell Gustafsson
3. Gustaf Söderling
4. Michael Lundh

3. Die Idee des Shampine-Testes

For example, such a user solving a set of equations describing a wind turbine got a stiffness indication from DE. The word “stiff” meant nothing to him, but the instruction to switch to another code did, and resulted in the effective solution of his problem. $\ldots$ The imperfection of the tests is that they may not detect stiffness; since they cost almost nothing, they do not harm and may prove useful.

L.F. Shampine (1977)1

1. Es soll hier kurz auf die grundlegende Idee des Shampine-Testes eingegangen werden. Die Stabilitätsbereiche von $P(EC)^iE$- und auch $P(EC)^i$-Verfahren für Picard Iteration zeigen deutlich, daß niedrigere Ordnungen einen wesentlich größeren Stabilitätsbereich aufweisen als höhere Ordnungen. Auf dieser einfachen Grundidee basiert letztlich die gesamte Überlegung. Thompson/Rodabaugh (1981) geben für die Adams-Formeln die folgenden Werte an, siehe Tabelle, die teilweise die Brauchbarkeit für schwach-steife Differentialgleichungen kennzeichnen. $P_i$ steht hier wie üblich für Prädiktor der Ordnung $i$, $E$ für Auswertung der Funktion und $C_j$ für Korrektur mit einer Formel der Ordnung $j$.

Mann erkennt deutlich (für die Adams-Formeln), daß die “Intervalllänge” des Stabilitätsgebietes schnell mit größer werdender Ordnung abnimmt. Ferner ist in diesem Falle die Kombination $P_kEC_{k+1}E$ am günstigsten. Wollte man eine noch weitere Vergrößerung des Stabilitätsintervalles vornehmen, so kann man durch Mittelung von Prädiktor und Korrektor für die Adams-Verfahren sogar eine bis zu zwölffache Verbesserung erzielen, siehe Thompson/Rodabaugh (1981). Man verwendet dann Verfahren der Form $C_k^\alpha = \alpha P_k + (1-\alpha) C_{k/k+1}$. Allerdings verändern sich dann auch die Fehlerkonstanten.

Bibliographisch:

1. D.J. Rodabaugh, S. Thompson: “A note on the relative efficiency of Adams methods” Computers & Mathematics with Applications, Vol. 7, Issue 5, pp. 401–403
2. Thompson, Skip
3. Rodabaugh, D.J.

2. Beginnt nun die Differentialgleichung steif zu werden, so erscheinen mehr und mehr die niedrigen Ordnungen, für die Schrittweiten- und Ordnungs-Steuerung “genauer” zu werden, während die höheren Ordnungen aufgrund der drohenden Instabilität zunehmend eine weniger glatte und stabile Lösung anzeigen. Bei nicht-steifen Regionen ist üblicherweise das Verhalten genau umgekehrt. Hier sind höhere Ordnungen bei höherer Genauigkeitsanforderung i.a. geeigneter.

Im Falle zunehmender Steifheit pendelt $H=h\lambda$ um den Rand des Stabilitätsbereiches. Diese Erscheinung lässt sich verbal wie folgt erklären. Eine Vergößerung der Schrittweite $h$ über den Rand des Stabilitätsbereiches hinaus, führt zu einem Ansteigen des lokalen Fehlers, welcher von der Fehlerkontrolle erkannt wird und durch ein Zurücknehmen der Schrittweite versucht wird entgegen zu treten und zwar so lange, bis $H$ wieder im Stabilitätsbereich liegt und sich damit der Ablauf wiederholt.

3. Diese verbale Beschreibung vereinfacht natürlich an manchen Stellen den wahren Sachverhalt. Zum einen wird man fragen, wie Steifheit auf den lokalen Fehler einwirkt und zum anderen wird eine Periodizität der Schwingung der Schrittweite um die Stabilitätsgrenze postuliert, aber nicht streng begründet. Denn, die Existenz und Einwirkung von Rückstellkräften belegt noch keine Schwingungs- oder Dämpfungserscheinungen. Den letzten Gedanken greift _{Hall, George}Hall (1985) und Hall (1986) auf. Dort wird ein Runge-Kutta-Verfahren untersucht auf Schwingungen um diese angedeutete Gleichgewichtslage. Die Analyse beruht auf der Untersuchung der nichtlinearen Rekurrenzgleichungen

wobei $E_i$ der vom Runge-Kutta-Verfahren geschätzte lokale Fehler ist und $H_n=\lambda h_n$. Diese Rekurrenzgleichung schwingt sich um einen Stabilisationspunkt ein, wenn sich eingeschleppte Störungen in erster Näherung wegdämpfen. Der Gleichgewichtspunkt liegt natürlich bei

Hier ist $\varepsilon$ die gewünschte Genauigkeitsanforderung, und $H_L$ liegt auf dem Rande des Stabilitätsgebietes. $\vartheta$ ist hier ein Sicherheitsfaktor, wie er in gängigen Programmen basierend auf Runge-Kutta-Verfahren üblicherweise benutzt wird. $S$ hängt vom Runge-Kutta Verfahren ab.

Mit dieser Analyse erhält man einen ersten Einblick über die Amplitude des Schwingens. Dämpfen sich die Störungen in erster Näherung nicht weg, ist also der Spektralradius einer gewissen Matrix größer als eins, so erwartet man heftige Ausschläge, und umgekehrt, falls sich Störungen schnell wegdämpfen, also der Spektralradius eins unterschreitet, so erwartet man wesentlich geringere Ausschläge um die Ruhelage. Diese Überlegungen wurden durch Hall (1985) und Hall (1986) empirisch bestätigt.

Sei $\rho$ der Spektralradius dieser hier nicht näher erörterten Matrix. Es wäre von großem Interesse die Eigenwerte der entsprechenden Matrix für die zyklischen Verfahren von J.M. Tendler genau zu kennen, insbesondere das Bild des Randes des Stabilitätsbereiches unter der Abbildung $\rho$, die sogenannte Hall-Kurve.

Bibliographisch:

1. Hall, George: “Equilibrium States of Runge-Kutta Shemes”, ACM TOMS, Vol 11, No 3, September 1985, pp.289–301
2. Hall, George: “Equilibrium States of Runge-Kutta Shemes: Part II”, ACM TOMS, Vol 12, No 3, September 1986, pp.183–192
3. George Hall

4. Der Shampine-Test kann von seiner grundsätzlichen Konzeption her ein unrichtiges Ergebnis liefern. Sei $A$ die Stabilitätsmenge für das Verfahren der höheren Ordnung, und sei $B$ die Stabilitätsmenge für das Verfahren der niedrigen Ordnung, und sei $C:=(A\setminus B)\cap\mathbb{C}^-$ die Differenzenmenge, also die Menge aller derjenigen Punkte aus $\mathbb{C}^-$, die zu $A$ aber nicht zu $B$ gehören. Wenn nun die Differenzenmenge $C$ nichtleer ist, so könnte die Schrittweiten- und Ordnungssteuerung eine Schrittweite $h$ vorschlagen, sodaß die Stauchung oder Streckung $H$ in $C$ liegt und damit das Verfahren der höheren Ordnung selbst im steifen Gebiete günstiger ist, als das Verfahren niedriger Ordnung.

Natürlich ist das einmalige Bestehen des Shampine-Testes zur globalen Klassifikation einer Differentialgleichung als steif oder nicht-steif, viel zu wenig. Ist der Shampine-Test jedoch 20-mal bis 50-mal hintereinander erfüllt, so ist dies ein sehr deutliches Zeichen für Steifheit. Mußte zudem ebenfalls die Ordnung auch noch runtergefahren werden und sind schon sehr “viele” Funktionsauswertungen bei niedrigen Ordnungen ($\le4$) durchgeführt worden, so hat man noch ein zusätzliches Anzeichen für Vorliegen von Steifheit.

5. Die gemachten Ausführungen lassen sich sehr gut durch Experimente bestätigen. Setzt man PMSSOC, so erhält man (analog wie bei \tracingpages=1 bei TeX) einen Ausdruck der Werte down, same, up, euler und bdf2, und man kann deutlich erkennen, daß wenn euler und bdf2 von unten gegen 1 rücken bei Vorliegen von Steifheit, daß dann auch die höheren Ordnungen rasch von oben gegen 1 streben und das Programm TENDLER bemüht ist die Ordnung zu senken. Zudem zeigt in dieser Situation das $QRP$, daß ein Schalten günstig ist. Den Ausdruck des $QRP_\nu$ und zahlreicher weiterer Größen, die während der Korrektoriteration anfallen, kann man sich durch Setzen von PMCVGCE ausgeben lassen.

4. Der Shampine-Test im Programm TENDLER

{Shampine, Lawrence F.}

The practical man frequently asks, “What is stiffnes?” There is no simple answer.

L.F. Shampine (1981a)

_{Shampine, Lawrence F.} Der Shampine-Test berechnet u.U. noch einen zusätzlichen Vergleichswert zu den Werten down, same und up. Die generierten Statistiken lassen darauf schliessen, daß dieser Test einen fast völlig vernachlässigen Einfluß auf das Schalten besitzt. Sehr wohl ist dieser Test dazu geeignet gewisse zusätzliche Informationen über die Differentialgleichungen zu erlangen. Dennoch, seine Berechnung kostet eine Normauswertung und ist damit recht aufwendig.

Die Berechnungen, die nun durchzuführen sind, gehen wie folgt vor sich. Der Wert bdf2 ist der Schrittweitenfaktor für die Ordnung 2, jedoch ohne einen Sicherheitsfaktor.

Für die Ordnung 1 wird der Wert euler berechnet nach der Vorschrift

Der Wert $\texttt{euler}$ ist der Schrittweitenfaktor für die Ordnung 1 und zwar wieder ohne Sicherheitsfaktoren.

Die Sicherheitsfaktoren wurden nicht berücksichtigt, um den Test restriktiv zu gestalten. Man beachte, daß ja die anderen Schrittweitenfaktoren wie down, same und up, sämtlich mit einem Sicherheitsfaktor belastet werden.

In einer früheren Version des Programmes TENDLER wurden beide Werte, somit bdf2 und auch euler immer ausgerechnet, falls dies möglich war, also bei den Ordnungen $p\ge4$. Da jedoch bei umfangreichen Tests immer

galt und daher die Berechnung von euler überflüssig war, wurde die Berechnungsvorschrift für den Shampine-Test dahingehend modifiziert, wie sie oben angegeben ist.

Der Shampine-Test gilt nun als bestanden, also die Variable shmptst ist somit TRUE, falls

Ein erstes Anzeichen für Steifheit liegt vor, falls der Shampine-Test nicht bestanden wurde. Geschieht dies “häufig”, oder wird über sehr lange Integrationsintervalle nur eine der niedrigen Ordnungen benutzt, so ist dies ein nachträgliches Anzeichen dafür, daß die zugrunde liegende Differentialgleichung in dem betrachteten Intervall wohl als steif zu bezeichnen ist.

5. Die Auswahl der Ordnung im Programm TENDLER

1. Grundsätzlich gilt für die Sicherheitsfaktoren $\sigma_{pj}$: Eine Bevorzugung der niedrigen Ordnung wird durch nahe bei 1 liegende Werte erzielt, eine Vermeidung der höheren Ordnungen durch nahe bei 0.5 liegende Werte. Einer Verringerung der Ordnung wird wesentlich größerer Zuspruch zuteil, als einer Ordnungserhöhung, weil den Verfahren niedrigerer Ordnung die folgenden beiden sehr erwünschten Eigenschaften zugute kommen:

1. bessere Stabilitätseigenschaften und
2. geringerer Rechenaufwand.

Dieser Effekt gilt insbesondere für die Schrittweiten- und Ordnungssteuerung im steifen Modus. Im nicht-steifen Modus ist der Trend zur Bevorzugung von Formeln niedriger Ordnung nicht so stark ausgeprägt.

Die Sicherheitsfaktoren $\sigma_{pj}$ für die Schrittweiten- und Ordnungssteuerung sind in den beiden Matrizen stsaffac[][] und nssaffac[][] abgelegt, und die Sicherheitsfaktoren für die Fehlerkontrolle sind im Vektor ecsaffac[] gespeichert. Man beachte, daß eine Veränderung der Ordnung, bei der Darstellung der Lösung in Form von rückwärtsgenommenen Differenzen, sehr einfach ist, im Gegensatz zu der Darstellung der Lösung in der Nordsieck-Darstellung. Diese Aussage gilt in entsprechender anderer Richtung für das Ändern der Schrittweite. In diesem Falle bietet die Nordsieck-Darstellung die größeren Vorteile.

2. Die eigentliche Ordnungsentscheidung sieht nun wie untenstehend aus.

Die vorgeschlagene Ordnung $p_{\rm neu}$ wird als neue Ordnung für den folgenden Zyklus nur dann ohne weitere Abfragen akzeptiert, falls der zu dieser Ordnung korrespondierende Schrittweitenfaktor factor, also einer der Größen down, same oder up, keine Schrittweitenverkleinerung vorschreibt, also $\texttt{factor}\ge1$ ist. Einer Schrittweitenverkleinerung, oder konstanter Schrittweite wird stets Rechnung getragen. Es wird an dieser Stelle angenommen, daß einer bevorstehenden Schrittzurückweisung durch die Fehlerkontrolle in den nächsten Stufen vorgebeugt werden kann.

Die generierten Statistiken deuten auf eine Bestätigung dieser Annahme hin: Die Schrittzurückweisungen sind rapide zurückgegangen. Bei einer Vergrößerung der Schrittweite wird die endgültige Akzeptanz davon abhängig gemacht, ob auch keine der möglichen Sperren, wie $\texttt{wait[}p\texttt{]}$ oder $\texttt{problem}$, aktiviert worden waren. Wird auch dieser letzte Test endgültig bestanden, so wird die vorgeschlagene Ordnung $p_{\rm neu}$ angenommen und hnext wird

Die weiter oben beschriebenen Sperrvariablen wait[] und problem, die schon mehrfach erwähnt wurden, sollen im folgenden genauer und ausführlicher erklärt werden.

6. Die beiden Sperren wait[] und problem im Programm TENDLER

Für die Sperre problem sind hier lediglich die folgenden Operationen zulässig:

Will man den Sperrmechanismus verfeinern oder verändern, so sind lediglich die drei obigen möglichen Klassenfunktionen zu modifizieren.

Für die Sperre wait[] sind die folgenden Funktionen die im Sinne von C++ einzig möglichen Klassenfunktionen:

1. SETWAIT(),
2. WAITCYCLE() und
3. die Abfrage WAITISOK()

Denkbar wäre die Änderung, daß man abgeht von der augenblicklich verwendeten dualwertigen Logik zu mehrfachwertiger Logik, oder sogar unscharfer Logik (fuzzy logic). Dies würde natürlich eine Veränderung der Argumentenlisten der beiden Funktionen ACTIVATE() und DEACTIVATE() erforderlich machen, die im Augenblick lediglich einstellig sind.

Benutzt werden diese Operationen von den Unterprogrammen bzw. Funktionen:

1. setconv(),
2. ecs(),
3. stssocs(),
4. nsssoc().

Die Funktion setconv() überwacht die erfolgreiche Konvergenz des Korrektors. Das Unterprogramm ecs() ist die Fehlerkontrolle (error control segment) und st/ns-ssocs() sind die Schrittweiten- und Ordnungssteuerungen (step-size and order control segment).

Ist die Ordnung $p>3$ und bleibt sie dort, so müssen mindestens ein oder zwei Zyklen mit konstanter Schrittweite durchgeführt werden. Die Fehlerkontrolle ecs() aktiviert mittels ACTIVATE() die Sperre problem und mittels SETWAIT() die Sperren wait[], bei einem Versagen des Fehlertests. Die Funktion setconv() schließlich aktiviert problem, falls wiederholte, langsame Konvergenz eine Refaktorisierung der Iterationsmatrix anzeigten. Eine dann möglicherweise direkt folgende Schrittweitenvergrößerung in der Schrittweiten- und Ordnungssteuerung wird durch das Aktivieren von problem unterbunden. Man behalte im Gesichtsfeld, daß die Iterationsmatrix $W=I-h\gamma J$, ja von der Schrittweite abhängig ist und somit eine Veränderung der Schrittweite $h$ die Iterationsmatrix “veralten” lässt.

Zusammenfassend hat man also die Abhängigkeiten:

und

7. Die Schrittweiten- und Ordnungssteuerung im Programm DE/STEP

Da die Theorie der Ordnungs- und Schrittweitenalgorithmen recht unvollständig ist, gibt es keine beste Vorgehensweise. Man muß die Fragmente der vorhandenen Theorie betrachten, die Folgerungen für den Rest des Codes, Experimente, die die Wirkung von Änderungen an den Algorithmen untersuchen, und schließlich die Erfahrung anderer Forscher.

L.F. Shampine, M.K. Gordon (1984)

1. In dem folgenden Abschnitt sollen kurz die grundlegenden Ideen und Algorithmen, die in dem Programm DE/STEP von L.F. Shampine und M.K. Gordon, während der Schrittweiten- und Ordnungssteuerung benutzt werden, erläutert werden, siehe Shampine/Gordon (1984). Zum Schluß werden Bemerkungen gegeben, zur Übernahme gewisser Teile aus diesem Segment, in andere Programme. Stellenweise erscheinen die in dem Programm DE/STEP zur Anwendung gekommenen Heuristiken genügend universell, sodaß Überlegungen über einen Einsatz dieser Strategien in anderen Programmen, lohnenswert sind. Da die Schrittweiten- und Ordnungssteuerung mit seinen Entscheidungen über die zu wählenden Schrittgrößen und Ordnungen, maßgeblichen Einfluß auf das Gesamtfunktionieren und die Gesamteffizienz des verwendeten Programmes haben, machen sich die Sorgfalt und Mühe, die man diesem Teil des Programmes angedeihen läßt, von selber mehr als bezahlt.

Das Programm DDASSL zur Lösung von differential-algebraischen Gleichungen der Form $F(t,y,\dot y)=0$, verwendet eine sehr ähnliche Startegie während der Schrittweiten- und Ordnungssteuerung, wie das Programm DE/STEP.

Insbesondere wird versucht die Schrittweite nur zu halbieren, zu belassen oder zu verdoppeln. Weiterhin sind die Darstellungen in beiden Programmen gleich: Beide Programme verwenden modifizierte rückwärtsgenommene Differenzen.

Dabei trifft das Programm DDASSL noch besondere Vorkehrungen, sodaß sich der führende Koeffizient der Formel nicht ändert. Diese Maßnahme ist bedeutsam für die Iterationsmatrix beim modifizierten Newton-Kantorovich Iterationsverfahren.

Anders als in den Programmen DIFSUB, EPISODE, GEAR, LSODE und LSODA wird hier zuerst eine geeignete Ordnung $p_{\rm neu}$ ausgewählt und erst anschliessend eine zu dieser Ordnung passende Schrittweite $h_{\rm neu}$. Das Programm DE/STEP arbeitet stets im $PECE$-Modus.

Es wird also kein Konvergenztest durchgeführt, und dies erspart zahlreiche Normauswertungen und zwar mindestens so viele, wie Schritte gemacht werden, häufig werden dadurch natürlich mehr Normauswertungen vermieden. Die Überlegungen bzgl. der Schrittweiten- und Ordnungssteuerung beruhen auf der Annahme, daß mit konstanter Schrittweite gerechnet wurde.

Biographisch:

1. Lawrence F. Shampine
2. Marylin K. Gordon

2. Der $PECE$-Zyklus wird nun wie folgt durchlaufen: Nach $P$ und $E$, wird geprüft, ob die Ordnung $p$ erniedrigt werden kann. Es wird an dieser Stelle nichts über eine Erhöhung der Ordnung oder einer Veränderung der Schrittweite entschieden, sondern lediglich eine Verminderung oder Beibehaltung der aktuellen Ordnung wird in Betracht gezogen. Die Regel zur Auswahl einer neuen Ordnung (höchstens vermindert oder belassen) ist:

Für eine Erniedrigung der Ordnung $p$ wird also verlangt, daß sowohl die um eins niedrigere Ordnung erkm1 günstiger ist, als auch zusätzlich, daß die um zwei tiefere Ordnung erkm2 ebenso noch besser ist, als die aktuelle Ordnung mit erk. Da bei der Ordnung $p=2$ eine solche Tendenz über zwei niedrigere Ordnungen nicht verfügbar ist, verlangt man, daß die um eins niedrige Ordnung erkm1 mindestens halb so klein, oder noch kleiner ist, als die aktuell benutzte Ordnung erk.

Nachdem nun die obigen Tests durchgeführt wurden und u.U. eine um eins kleinere Ordnung $p$ vorschlugen, wird der Fehlertest vollzogen. Man befindet sich hier immer noch hinter $P$ und $E$, aber noch vor $C$ und $E$. Der Fehlertest entscheidet, ob der Schritt wiederholt werden muß. Außer in Notfällen, wird hier lediglich die Schrittweite $h$ halbiert. Ein Notfall liegt vor, bei wiederholtem Nichtbestehen des Fehlertestes. Wird der Fehlertest nicht bestanden, so werden die Stufen $C$ und $E$ nicht durchlaufen.

Ist der augenblickliche Schritt durch die Fehlerkontrolle akzeptiert worden, so folgen nun $C$ und $E$; es wird nun also korrigiert und die Schlußauswertung der Funktion $f$ wird durchgeführt. Danach wird schnell geprüft, ob schon eine Entscheidung über eine Ordnungsreduktion gefällt wurde---w.o. angegeben. Die endgültige Entscheidung für eine neue Ordnung für den folgenden Schritt, geschieht nach der Regel

Die Ordnung $p$ wird also endgültig erniedrigt, falls dies vorher schon bei der Bestimmung von $\hat p_{\rm neu}$ feststand, oder falls die um eins niedrigere Ordnung erkm1 sowohl besser ist, als die aktuelle Ordnung erk, als auch zugleich besser ist als die um eins höhere Ordnung erkp1. Die Ordnung $p$ wird erhöht, falls die maximal mögliche Ordnung 12 noch nicht erreicht wurde und falls eine Ungleichungskette der Form

vorliegt. Außerdem müssen dann mindestens $(k+1)$ Schritte mit konstanter Schrittweite durchgeführt worden sein. Die Zählung der Schrittweite mit konstantem Betrag, wird mit Hilfe des Zählers ns vorgenommen.

Ist die Ordnung vorher gleich eins, so kann eine Tendenz über mehrere Ordnungen nicht überprüft werden und es wird dann verlangt, daß die um eins höhere Ordnung mindestens halb so klein oder noch kleiner ausfällt.

3. Ist nun die neue Ordnung $p$ aufgrund der obigen Überlegungen bestimmt worden, so wird jetzt die zu dieser Ordnung $p$ passende Schrittweite $h$ berechnet. Diese Berechnungen werden also ganz zum Schluß durchgeführt, also nach $PECE$. Als Hilfsgröße sei für das weitere erklärt der Schrittweitenfaktor $r$ mit

Die neue Schrittweite $h_{\rm neu}$ wird nun wie folgt berechnet:

Die Startphase wird verlassen wenn ein Schritt zurückgewiesen wurde, die Ordnung erniedrigt wurde, oder die maximal mögliche Ordnung erreicht wurde. Zur Vermeidung des Rechnens an der Grenze der Maschinengenauigkeit wird tatsächlich weitergerechnet mit der Schrittweite

$u$ ist hierbei die kleinstmögliche Zahl, die auf dem Rechner noch dargestellt werden kann, sodaß $1+u>1$.

DE/STEP ist eines der Programme, welches dem Rechnen nahe der Grenzgenauigkeit große Aufmerksamkeit widmet.

Kann die Schritweite im nächsten Schritt verdoppelt werden, fällt also der Schrittweitenfaktor $r$ größer als zwei aus, so wird die Schrittweite $h$ auch tatsächlich verdoppelt, jedoch sonst nicht weiter vergrößert. Man beachte, daß das Programm DE/STEP die Schrittweiten höchstens verzweifachen kann.

Gilt nicht $r>2$ und ist $\texttt{erk}\le{1\over2}\varepsilon$, man zielt hier also auf doppelt so genaue Ergebnisse, so wird die Schrittweite $h$ derart belassen, wie sie ist. Im anderen Falle wird sie schlimmstenfalls halbiert, aber mindestens jedoch um 10% herabgesetzt.

Die Schrittweitenänderungen sind in diesem Lichte recht dual: man verdoppelt, belässt oder halbiert die Schrittweite. Es sei nocheinmal betont, daß wenn die Schrittweite $h$ nicht verdoppelt wird, sie überhaupt in keiner anderen Weise vergrößert wird. Diese Entscheidung begünstigt konstante Schrittweitenfolgen, welche für lineare Mehrschrittverfahren, in der Darstellung mit rückwärtsgenommenen Differenzen, besonders vorteilhaft ist. Sehr wohl wird jedoch relativ häufig die Ordnung $p$ verändert. Ein Wechsel der Ordnung ist für ein lineares Mehrschrittverfahren basierend auf einer Differenzendarstellung äußerst einfach und äußerst rechengünstig.

4. Wie schon mehrfach geschehen, sei hier abermals darauf hingewiesen, daß die zahlreich auftretenden Minima, Maxima und Sicherheitsfaktoren wie $\varepsilon/2$, $9\over10$ und dergleichen, nicht etwa einem Beweise entspringen, sondern heuristischen Ursprunges sind. Diese Überlegungen basieren auf Erfahrungen und Plausibilitätsüberlegungen. Ihre Bestätigung erhalten sie durch gutes, sicheres und effizientes Funktionieren des resultierenden Programmes.

Es handelt sich allerdings nicht um ad hoc Maßnahmen, sondern um Überlegungen, zu denen es triftige Hinweise gibt. An sehr speziellen Testproblemen kann ihre Wirksamkeit dargelegt werden.

Wollte man die Schrittweiten- und Ordnungssteuerungs-Strategien, wie sie in dem Programm DE/STEP verwendet werden, w.o. beschrieben, übernehmen als Schrittweiten- und Ordnungssteuerungsstrategie für andere Programme, so ist zu beachten, daß eine so duale Schrittweitenauswahl, die Iterationsmatrix $W=I-\gamma J$ ebenfalls recht abrupt ändert. Es ist also so, daß eine glattere und langsamere Anpassung der Schrittweite $h$ an den Verlauf der Lösung, bei steifen Differentialgleichungen, bei denen die impliziten Formeln i.a. mit Hilfe eines modifizierten Newton-Kantorovich Iterationsverfahren gelöst werden, Vorteile haben kann, mit Hinblick auf die Iterationsmatrix $W$. Darum erscheint eine Verwendung der Nordsieck-Darstellung, wie sie bei dem Programm LSODA verwendet wird, eine günstige Wahl zu sein.

Zumindestens wäre der Hindmarsh-Test entscheidend abzuändern, wenn man eine Schrittweiten- und Ordnungssteuerungsstrategie wie in dem Programm DE/STEP, verwenden wollte. Allerdings zeigt das gute Funktionieren des Programmes DDASSL, daß der Aspekt der langsamen Änderung der Iterationsmatrix bzgl. der Schrittweite $h$, nicht von alles entscheidender Bedeutung ist.

Die anderen Überlegungen, die in dem Programm DE/STEP zum Tragen kommen, also Bestimmung der Ordnung vor der Bestimmung der Schrittweite und die Notwendigkeit einer Tendenz im Verhalten von Termen, sind überlegenswerte Ansatzpunkte auch für andere Programme. Das Programm TENDLER übernahm die Gedanken der Tendenz genau von dem Programm DE/STEP.

Auch die Möglichkeit der Verkleinerung der Schrittweite, selbst nach erfolgreichem Schritt, haben Parallelen zwischen dem Programm DE/STEP und TENDLER. Die Programme STINT, LSODE und LSODA können die Schrittweite nur dann überhaupt verkleinern, bei einem Fehlschlag oder bei mehrfachem Konvergenzversagen.

5. Um auch einen Überblick über die Größe der Argumentenliste des Programmes DE/STEP zu erhalten, seien diese Werte hier vermerkt:

Das Programm DE/STEP setzt sich also aus lediglich 4 Unterprogrammen zusammen, von denen das Unterprogramm de die Hauptlast der algorithmischen Leistung vollbringt. Wie allgemein üblich handelt sich bei allen hier überhaupt erwähnten Programmen um Unterprogrammsammlungen und nicht etwa um selbstständige, autonom ablauffähige Programme. Erst durch ein vom Benutzer für sein spezielles Problem angepasstes Aufrufprogramm, werden die erwähnten Unterprogrammsammlungen genutzt.