25th April 2026, 14 min read

TENDLER: 8. Ausblick

Fortsetzung der TENDLER Programmbeschreibung.

Inhalt.

1. Erweiterungen des Programmes TENDLER
2. Erweiterungen des Verständnisses einiger Programmsegmente
3. Literatur

Perhaps we will even one day find Pulitzer prizes awarded to computer programs.

Donald Ervin Knuth (1984)

Finally, we recognize that the last word on ODE methods has not been said.

George Dennis Byrne, Alan C. Hindmarsh (1987)

Hier sollen einige wenige Stichwörter geliefert werden für Erweiterungen des Programmes TENDLER, als auch für die Theorie, die hinter zahlreichen Entscheidungen steckt. Obwohl das Programm TENDLER zahlreiche Anwender in seiner Jetztform zufrieden stellen sollte, bleibt noch ein gewisser Freiraum für Ergänzungen doch erhalten. Einige von den möglichen Erweiterungen wurden schon bei der Entwicklung des Programmes TENDLER mit einbezogen. Sie sind damit relativ einfach nachzurüsten, während andere Ergänzungen hingegen, größere Änderungen und Modifikationen im Programm TENDLER nötigenfalls erforderlich machen.

1. Erweiterungen des Programmes TENDLER

The first user community sees the code as a means to an end while the other see it as an end in itself.

Peter E. Tischer, Gopal K. Gupta, A.J. Maeder (1986)

1. Augenblicklich verwendet das Programm TENDLER für die Schrittweiten- und Ordnungssteuerung, die Fehlerkontrolle und für das Prädiktorsegment, rückwärtsgenommene Differenzen und zwar in ihrer klassischen Form. Das Programm DE/STEP verwendet modifizierte rückwärtsgenommene Differenzen. Es werden also nicht die Werte $\nabla^0f_{n+1},\ldots,\nabla^{p+2}f_{n+1}$ abgespeichert, sondern

$$ h_{n+1}(h_{n+1}+h_n)\cdots(h_{n+1}+h_n+\ldots+h_{n+2-i}) f[t_n,t_{n-1},\ldots,t_{n-i+1}],\qquad i\gt 1. $$

Im Spezialfall äquidistanter Schrittweiten geht diese Darstellungsform über in die klassische Darstellung. Im Programm DE/STEP lässt sich für die Koeffizienten der Verfahren eine vergleichsweise einfache Rekursion herleiten, mit der man diese berechnen kann und dies selbst im Falle nicht-äquidistanter Knoten. Diese Darstellung hat Vorteile bei Schrittweitenwechseln.

Auch schon Tendler (1973) in seiner Dissertation weist darauf hin, daß dies eine möglicherweise denkenswürdige Erweiterung wäre. Für die zyklischen Formeln von Tendler ist jedoch auf Anhieb keine einfache und sehr schnell berechenbare Rekursion für die Koeffizienten bekannt. Im Gegensatz zu einstufigen, linearen Mehrschrittverfahren, wo die Koeffizienten des Verfahrens nach ausreichend vielen Schritten mit konstanter Schrittweite sich nicht mehr ändern, sind die Stufen bei zyklischen linearen Mehrschrittverfahren häufig nicht gleich. Demzufolge ändern sich also auch u.U. die in dem Programm verwendeten Koeffizienten. Nur wenn ausreichend viele Schritte mit konstanter Schrittweite durchgeführt wurden und wenn gleichzeitig immer der Zyklus gleicher Ordnung benutzt wurde, stellt sich eine gleichförmige Variation der Koeffizienten ein.

Man kann hier natürlich noch den besonders günstigen Sonderfall betrachten, bei dem die Verfahren niedriger Ordnung in den Verfahren höherer Ordnung eingebettet sind, wie dies beispielsweise bei den Adams-Moulton-Verfahren in der Darstellung rückwärtsgenommener Differenzen der Fall ist. Es scheint im Augenblick kein Programm basierend auf zyklischen Formeln zu geben, welches modifizierte rückwärtsgenommene Differenzen verwendet. Die Programme DIFJMT, STINT, TENDLER und ODIOUS benutzen sämtlich die klassischen rückwärtsgenommenen Differenzen.

Der Aufwand beim Schrittweitenwechsel (nicht beim Ordnungswechsel) beträgt ${\cal O}(3/2{\mskip 3mu}p^2)$ Vektoroperationen der Form

$x\gets x+y$ und
$x\gets x+ay$, $x,y\in\mathbb{R}^n$, und $a\in\mathbb{R}$.

Dieser Aufwand leitet sich daraus ab, daß $(p-1)$-mal interpoliert werden muß, mit je $(p+1)$ Vektoroperationen. Anschließend ist die Differenzentabelle neu aufzubauen, was ${1\over2}p(p+1)$ Vektoroperationen benötigt und schließlich muß $z=hf$ angepasst werden. Grob gesprochen kann man damit sagen, daß jeder der $p$ gespeicherten Werte, einen Aufwand von ${3\over2}p$ Vektoroperationen zur Anpassung an ein neues Gitter benötigt.

Dieser Aufwand wird umso höher, desto größer die Ordnung ist, was erneut eine Favorisierung der niedrigen Ordnungen nahe legt. Ist $p\le\ell$, d.h. ist die Ordnung $p$ kleiner oder gleich der Anzahl der Stufen $\ell$ innerhalb eines Zykluses, so entspricht der Aufwand für das Wechseln der Schrittweite in etwa dem Aufwand für das Aktualisieren der Differenzentabelle während des Durchlaufens eines gesamten Zykluses.

2. Für lineare Differentialgleichungen der Form $\dot y=Ay+f(t)$, scheint es angebracht zu sein, einen eigenen Treiber bereitzustellen. Differentialgleichungen dieser Form treten vergleichsweise häufig in den Anwendungen auf, und für einen gesonderten Treiber ist es sehr einfach, die durch die speziellen Eigenschaften sich darbietenden Einsparungen und Effizienzsteigerungen zu nutzen, während umgekehrt für einen völlig ungeübten Benutzer diese Einsparungsmöglichkeiten nicht auf der Hand liegen. Zum einen beträfe eine solche Erweiterung die Behandlung der Jacobimatrix $J$, welche hier konstant ist, und zum anderen beträfe dies die Tatsache, daß bei $P(EC)^i{E}$-Verfahren, die Differentialgleichung u.U. mehrfach an einem Gitterpunkt ausgewertet wird, jedoch nur alleine sich $y$ ändert, $t$ bleibt dabei fest. Diese Bemerkung gilt grundsätzlich. Der Benutzer sollte sich in einer static Variablen stets merken, ob die Funktion beim selben $t$ ausgewertet wird oder nicht. Diese Überprüfungen lohnen sich nur dann, wenn $f$, wie oben, “aufwendig” ist.

3. Neben einem symbolic-dump, kann es gelegentlich nützlich sein, von sämtlich wichtigen Größen im diffeqdescriptor eine Ausgabe in eine Datei zu erzeugen. Weiterhin wären dann Hilfsroutinen bereitzustellen, die es erlauben, dort weiter zu rechnen, wo vorher aufgehört wurde, z.B. mit Routinen ddsave() und ddrestore(). In größeren Rechenzentren, mit seinen zahlreichen up- und down-loads, wäre eine solche Erweiterung von Nutzen. Zugleich wäre zu überprüfen, ob damit einhergehend, diese Ausgabe derart gestaltet wird, daß folgendes Szenario leicht möglich wird: Man erledigt auf einem Vorrechner geeignete Vorbereitungsrechnungen und Nachberechnungen (wie z.B. Plots, etc.) und überlässt die eigentliche Rechenarbeit (number crunching) einem nachgeschalteten größeren Rechner. Die verwendeten Dateiformate hätten dann auf beiden Rechnern gleichartig auszusehen. Bei dem Programm LSODE und dessen Modifikationen muß man dies z.Z. ebenfalls selber programmieren.

4. Die Lösung großdimensionaler Differentialgleichungen erzwingt, daß die Routinen für die $LU$-Zerlegung angepasst werden müssen, um die i.d.R. vorliegenden Bandstrukturen auszunutzen. Dies wäre eine sehr vordringliche Aufgabe. In diesen Zusammenhang böten sich auch zahlreiche Erweiterungen an, wie z.B. Enright's update technique, auf die schon an anderer Stelle aufmerksam gemacht wurde. Diese Neuberechnungsstrategie für die Iterationsmatrix $W$ liesse sich im Vergleich zur Änderung der rückwärtsgenommenen Differenzen in modifizierte Differenzen, sehr leicht einbauen, da hier schon entsprechende Vorbereitungen vorhanden sind.

5. Ein äusserst nützlicher Zusatz wäre die Möglichkeit, Probleme mit impliziter Zeitvorgabe zu lösen, auf die ebenfalls schon an mehrfacher Stelle hingewiesen wurde. Aufgrund der Darstellung der Lösung und der Wechselwirkung mit anderen Teilen eines Differentialgleichungslösers, lässt sich keiner der gängigen Unterprogramme zur Lösung von Zeitvorgabeproblemen unmodifiziert übernehmen, sondern setzt gewisse Änderungen voraus. Da die rückwärtsgenommenen Differenzen hoher Ordnung schon berechnet sind, empfiehlt es sich diese auch zu nutzen. Die mit einer solchen Erweiterung verbundenen Schwierigkeiten werden nicht selten gründlich unterschätzt. Shampine/Hiebert (1980)1 äussern sich wie folgt hierzu:

Thus one can normally detect only an odd number of roots (counting multiplicities) in the course of a step. This is a fundamental limitation. It is popularly assumed that the step selected by the ODE solver will be appropriate for “seeing” roots. This is not true. $\ldots$
There has been relatively little attention given to the task described. As a result few codes provide such a capability, and the literature is sparse. $\ldots$
In our observation the difficulty of the task has been grossly underestimated by those wanting such a capability. There are also some significant software difficulties which are not generally appreciated.

Hairer/Wanner/Nørsett (1987) schlagen hier die Verwendung des Newton-Verfahrens vor, um innerhalb eines Schrittes nach einer Nullstelle zu suchen. Einfachste Beispiele zeigen, daß man hier zahlreiche Nullstellen ungerader Ordnung nicht finden kann, zumal man stets überprüfen muß, ob die Ableitung nicht gerade einen Nulldurchgang hat. L.R. Petzold in dem Programm LSODAR verwendet eine Kombination von Bisektion und regula falsi, wobei zahlreiche Sonderfälle separat behandelt werden, wie z.B. Nullstellen zu Beginn der Integration und dergleichen mehr.

6. Für viele kleindimensionale Differentialgleichungen ist die Rechenzeit für die Übersetzung und die Bindezeit länger als die eigentliche Rechenzeit für das Programm TENDLER. Die Funktionsauswertungen bestehen aus so wenigen arithmetischen Operationen, daß auch die Anzahl der Funktionsauswertungen beinahe völlig unerheblich ist. Gerade bei diesen Differentialgleichungen sind Runge-Kutta-Verfahren den Programmen basierend auf linearen Mehrschrittverfahren in der Rechenzeit überlegen, jedoch sind diese ohnehin unbedeutend. Dennoch ist der edit-compile-load Zyklus für den Benutzer zeitmässig i.d.R. nicht vernachlässigbar. Hier bietet es sich an, den Funktionsausdruck in symbolischer Form zu akzeptieren und den dann in postfix-Notation umgeformten Ausdruck weiter zu gebrauchen und nur mit diesem zu rechnen. Das Programm würde dann beispielsweise aus einer Datei die Funktion und wichtige Steuerungsgrößen, wie z.B. Anfangszeit, Endzeit, Statistikmodus, Druckmodus, Genauigkeitswunsch und dergleichen, lesen und sofort mit der Ausführung beginnen. Gerade für ein Programm basierend auf Mehrschrittverfahren würden die zusätzlichen Kosten bei den Funktionsauswertungen nicht ins Gewicht fallen. Weiter könnte man daran denken, die Parameterübergabe zahlreicher Funktionen nicht zu basieren auf numerischen Werten, sondern auf Zeichenketten\ifabzwick. \else, so wie jetzt schon bei der Funktion ssetmodes(). Desweiteren könnte eine weitere Benutzersicherung der Gestalt eingebaut werden, daß in dd eine calling-sequence-number integriert wird, zur Kontrolle der Aufrufreihenfolge.

7. Das Programm DE/STEP verwendet als Prädiktor ein Verfahren, welches eine um eins niedrigere Konsistenzordung hat als der Korrektor. Dies wird natürlich besonders nahe gelegt durch die unterschiedliche Anzahl der Startwerte der Adams-Bashforth und der Adams-Moulton Verfahren gleicher Konvergenzordnung. Jedoch zeigen die Stabilitätsbereiche insbesondere für kleinere Ordnungen hier Vorteile. Es scheint nicht bekannt zu sein, ob dies auch für die Verfahren von Tendler gilt.

Wäre der Stabilitätsbereich tatsächlich bedeutsam größer, so wäre es u.U. angeraten, die Prädiktorformeln auszutauschen. Wichtig ist der Stabilitätsbereich insbesondere für die Picard-Iteration. Die Stabilitätsbereiche ohne Berücksichtigung der Prädiktorformeln gehen davon aus, daß die Formeln ausreichend genau aufgelöst werden. Bei einem modifizierten Newton-Kantorovich Iterationsverfahren muß dies nicht immer erfüllt sein. Ein Einbau neuer Prädiktorformeln ist einfach möglich, wie bei der Beschreibung des Prädiktorsegments deutlich wurde. Genauso einfach und unkompliziert ist auch der Austausch von Korrektorformeln. Es wurde auch darauf aufmerksam gemacht, daß zyklische Formeln mit 4 oder sogar 5 Stufen nicht durchweg unbrauchbar sein müssen.

8. Die Schrittweite im Programm DE/STEP wird im großen und ganzen nur halbiert oder verdoppelt, Wie schon bei der Beschreibung der Schrittweiten- und Ordnungssteuerung des Programmes DE/STEP deutlich gemacht wurde, ist es u.U. sehr vorteilhaft die Schrittweite nur zu verdoppeln, zu halbieren, oder zu belassen, von den Sonderfällen bei mehrfachen Mißlingen eines Schrittes einmal abgesehen.

Es wäre von Interesse, dieses Vorgehen in das Programm TENDLER zu übernehmen. Die mögliche Gefahr dabei ist, daß sich die Iterationsmatrix $W$ recht abrupt ändert, bei einer solch dualen Änderung der Schrittweite. Demzufolge wäre der Hindmarsh-Test entsprechend zu modifizieren um zu vermeiden, daß jede Änderung der Schrittweite eine Refaktorisierung der Iterationsmatrix $W$ nach sich zieht. Hierzu allerdings müsste man wissen, ob die Konvergenzeigenschaften der Korrektoriteration unter einer stark veränderten Iterationsmatrix ungünstig beeinflusst werden, oder ob durch eine sonstwie höhere Genauigkeit der gespeicherten Werte hier nicht eine Kompensation eintritt.

Zu berücksichtigen ist jedoch, daß das Programm DDASSL durchaus mit dieser dualen Änderung der Schrittweite erfolgreich zu arbeiten in der Lage ist, wie z.B. Tests in Deuflhard/Hairer/Zugck (1987) deutlich machten. Allerdings macht der letztgenannte Artikel keine Aussagen über den globalen Fehler während der Tests.

Shampine (1980)2 §4, liefert erste Analysen bzgl. der Abhängigkeit der Korrektoriteration und der Abhängigkeit der Konvergenzeigenschaften von der Iterationsmatrix. Es bietet sich in diesem Zusammenhang an, einen Vektor von Hindmarsh-Funktionen zu benutzen, so, wie dies augenblicklich auch schon bei der Korrektoriteration gehandhabt wird. Zusätzlich wäre die Verwaltung der Ringliste zu modifizieren, da augenblicklich in ihr für so viele vergangene $y$-Werte und $z$-Werte kein Platz mehr frei ist.

Überlegenswert wäre hier, ob man gewisse Ringelemente, nämlich das letzte, nicht einer Sonderbehandlung unterzieht und nur die nachfolgenden Elemente überschreibt. Eine Verdopplung der Schrittweite ist bei den höheren Ordnungen nur möglich, wenn mehr als nur der letzte Zyklus gespeichert werden. Eine Halbierung der Schrittweite hat den positiven Nebeneffekt, daß der Aufwand für die Interpolation der Zwischergitterpunkte fast halbiert wird. Bei einer Vergrößerung der Schrittweite um das doppelte, braucht sogar überhaupt nicht interpoliert werden.

9. Einbau einer Option, die es erlaubt, stets die maximal mögliche Grenzgenauigkeit als Toleranz zu benutzen, z.B. eps=0. Gleichfalls Einbau einer Option, die beliebig viele Funktionsauswertungen gestattet, ohne, daß das Programm TENDLER mit einer Fehlermeldung die Integration abbricht, z.B. maxnfe=0.

2. Erweiterungen des Verständnisses einiger Programmsegmente

Diese Codes lassen noch viele Wünsche offen, und es werden noch eine Vielzahl vollkommen verschiedener Ansätze untersucht in der Hoffnung, einen Test zu finden, der befriedigender arbeitet.

Lawrence Fred Shampine, Marylin K. Gordon (1984)

1. Das Verhältnis von $\tau$ zu $\varepsilon$ wurde schon bei der Beschreibung des Konvergenztestes während der Korrektoriteration, als sehr einflußreich auf die Gesamteffizienz charakterisiert. Es wäre in hohem Maße wünschenswert hier genauere Analysen zur Verfügung zu haben. Diese wären selbst dann von Interesse, wenn sie lediglich für die Testgleichung $\dot y=\lambda y$ gelten würden.

2. Die Stabilitätsbereiche unter Berücksichtigung einer “zu alten” Iterationsmatrix, würden großen Aufschluß verschaffen, bei dem Verhalten des Lösers bei steifen Gleichungen. \ifabzwick\else Von Interesse wäre gleichfalls, wie schnell die Stabilitätsbereiche sich vergößern, bei sich veränderndem $r$ bei ${\cal Y}_r$.

3. Von großem Interesse ist weiter, ob die beliebige zyklische Kombination von $A$-stabilen Verfahren wiederum $A$ stabil ist. Bei $D$-Stabilität muß dies nicht grundsätzlich der Fall sein. Wichtiger ist aber noch, ob die Kombination $A_\infty^0[\alpha]$-stabiler, bzw. $S_\infty^0[\delta]$-stabiler Formeln der speziellen Gestalt wie bei den Formeln von Tendler, erneut noch $A_\infty^0[\alpha]$- bzw.
$S_\infty^0[\delta]$-stabil ist. Für die Ordnungssteuerung wäre diese Information hilfreich.

4. Wollte man tatsächlich neue zyklische, steif-stabile Formeln generieren, so könnte man auch daran denken, daß man innerhalb des Zykluses die Schrittanzahl der Stufen sukzessive erhöht. Tischer (1983)2 und Tischer/Sacks-Davis (1983)3 tuen dies. Allerdings wurde schon darauf aufmerksam gemacht, daß sich der Arbeitsaufwand deutlich erhöhen kann.

3. Literatur

Byrne, George Dennis
Byrne, George D. und Hindmarsh, Alan C.: “Stiff ODE Solvers: A Review of Current and Coming Attractions”, Journal of Computational Physics, Vol 70, No 1, May 1987, pp.5—62
Deuflhard, Peter Jochen (1944–2019)
Deuflhard, P. und Hairer, Ernst und Zugck, J.: “One-step and Extrapolation Methods for Differential-Algebraic Systems, Numerische Mathematik, Vol 51, 1987, pp.501–516
Enright, Wayne H.: “Improving the Efficiency of Matrix Operations in the Numerical Solution of Stiff Ordinary Differential Equations”, ACM TOMS, Vol 4, No 2, June 1978, pp.127–136
Gordon, Marylin K.
Hairer, Ernst (*1949)
Hairer, E., Wanner, G., Nørsett, S.: Solving Ordinary Differential Equations I – Nonstiff Problems. Springer Berlin, Heidelberg (2008), https://doi.org/10.1007/978-3-540-78862-1
Hindmarsh, Alan C. (*1942)
Knuth, Donald Ervin (*1938)
Knuth, Donald Ervin: “Literate Programming”, The Computer Journal, Vol 27, No 2, 1984, pp.97–111
Nørsett, Syvert Paul (1944–2025)
Pulitzer, Joseph (1847–1911)
Shampine, Lawrence Fred
Shampine, Lawrence Fred + Gordon, Marylin K.: Computer Solution of Ordinary Differential Equations: The Initial Value Problem, Freeman & Co Ltd, 1975
Shampine, Lawrence Fred und Hiebert, Kathie L.: “Detecting Stiffness with the Fehlberg (4,5) Formulas”, Computer and Mathematics with Applications, Vol 3, No 1, 1977, pp.41–46
Shampine, Lawrence Fred und Hiebert, Kathie L.: “Implicitly Defined Output Points For Solutions of ODEs”, Sandia report sand80-0180, Sandia National Laboratories, Albuquerque, New Mexiko (USA), February 1980, 24 S.
Shampine, Lawrence Fred: “Implementation of Implicit Formulas for the Solution of ODEs”, SIAM Journal on Scientific and Statistical Computing, Vol 1, No 1, March 1980, pp.103–118
STINT, source
Tendler, Joel Marvin (1943–2022)
Tendler, J.M.: A Stiffly Stable Integration Process Using Cyclic Composite Methods. Ph.D. thesis, Syracuse University, Syracuse, New York (1973)
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
Tischer, Peter E. und Gupta, Gopal K. und Maeder, A.J.: “A Software Engineering Approach for ODE Solvers”, in “Computational Techniques and Applications”, CTAC-85, Editor John Noye and Robert May, North Holland, Amsterdam New York Oxford Tokyo, 1986, pp.299–308
Wanner, Gerhard (*1942)