7th April 2026, 36 min read

TENDLER: 7. Start, Statistiken und Schalten

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
6. TENDLER: 6. Die Schrittweiten- und Ordnungssteuerung
7. TENDLER: 7. Start, Statistiken und Schalten

Inhalt.

• 1. Optionale Voreinstellungen
• 2. Die Startphase im Programm TENDLER
• 3. Die Statistiken im Programm TENDLER
• 4. Der Quotient Residuum Pseudoresiduum
• 5. Das Schalten im Programm TENDLER
• 6. Bewertung des Programmes TENDLER
• 7. Rechnung nahe der Maschinengenauigkeit

Mehrschrittverfahren benötigen, wie der Name andeutet, anders als Einschrittverfahren, wie z.B. Runge-Kutta-Verfahren, geeignete Startwerte. Diese Startwerte können natürlich nicht mit dem Mehrschrittverfahren selber bestimmt werden, sondern müssen anderweitig besorgt werden.

Eine Möglichkeit bestände darin, mit einem Runge-Kutta-Verfahren so lange zu rechnen, bis genügend viele Startwerte vorhanden sind und dann auf das Mehrschrittverfahren überzuwechseln. Bei steifen Differentialgleichungen sind jedoch blockimplizite oder semiimplizite Runge-Kutta-Verfahren vonnöten, welche bei sehr großdimensionalen Problemen den linearen Mehrschrittverfahren i.a. rechenzeitmässig unterlegen sind, siehe Enright/Hull/Lindberg (1975). Nebenher ist die Programmierung eines blockimpliziten Verfahrens für steife Gleichungen keine einfache Aufgabe. Einfacher ist, wenn mit einem Mehrschrittverfahren gestartet wird, welches einen Startwert weniger benötigt und man diesen Gedanken rekursiv fortsetzt. Man startet also mit dem impliziten Euler-Verfahren und steigert sich dann bis zum gewünschten Mehrschrittverfahren hoch. Während der Startprozedur braucht man die Schrittweite jedoch nicht konstant zu lassen, sondern wenn es geboten erscheint, verändert man sie. Da die Schrittweiten- und Ordnungssteuerung ehedem schon programmiert wurde, kann man sie auch gleich mitbenutzen. Diese Startstrategie wird nicht nur in dem Programm TENDLER so vollzogen, sondern auch in vielen weiteren gängigen Programmen basierend auf linearen Mehrschrittverfahren, wie z.B. DE/STEP, EPISODE, GEAR, LSODE und zahlreichen anderen Lösern.

Eine Besonderheit des Programmes TENDLER ist, daß jedes Segment mit statistischen Zählern versehen wurde. Diese Zähler messen mehr als nur die Anzahl der Funktionsauswertungen, der Jacobimatrixauswertungen, $LU$-Zerlegungen und Schritte, wie man es beispielsweise auch bei dem Programm LSODE findet. Die Art und Weise wie das Problem vom Programm TENDLER gelöst wird, erfährt eine genauere Durchleuchtung. Diese statistischen Daten sind für den Anwender aber auch für den Programmierer von Interesse. Wegen der zahlreichen Einstellungsmöglichkeiten, die das Programm TENDLER bietet, können, wenn man erst einmal einen Überblick über die Art und Weise der Lösungsgenerierung besitzt, diese Optionen vom Benutzer umsichtiger und gezielter eingesetzt werden. Der Programmierer kann hier Programmfragmente entdecken, welcher Beachtung bedürfen und Programmfragmente, welcher keiner besonderen Beachtung bedürfen. Z.B. wäre die Aussage, daß der Shampine-Test niemals anzieht, für die Wartung dieses Programmsegmentes von besonderem Interesse. Oder die Aussage wie sich die Schritte auf die Ordnungen verteilen, oder welche Stufe innerhalb eines Zykluses die meisten Fehlertest- oder Konvergenzversagen auf sich vereinigt. Aussagen dieser Art lassen sich mit den im Programm TENDLER eingebauten Statistiken machen. Sie wurden schon mehrfach erwähnt.

Beim Start werden zahlreiche wichtige Integrationsgrößen voreingestellt. Viele dieser Größen können vom Benutzer jederzeit verändert werden, oder gleich zu Anfang anders besetzt werden. I.d.R. besteht kein Anlaß hier Abweichungen vom Üblichen vorzunehmen, zumal dies u.U. ein erweiteres Verständnis voraussetzt. Es können jedoch besondere Situationen hier Anlaß geben um gewisse Voreinstellungen anders zu setzen, z.B. wenn man seine eigenen Treiber schreiben will.

1. Optionale Voreinstellungen

1. Mit dem Unterprogramm setstep() kann die minimale Schrittweite $h_{\min}$, die Anfangsschrittweite $h_0$ und die Maximalschrittweite $h_{\max}$ eingestellt werden, jedoch muß dies vom Benutzer nicht extra durchgeführt werden. Ruft er dieses Unterprogramm nicht auf, so werden Voreinstellungen wirksam. Die Einstellung der minimalen, maximalen und anfänglichen Schrittweite ist natürlich auch anders möglich und muß nicht notwendig über dieses Unterprogramm erfolgen, jedoch stellt das Unterprogramm setstep() eine komfortable Benutzerschnittstelle dar, die zudem die Werte einem rigorosen Test unterzieht.

Ähnlich verhält es sich mit dem Unterprogramm setmodes(), welches für die Einstellung dient

1. der Art und Weise der Benutzung der Jacobimatrix $J$,
2. der Iterationsart,
3. der Druckkontrolle und
4. der Art der Fehlerkontrolle.

Das Unterprogramm setaccuracy() bietet eine einfache Schnittstelle zur Einstellung der Toleranz, der Maximalordnung und der maximal erlaubten Anzahl der Funktionsauswertungen.

2. Die Unterprogramme setaccuracy(), setstep() und setmodes() müssen nicht aufgerufen werden. Alle Werte sind mit häufig auftretenden und der Differentialgleichung angepaßten Werten vorbesetzt. Beispielsweise berücksichtigt die voreingestellte Toleranz die Dimension der Differentialgleichung. Weiter ist es nicht nötig alle Parameter der Unterprogramme zu spezifizieren.

Der Wert AUTO_SELECT sorgt dann dafür, daß die Voreinstellungen aktiv werden. Zulässig ist beispielsweise der Aufruf $$ setmodes (dd, AUTO_SELECT, SWITCHING, PMSSOC|PMDUMP, AUTO_SELECT), $$ bei dem also der Wert der Variablen jacmode absichtlich unspezifiziert bleibt.

2. Die Startphase im Programm TENDLER

Ganz zu Beginn der Integration startet das Programm TENDLER mit der niedrigst möglichen Ordnung: mit der allerersten. Hat der Benutzer gute Kenntnisse über den Verlauf der Lösung “vor” dem Integrationsanfang, so kann das Programm TENDLER auch von einer anderen Ordnung aus gestartet werden. Dieser Quereinstieg ist jedoch, abgesehen von einer Fortsetzung einer schon früher durchgeführten Integration, nicht die Regel. In dieser Form wird er daher auch gegenwärtig nicht gesondert unterstützt.

Der übliche Start verläuft wie folgt. Die Ordnung $p$ wird auf $p=1$ gesetzt, die Statistikzähler werden ggf. gelöscht, der Gewichtsvektor bei der Normberechnung wird initialisiert, und die Startschrittweite $h_0$ wird berechnet nach der Formel

Hierbei bezeichnet $\left\Vert\cdot\right\Vert$ die zu $\left\Vert\cdot\right\Vert_w$ korrespondierende ungewichtete Norm, die vom Benutzer beliebig geändert werden kann und sogar zur sonst verwendeten Norm keinerlei Bezug haben muß. Die beiden reellen Werte $a,b$ sind die Integrationsgrenzen: Integriert wird von $a$ nach $b$. Entsprechend wird direction gesetzt.

In die Bestimmung der Anfangsschrittweite gehen die Integrationsgrenzen ein. In die Schaltentscheidung gehen sie nicht ein. Für die Berechnung einer geeigneten Anfangsschrittweite gibt es sehr ausgefeilte Techniken und Strategien, siehe Gladwell/Shampine/Brankin (1987). Die im Programm TENDLER verwendete Rechenvorschrift für die Bestimmung einer Anfangsschrittweite, w.o. beschrieben, ist nicht besonders ausgefeilt, dafür aber sehr einfach zu errechnen.

Die minimale Schrittweite $h_{\min}$ wird, wenn dies nicht vom Benutzer anders gesetzt wird, vorbesetzt mit dem Wert

und die maximal mögliche Schrittweite wird gesetzt zu

falls dies möglich ist, also falls $\mathopen|b-a\mathclose| > 8u\max(\mathopen|a\mathclose|,\mathopen|b\mathclose|)$. Hierbei ist $u$ die kleinste positive, auf der Rechenanlage darstellbare Zahl, sodaß gilt $1+u>1$.

Eine Fülle von Parametern sind weiterhin ganz am Anfang der Integration zu setzen. Hierzu gehören:

1. hused,
2. nqused,
3. inhibitJ,
4. itmode,
5. mnicvr,
6. fixcvr,
7. canstc,
8. cslow,
9. cquick,
10. atlst2,
11. shmptst,
12. h,
13. problem,
14. wait,
15. hfr,
16. etc.

Bibliographisch:

1. Gladwell, Ian + Shampine, Lawrence F. + Brankin, R.W.: “Automatic selection of the initial step size for an ODE solver”, Journal of Computational and Applied Mathematics, Vol 18, 1987, pp.175–192

3. Die Statistiken im Programm TENDLER

At first we were reluctant to incorporate such a capability, but NASA wanted the flexibility. $\ldots$ Our first reaction was to want a sign that would come up and say “your warranty is void” when they used the capability. $\ldots$ On one flight an overflow occurred in a counter that was keeping track of the number of Reaction Control Jet firings.

Tony (Anthony J.) Macina (1984)

Insgesamt 79 Größen, nicht berücksichtigt die Statistikvektoren, dienen der Protokollierung und Erfassung der Effizienz und Nutzungsart des Programmes TENDLER. Sie dienen weiterhin zur Erkennung von Effizienzreserven und zur Vermeidung der Programmieroptimierung an Stellen, die von geringer Bedeutung sind. Hierzu dient u.a. auch der Profiler. Profiler und Statistiken ergänzen sich gegenseitig. Während die Statistiken problemorientiert sind, ist der Profiler natürlich programmorientiert. Die statistischen Zähler sind dem Benutzer jederzeit direkt zugänglich. Zur Aufbereitung der statistischen Zähler dient das Unterprogramm prtstat(). Dieses Unterprogramm druckt in Abhängigkeit der Druckeinstellung nur diejenigen Werte, die maßgeblich von Interesse sind. Zudem werden die statistischen Daten miteinander in Beziehung gesetzt um Querverbindungen aufzuzeigen. Überdies werden auffallende Ergebnisse besonders hervorgehoben. Andererseits werden u.U. andere Daten gar nicht ausgegeben, was beispielsweise bei zweifachen Korrektorversagen schnell einleuchtet, was ja eben nicht immer auftritt. Sonderbehandlung findet auch der Integrationsanfang. Zur Interdependenz verschiedener Zähler: Beispielsweise stehen die Anzahl der erfolgreichen Schritte mit der Anzahl der zurückgewiesenen Schritte, dem Konvergenzversagen und dem Fehlertestversagen in enger Wechselwirkung. Für jeden Punkt existieren eine Reihe von Zählern. Beim Konvergenzversagen sind es über 10 Stück.

Gemeinsam ist den Routinen für die Statistik, daß sie vergleichsweise aufwendig zu programmieren sind, aber, wie profile-listings zeigen, nur marginalen Einfluß auf die Gesamtrechenzeit haben. Z.B. ist das Unterprogramm prtstat() das längste Unterprogramm überhaupt im Programm TENDLER: es ist über 500 Zeilen lang. Dies ist insoweit besonders hervorhebenswert, weil anders als man anfangs denken könnte, es mit Hochzählung von Registern nicht getan ist. Beispielsweise muß für die Variation der Schrittweite mehrfache Genauigkeit nachgebildet werden, da sonst leicht der Bereich der darstellbaren Gleitkommazahlen überschritten (overflow) oder unterschritten (underflow) wird. Ferner sind gewisse Register herunterzuzählen, falls andere Register hochgezählt werden, beispielsweise mni1diverg und mni2diverg. Dies liegt u.a. daran, weil prtstat() nichts darüber weiß, von wo es, wie, aufgerufen wird. Bei nicht so großdimensionalen Differentialgleichungen zieht jedoch die Funktion prtstat() aufgrund der zahlreichen Druck- und Formatieranweisungen prozentual einen größeren Teil der Gesamtrechenzeit auf sich, allerdings weniger als die Normauswertungen und weniger als die Funktionsauswertungen. Bei sehr kleindimensionalen Differentialgleichungen, oder sehr einfach auszuwertenden Funktionen kann sich dieses Verhältnis stark zu ungunsten der Funktion prtstat() verschieben.

Wichtig ist zu vermerken, daß prtstat() die schon generierten Statistiken im großen und ganzen optisch aufbereitet und dann nur ausgibt. Diese Funktionalität kann durch Setzen von PMSTAT in prtmode ausgeschaltet werden. Es werden dann keine Statistiken ausgedruckt. Erzeugt werden die Statistiken allerdings weiterhin noch. Die Erzeugung kann nicht unterdrückt werden. Die Unterdrückung der Erzeugung ist allerdings auch gar nicht nötig, da die Rechenzeiten hierfür deutlich unter 1% liegen, bei typischen Differentialgleichungsproblemen. Auch dies läßt sich mit einem Profiler belegen.

4. Der Quotient Residuum Pseudoresiduum

Eine der weiteren besonderen Auffälligkeiten des Programmes TENDLER ist die Fähigkeit, automatisch zu erkennen, ob die Differentialgleichung an der aktuellen Stelle nun als steif oder als nicht-steif zu betrachten ist und um dann anschliessend mit dem für diesen Falle geeigneten Verfahren fortzufahren.

U.U. unabhängig von ihrer ursprünglichen Bedeutung seien die folgenden Begriffe eingeführt.

1. Definition: (1) Sei $g\colon\mathbb{C}^n\to\mathbb{C}^n$ eine stetige Funktion und $a\in\mathbb{C}^n$ beliebig. Dann heißt die Folge $(x^i)_{i=0}^\infty$ gebildet gemäß

die zu $[g,a]$ gehörende Picard Folge oder auch Picard Iteration.

(2) Der Wert $s^i=x^{i+1}-x^i$ heißt $i$-tes Pseudoresiduum zur $[g,a]$ gehörenden Picard Folge. Sind Mißverständnisse ausgeschlossen, so heißt auch die Zahl $\mathopen|s^i\mathclose| = \mathopen|x^{i+1}-x^i\mathclose|$ ebenfalls Pseudoresiduum.

Das Pseudoresiduum ist also nichts anderes als die Differenz aufeinanderfolgender Folgenglieder. Die Betonung auf $g\colon\mathbb{C}^n\to\mathbb{C}^n$ ist kein Zufall, da Programme wie TENDLER, LSODE, etc., z.Z. während der Korrektoriteration nicht überprüfen, ob sie den Definitionsbereich nicht schon verlassen haben. Im Programm TENDLER ist diese Eigenschaft insoweit schon berücksichtigt, als daß die Auswertevorschrift für die Funktion $f$ der Differentialgleichung $\dot y=f(t,y)$ eine Statusanzeige zurückliefern kann, ob der Definitionsbereich verlassen wurde, oder ob der zulässige Darstellungsbereich für Gleitkommazahlen überschritten wurde (under- bzw. overflow). In einer früheren Version des Programmes TENDLER wurde hiervon auch tatsächlich Gebrauch gemacht. Jedoch in der gegenwärtigen Fassung wurden aus Effizienz- und Übersichtlichkeitsgründen diese zusätzlichen Überprüfungen wieder entfernt.

2. Definition: (1) Sei $f\colon\mathbb{C}^n\to\mathbb{C}^n$ eine stetige Funktion, $(W_i)_{i=0}^\infty $ eine Folge invertierbarer Matrizen und $a\in\mathbb{C}^n$ beliebig. Dann heißt die Folge $(x^i)_{i=0}^\infty$ gebildet gemäß

die zu $[f,(W_i)_{i=0}^\infty,a]$ gehörende Newton-artige Folge oder auch Newton-artige Iteration.

(2) Gilt $W_i\equiv W$, so heißt $(x^i)_{i=0}^\infty$ w.o., die zu $[f,W,a]$ gehörende Newton-Kantorovich-artige Folge oder Newton-Kantorovich-artige Iteration.

(3) Der Wert $r^i=f(x^i)$ einer Newton-artigen oder Newton-Kantorovich-artigen Iteration heißt $i$-tes Residuum. Falls Zweifel nicht bestehen, nennt man auch $\mathopen|r^i\mathclose|=\mathopen|f(x^i)\mathclose|$ Residuum.

(4) Der Wert $s^i=x^{i+1}-x^i$ einer Newton-artigen oder Newton-Kantorovich-artigen Iteration heißt Pseudoresiduum, analog aber ebenfalls $\mathopen|s^i\mathclose|=\mathopen|x^{i+1}-x^i\mathclose|$.

(5) Die zu $[f,(W_i)_{i=0}^\infty,a]$ gehörende Newton-artige Folge $(x^i)_{i=0}^\infty$ hat den Wert

und den Wert

QRP erinnert an: Quotient Residuum Pseudoresiduum. Genauer: zu $[f,(W_i)_{i=0}^\infty,a]$ gehörende $\text{QRP}_i$ bzw. $\text{QRP}$.

Die Werte $\text{QRP}_i$ kann man auch schreiben als

Es gilt $x^{i+1}-x^i\ne0 \iff f(x^i)\ne0$, da $\det W_i \ne 0$.

3. Lemma: Für jede Matrixnorm $|\cdot|$ verträglich zur Vektornorm $|\cdot|$ gilt

Beweis: Die zweite Ungleichung ist klar wegen $\mathopen|Ax\mathclose| \le \mathopen|A\mathclose| \mathopen|x\mathclose|$. Sei $x\ne0$ und $Ax\ne0$, wegen $\det A\ne0$. Dann ist

    ☐

4. Allgemein gilt

5. Aus dem Lemma folgt für $\text{QRP}_i \ne 0$ sofort

und genauso

Sei $W_i=A_i+B_i$. Dann gilt für $x\ne0$,

Ist $x=x^{i+1}-x^i\ne0$ und $\det A_i\ne0$, so gilt

Im besonders wichtigen Spezialfall $A_i\equiv I$ folgt also $\text{QRP}_i \ge 1 - \mathopen|B_i\mathclose|$.

5. Das Schalten im Programm TENDLER

Unfortunately, a change to an existing code which is simple in concept is often suprisingly difficult in practise. Providing an alternative iteration method is only part of the task. To properly solve nonstiff problems, some basic tactics, such as those for adjustment of step size, must be different from those for solving stiff problems. Providing an alternative set of tactics implies a substantial software development effort.

L.F. Shampine (1981)

1. Die obigen Überlegungen bzgl. des QRP werden nun angewendet bei der Schaltentscheidung, nämlich dem Wechsel der Iterationsart. Das Programm TENDLER wechselt beim Übergang vom steifen in den nicht-steifen Teil nicht die Formeln, sondern die Art, wie die Implizitheit der Formeln aufgelöst wird. Das Schalten ist also ein Wechseln von der modifizierten Newton-Kantorovich Iteration zur Picard Iteration.

Zu lösen ist

Der Subskript $m+1$ falle im weiteren zur Schreibvereinfachung weg. Die Picard Iteration zur obigen Gleichung lautet

wobei $y^0$, der Startwert der Iteration, durch einen getrennten Rechengang ermittelt wird. Die Newton-Kantorovich Iteration für obige implizite Gleichung lautet

mit

Hinreichende Konvergenzbedingung für die Picard Iteration ist $\mathopen|h\gamma J\mathclose| < 1$. Ist diese Bedingung erfüllt, so ist $W$ invertierbar und der Wert $\text{QRP}_\nu \ne 0$ der zu $[g,W,y^0]$ gehörenden Newton-Kantorovich-artigen Iteration erfüllt

Es ist $(*)$ eine notwendige Bedingung für $\mathopen|h\gamma J\mathclose|<1$. Allerdings folgt aus $(*)$ nicht immer $\mathopen|h\gamma J\mathclose|<1$, wie man z.B. erkennt anhand von

2. Im Programm TENDLER wird verlangt, daß $\text{QRP}_\nu\le1$. Getestet wird dies in der Form

Bei der eigentlichen Programmierung wird die obige Quotientenbildung umgangen. Der Zähler entspricht aber nun formal genau dem Pseudoresiduum der Picard Iteration, welche sich ja schreiben läßt als

Nocheinmal sei betont, daß $\text{QRP}_\nu\le2$ nur eine notwendige Bedingung ist für $\mathopen|h\gamma J\mathclose| < 1$. Wenn nun trotzdem die “unerlaubte” Rückrichtung benutzt wird, muß diese durch weitere Bedingungen zumindestens erhärtet werden. Dies ist im Programm TENDLER der Fall, allerdings nicht beim Programm NTI/SIMPLE.

Hinzugezogen zur Schaltentscheidung werden:

1. wie lange schon nicht geschaltet wurde,
2. der Shampine-Test,
3. die Sperrvariable problem und
4. die Scheinkonvergenzrate der Picard Iteration.

Obwohl also die Schrittweiten- und Ordnungssteuerung keine Schaltentscheidung alleine herbeiführen kann, beeinflußt sie diese jedoch.

Für die Iterierten

berechnet man den Wert $\text{QRP}_\nu$ über

und für die Iterierten

würde man sinngemäß rechnen

3. Der Indikator alwnst (always nonstiff) wird während der Newton-Kantorovich-Iteration auf

gesetzt. Vor jedem Zyklus wird diese Variable initialisiert mit

Hierbei ist $f_\nu$ die Scheinkonvergenzrate der Picard Iteration. Als Startwert für die Schein-Picard-Iteration nimmt man den Wert, der von der modifizierten Newton-Kantorovich Iteration geliefert wird. $\delta_\nu$ ist die Norm des Pseudoresiduums bei der Newton-Kantorovich Iteration, und $r_\nu$ ist die Norm des Pseudoresiduums der Schein Picard Iteration. Auch hier wird wieder der Startwert für die Schein Picard Iteration von der Newton-Kantorovich Iteration geliefert. Ein Wechsel vom steifen Teil in den nicht-steifen Teil kann als ein Wechsel von Integrationsformeln gedeutet werden, dessen Stabilität getrennt untersucht werden muß. Man denke in diesem Falle daran, daß die Stabilitätbereiche von Picard $P(EC)^i\{E\}$-Verfahren abhängig sind von der Anzahl der durchgeführten Iterationen. Aufgrund dieser Überlegung sollte ein Wechsel von steif nach nicht-steif nur durchgeführt werden, wenn keinerlei Bedenken sonstiger Art hiergegen einzuwenden wären, wie z.B. wiederholt aufgetretene Konvergenzschwierigkeiten, Minuspunkte, Fehlertestversagen, u.ä. Gleichzeitig wird verlangt, daß die Schein Picard Iteration eine günstig kleine Konvergenzrate aufweist, in der Hoffnung, daß die wirkliche Picard Iteration sich wohl ähnlich verhält.

Die eigentliche Schaltentscheidung für den nächsten Zyklus läßt sich programmnah wie untenstehend beschreiben. Diese Entscheidung wird ganz zum Schluß des Zykluses, also nach der Schrittweiten- und Ordnungssteuerung durchgeführt.

switch (itmode) {
    case MNI:      ...
    case MNIFLOAT:
        ...
        if (alwnst && "Schalten erlaubt") {
            // zähle Anzahl der offensichtlich nicht-steifen Zyklen
            if ("mehr als" MAXWAIT "offensichtlich nicht-steife Zyklen) {
                // wechsele in den nicht-steifen Modus über
                "nicht-steifer Zyklus registriert"
            }
        } else
            // steifer Zyklus registriert
        break;
    case FIXPECE:  ...
    case FIX:      ...
    case FIXFLOAT: ...
    default: ...
}

Geschaltet wird von steif nach nicht-steif also wenn:

1. alwnst erfüllt ist, was heißt, daß bisher nicht-steifes Verhalten während der und nach der Korrektoriteration aufgetreten ist und das Schalten überhaupt erlaubt ist;
2. die Anzahl der Zyklen, bei denen die Bedingung (1) erfüllt ist, mehr als MAXWAIT-mal aufgetreten ist.

Geschaltet von nicht-steif nach steif wird, falls

1. `itmode = FIXFLOAT ist und gleichzeitig
2. Divergenz der Picard-Iteration vom Konvergenztest gemeldet wurde.

Der Konvergenztest der Picard-Iteration enthält zusätzliche Maßnahmen, sodaß nach einer gewissen Anzahl von Stufen, mindestens eine Iteration mit zwei oder mehr Iterationsschritten durchgeführt wird. Dies dient dazu, daß die Konvergenzrate erneut berechnet werden kann und nicht lediglich geschätzt werden muß.

4. Zur Schaltentscheidung könnte man auch direkt die Norm der Jacobimatrix, also $\mathopen|J\mathclose|=\mathopen|f_y\mathclose|$, heranziehen, wie dies z.B. in dem Programm LSODA getan wird. Im Programm TENDLER wird dies nicht gemacht. Dies hat maßgeblich zwei Gründe.

(1) $J$ wird über mehrere Schritte, nicht selten über 20–40 Schritte, konstant gelassen, sodaß die Aktualität von $J$ für eine Momentanentscheidung (Schalten oder nicht Schalten) in Frage gestellt ist. Obwohl bei steifen Differentialgleichungen es eine natürliche Annahme ist, daß sich $J$ nur langsam ändert, muß dies nicht in gleichem Maße für $\mathopen|J\mathclose|$ gelten (trotz der Stetigkeit von $|\cdot|$) und “bei der Schwelle” von steif nach nicht-steif u.U. schon gar nicht. Bei den Programmen LSODE, LSODA und LSODAR wird $J$ genau dann neu berechnet, wenn $W=I-h\gamma J$ neu berechnet wird, sodaß in diesen Programmen $J$ i.a. aktueller ist. Das Programm VODE, siehe Hindmarsh/Brown/Byrne (1989), kann wie das Programm TENDLER, die Jacobimatrix $J$ unabhängig von $W$ konstant lassen. Z.Z. scheint es auch keine schaltfähige Version des Programmes VODE zu geben, insbesondere keine Version, die $\mathopen|J\mathclose|$ zum Schalten verwendet.

(2) Eine Berechnung von $\mathopen|J\mathclose|$ ist ca. $n$-mal so aufwendig, wie die Berechnung von $\mathopen|f(t,y^\nu)\mathclose|$. Das Pseudoresiduum $\mathopen|y^{\nu+1}-y^\nu\mathclose|$ muß ohnehin für die Konvergenzrate und insbesondere für den Konvergenztest berechnet werden. Obwohl bei großdimensionalen Differentialgleichungsproblemen $J$ i.d.R. Bandstrukturen oder allgemeiner dünn besetzt ist, also das “$n$-mal” abzuschwächen wäre, ist die Berechnung von $\mathopen|J\mathclose|$ stets um ein Vielfaches aufwendiger als die Berechnung von $\mathopen|f(t,y^\nu)\mathclose|$.

Die Gründe (1) und (2) favorisieren die Benutzung des $\text{QRP}_\nu$.

5. Die Konvergenzrate $f_\nu$ bei der Picard Iteration ist gleichzeitig eine untere Schranke für die Lipschitzkonstante $L$ der Funktion $f(t,y)$. Aus $$ y^{\nu+1} = h\gamma f(y^\nu) + \psi, \qquad y^\nu = h\gamma f(y^{\nu-1}) + \psi, $$ folgt sofort durch Subtraktion $$ f_\nu = {\mathopen|y^{\nu+1}-y^\nu\mathclose| \over \mathopen|y^\nu-y^{\nu-1}\mathclose|} \le \mathopen|h\gamma\mathclose| L \qquad (\nu\ge1). $$ Im Programm TENDLER wird während der Picard Iteration geprüft, ob $$ f_\nu \le 1/5 \qquad (\nu>0). $$ Dies ist implizit ein Test auf die Größe der Lipschitzkonstanten $L$, in Abhängigkeit von der Schrittweitengröße $\mathopen|h\mathclose|$. Wird dieser Test nicht bestanden, so wird Divergenz angenommen und somit, falls erlaubt, geschaltet.

Die Schaltentscheidung im Programm TENDLER benötigt während der Picard Iteration auf jeden Fall einen Konvergenztest, da andernfalls nicht in den steifen Modus übergewechselt werden kann. Die Unterprogramme fixpece() und fix() sind damit nicht benutzbar. Dies wird vom Programm TENDLER überprüft. D.h. wählt der Benutzer ausdrücklich diese Iterationsarten an, so kann er nicht gleichzeitig eine Schaltfähigkeit wollen.

Die bisherigen Erfahrungen mit der Schaltentscheidung über den Wert $\text{QRP}_\nu$ einer Newton-Kantorovich-artigen Iteration stellen praktisch im großen Maße zufrieden, obwohl theoretische Fragen im Prinzipiellen noch nicht ganz überwunden sind.

1. Zum einen findet dies seine Erklärung in der Tatsache, daß die “verbotene” Rückrichtung benutzt wird (so verboten wie: $f'(x_0)=0$ $\Rightarrow$ $f(x_0)$ extremal).
2. Zum anderen darin, daß in der Iterationsmatrix $W$, als Bestandteil des $\text{QRP}_\nu$, nicht nur ein zu altes $h\gamma$ einfließt, sondern auch ein zu altes $J$, also insgesamt zwei zu alte Werte. $W$ wird ja nicht immer dann faktorisiert, wenn $h$, oder $\gamma$, oder $J$ sich ändern.

Bibliographisch:

1. Shampine, Lawrence F.: “Type-Insensitive ODE Codes Based on Implicit $A$-Stable Formulas”, Mathematics of Computation, Vol 36, No 154, April 1981, pp.499–510
2. Peter N. Brown, George D. Byrne, and Alan C. Hindmarsh: “VODE: A Variable-Coefficient ODE Solver”, SIAM Journal on Scientific and Statistical Computing, Vol. 10, No 5, September 1989, pp.1038–1051

6. Bewertung des Programmes TENDLER

Die bisherigen Erfahrungen mit dem Programm TENDLER sind als durchaus erfolgversprechend einzustufen. Z.B. werden die Orbit Gleichungen ND sämtlich korrekt als nicht-steif erkannt, auch das Erde-Mond-Satellit Problem wird korrekt als nicht-steife Gleichung erkannt und demenstprechend gelöst.

Bei der van der Polschen Gleichung erfolgt ein korrektes Pendeln zwischen Newton-Kantorovich-artiger Iteration und Picard Iteration, also zwischen mnifloat() und fixfloat().

Dennoch ist das Vorhandensein von Effizienzreserven nicht ausgeschlossen. Der Einbau einer getrennten Schrittweiten- und Ordnungssteuerung erfolgte, als deutlich wurde, daß für den nicht-steifen Modus die Steuerung für den steifen Modus wohl zu vorsichtig und zögerlich die Schrittweite und Ordnung erhöhte. Durch diese zwei unterschiedlichen Schrittweiten- und Ordnungssteuerungen konnten deutliche Effizienzsteigerungen festgestellt werden, auch was den globalen Fehler anbetraf.

Die strengen Forderungen an die Schrittweiten- und Ordnungssteuerung für den steifen Teil wurden durch die Ergebnisse von Gaffney (1984) nahe gelegt, als auch durch entsprechende Protokollierungen mit dem Programm TENDLER.

Bei steifen Gleichungen, z.B. A-Gruppe und B-Gruppe aus STDTST, zeitigte die Schaltautomatik im Programm TENDLER klar günstigere Werte, als wenn man die Iterationsart auf mni() oder mnifloat() fixierte. Dies liegt maßgeblich an der Startphase, wo die Differentialgleichungen i.d.R. noch nicht genügend steif sind, um einen steifen Integrationsmodus zu erzwingen. Man könnte sogar soweit gehen und sagen, daß die Gewinne, die aus der Startphase stammen, so groß sind, daß nur aufgrund alleine der Startphase, sich eine Schaltautomatik auf jeden Fall anbietet. Dies liegt daran, daß während der Startphase die Schrittweite stark progressiv zunimmt, mit entsprechenden Konsequenzen für die Iterationsmatrix, welche ja von der Schrittweite abhängt. Bei den steifen Gleichungen aus STDTST kommt hinzu, daß diese Gleichungen zu Anfang auch tatsächlich nicht-steif sind, z.B. läßt sich das Problem B5 über einen anfänglich erstaunlich weiten Bereich schnell und genau mit Picard Iteration lösen (transiente Phase).

Nicht-steife Gleichungen auch komplizierterer Art konnten mit der Schaltautomatik wesentlich günstiger gelöst werden, als wenn man den steifen Modus ausdrücklich vorgibt. Unter den Alternativen FIXFLOAT, FIXPECE und SWITCHING konnte in diesem Falle kein eindeutiger Favorit ausgesondert werden. Jedoch haben die separaten Iterationsarten fixpi() und fpece() sich als durchaus zweckmässig erwiesen, wenn man im voraus um den Typus der Differentialgleichung weiß, also steif oder nicht-steif.

Vergleicht man jedoch die derzeitige Version des Programmes TENDLER mit den veröffentlichten Daten von DE/STEP, siehe Shampine/Gordon (1984), so ist das Programm DE/STEP bzgl. Funktionsauswertungen und globalem Fehler das eindeutig überlegenere Programm für nicht-steife Gleichungen.

Z.B. für die Differentialgleichung ND$i$ mit verschiedenen Ellipsenexzentrizitäten benötigt das Programm TENDLER zwei- bis dreimal soviele Funktionsauswertungen bei 100–1000 mal so großen globalen Fehlern. Für diese Unterschiede könnten mehrere Gründe eine Rolle spielen.

1. Die Fehlerkonstanten der zyklischen Formeln von Tendler sind größer als die Fehlerkonstanten für die Adams-Formeln.
2. Im Programm TENDLER gibt es kein Verfahren mit der Ordnung größer als sieben, während beim Programm DE/STEP Verfahren bis zur Ordnung zwölf verwendet werden können. Dies stellt wahrscheinlich mit einer der Hauptgründe dar, denn das Programm DE/STEP benutzt tatsächlich zur Integration der ND Gleichung fast ausschließlich Ordnungen größer als sieben, siehe Shampine/Gordon (1984).
3. Die Darstellung der Lösung und die damit verbundenen Probleme beim Schrittweitenwechsel unterscheiden sich.
4. Die Strategien im Programm TENDLER sind weiterhin auf Robustheit und Vorsicht ausgelegt.

Die Eigenschaft, daß bei den ND Gleichungen das Programm TENDLER stark von der exakten Lösung abdriftete, ist keine Eigenheit des Programmes TENDLER. Sowohl das Programme LSODA, als auch das Programm EPISODE arbeiten bei dieser Gruppe von Differentialgleichungen, insbesondere bei stärkeren Exzentrizitäten, vergleichsweise ungenau. Auch hier ist 100- bis 1000-faches Überschreiten der geforderten Genauigkeit zu beobachten, siehe Hairer/Wanner/Nørsett (1987) und siehe Bruder/Strehmel/Weiner (1988). Das Programm RKF4RW arbeitet bei der Gleichung ND noch ungenauer als die Programme PAI4 und LSODA.

Vergleiche bzgl. der Rechengeschwindigkeit liegen allerdings noch nicht vor. Dies hat mehrere Gründe. Dies hängt zum einen wesentlich mit dem “versioning” zusammen. Das Programm TENDLER ist kein unveränderliches Programm. Im Gegenteil. Verschiedene Versionen unterscheiden sich innerlich nicht unerheblich. Vorüberlegungen, Erfahrungen, Testläufe und Literaturberichte flossen beispielsweise in den Konvergenztest ein, jedoch nicht sämtlich auf einmal. Aber auch der Einbau einer zweiten Schrittweiten- und Ordnungssteuerung geschah erst nach Probeläufen mit einer Vorversion, die einen solchen Einbau als lohnenswert andeuten ließ. Auch die Auswertung der generierten Statistiken lieferte zahlreiche Hinweise für mögliche Verbesserungen. Durch die profiling Ergebnisse rückte die Normroutine und der Prädiktor in das engere Blickfeld. Durch die Protokolliermöglichkeiten wurden die Fehlerkontrolle, der Schaltteil, die Schrittweiten- und Ordnungssteuerung und der Konvergenztest näher beleuchtet.

Dennoch gibt es Anhaltspunkte für eine grobe Einordnung. Bei der van der Polschen Gleichung $$ \ddot y = 100(1-y^2)\dot y-y, \qquad y(0)=2, \quad \dot y(0)=0, $$ auf dem Intervall $[0,250]$ erfolgt ein periodisches Schalten zwischen Picard- und modifizierter Newton-Kantorovich Iteration. Ein Vergleich mit den Ergebnissen von LSODA zeigt, daß das Programm TENDLER leicht weniger Funktionsauswertungen ($-10%$), weniger als halb soviele Jacobimatrixauswertungen, aber um ein Drittel mehr $LU$-Zerlegungen benötigte. Beide Programme wurden mit einer Genauigkeitsanforderung von $10^{-6}$ konfrontiert. Über globale Fehler sind für beide Programme z.Z. für diese Differentialgleichung keine Daten vorhanden. Interessant ist bei der van der Polschen Gleichung noch, daß im Programm TENDLER die Normroutine ca. 10%, die Rücksubstitutionen 2% und die Funktionsauswertungen unter 2% der gesamten Rechenzeit beanspruchen. Man kann beobachten, daß das Programm LSODA eine eher “unruhige” Schaltung verwendet, welche gelegentlich zwischen den beiden Integrationsverfahren öfters hin- und her wechselt. Hingegen beim Programm TENDLER ist in diesem Falle und in weiteren Fällen eine “zitterfreie” Schaltung zu beobachten.

Es wäre beim Programm TENDLER möglich, die Schaltfähigkeit früher anziehen zu lassen, z.B. durch $\text{QRP}\le11/10$, jedoch beim Problem E4 treten dann eher unerwünschte Effekte ein. Tatsächlich zeigt das Programm LSODA bei der Differentialgleichung E4 leichte Schwierigkeiten, insbesondere bei höheren Genauigkeiten. Das Programm TENDLER löst E4 für alle Genauigkeiten genau und effizient, insbesondere schneller als das Programm STINT.

Die Schaltautomatik erkennt bei E4 korrekt nur die Startphase als nicht-steif und den Rest als steif.

Müßte man alle Ergebnisse in einem einzigen Satze zusammenfassen, so wäre zu antworten, daß eine Schaltautomatik für zyklische Mehrschrittverfahren zum Wechseln zwischen Newton und Picard Iteration, eine sinnvolle Erweiterung darstellt.

Weiß man von einer Differentialgleichung im voraus, daß sie ausschließlich nicht-steif ist, so ist das Programm DE/STEP ggü. dem Programm TENDLER vorzuziehen. Im Falle steifer, oder wechselhaft steifer Gleichungen ist das Programm TENDLER dem Programm DE/STEP klar und deutlich überlegen. Für steife Gleichungen ist das Programm DE/STEP allerdings auch gar nicht konzipiert.

Die Tatsache, daß Formeln für steife Gleichungen auf nicht-steife Probleme angesetzt, i.a. weniger effizient sind als Formeln für nicht-steife Gleichungen, ist bekannt, selbst dann, wenn man die Formeln für steife Gleichungen im nicht-steifen Modus (z.B. Picard $P(EC)^i{E}$ Modus) betreibt.

Z.B. ist das Programm LSODE mit der Einstellung mf=20 angesetzt auf nicht-steife Gleichungen, i.a. weniger effizient als die Einstellung mf=10, also BDF Picard $P(EC)^i$ Modus vers. Adams Picard $P(EC)^i$ Modus, $i\in{1,2,3}$.

Die zyklischen Formeln von Tendler (1973) sind gezielt für steife Differentialgleichungen entwickelt worden. Ihr Einsatz für nicht-steife Gleichungen war nicht speziell vorgesehen. Die Ergebnisse mit dem Programm TENDLER zeigen, daß die Formeln von Tendler (1973) durchaus für nicht-steife Gleichungen geeignet sind.

Es gibt jedoch für rein nicht-steife Gleichungen ggü. den Tendlerschen Formeln bessere Alternativen, nämlich z.B. die im Programm DE/STEP verwendeten Adams-Formeln. Besser hier im Sinne von kleineren Fehlerkonstanten, parasitäre Wurzeln bei Null und Schrittweitenwechseleigenschaften.

Bibliographisch:

1. Bruder, Jürgen und Strehmel, Karl und Weiner, Rüdiger: “Partitioned Adaptive Runge-Kutta Methods for the Solution of Nonstiff and Stiff Systems”, Numerische Mathematik, Vol 52, Fasc 6, 1988, pp.621–638
2. Petzold, Linda Ruth: “Automatic Selection of Methods for Solving Stiff and Non-Stiff Systems of Ordinary Differential Equations”, SIAM Journal on Scientific and Statistical Computing, Vol 4, No 1, March 1983, pp.136–148

7. Rechnung nahe der Maschinengenauigkeit

1. Ein wesentlicher Aspekt bei der Integration gewöhnlicher Differentialgleichungen ist das Verhältnis von geforderter Genauigkeit zu verfügbarer Rechengenauigkeit auf der zugrundeliegenden Rechenanlage. Ab einer gewissen Stelle führt ein weiteres Verkleinern der Schrittweite prinzipiell nicht mehr zu einer Erhöhung der Genauigkeit. Während dieser Effekt sonst im Zusammenhang steht mit dem Werteverhalten höherer Ableitungen der zu integrierenden Funktion, ist dieser Effekt eine Folge der Endlichkeit der Rechengenauigkeit. Der Grenzübergang $h\to0$ kann auf der Rechenanlage nicht vollzogen werden. Dieser an sich selbstverständliche Sachverhalt hat aber erhebliche Konsequenzen für einen Löser. Gewisse sehr kleine Fehlertoleranzen $\varepsilon$ müssen als unzulässige Eingaben abgefangen und zurückgewiesen werden.

Gewisse Teile des Programmes TENDLER berechnen Ausdrücke der Form $$ y\gets y+h\Sigma, $$ wobei $\Sigma$ hier für eine hier nicht näher interessierende Summe steht. In Wirklichkeit wird aber berechnet $$ y\gets y\oplus h\Sigma = \left(y+h\Sigma\right)(1+\tau)\approx (1+\tau)y, \qquad \mathopen|h\Sigma\mathclose| > \hbox{klein}. $$ Im Falle sehr hoher Genauigkeitsanforderungen kann man erwarten, daß $\mathopen|y\mathclose| \gg \mathopen|h\Sigma\mathclose|$. Es ist naheliegend zu fordern, daß $$ \varepsilon > u y, \qquad 1\oplus u > 1, $$ wobei $u$ die kleinste positive Zahl ist, sodaß auf der Rechenanlage gilt $1\oplus u>1$.

Um zusätzliche Sicherheit zu erlangen, wird im Programm TENDLER die Genauigkeitsanforderung weiter restringiert und zwar muß gelten $$ \varepsilon \ge 40u \left\Vert y\right\Vert_w,\quad \left\Vert y\right\Vert_w = \texttt{wnorm}(n,y,w) = \sqrt{{1\over n}\sum_{1\le i\le n}\left(y_i\over w_i\right)^2}. $$ Es ist für den Benutzer möglich diese Restriktion zu umgehen, indem er sein eigenes Überprüfprogramm chkinp() schreibt. Dies könnte in gewissen, sehr seltenen Fällen von Nutzen sein, wenn sich herausstellen sollte, daß die Überprüfungen zuviel Rechenarbeit kosten. Das Steuerprogramm tendler() vollzieht selber keine Überprüfungen, sondern delegiert alle diese Überprüfungen weiter an das Unterprogramm chkinp(). Das Unterprogramm chkinp() berechnet auch selbständig eine für die Differentialgleichung passende Maximalordnung und Genauigkeitsanforderung, wenn diese Werte mit AUTO_SELECT vorbesetzt sind.

Man wird erwarten, daß sehr nahe bei der Grenze der Rechengenauigkeit der Rechenanlage, also bei circa $$ \varepsilon\approx 40u \left\Vert y\right\Vert_w, $$ weitere Effekte eine Rolle spielen.

In der augenblicklichen Version des Programmes TENDLER werden hier keine extra Algorithmen verwendet. Allerdings nehmen gewisse Tests, wie z.B. die Konvergenztests, den Effekt des Rechnens nahe der Maschinengenauigkeit in ihre Entscheidungen mit auf. Das Programm DE/STEP von Shampine/Gordon (1984) verwendet sogar zusätzlich getrennte Algorithmen für die Berechnung des Prädiktors und des Korrektors, nämlich während der $P$- und $CE$-Phase beim $PECE$-Durchlauf.

2. Die grundlegende Idee ist hierbei die folgende. Sie soll anhand eines einfachen, ausführlich zu besprechenden Beispiels verdeutlicht werden. Der Wert der Riemannschen Zetafunktion (Riemann, Georg Friedrich Bernhard (1826–1866)) an der Stelle 4 sei zu bestimmen. Es ist $$ \zeta(4) = \sum_{j=1}^\infty {1\over j^4} = {\pi^2\over90} \approx 1.082 323 233 711 138. $$ Zur Schätzung, wieviele Glieder $n$ der Summe berücksichtigt werden müssen, um den Wert auf $\pm5\cdot10^{-11}$ (gerundete zehnstellige Genauigkeit) genau zu ermitteln, rechnet man $$ 0 \le \zeta(4) - \sum_{j=1}^n {1\over j^4} = \sum_{j=n+1}^\infty {1\over j^4} \le \int_n^\infty{dx\over x^4}={1\over3n^3}, $$ also $$ {1\over3n^3} < 5\cdot10^{-11},\qquad\hbox{damit}\qquad n>1883. $$

Die Summe wird nun mit drei verschiedenen Algorithmen $A$, $B$ und $C$ berechnet. Die Variable sum sei zu Anfang immer gleich 0.

Der Algorithmus $A$ laute

j=1,2,...,n
    sum += 1 / j^4

Der Algorithmus $B$ sei:

j=n,n-1,...,2,1
    sum += 1 / j^4

Schließlich der letzte Algorithmus $C$ sei:

roundoff_error = 0
j=1,2,...,n
    next_term = (1 / j^4) + roundoff_error
    temp_sum  = sum + next_term
    roundoff_error = next_term - (temp_sum - sum)
    sum       = temp_sum

Das erwartete Ergebnis ist: $$ {\pi^2\over90}\approx1.082\ 323\ 233\ 711\ 138 $$

Man erhält für unterschiedliche $n$ die Tabelle

Der Algorithmus $A$ summiert in aufsteigender Reihenfolge und der Algorithmus $B$ in umgekehrter Reihenfolge. In diesem speziellen Falle ist das Vorgehen wie bei Algorithmus $B$ optimal, da die Reihe monoton konvergent ist und beim numerischen Rechnen stets darauf zu achten ist, daß die Summanden in aufsteigender Reihenfolge ihrer Beträge zu summieren sind. Bei nicht monoton konvergenten Reihen ist diese Bedingung aber schwierig einzuhalten, wenn man nicht der Rechnung einen Sortierprozeß voran stellen will. Diese Sortierung wird auch stellenweise tatsächlich durchgeführt. Insbesondere wird vorgeschlagen dies nicht zu programmieren, sondern dies eine Ebene darunter direkt auszuführen, also Mikroprogrammierung oder u.U. direkte Implementierung auf der Gatterebene. Gerade bei Matrix-Vektor-Multiplikationen wird dieses Vorgehen teilweise angewandt, man vgl. hierzu die Arbeiten einführenden Charakters, siehe Bohlender/Grüner/Gudenberg (1982), siehe Rump (1980), sowie Rump (1983).

Der letzte Algortihmus $C$ verlangt zusätzliche Erklärungen. Es ist $$ \texttt{next_term} = (1/j^4) + \texttt{roundoff_error} $$ hierbei der nächste Summand für temp_sum, allerdings wurde der Wert next_term korrigiert und zwar mit roundoff_error. Die Summe $$ \texttt{temp_sum} = \texttt{sum} + \texttt{next_term} $$ verliert an Genauigkeit, eben wegen Rundung. Das was man jedoch an Genauigkeit verliert, holt man sich nun zurück. Ist jetzt $\texttt{sum}\gg\texttt{next_term}$, so ist der einzige Teil, der von next_term zu temp_sum hinzugeliefert wird, genau $\texttt{temp_sum}-\texttt{sum}$. Der Rest, also das was durch Rundung verloren gegangen ist, ist genau $$ \texttt{next_term} - \bigl(\texttt{temp_sum} - \texttt{sum}\bigr). $$ Also speichert man genau diese Größe für den nächsten Additionsschritt, $$ \texttt{roundoff_error} = \texttt{next_term} - \bigl(\texttt{temp_sum} - \texttt{sum}\bigr). $$ Schließlich setzt man sum = temp_sum.

Man vgl. zu diesen Bemerkungen Flanders (1987) aber auch Butcher (2016). Eine Beschreibung dieses Rechenkniffs findet man ebenfalls bei Kahan (1965).

Dieser Rechenkniff funktioniert nicht, wenn anstatt gerundeter Maschinenarithmetik nur abgeschnittene Arithmetik benutzt wird und wenn vor der Normalisierung nicht gerundet oder abgeschnitten wird. Auf Rechenanlagen vom Typ IBM 650, IBM 1620, Univac 1107 und CDC 3600 funktioniert dieser Rechenkniff also nicht. (Die genannten Maschinen sind heute veraltet.) Knuth (1971) verwendet einen axiomatischen Ansatz zur Erklärung des obigen Sachverhalts. Der Rechenkniff hat dann seine Erklärung in den beiden folgenden Aussagen.

3. Hilfssatz: Seien $u$ und $v$ normalisierte Gleitkommazahlen. Dann gilt $$ \bigl((u\oplus v)\ominus u\bigr)+\bigl((u\oplus v)\ominus\bigl(( u\oplus v)\ominus u)\bigr)=u\oplus v, $$ falls in keiner der obigen Ausdrücke ein Exponentenüberlauf oder Exponentenunterlauf entsteht.

Die Aussage des Hilfssatzes lautet umgeschrieben wie folgt, wobei zur einfacheren Schreibweise die Abkürzungen eingeführt werden $$ \eqalign{u'&:=(u\oplus v)\ominus v,\cr u''&:=(u\oplus v)\ominus v',\cr} \qquad \eqalign{v'&:=(u\oplus v)\ominus u,\cr v''&:=(u\oplus v)\ominus u'.\cr} \tag{*} $$ Hiermit gilt $$ u\oplus v = u'+v'' = u''+v', $$ was eine stärkere Aussage ist als $$ u\oplus v = u'\oplus v'' = u''\oplus v'. $$

4. Satz: Unter der Voraussetzung des Hilfssatzes und $(*)$ gilt $$ u+v = (u\oplus v)+\bigl((u\ominus u')\oplus(v\ominus v'')\bigr). $$

Beweis: siehe Knuth (1971), §4.2.2.     ☐

Dieser Satz liefert eine explizite Formel für die Differenz zwischen $u+v$ und $u\oplus v$ in Form berechenbarer Größen. Es ist somit möglich die Genauigkeit zu steigern, indem man die “Korrekturterme” $(u\ominus u')\oplus(v\ominus v'')$ akkumuliert. Der Beweis des obigen Satzes beruht auf der Tatsache, daß gilt $$ u\oplus v = \mathop{\rm round}(u+v,p) , $$ wobei für $x\in\mathbb{R}$ und $p\in\mathbb{N}$, $$ \eqalign{ \mathop{\rm round}(x,p) &:= \hbox{$x$ auf $p$ Stellen gerundet}\cr &:=\cases{ b^{e-p}\lfloor b^{p-e}x+{1\over2}\rfloor,&falls $b^{e-1}\le x<b^e$;\cr 0,&falls $x=0$;\cr b^{e-p}\lceil b^{p-e}x-{1\over2}\rceil,&falls $b^{e-1}\le-x<b^e$.\cr } } $$ Gilt dies nicht, so kann man den obigen Satz nicht mehr ohne weiteres heranziehen um eine Steigerung der Rechengenauigkeit zu erzielen.

Im Rahmen der Maschinengenauigkeit kann es durchaus vorkommen, daß gilt $2(u^2\oplus v^2) < (u\oplus v)^2$. Bei der Berechnung der Standardabweichung nach der Formel $$ \sigma = {1\over n}\sqrt{n\sum_{i=1}^n x_i^2-\left(\sum_{i=1}^n x_i\right)^2} $$ kann es somit passieren, daß man die Wurzel aus einer negativen Zahl zieht.

5. Nun darf die Effektivität dieses “Doppelrundungstrickes” nicht darüber hinweg täuschen, daß dieser Genauigkeitsgewinn fünfmal so aufwendig ist, wobei zusätzlich noch ein extra Vektor an Speicherplatz hinzukommt. sum und roundoff_error sind Vektoren. Zu Beginn der Summation müssen alle Komponenten von roundoff_error[] gelöscht werden. In dem Program DE/STEP von Shampine/Gordon (1984) wird dieser Aufwand auch tatsächlich getrieben, allerdings nur dann, falls $$ {1\over2}\varepsilon \le 100\cdot 2u \left\Vert y\right\Vert_w, \qquad \left\Vert y\right\Vert_w = \sqrt{\sum_{i=1}^n\left(y_i\over w_i\right)^2}. $$ In Abhängigkeit dieses Bedingung wird die logische Variable nornd (no round) gesetzt, die dann jedes Mal beim Prädiktorschritt und beim Korrektorschritt abgefragt wird:

if (nornd)
    // rechne wie üblich
else
    // tribe mehr Aufwand

An Speicherplatz werden hierfür $\phi_{15}$ und $\phi_{16}$, zwei Vektoren, zur Hilfe genommen. Im Programm TENDLER bräuchte man keinen zusätzlichen Speicherplatz anzufordern. Zur Realisierung dieses “Rechentricks” wäre in den nachfolgenden Blättern, bzw. in den nachfolgenden Spalten, der ringzyklischen Liste, ausreichend Platz. Die theoretische Begründung für die Benutzung dieses “Tricks” gelten nicht unverändert bei steifen Differentialgleichungen, sodaß diese Option vorerst nicht in das Programm TENDLER eingebaut wurde.

Shampine/Gordon (1984) zeigen anhand einer Beispieldifferentialgleichung (Orbitgleichung ND, mit Ellipsenexzentrizität $\varepsilon=0.6$) die Auswirkung dieser Steuerung mit Hilfe der Variablen nornd. Im folgenden bezeichne $\varepsilon$ wieder die vom Benutzer geforderte Genauigkeit. nfe gibt die Anzahl der Funktionsauswertungen an, und $\Delta$ ist der größte gemessene globale Fehler.

Die Ergebnisse einmal mit der eingeschalteten Steuerung sind:

Ohne die Steuerung lauten die Resultate:

Ein Stern neben der geforderten Genauigkeitsschranke zeigt an, daß die Fehlerschranke vom Programm DE/STEP vergößert werden mußte, um die Integration abzuschliessen. Man erkennt, daß die globalen Fehler ca. um eine Größenordnung genauer sind bei Benutzung der Steuerung. Dabei sind die benötigten Funktionsauswertungen fast völlig gleich geblieben.

Bibliographisch:

1. Bohlender, G., Grüner, K., Gudenberg Wolf von, J.: “Realisierung einer optimalen Arithmetik”, Elektronische Rechenanlagen 24 (1982), Heft 2, pp.68–72
2. Butcher, John Charles.: Numerical Methods for Ordinary Differential Equations. John Wiley & Sons, Chichester (2016)
3. Flanders, Harley: “Controlling Roundoff Error in Sums”, The College Mathematics Journal, Vol 18, No 2, March 1987, pp.153–156
4. Kahan, William M.: “Further Remarks on Reducing Truncation Errors”, CACM, Vol 8, No 1, January 1965, page 40
5. Knuth, Donald Ervin: The Art of Computer Programming, Volume 2 / Seminumerical Algorithms, Addison-Wesley Publishing Company, Reading (Massachusetts) Menlo Park (California) London Don Mills (Ontario)
6. Rump, S.M.: “Notiz zur Genauigkeit der Arithmetik in Rechenanlagen”, Elektronische Rechenanlagen 22 (1980), Heft 5, pp.243–244
7. Rump, S.M.: “Mathematik auf dem Rechner—Computer Algebra und neue Numerik”, Elektronische Rechenanlagen 25 (1983), Heft 6, pp.126–132