, 62 min read

TENDLER: 4. Die Korrektoriteration

In diesem Abschnitt wird eines der mitaus wichtigsten Kernstücke eines jeden Lösers vorgestellt: die Korrektoriteration. Die Effizienz dieses Segments wirkt maßgeblich zurück auf die anderen Segmente des Programmes. Nun gilt dies natürlich auch für andere Segmente, jedoch ist hier bei der Korrektoriteration die Hebelwirkung besonders ausgeprägt. Auch für das Schalten spielt sie im Programm TENDLER eine zentrale Rolle.

Bei den zyklischen Formeln von Tendler (1973), sind sämtliche Stufen aller Zyklen echt implizite lineare k-Schrittverfahren, für $k=1,\ldots,7$, analog wie die BDF$i$, für $i=1,\ldots,6$. Zur “Beseitigung” der Implizitheit ist ein irgendwie geartetes Iterationsverfahren nötig, da eine geschlossene Auflösung in symbolischer Form bei allgemeinem $f$, in voller Allgemeinheit nicht möglich ist und häufig auch gar nicht zweckmässig wäre. Bei speziellen Differentialgleichungen, z.B. bei linearen Differentialgleichungen gelingt eine Auflösung. Gerade im Sonderfall linearer Differentialgleichungen, welche eine für die Anwendungen besonders wichtige Rolle spielen, wäre es möglich, aus dem a priori Wissen über die Differentialgleichung, Effizienzsteigerungen für einen Differentialgleichungslöser zu erhalten. Das Programm TENDLER in seiner jetztigen Form nutzt wie viele andere Löser, diese möglichen Vorteile nicht voll aus.

In jedem Zeitschritt wird ein Programm basierend auf impliziten linearen Mehrschrittverfahren (und auch viele weitere Verfahren) mit einer i.a. nichtlinearen Gleichung der Form $y_{m+1}=h\gamma f(y_{m+1})+\psi$, und $y_{m+1}=?$ konfrontiert, $h$, $\gamma$, $\psi$ fest. Werden 1000 Schritte durchgeführt, so ist also auch 1000-mal eine solche Gleichung aufzulösen. Die Anzahl kann sich noch weiter erhöhen, falls Schritte aufgrund von Fehlerversagen zurückgewiesen werden. Man könnte jetzt daran denken einen allgemeinen Gleichungslöser aufzurufen, wie ihn zahlreiche Programmsammlungen, wie etwa NAG, CERN, Harwell, TOMS, $\ldots$, bereitstellen. Bekannte Differentialgleichungslöser, wie LSODE, DE/STEP, $\ldots$, tun dies nicht. Es erweist sich als erheblich effizienter, die Gleichungsauflösung selbst in die Hand zu nehmen.

Bei steifen Differentialgleichungen geschieht die Gleichungsauflösung in vielen bekannten Programmen mit Hilfe eines Newton-Verfahrens, welches man zur Effizienzsteigerung des Differentialgleichungslösers (nicht notwendig zur Steigerung der Effizienz des Gleichungslösers) geeignet modifiziert. Bei nicht-steifen Differentialgleichungen geschieht die Gleichungsauflösung in fast aller Regel mit Hilfe von Picard Iteration (Fixpunktiteration), manchmal wird jedoch auch hier eine Variante eines Newton-Jacobi Iterationsverfahrens benutzt. Die Anzahl der Iterationen bei diesen Iterationsverfahren liegt fast immer deutlich unter 10, allerdings gibt es hier wenige Ausnahmen. Vorgeschlagen werden jedoch auch Mehrgitterverfahren, siehe Sommeijer/Houwen (1984).

Die bei der Newton Iteration anfallenden linearen Gleichungssysteme stellen erneut eine Arbeitsbelastung dar, für einen Löser konzipiert für steife Differentialgleichungen. Für die Lösung dieser linearen Gleichungssysteme werden eine Reihe von Lösungsverfahren angewendet: volle $LU$-Zerlegung (sehr häufig), Band-$LU$-Zerlegung (häufig), Gauß-Seidel Iteration (gelegentlich), $LU$-Zerlegung für dünn besetzte Matrizen und sogar Krylow-Unterraumiteration, siehe Gear/Saad (1983). Die augenblickliche Version des Programmes TENDLER verwendet an dieser Stelle lediglich eine volle $LU$-Zerlegung. Diese Beschränkung ist nicht durch Formeln oder Programmiermethodik erzwungen, und es gibt kein Hindernis hier einen Eingriff vorzunehmen. Im Gegenteil, dieser Eingriff ist an geeigneten Stellen schon mitbedacht und vorbereitet, siehe wprepare(). Durch die Beschränkung auf eine volle $LU$-Zerlegung wird die Anwendbarkeit auf steife Differentialgleichungen großer Dimension, wie sie beispielsweise bei Semidiskretisierungen partieller Differentialgleichungen auftauchen, erheblich eingeschränkt, eigentlich sogar fast gänzlich verunmöglicht. Dies stellt eine echte Beschränkung dar. Will man dennoch solche großdimensionalen Gleichungen lösen, so muß man beispielsweise auf das Programm LSODE oder einer seiner Varianten, zurückgreifen, siehe Hindmarsh (1983). Das Programm LSODE ist auch in der Unterprogrammsammlung NAG in leicht veränderter Form enthalten, siehe Gladwell/Berzins/Brankin (1988).

Zuerst wird schnell der Zusammenhang zwischen Picard und Newton Iteration hergestellt. Beide Iterationsarten sind bei Anwendung auf die Keplersche Gleichung sehr ähnlich. Anschliessend werden beide Iterationen einzeln vorgestellt. Neben dem eigentlichen Konvergenztest gibt es weitere Segmente, welche letztlich die einzige Aufgabe haben, den Konvergenztest zu ergänzen.

Zum Abschluß wird verglichen, wie andere bekannte Löser für steife Differentialgleichungen, basierend auf Mehrschrittverfahren, hier die Programme STINT, DIFSUB, GEAR, EPISODE, LSODE und DSTIFF, den Konvergenztest gestalten. Der Konvergenztest hat auch bei diesen Programmen entscheidenden Einfluß auf die Effizienz des gesamten Restprogrammes.

Bibliographisch:

  1. Sommeijer, B.P. und van der Houwen, P.J.: “Algorithm 621. Software with Low Storage Requirements for Two-Dimensional, Nonlinear, Parabolic Differential Equations”, ACM TOMS, Vol 10, No 4, December 1984, pp.378–396
  2. Hindmarsh, Alan C.: “ODEPACK, A Systematized Collection of ODE Solvers”, in Robert S. Stepleman, M. Carver, R. Peskin, William F. Ames and Robert Vichnevetsky (Editors): “Scientific Computing—Applications of Mathematics and Computing to the Physical Sciences”, IMACS Transactions on Scientific Computation, Volume 1, Noth-Holland Publishing Company, Amsterdam New York Oxford, 1983, pp.55–64
  3. Gear, Charles William und Saad, Youcef: “Iterative Solution of Linear Equations in ODE Codes”, SIAM Journal on Scientific and Statistical Computing, Vol 4, No 4, December 1983, pp.583–601
  4. Gladwell, Ian und Berzins, M. und Brankin, R.W.: “The Stiff Integrators in the NAG Library”, ACM SIGNUM Newsletters, Vol 3, No 2, April 1988, pp.16–23

1. Zusammenhang zwischen Picard und Newton Iteration

A choice of the basic variable has to be made. The literature is confusing as to what the possibilities are and the consequences of the choice. These matters are clarified.

L.F. Shampine (1980)

Zur Lösung allgemeiner, nichtlinearer Gleichungssysteme im $\mathbb{R}^n$ existieren eine Fülle von Verfahren. Zwischen zwei dieser Verfahren bestehen enge Verwandtschaftsbeziehungen: zwischen der Picard Iteration und der Newton Iteration.

Ist

$$ x = \varphi(x), \qquad\hbox{für}\quad x\in\mathbb{C}^n, $$

so iteriert man mit Hilfe des Newton-Verfahrens nach der Vorschrift

$$ x^{\nu+1} = x^\nu - (I-\varphi')^{-1}(x^\nu - \varphi(x^\nu)). $$

Setzt man $\varphi'=0$, so erhält man direkt die Rechenvorschrift für die Picard Iteration

$$ x^{\nu+1} = \varphi(x^\nu), \qquad\hbox{oder}\qquad x^{\nu+1} = x^\nu+(\varphi(x^\nu)-x^\nu), $$

welches auch als Residueniteration bezeichnet wird. Andere Bezeichnungen sind Fixpunktiteration oder repetierende Substitution.

Für die Keplersche Gleichung

$$ y = \varepsilon f(y)+\psi =: \varphi(y), \qquad y\in\mathbb{R}^n, $$

braucht man keine weitere Erklärung, da sie in der obigen allgemeinen Form der Fixpunktschreibweise enthalten ist. Diese Erklärungen seien jetzt schon gegeben, obwohl trivial, da es dann später, in unübersichtlicherem Zusammenhang, leichter ist den einfachen Kern herauszulesen. Für Spezialfälle der eindimensionalen Keplerschen Gleichung existieren speziell angepasste Lösungsverfahren. Einen Überblick findet man beispielsweise bei Dandby (1987), und ein spezielles Verfahren für den eindimensionalen Fall wird beschrieben bei Mikkola (1987). Im Programm TENDLER wird zur Lösung der mehrdimensionalen Keplerschen Gleichung ein modifiziertes Newton-Kantorovich Iterationsverfahren benutzt.

Bibliographisch:

  1. Shampine, Lawrence: “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
  2. Picard, Charles Emile (1856–1941)
  3. Newton, Sir Isaac (1643–1727)
  4. Dandby, J.M.A.: “The Solution of Kepler's Equation, III”, Celestial Mechanics, Vol 40, 1987, pp.303–334
  5. Mikkola, Seppo: “A Cubic Approximation for Kepler's Equation”, Celestial Mechanics, Vol 40, 1987, pp.329–334

2. Newton Iteration im Programm TENDLER

Therefore the usual iteration mode will reintroduce the numerical instability for insufficiently small $h$ which is our objective to avoid.

Irwin W. Sandberg, H. Shichman (1968)

Neben der Tatsache, daß bei steifen Differentialgleichungen häufig die benutzten Formeln für große Schrittweiten instabil werden, tritt jetzt noch eine weitere Schwierigkeit hinzu. Es ist nicht länger möglich durch einfache Fixpunktiteration, die impliziten Formeln nach den $y_{m\ell+i}$ aufzulösen, weil dies erneut eine erhebliche Einschränkung der Schrittweite $h$ bedeuten würde, ganz entgegen dem Verhalten der exakten Lösung der Differentialgleichung, welche sehr wohl große Schrittweiten geradezu anbieten würde. Es besteht natürlich ein enger Zusammenhang zwischen Instabilität und Nichtkonvergenz der Picard Iteration.

Nun würde einfache Fixpunktiteration stets konvergieren, bei Erfülltsein von

$$ \mathopen\|h {\beta_{ii} \over \alpha_{ii}} J \mathclose\| \lt 1, \qquad \hbox{$J=f_y$ Jacobimatrix von $f$.} $$

Dies ist immer erreichbar für genügend kleine Schrittweite $h$. Bei steifen Poblemem ist man jedoch gerade daran nicht interessiert. Zur Vermeidung dieser Problematik könnte man deswegen das Newton Iterations-Verfahren verwenden, mit der Vorschrift

$$ y_{m\ell+i}^{\nu+1} = y_{m\ell+i}^\nu - W_\nu^{-1} \bigl[ y_{m\ell+i}^\nu - h {\beta_{ii} \over \alpha_{ii} } \dot y_{m\ell+i}^\nu + \underbrace{ {1\over\alpha_{ii}} \sum_{j=-k+i}^{i-1}[\alpha_{ij}y_{m\ell+j}-h\beta_{ij}\dot y_{m\ell+j}] }_{\hbox{${}=:(-\psi)$ konstant}} \bigr], $$

mit

$$ W_\nu := I - h{\beta_{ii}\over\alpha_{ii}} f_y(t_{m\ell+i},y^\nu_{m\ell+i}) \in \mathbb{R}^{n\times n}, $$

um die gesuchten $y_{m\ell+i} \in \mathbb{R}^n$ zu bestimmen.

1. Um den gewaltigen Rechenaufwand, der mit diesem Verfahren verbunden ist, einzugrenzen, wird ein modifiziertes Newton-Kantorovich Verfahren verwendet. Die Modifizierung des Newton-Kantorovich-Verfahrens besteht darin, daß die Iterationsmatrix $W$, hier mit

$$ W = (I_n - h{\beta_{ii}\over\alpha_{ii}} J),\qquad J = f_y(t_{m\ell},y_{m\ell}), $$

nicht in jedem Integrationschritt ausgewertet wird.

Der Unterschied des Kantorovich-Verfahrens gegenüber dem Newton-Verfahren (auch Newton-Raphson Iteration) bestand ja darin, die Iterationsmatrix $W_\nu$ nur am Anfang der Iteration zu berechnen und damit nur eine einzige $LU$-Zerlegung zu benötigen, im Gegensatz zum Newton-Verfahren, welches pro Iterationsschritt eine $LU$-Zerlegung benötigt.

Man erkennt, daß die Jacobimatrix $J$ ausgewertet wird an der Stelle $(t_{m\ell},y_{m\ell})$, also für die letzte Stufe des letzten, gerade erfolgreich abgeschlossenen Zyklus (vom Start einmal abgesehen). Dies ist eine bewußt vorsichtige Strategie.

Andere Strategien sind denkbar und werden auch tatsächlich in existierenden Programmen angewendet. Das Programm LSODE beispielsweise wertet die Jacobimatrix $J$ mit dem gerade aktuellen Prädiktorwert aus. Da der Schritt durch anhaltende Divergenz, durch Nichtinvertierbarkeit von $W$ und insbesondere durch Fehlertestversagen zurückgewiesen werden kann, ist $(t_{m\ell},y_{m\ell})$ eine geeignete Wahl.

Bzgl. der Variabilität von $J$, beachte man auch die nachfolgenden Ausführungen von Shampine (1980).

Die Iterationsmatrix $W$ wird in dem Programm TENDLER nur dann berechnet, falls

  1. das modifizierte Newton-Kantorovich Verfahren selbst nach “ausreichend vielen” (häufig drei) Iterationen den Konvergenztest nicht bestand,
  2. u.U. nachdem die Fehlerkontrolle nicht passiert wurde, oder
  3. die Schrittweite mal $\beta_{ii}/\alpha_{ii}$ sich um mehr als 40% seit der letzten Berechnung geändert hat.

Man beachte, daß die Iterationsmatrix eine Funktion aus dem Produkt von $h$ und $\beta_{ii} / \alpha_{ii}$ ist. Es wird sogar eigentlich nur eine “mittlere” Iterationsmatrix genommen, die sich aus dem Mittel der verschiedenen $\beta_{ii} / \alpha_{ii}$ der einzelnen Stufen zusammensetzt. Die Konvergenzeigenschaften werden dadurch nicht wesentlich verändert.

Es sei wieder $\gamma = \beta_{ii} / \alpha_{ii}$. Deutlich gemacht werden sollte folgendes: Das unmodifizierte Newton-Verfahren wird auf mehrere Art und Weise verändert, mit dem Ziel der Aufwandsreduktion. Dreh- und Angelpunkt ist hierbei die Iterationsmatrix $W$. Folgende Veränderungen werden durchgeführt:

  1. Die Iterationsmatrix $W$ wird nicht in jedem Iterationsschritt neu ausgewertet.
  2. Die Iterationsmatrix $W$ wird nicht in jedem Integrationsschritt neu ausgewertet.
  3. Wenn $W$ neu ausgewertet wird, so wird nicht notwendigerweise auch gleichzeitig die Jacobimatrix $J$ mit ausgewertet. Es wird häufig (typisch: alle 20–40 Schritte) ein und dieselbe Jacobimatrix $J$ verwendet, allerdings mit aktualisierter Schrittweite $h$ und formelabhängigen Größen $\beta_{ii}$ und $\alpha_{ii}$.

2. Die Beibehaltung der Iterationsmatrix $W$ über mehrere Integartionsschritte gliedert sich erneut auf in

  1. $h\gamma$ ändert sich,
  2. $h\gamma$ ändert sich nicht.

Hier sprechen Tischer/Gupta (1983)1:

  1. im ersten Fall von einem modifizierten Newton-Verfahren und
  2. im zweiten Falle von einem simplifizierten Newton-Verfahren.

Bibliographisch:

  1. Tischer, Peter E. und Gupta, Gopal K.: “A Cyclic Method Stiff ODE Solver”, Technical Report 38, Department of Computer Science, Monash University, Clayton, Victoria, Australia, June 1983, i+45 S.

Der erste Fall tritt beispielsweise auf bei der Änderung der Schrittweite $h$, aber auch alleine schon das Durchlaufen der Tendlerschen Zyklen mit seinen i.d.R. nichtgleichen Stufen. Allerdings liegen diese letzten Variationen immer unter 20%. Von Wichtigkeit werden die Überlegungen von Shampine (1980), S.111:

The fragments of theory available for stiff problems give a great deal of attention to model problems. For the typical models to have any validity, it is necessary that the solution of the differential equation be slowly varying and the Jacobian be approximately constant in the neighborhood of the solution. We did not make it clear enough in $\ldots$ that there are important practical reasons for supposing this situation. If the problem is stiff, we hope to be able to use “large” step sizes $h$. According to

$$ y^*-y^{m+1} = (I-h\gamma J)^{-1}h\gamma({\cal J}_m-J)(y^*-y^m), $$

this is not likely to be possible unless $J$ is close to ${\cal J}_m$ for each $m$. This requires that the Jacobian be nearly constant near $y^*$. Using a factorization of $I-h\gamma J$ for several steps is not likely to be possible unless the solution changes slowly and the Jacobian changes slowly along the solution.

Weiter

On general grounds we can argue that most codes form a new $J$ far too frequently. For the common case of Jacobians which are fairly expensive, we would like to avoid this wasted effort. Sometimes Jacobians are inexpensive, e.g., they are “free” for linear problems. Avoiding the evaluation of Jacobian in such a case is unnecessary, but should do no harm. In the absence of information about the relative cost of evaluating a Jacobian, we shall presume that it is substantial and seek to minimize it. $\ldots$ If $J$ is a Jacobian approximation based on a solution of about the current size, the need for a new iteration matrix is likely to be due almost solely to changing $h\gamma$, in part because we have seen that $J$ does not need to be very accurate. In either case the arguments suggest that a good tactic is to keep a copy of $J$, and when forming a new iteration matrix, to first try the old $J$. The substantial increase of the storage required by the code is probably why few codes do this. $\ldots$ We have argued that the typical code forms a new approximate Jacobian too often. We advocate keeping a copy of $J$ and judiciously reusing it. This does require extra storage, but we think it worthwhile.

Das Programm TENDLER hält sich vollständig an diese Empfehlungen.

Die Iterationsmatrix $W$ wird noch zusätzlich dahingehend modifiziert, daß die Jacobimatrix nur dann ausgewertet wird, falls

  1. die Integration startet,
  2. $h\gamma$ sich um mehr als 100% geändert hat,
  3. entsprechend viele Minuspunkte durch langsame Konvergenz angesammelt wurden.

Durch die getrennte Behandlung von $W$ und $J$ sind zwei Matrizen zu speichern (dies ist ein Nachtei), jedoch kann die Anzahl der Jacobimatrixauswertungen drastisch reduziert werden (dies ist ein Vorteil). In diesem Zusammenhang schreibt Enright (1978):

One unexpected observation, which is illustrated in the results of problem D4, is that on the nonlinear problems tested, where the Jacobian often changed significantly during integration, convergence difficulty was seldom encountered by MOD-EPISODE. As a result, Jacobian evaluations $\ldots$ are rare for MOD-EPISODE.

In dem Programm MOD-EPISODE wird Enright's update technique, siehe Enright (1978), verwendet. Mit Hilfe dieser Auffrischungstechnik für die Iterationsmatrix $W$, ist es möglich $W$ zu aktualisieren, ohne $W$ zu refaktorisieren, vorausgesetzt, daß nur $h\gamma$ aktualisiert werden muß und nicht etwa $J$. Im Gegenzug allerdings nehmen die Kosten für die Lösung eines Systems $Ax=b$ mit faktorisierter Matrix $A$, leicht zu, und zwar um ${1\over2}n^2$ Operationen. Gleichzeitig ist noch ein zusätzlicher Pivotvektor vonnöten. Diese Ergebnisse von Enright (1978) verdeutlichen besonders hervorhebenswert, daß für die Konvergenz des also mehrfach modifizierten Newton-Verfahrens, der Wert $h\gamma$ von entscheidender Bedeutung ist, nicht jedoch so sehr von $J$.

Zu ähnlichen Schlüssen gelangen Karasalo/Kurylo (1981), Kopie. Das Nachfolgerprogramm von LSODE (bzw. in gewisser Hinsicht mehr das Nachfolgerprogramm von EPISODE) von Hindmarsh/Brown/Byrne (1989), nämlich VODE, verwendet ebenfalls, wie das Programm TENDLER, zwei getrennte Matrizen $W$ und $J$.

Erkennt man nun, daß an $J$ nicht hohe Genauigkeitsansprüche gestellt werden müssen, werden zwei Wege deutlich, um weitere Ersparnisse zu erzielen:

  1. Gewisse Elemente von $J$ können zu Null gesetzt werden, falls sie nicht zu “groß” sind. Damit ist es beispielsweise möglich, Bandstrukturen in $W$ zu erzielen und hierfür geeignete Routinen der linearen Algebra zu verwenden, z.B. Bandelimination.
  2. Wird $J$ durch numerische Differentiation bestimmt, was bei komplizierten Problemen eine große Arbeitserleichterung für den Anwender bedeutet (man denke beispielsweise an NC5), so ist es u.U. für den Benutzer einfacher, ein ungenaues Schema zu verwenden.

Üblich bei gängigen Differentialgleichungslösern ist eine Approximation der Jacobimatrix durch

$$ {\partial f\over\partial_i e_i} \approx {f(y+\sigma e_i)-f(y)\over\sigma_i}, \qquad e_i = \mathop{\rm col}_{\nu=1}^n\delta_{i\nu} = \pmatrix{0\cr\vdots\cr1\cr\vdots\cr0\cr}, \qquad \hbox{$\sigma_i$ klein.} $$

Die Programme unterscheiden sich in der Wahl von $\sigma$.

3. Auch Tischer/Gupta (1983)1, siehe oben, analysieren die Korrektoriteration, ausgehend von den Überlegungen in Shampine (1980). Eine weitere Modifikation des Newton-Verfahrens, wie sie bei der allgemeinen Lösung von nichtlinearen Gleichungen benutzt wird, nämlich die Dämpfung, findet in dem Programm TENDLER keinerlei Anwendung.

Dies hat mehrere Gründe.

  1. Zum einen befindet sich das Programm TENDLER zu einem maßgeblichen Teil im $P(EC)$-Modus, es wird also nur ein einziges Mal pro Schritt iteriert. Dadurch kann eine Dämpfungsstrategie, welche auf dem Verhalten vergangener Iterierten beruht, ersteinmal schon gar nicht mehr eingesetzt werden.
  2. Zweitens, wenn nun in seltenen Fällen doch mehr als zweimal iteriert wird, so bedeutet dies, daß i.d.R. die Schrittweite $h$ zu groß ist oder die Iterationsmatrix $W$ zu alt ist, $J$ zu alt, o.ä. und nicht, daß das Newton-Verfahren als solches Konvergenzschwierigkeiten hat. Dies erkennt man daran, wenn man im Programm TENDLER eine entsprechende Option setzt und sich somit entsprechende Werte, wie Konvergenzrate u.ä., vor und nach Versagen ausgeben lässt.
  3. Drittens zeitigten die Testergebnisse von Tischer/Gupta (1983)1, S.23, daß das Programm ODIOUS mit gedämpften Newton-Verfahren nicht selten einen guten Prozentpunkt schlechter arbeitet (Funktionsauswertungen, Schritte), als das ungedämpfte Verfahren. Hierbei gilt zu bemerken, daß das Programm ODIOUS hauptsächlich im $P(EC)(EC)$-Modus arbeitet, sodaß der Steuerung des Dämpfungsfaktors noch vergleichsweise “viele” Iterierte zur Verfügung stehen.

Umgekehrt zeigen Ergebnisse anhand des Programmes STRIDE, daß das Newton-Verfahren mit Dämpfungsstrategie Ersparnisse im Bereich 8–22% Funktionsauswertungen und 18–22% Jacobimatrixauswertungen liefern kann. Dies verdeutlicht, daß die Newton Iteration bei impliziten Runge-Kutta-Verfahren ein anderes Gesamtgewicht auf das Restprogramm ausübt, als anders etwa bei linearen Mehrschrittverfahren und zyklischen Verfahren.

Erneut verdeutlicht dies, daß Änderungen in einem Programm nicht die gleichen Änderungen in anderen Progammen direkt nahe legen.

Nun ist es natürlich möglich eine Dämpfungsstrategie auf anderen Kriterien beruhen zu lassen, als auf dem Konvergenzverhalten vergangener Iterierter. Shampine (1980) weist darauf hin, daß die Iterationsmatrix $W$ “steife Komponenten” und “nicht-steife Komponenten” unterschiedlich bewertet. Dies hängt u.a. an den Eigenrichtungen von $W$. Es ist also nicht so, daß sich eine Dämpfungstechnik in das Programm TENDLER grundsätzlich nicht einbauen lässt. Das Gegenteil ist der Fall. Aufgrund der weitreichenden Kapselung von Aufgaben würde sich dies als eine einfache Aufgabe darstellen. Dennoch gab es, s.o., mehrere Gründe eine Dämpfungstechnik nicht einzubauen.

4. Die Ableitung für die endgültige Programmierung der modifizierten Newton-Kantorovich Iteration verläuft nun wie folgend. Die i.d.R. nicht-linearen Gleichungen der Form

$$ y = h\gamma f(y)+\psi, $$

werden mit Hilfe des Newton-Kantorovich Iterationsverfahrens gelöst der Form

$$ \eqalign{ y^{\nu+1} &= y^\nu - W^{-1}(y^\nu - h\gamma f(y^\nu)-\psi)\cr &= y^\nu + \gamma W^{-1}(hf(y^\nu) - z^\nu).\cr } $$

Der Subskript $m\ell+i$ wurde hier zur Schreibvereinfachung weggelassen. Für die Variable $z=hf(y)$ ist die Ableitung völlig analog. Wegen $z = hf(\gamma z+\psi)$ ist

$$ \eqalign{ z^{\nu+1} &= z^\nu - W^{-1}(z^\nu - hf(y^\nu))\cr &= z^\nu + W^{-1}(hf(y^\nu) - z^\nu).\cr } $$

5. Im allgemeineren Falle der Differentialgleichung $E\dot y=f(y)$ erhält man die verallgemeinerte Keplersche Gleichung

$$ Ey = h\gamma f(y)+\psi = \gamma z+\psi, $$

mit der Abkürzung $z=hf(y)$ dann entsprechend

$$ z = {1\over\gamma}(Ey-\psi). $$

Die Iterationsmatrix $W$ ist nun $W=E-h\gamma J$ und man erhält

$$ \eqalign{ y^{\nu+1} &= y^\nu - W^{-1}(Ey^\nu - h\gamma f(y^\nu) - \psi)\cr &= y^\nu + \gamma W^{-1}(hf(y^\nu) - z^\nu), \qquad\qquad |{}\cdot E |{}-\psi |{}\cdot{1\over\gamma}\cr } $$

multipliziert man von links mit der Matrix $E$, zieht man danach $\psi$ ab und teilt schließlich durch $\gamma$, so erhält man den gewünschten Ausdruck für $z^{\nu+1}$ zu:

$$ z^{\nu+1} = z^\nu + EW^{-1}(hf(y^\nu)-z^\nu) . $$

Biographisch:

  1. Hindmarsh, Alan C.
  2. Brown, Peter N.
  3. Byrne, George D.
  4. Karasalo, Ilkka
  5. Kurylo, John

3. Die Picard Iteration im Programm TENDLER

Hier reduziert sich $W$ auf $W=I$ wegen $J=0$, und es ergibt sich

$$ \eqalign{ y^{\nu+1} &= y^\nu + \gamma(hf(y^\nu)-z^\nu),\cr z^{\nu+1} &= hf(y^\nu).\cr } $$

Die Programmierung der Picard Iteration ist nun wie folgt:

$$ \eqalign{ \hat z &\gets hf(y^\nu),\cr \Delta &\gets \hat z - z^\nu,\cr y^{\nu+1} &\gets y^\nu + \gamma\Delta,\cr z^{\nu+1} &\gets \hat z.\cr } $$

Die nullten Iterierten $y^0$ und $z^0$ werden vom Prädiktor geliefert. Bei Picard Iteration mit fest vorgegebener Anzahl der Iterationen, kann man bei $2k+1$ bzw. $2k$ Iterationen, $k$ Vektorkopieranweisungen durch geeignete Programmierung einsparen.

4. Der Konvergenztest im Programm TENDLER

There is already a problem, however, which you will have noticed if you are the kind of person who does not divide blindly. Even for $h\ne0$ we might have $g(a+h)-g(a)=0$, making the division and multiplication by $g(a+h)-g(a)$ meaningless.

Michael Spivak (1967)

Ein Konvergenztest in einem Differentialgleichungslöser soll feststellen, wann genügend viele Iterationen ausgeführt wurden. Es handelt sich im eigentlichen Sinne des Wortes nicht um einen Konvergenztest, sondern um einen Test auf Abbruch der Iteration. Da prinzipiell nur endlich viele Iterationen durchgeführt werden, ist es auch schon aus diesem Grunde möglich eine Funktion anzugeben, sodaß keine Konvergenz vorliegen braucht. Die Zahl der Iterationen ist i.a. sehr gering. Sie liegt i.a. zwischen 1 und 3. Bei den Programmen LSODE, EPISODE, VODE, GEAR, DIFSUB, etc., liegt die Anzahl der Iterationen in ebensolchen Bereichen. Beim Programm TENDLER kann man die maximale Anzahl der Iterationen frei wählen. Erfahrungen lassen aber Werte größer als 5 nicht als besonders effizient erscheinen. Dies gilt allerdings ersteinmal nur für das Programm TENDLER.

1. Es sei $\delta_\nu=\mathopen|y_N^{\nu+1}-y_N^\nu\mathclose|$ die Norm des Pseudoresiduums bei der Newton-Kantorovich Iteration, und $r_\nu=\mathopen|y_P^{\nu+1}-y_P^\nu\mathclose|$ sei die Norm des Pseudoresiduums bei der Picard Iteration. Sei $k_\nu$ die Konvergenzrate des modifizierten Newton-Kantorovich Iterationsverfahrens, und sei $f_\nu$ die Konvergenzrate der Picard Iteration. $\hat k_\nu$, bzw. $\hat f_\nu$, kennzeichnen die Konvergenzraten der entsprechend vorherigen Iteration, entsprechend sei $\hat k_\infty$, bzw. $\hat f_\infty$ die letzte Konvergenzrate von der vorherigen Iteration (vorheriger Zeitschritt). Vor jedem Konvergenztest werden die beiden Werte $k_\nu$ und $f_\nu$ wie folgt aktualisiert

$$ k_\nu\gets\cases{ \displaystyle\max\left({\delta_\nu\over\delta_{\nu-1}}, {3\over20}k_{\nu-1}\right), & für $\nu\gt 1$,\cr \displaystyle\min\left({1\over10},{61\over50}\hat k_{\infty}\right), & für $\nu=1$\cr } $$

und

$$ f_\nu\gets\cases{ 1/10, & falls $100r_\nu\le\tau$ und $\nu\gt 1$,\cr \displaystyle\min\left({r_\nu\over r_{\nu-1}}, {3\over20}f_{\nu-1}\right), & falls $100r_\nu\gt \tau$ und $\nu\gt 1$,\cr \hat f_\infty, & für $\nu=1$.\cr } $$

D.h. anfangs steigt die Konvergenzrate für das modifizierte Newton-Kantorovich Iterationsverfahren nicht über $1/10$, aber verschlechtert sich jedesmal um 22%, falls in jedem Schritt nur einmal iteriert werden sollte, also stets der $P(EC)$-Modus benutzt wird. Diese in geometrischer Progression ansteigende Verschlechterung dient dazu, daß “gelegentlich” auch der $P(EC)^i$-Modus $(i>1)$ verwendet wird und somit eine “echte” Konvergenzrate berechnet werden kann und nicht lediglich eine Schätzung benutzt werden muß. Bei der Picard Iteration wird die logische Variable atlst2 (at least two) dazu eingesetzt, daß u.U. zweimal iteriert werden muß. Der Hintergrund für die Minima und Maxima liegt in der Vermeidung von Ausreissern. Z.B. dient $\max(\delta_\nu/\delta_{\nu-1},3k_{\nu-1}/20)$ dazu, daß die Konvergenzrate $k_\nu$ für die Newton-Kantorovich Iteration nicht unterhalb 15% ihres Vorwertes fällt. Vermeidet man aus vermeintlichen Effizienzgründen diese Minima und Maxima, so kann man tatsächlich beobachten, daß ein einziger Ausreisser, die Konvergenzrate über zahlreiche Schritte verschlechtern kann, und dies obwohl schon der nächste Zeitschritt wieder brauchbarere Werte liefert.

Diesen Effekt kann man durch Testläufe gut beobachten. In dem Programm TENDLER gibt es eine Option sich diese Werte ausdrucken zu lassen (analog wie \tracingpages in TeX), sodaß man diesen Effekt gut beobachten kann. Diese Option ist aktivierbar ohne Neuübersetzung.

$f_0$ wird nicht gesondert behandelt, da bei einmaligem Iterieren keine Schaltentscheidung möglich ist. Bei zwei Iterationen stehen aber dann genügend viele benötigte Werte zur Verfügung zur Berechnung einer Konvergenzrate. Ganz zu Beginn der Integration werden $k_\nu$ und $f_\nu$ initialisiert mit

$$ k_0\gets{1\over2}, \qquad f_0\gets{1\over2}. $$

Das Programm LSODE benutzt hier den Wert $0.7$ und die Programme STINT, GEAR, EPISODE und zahlreiche weitere, verwenden hier den Wert $1$.

Die Überprüfung, ob $100r_\nu \le \tau$ ist, dient dazu, Besonderheiten durch übermässig kleine Pseudoresiduen nicht überreichlich Gewicht beizumessen und insbesondere eine Division durch Null zu vermeiden. Auch bei Rechnung nahe der Grenze der Maschinengenauigkeit greift dieser Zusatztest. Würde man diesen einfach durchzuführenden und schnell berechneten Test unterlassen, so könnte es passieren, daß eine “zufällig” gebildete große oder kleine Konvergenzrate gebildet würde und dem Konvergenztest seine Kontinuierlichkeit entzieht. Das Setzen auf den Wert $1/10$ besitzt eine gewisse Willkürlichkeit. Da jedoch $100r_\nu\le\tau$ erfüllt ist, liegt es nahe anzunehmen, daß die Iteration konvergiert. Diesen Sachverhalt sollte die Konvergenzrate in angemessener Weise widerspiegeln. Hierbei ist der Wert $1/10$ sogar noch vergleichsweise hoch. Typische Konvergenzraten liegen nicht selten eine Zehnerpotenz darunter. Es ist häufig zu beobachten, daß wenn $100r_\nu\le\tau$ erfüllt ist, daß dann sehr schnelle Konvergenz vorliegt und wenn es dann doch zu zweimaligen Iterieren kommt, die vergleichsweise hohe Konvergenzrate von $1/10$ heruntergezogen wird. Aus Gründen, die mit dem Schalten zusammenhängen, werden “gelegentlich” 2 Iterationen ausdrücklich gefordert.

Bei unkritischem Betrachten der Überlegungen, wann eine Iteration abzubrechen ist, können Mißdeutungen auftreten. Es sei deutlich darauf hingewiesen, daß die obigen Plausibilitätsüberlegungen keinerlei Beweiskraft besitzen, da z.B. durch eine kleine Schrittweite $\mathopen|h\mathclose|$ das Pseudoresiduum nichts über Konvergenz aussagen muß---man denke an die Reihe $\sum(1/n)$. Dennoch zeigt sich, daß im Verein mit den restlichen Segmenten des Programmes TENDLER der Test so arbeitet, wie er geplant wurde. Es sei angemerkt, daß in anderen Programmen, wie z.B. LSODAR, nicht rigoroser vorgegangen wird. Hierbei ist es so, daß das Programm LSODAR sehr viel mehr Aufwand treibt beim Konvergenztest, wie die Programme DIFSUB, GEAR, EPSISODE, STINT und weitere. Allerdings handelt es sich bei den zuletzt genannten Programmen auch um “ältere” Programme. Beim Konvergenztest für das modifizierte Newton-Kantorovich Iterationsverfahren braucht kein Test der Form $100\delta_{\nu-1}<\tau$ verwendet werden, da andernfalls schon vorher Konvergenz gemeldet worden wäre. Bei Picard Iteration kann Konvergenz angezeigt sein und trotzdem wird weiter iteriert aufgrund der Variablen atlst2. Ein Iterationsschritt des modifzierten Newton-Kantorovich Iterationsverfahrens ist für $n>1$ aufwendiger als ein Iterationsschritt mit der Picard Iteration.

Der eigentliche Test ist nun

$$ \delta_\nu k_\nu\le\tau, \qquad\hbox{bzw.}\qquad \delta_\nu f_\nu\le\tau. $$

Danach wird geprüft, ob die Iteration abzubrechen ist, falls einer der beiden Bedingungen

$$ k_\nu\le{47\over100}, \qquad\hbox{bzw.}\qquad f_\nu\le{1\over5} \qquad(\nu\gt 1) $$

verletzt ist. Dieser Test wird also nur durchgeführt, wenn eine “echte” Konvergenzrate gebildet wurde. Im Falle von Konvergenz wird bei der modifizierten Newton-Kantorovich Iteration die Funktion setconv() aufgerufen, welche Einfluß nimmt auf die Statusvariable iweval. Die Funktion setconv() ergreift besondere Maßnahmen bei Erkennung “langsamer” Konvergenz. Bei Picard Iteration wird sofort zum Fehlerkontroll-Segment verzweigt. Die Funktion setconv() wird bei Picard Iteration nicht aufgerufen.

2. Der Wert $\tau$ wird vor jedem Zyklus besetzt mit $\tau\gets\varepsilon/16$. Die Wahl dieser Größe ist kritisch. Der Wert $\varepsilon/10$ ist schon zu “groß”, der Wert $\varepsilon/22$ wiederum zu “klein”. Dies soll etwas genauer erklärt werden. Manchem Programm kann seine Leistungsfähigkeit, seine Effizienz geraubt werden durch eine zu kleine Wahl von $\tau$. Bei dem Programm STINT liegt die Anzahl der Korrektoriterationen pro Schritt bei 1.8 (grob gemittelt). Natürlich hängt dieses Verhältnis ab von der zugrundeliegenden Differentialgleichung, von der geforderten Genauigkeit, von der verwendeten Rechenanlage und weiterem. Bei einer früheren Version des Programmes ODIOUS lag dieser Mittelwert ebenfalls bei 1.8 (erneut grob gemittelt, unter erneuter Mißachtung von sonstigen Abhängigkeiten, w.o.). Bei dieser früheren Version vom Programm ODIOUS war $\tau=\varepsilon/100$. Bei dem Programm LSODE liegt das Verhältnis von Iterationen zu Schritten bei ca. 1.2–1.3. Nun ist aber bei den zyklischen Formeln von Tendler (1973) wenig Grund ersichtlich, warum diese Werte derart differieren sollen, zumal die insgesamt 24 Stufen zu 62.5% aus den BDF bestehen, mithin die Formeln von J.M. Tendler eine matrixwertige Verallgemeinerung der BDF darstellen.

Durch Probeläufe mit dem Programm TENDLER wurden nun festgestellt, daß ohne Einbußen an Genauigkeit, ein Übergang von einem quasi $P(EC)(EC)$-Betrieb auf einen $P(EC)$-Betrieb übergewechselt werden kann, ohne Einbrüche an Schrittzahl, Jacobimatrixauswertungen, Fehlerversagen und dergleichen feststellen zu müssen.

In der Wahl des Wertes $\tau$ differieren viele Programme und dies selbst dann, wenn diese Programme auf den gleichen Formeln basieren. Manche Programme bringen $\tau$ mit dem Diskretisierungsfehler in Verbindung, manche mit der Ordnung des Verfahrens, manche mit beiden, manche mit gar nichts.

Für alle diese Auswahlen gibt es Erklärungen, Plausibilitätsüberlegungen, aber keinen mathematischen Beweis, der aussagt, daß nur diese oder jene Form allen anderen Formen grundsätzlich überlegen ist. Allen haftet letztendlich die gleiche Willkürlichkeit an, wie eben der Wahl $\tau=\varepsilon/16$. Der Wert $\tau=\varepsilon/16$ im Programm TENDLER darf nicht notwendig zur Annahme verleiten, daß zwischen $\tau$ und $\varepsilon$ immer ein solch einfacher funktionaler Zusammenhang besteht. Ebenso denkbar ist, daß Änderungen anderer Segmente letztlich wieder andere Konstellationen günstiger erscheinen lassen. Shampine (1979) schreibt hierzu:

An unsolved basic problem with the acceptance test is the relationship between an appropriate tolerance $\tau$ for $\mathopen|\varrho\mathclose|$ or $\left|y^m-y^*\right|$ and the tolerance $\varepsilon$ for the local truncation error. As far as we know, the relation was selected by experiment in all the current codes. There are two issues here. One is the effect on the remainder of the integration due to computing $y^*$ inaccurately. The other is the effect of inaccurate solution values on other algorithms in the code, such as that for estimating the local truncation error. $\ldots$ Thus it is a waste of effort to take $\tau$ very small compared to $\varepsilon$.

An anderer Stelle, siehe Shampine (1980), hierzu weiter:

We have not said anything about the relationship between the tolerance $\tau$ in $r/(1-r)\left|y^{m+1}-y^m\right|$ and the accuracy $\varepsilon$ desired of the solution of the differential equation. This is a research question whichs needs attention. It is clear that $\tau$ must be smaller than $\varepsilon$. A small $\tau$ is one way to compensate for an inaccurate estimate of $r$ used in the convergence test. However the smaller $\tau$ is made, the more it costs to compute $y^*$. Experiment says that $\tau$ a great deal smaller than $\varepsilon$ does not improve the solution of the differential equation. In current codes the relationship appears to be ad hoc. The choice varies wildly, but $\tau=0.1\varepsilon$ is representative.

Die Simplifizierung der Tests auf $\delta_\nu k_\nu$ bzw. $r_\nu f_\nu$, anstatt $k_\nu/(1-k_\nu)\delta_\nu$ bzw. $f_\nu/(1-f_\nu)r_\nu$, hat seinen Grund in der Begrenzung der Konvergenzraten. Hierbei ist $1-k_\nu\approx1$ bzw. $1-f_\nu\approx1$. Wäre dies nicht der Fall, so wäre im Programm TENDLER ohnehin schon vorher Divergenz gemeldet worden. Man hätte aber auch sonst allen Grund Divergenz anzunehmen.

Man könnte nun angesichts der Wichtigkeit dieses Wertes $\tau$ sich der Mühe unterziehen, mit einem Ansatz der Form

$$ \tau = {a\varepsilon + b\over cp^2+dp+e}, $$

beispielsweise, ggf. unter Berücksichtigung der Beträge der Fehlerkonstanten o.ä., und durch eine Batterie von Testdifferentialgleichungen mit Hilfe eines Programmes zur Lösung von Optimierungsaufgaben, unter Nebenbedingungen ggf., die günstigsten Werte $a$, $b$, $c$, $d$ und $e$ zu ermitteln versuchen. Z.B. hängt in den Programmen DIFSUB und STINT der Wert $\tau$ von der Ordnung $p$ ab. Jedoch muß man berücksichtigen, daß auch die Wahl von Sicherheitsfaktoren, wie sie bei der Fehlerkontrolle und der Schrittweiten- und Ordnungssteuerung auftreten, Rückwirkungen haben. Beispielsweise konnte folgendes beobachtet werden: Ist $\tau$ zu groß, so wird der Fehlerkontrolle fast gänzlich alleine die Entscheidung über die Güte einer Näherung aufgebürdet, desto häufiger versagt sie. Die Korrektoriteration hatte nicht genügend Zeit, die “Steifheit wegzudämpfen”; für eine Präzisierung des letzten siehe Shampine (1980).

Diesen Effekt kann man teilweise entgegenwirken, indem man die Fehlerkontrolle sehr restriktiv macht, also schon sehr früh versagen lässt, bevor dies irgend etwas mit Steifheit zu tun haben könnte. Durch Wahl geeigneter Sicherheitsfaktoren ist dies leicht zu bewerkstelligen. Man befindet sich in der gleichen mißlichen Situation, daß eine Minimierung der Funktion $f(x_1,x_2,\ldots,x_n)$ in der Variablen $x_1$, noch lange keine Minimierung der Funktion $f$ schlechthin zu bedeuten hat, d.h. eine Minimierung in allen Variablen. An dem letzten ist man ja eigentlich interessiert, nämlich dem Maximum an Effektivität eines Differentialgleichungslösers.

Eine weitere Bestätigung für die Wahl des $P(EC)$-Modus als dominierende Betriebsart, erhält man, wenn man das Programm TENDLER im Newton-Kantorovich- bzw. Picard-Modus, mit fest vorgegebener Anzahl der Iterationen, benutzt. Hier zeigen Testläufe, daß die globalen Fehler sich nicht merklich ändern, allein einzig das Verhältnis von Funktionsauswertungen zu Schritte, hat fast genau den Wert erreicht, denn er auch theoretisch erreichen müßte, also $i$, für die $P(EC)^i$-Betriebsart, und 2, für die $P(EC)E$-Betriebsart.

Diese empirisch gewonnene Erkenntnis deutet an, daß der Sachverhalt, der vom Startfehlersatz von Liniger (1971) ausgedrückt wird, allgemeiner gilt, als der ihn tragende Beweis. Der Beweis des Startfehlersatzes von Liniger (1971) ist an Bedingungen geknüpft, die im Falle steifer Differentialgleichungen nicht zutreffen.

Dieser Sachverhalt ist jedoch aus anderem Zusammenhang gut bekannt. Die Formeln, die Tischer (1983)2 (PhD, siehe unten) und Tischer/Sacks-Davis (1983)3 vorstellen, sind beweistechnisch geknebelt an die Bedingung, daß $\left|hDK_2\right|<1$ ist [$D=\sup\left|A^i\right|$, $K_2$ Lipschitzkonstante bzgl. $v$ bei $\varphi(u,v)$]. Die Tatsache, daß das Programm ODIOUS, welches auf diesen zyklischen Formeln von Tischer (1983)2 und Tischer/Sacks-Davis (1983)3 basiert, in seiner Leistungsfähigkeit bei den DETEST-Problemen dem Programm LSODE in fast nichts nachsteht, ja sogar bei der besonders steifen Testdifferentialgleichung B5, seine gänzliche Überlegenheit zeigt, siehe Tischer/Gupta (1985)2, macht deutlich, daß die Aussagen, wie sie die derzeitige Theorie macht, allgemeiner gelten, als die sie stützenden Beweise.

Bibliographisch:

  1. 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.
  2. Tischer, Peter E. und Gupta, Gopal K.: “A Cyclic Method Stiff ODE Solver”, Technical Report 38, Department of Computer Science, Monash University, Clayton, Victoria, Australia, June 1983, i+45 S.
  3. 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
  4. Tischer, Peter E. und Gupta, Gopal K.: “An Evaluation of Some New Cyclic Linear Multistep Formulas for Stiff ODEs”, ACM TOMS, Vol 11, No 3, September 1985, pp.263–270

3. Der Vorteil der Benutzung von MNI, FIX und FIXPECE, also den Iterationsarten mit fest vorgegebener Anzahl an Iterationen, liegt darin, daß kein Konvergenztest durchgeführt wird, keine Konvergenzraten gebildet werden (somit keine Schaltmöglichkeit besteht), insbesondere keine Normauswertung während der Korrektoriteration erfolgt. Nun ist der Einfluß der Normberechnungen erwartungsgemäß zwar nicht alles dominierend, jedoch stellt er einen bei flüchtiger Betrachtung häufig unterschätzten Beitrag dar. %bisweilen Der Anteil an der Gesamtrechenzeit variiert selbstverständlich. Er hängt vom Differentialgleichungsproblem ab und auch von der Wahl der Norm schlechthin. Z.B. ist $\left|{}\cdot{}\right|_\infty$ i.a. signifikant billiger zu berechnen als $\left|{}\cdot{}\right|_2$. Bei kleindimensionalen Differentialgleichungsproblemen beherrscht das Drucken oder Plotten die Gesamtrechenzeit.

Nicht untypisch ist jedoch, daß die Funktionsauswertungen ca. 5–10% bei einfach auszuwertenden Funktionen und 10–20% bei aufwendigeren Funktionen ausmachen und die Rücksubstitutionen und Normberechnungen jeweils ca. 10%. Dies wird anfänglich möglicherweise überraschen. Bedenkt man jedoch, wie häufig Normberechnungen notwendig sind, so werden die Zahlen verständlich. Nocheinmal sei betont, daß die Werte 20% und 10% nicht feststehende Zahlen sind, sondern variieren, u.U. nicht unerheblich, jedoch einen typischen Bereich abstecken. Z.B. bei den Problemen B2 bis B5 benötigen die Normberechnungen insgesamt ca. 10%, die Rücksubstitutionen ca. 8% und die Funktionsauswertungen ca. 3% der Gesamtrechenzeit. Bei dieser Schar von Differentialgleichungen ist die Auswertung der Funktion sehr einfach. Es werden in diesem Falle ca. zwei bis dreimal soviele Normauswertungen durchgeführt wie Funktionswerte bestimmt werden.

Bei den zwei Iterationsarten MNI und FIX kann selbstverständlich nach jedem Zyklus die Anzahl der Iterationsschritte geändert werden, genauso, wie prinzipiell nach jedem Zyklus der Benutzer jede beliebige Iterationsart anwählen kann, also MNI, MNIFLOAT (Newton-Kantorovich mit Konvergenztest), FIX, FIXFLOAT (Picard mit Konvergenztest) und FIXPECE. Natürlich kann auch die maximale Anzahl der Iterationen für MNIFLOAT und FIXFLOAT vom Benutzer auf jeden beliegen Wert größer Null gesetzt werden. Vom Programm TENDLER werden hier gewisse Voreinstellungen wirksam.

4. Beispiel: Bezüglich der Konvergenzraten beachte man, daß es sehr wohl sein kann, daß

$$ \mathopen|y^2-y^1\mathclose| \gg \mathopen|y^1-y^0\mathclose| $$

und dies obwohl die Folge der Iterierten konvergiert, sogar sehr schnell. Diese grundsätzliche Problematik sieht man beispielsweise für

$$ y^i = \cases{1/2^i, & für $i\mathbin\%13=0$,\cr 1/10^i, & für $i\mathbin\%13\ne0$,\cr } \qquad (i\ge0). $$

Darüberhinaus erkennt man sogar, daß unendlich häufig gilt

$$ \left|y^{13\nu}-y^{13\nu-1}\right| \gg \left|y^{13\nu-1}-y^{13\nu-2}\right|, \qquad (\nu\gt 0). $$

Den Quotienten $\mathopen|y^{\nu+1}-y^\nu\mathclose| / \mathopen|y^\nu-y^{\nu-1}\mathclose|$ kann man natürlich beliebig groß gestalten, durch passende Änderungen. Es muß also zuerst geprüft werden, ob der Betrag des Residuums oder der Betrag des Pseudoresiduums ausreichend klein ist. Erst danach ist zu überprüfen, ob die Konvergenzrate Anlaß zur Besorgnis gibt. Auf der anderen Seite gibt Shampine (1980) ein Beispiel einer steifen Differentialgleichung an, wo die Dinge genau umgekehrt liegen, also erst das Pseudoresiduum und danach das Residuum getestet werden sollte. Den Term $1/2^i$ kann man auffassen als entstandenen Rundungsfehler. Die Zahl 13 ist hier natürlich völlig willkürlich gewählt.

In dem obigen Beispiel ist $(y^\nu)_{\nu=0}^\infty$ keine kontrahierende Folge. Jedoch ist man nicht an Kontraktion per se interessiert, sondern an Konvergenz. Kontraktion und gewisse Selbstabbildungseigenschaften sind lediglich eine hinreichende Bedingung für Konvergenz; notwendig sind diese Bedingungen für Konvergenz nicht.

5. Der Divergenzfall: Wird nun der oben beschriebene Konvergenztest nicht bestanden, so kommen folgende Strategien zur Anwendung. Eine geschickte Wahl dieser Strategie ist wichtig für die Effizienz des Programmes und das Weiterrechnen nach einem Fehlschlag.

Trotz des Fehlens einer abgeschlossenen mathematischen Theorie, verlangen die Heuristiken große Aufmerksamkeit.

Den obigen Fall muß man eher als sehr große Ausnahme werten. Häufiger ist der nachstehende Fall, obwohl natürlich auch dieser Fall eher selten auftritt.

Die maximale Anzahl der Iterationen ist in den einzelnen Iterationsarten begrenzt durch die Variablen mnimaxit und fixmaxit. Es ist auch möglich, daß der Benutzer die Anzahl der Iterationen fest vorgibt. Dann werden jedoch andere Unterprogramme aufgerufen. Um die ohnehin komplexen Konvergenztestuntersuchungen nicht mit weiteren Abfragen zu überfrachten, werden für sich getrennte Unterprogramme aufgerufen, die zudem die speziellen Bedürfnisse der Iterationsart berücksichtigen und in dieser Hinsicht optimiert sind. In der Funktion setiweval() wird der Hindmarsh-Test vollzogen.

Für das modifizierte Newton-Kantorovich Iterationsverfahren sieht die grundsätzliche Steuerlogik wie folgt aus:

6. Singularität der Iterationsmatrix: Im Falle einer nicht-invertierbaren Iterationsmatrix $W$ wird einfach weiterverfahren, und zwar so als ob Divergenz gemeldet worden wäre. Zwar wird von mnifloat() als auch mni() ein spezielles Zeichen gesetzt, welches Singularität der Iterationsmatrix $W$ anzeigt, doch wird dieses augenblicklich nicht gesondert verarbeitet. Prinzipiell besteht zwischen Divergenz und Nichtinvertierbarkeit ein großer Unterschied. Sei $n$ die Dimension der zu lösenden Differentialgleichung. Es kann nun nur maximal $n$ verschiedene Eigenwerte des Eigenwertproblems

$$ J{\bf x} = \lambda{\bf x}, \qquad\lambda = {1\over h\gamma} $$

geben. Man beachte die Form der Iterationsmatrix $W=I-\lambda^{-1}J$. Kennt man das Spektrum der Jacobimatrix $J$ im voraus, so kann man sofort die “verbotenen” Werte $h_i\gamma_i=\lambda_i^{-1}$ angeben, für $i=1,\ldots,n$. I.d.R. kennt man natürlich das Spektrum von $J$ nicht, zumal $J$ von der Zeit $t$ abhängen kann.

Weiterhin sieht man, daß die Drittelung der Schrittweite sofort erneut zu einer nicht-invertierbaren Iterationsmatrix $W$ führen kann, wie man anhand der speziellen Jacobimatrix $J=\bigl({1\atop0}{0\atop3}\bigr)$ und einer Schrittweitendrittelung ersehen kann.

Umgekehrt ist es gut möglich, daß eine Vergrößerung der Schrittweite $h$ eine viel angemessenere Reaktion auf Nichtinvertierbarkeit darstellen könnte. Es gibt kein Programm, welches dies tut, dennoch wäre ein solches Verhalten während der Startphase nicht abwegig, wenn man gewisse Zusatzinformationen der Schrittweiten- und Ordnungssteuerung mit berücksichtigen würde, was selbstverständlich einen erhöhten programmiermässigen Aufwand darstellen würde. Aus Sicherheitsgründen verkleinert man natürlich lieber die Schrittweite $h$, vergößern kann man sie dann anschliessend immer noch.

Vorstellbar ist dabei jedoch das folgende Szenarium mit der Jacobimatrix $J=\bigl({1\atop0}{0\atop5}\bigr)$. Sei $h=1$, $W$ sei nicht invertierbar, $h$ wird gedrittelt, in der Schrittweiten- und Ordnungssteuerung dann jedoch später wieder verdreifacht, $\ldots$ . Man erkennt zusätzlich, daß eine Veränderung der Ordnung $p$ alleine ebenfalls schon eine mögliche Quelle für Nichtinvertierbarkeit der Iterationsmatrix $W$ darstellt, falls $W$ refaktorisiert wird. Die hier beinahe sehr in das Detail gehenden Überlegungen spielen allerdings bei praktischen Differentialgleichungen nicht die Rolle, die man vermuten könnte. Statistiken zeigen, daß Nichtinvertierbarkeit der Iterationsmatrix $W$ so gut wie gar nicht vorkommt: es wurde bisher noch nie beobachtet. Dies liegt jedoch z.T. daran, daß die $LU$-Zerlegung im Programm TENDLER keine scharfe Rangentscheidung trifft, sondern einfach u.U. mit einer sehr schlecht konditionierten Matrix weiterrechnet. Hierbei stellt man sich vor, daß die Kondition einer Matrix ihre “Nähe” zur Nichtinvertierbarkeit beschreibt. Diese Deutung kann man streng präzisieren, obwohl sie teilweise irreführend sein kann, siehe Wilkinson, J.H., Wilkinson (1965). Eine Rangentscheidung ist aber auch als solches ein numerisch diffiziles Problem, da die Menge der nichtinvertierbaren Matrizen eine abgeschlossene Nullmenge in der Menge aller Matrizen ist.

7. Es ist wichtig zu vermerken, daß erst versucht wird bei Konvergenzversagen die Iteration bei derselben Stufe neu zu starten und zwar mit aufgefrischter Iterationsmatrix $W$. Erst wenn dies nicht mehr möglich ist, d.h. wenn $W$ schon frisch berechnet wurde, erst dann die Schrittweite $h$ oder die Ordnung $p$ verändert werden. Im letzteren Falle wird der gesamte Zyklus wiederholt und nicht nur die betreffende Stufe. Programmiermässig wäre es genauso gut möglich nur an der besagten Stufe neu zu beginnen, jedoch hätte man dann den angefangenen, nicht zu Ende geführten Zyklus, als Zyklus für sich zu betrachten. Dieser Teilzyklus hat natürlich ein anderes Stabilitätverhalten, insbesondere ist nun überhaupt nicht mehr garantiert, daß dieser Teilzyklus $A_\infty[\alpha]$- bzw. $S_\infty[\delta]$-stabil für gewisse $\alpha$ und $\delta$ ist. Beispielsweise sind die ersten beiden Stufen des Vefahrens 7.ter Ordnung die BDF7, welche nicht $D$-stabil sind und schon erst recht nicht $A_\infty[\alpha]$-stabil. Allerdings bedeutet ein Wechsel der Schrittweite $h$ oder der Ordnung $p$ ohnehin einen Wechsel zu einem mehrstufigen Verfahren, nicht notwendig zyklisch.

Die Programme DIFJMT, STINT und ODIOUS sind ebenfalls, wie das Programm TENDLER, vollständig zyklusweise orientiert. Dennoch ist es wichtig den folgenden Sachverhalt deutlich zu machen: Das Wiederholen eines Zykluses, inklusive der dazu benötigten Interpolation, ist nicht aufwendig, vorausgesetzt, daß dies nicht ständig passiert. Aufwendiger ist, daß eine starke Pro- oder Degression der Schrittweite $h$ zu verändertem Konvergenzverhalten der Korrektoriteration führt, insbesondere ist u.U. die Iterationsmatrix $W$ erneut zu faktorisieren. Ein Wiederholen einer einzigen Stufe innerhalb der Formeln von Tendler (1973) setzt eine Differenzentabelle, oder eine sonstige extra Verwaltung für die skalierten Ableitungswerte $z_{m\ell+i}=hf_{m\ell+i}$, voraus. Dies geschieht weder in den Programmen DIFJMT und STINT, noch in dem Programm TENDLER. Denn dies würde zusätzlichen Verwaltungsaufwand erfordern und zusätzlichen Speicherplatz beanspruchen, von grundsätzlichen theoretischen Bedenken einmal ganz abgesehen, welche gegen ein solches Vorgehen sprechen. Die skalierten Ableitungswerte $z_{m\ell+i}=hf_{m\ell+i}$ werden innerhalb eines Zykluses lediglich als Hilfsgrößen betrachtet und bei Beginn eines neuen Zykluses einfach überschrieben.

Das Programm TENDLER macht von den $z_{m\ell+i}$ bzgl. der Interpolation keinerlei Gebrauch. Das Programm STINT verwendet die Werte $z_{m\ell+i}$ bei der Interpolation, aus dem einfachen Grund, weil eine entsprechend hohe Differenz $\nabla^{k+1}y$ überschrieben wurde und nun restauriert werden muß. Dies geschieht unter Zuhilfenahme der Beziehung

$$ \nabla^py\approx p\cdot\left(hf-\sum_{i=1}^{p-1}{1\over i}\nabla^iy\right), $$

wegen

$$ {1\over h}\sum_{i=1}^\infty {1\over i} \nabla^i y = \dot y. $$

Es sei $E$ der Shiftoperator und $D$ der Differentiationsoperator. Die obige Gleichung ist richtig, wegen

$$ E = e^{hD} = {1\over1-\nabla}, \qquad\hbox{also}\qquad hD = -\ln(1-\nabla) = \nabla+{\nabla^2\over2}+{\nabla^3\over3}+\cdots{} , $$

bei rein formaler Rechnung, die sich aber streng begründen lässt. Die obige Approximation ist in dem Programm TENDLER vollständig überflüssig, aufgrund einer verschiedenen Organisation der Differenzentabelle der einzelnen Stufen. In dem Programm STINT ist diese Rechnung nötig bei einem zurückgewiesenen Zyklus, falls kurz vorher die Schrittweiten- und Ordnungssteuerung, oder der Benutzer, eine Veränderung der Schrittweite erzwangen, insbesondere falls vorher die Ordnung höher war.

Weiterhin ist das Durchlaufen dreier Stufen in einem Zyklus vom Wesen her u.U. einfacher und effizienter, als die dreifache Ausführung eines einstufigen Verfahrens. Dies wurde auch durch empirische Untersuchungen z.T. erhärtet, siehe Gaffney (1984). Schließlich sind Zurückweisungen ohnehin die Ausnahme, denn es ist eine der Hauptzwecke der Schrittweiten- und Ordnungssteuerung genau dies zu erreichen. Dennoch ist i.d.R. der benötigte Speicherplatz für zyklische lineare Mehrschrittverfahren in linearen Termen der Dimension $n$ der Differentialgleichung höher, als bei einstufigen linearen Mehrschrittverfahren. Dies natürlich i.d.R. nur, falls die Schrittweite geändert werden soll, sonst gilt dies so ohne weiters nicht.

5. Verhinderung langsamer Konvergenz in TENDLER

I appeal to my colleagues to remember that in any field of numerical computation, technical papers are often less influential than computer programs!

Ll.N. Trefethen (1985)

Nach erfolgreichem Abschluß der Korrektoriteration, wird die Funktion setconv() aufgerufen, welche anhand der Anzahl der benötigten Iterationen und aufgrund vergangener Konvergenzbeobachtungen, mitentscheidet, ob und wie, eine Refaktorisierung der Iterationsmatrix $W=I-h\gamma J$ vorzunehmen ist. Dies geschieht in Form von Ansammeln von Minuspunkten, man vgl. hier auch Knuth (1986), PDF.

Gelingt die Konvergenz des Korrektors gleich im ersten Iterationschritt, so werden alle bisher angesammelten Minuspunkte vergessen, da jetzt angenommen wird, daß die Konvergenz genügend schnell stattfindet. Bei drei Iterationen bis zur Konvergenz, oder aber bei “großer” Konvergenzrate, werden “viele” Minuspunkte aufgerechnet, andernfalls werden “wenige” Minuspunkte aufgesammelt.

Überschreitet nun die Summe der Minuspunkte eine gewisse Grenze an maximal möglichen Minuspunkten, nämlich MAXSLOW, so wird die Sperrvariable problem aktiviert und die Jacobimatrix $J$ wird nun berechnet und anschliessend die Iterationsmatrix $W=I-h\gamma J$ refaktorisiert. Nebenläufig sei erwähnt, daß die Funktion setconv() selber keine Berechnungen oder Faktorisierungen durchführt, sondern nur entsprechende Markierungen setzt.

Wie üblich wird für statistische Zwecke gezählt, wie häufig durch diese Maßnahme Einfluß auf die Iterationsmatrix $W$ und damit auf die Konvergenz des modifizierten Newton-Kantorovich Verfahrens, genommen wird. Gleichzeitig kann damit erkannt werden, wie wirkungsvoll dieser Teil der Konvergenzüberwachung arbeitet. Derzeitige Erfahrungen lassen erkennen, daß dieser Test bei vernachlässigbaren Kosten, eine lohnenswerte Veränderung darstellt.

Die augenblickliche Gewichtung der Minuspunkte ist wie folgt: 9 Minuspunkte für dreifache Iteration, 9 Minuspunkte für zweifache Iteration mit einer Konvergenzrate größer 0.09 und sonst bei lediglich zweifacher Iteration ein einziger Minuspunkt. Ein Eingriff, also ein Verändern der Sperre problem und der Statusvariablen iweval, geschieht bei insgesamt 18 Minuspunkten.

6. Der Hindmarsh-Test im Programm TENDLER

Hindmarsh, Alan C. Dieser Test ist natürlich nicht von Hindmarsh originär neu erfunden worden, jedoch um einen Namen zu haben, sei dieser Test im weiteren so genannt. Hindmarsh hat dennoch diesen Test in vielen Punkten verfeinert und in dem Programm GEAR erfolgreich angewendet. Um auch ein theoretisches Verständnis für diesen Test zu erhalten, folgen nun einige Bemerkungen.

Es soll hier auch um die Beantwortung der folgenden beiden Fragen gehen:

  1. Gegen welchen Wert konvergiert das modifizierte Newton-Kantorovich-Verfahren, wenn man unendlich lange iterieren würde, bei genügend kleiner Differenz $\mathopen|H-\tilde H\mathclose|$ und
  2. in welcher Beziehung steht dieser Grenzwert mit dem Werten aus vorhergehenden Schritten?

Wendet man sich dem skalaren und autonomen Testproblem $\dot y=\lambda y$ zu und benutzt man als Prädiktor das explizite Euler-Verfahren

$$ y_{n+1} = y_n + hf_n $$

und als Korrektor das implizite Euler-Verfahren

$$ y_{n+1} = y_n + hf_{n+1}, $$

und behandelt man das dabei auftretende Nullstellenproblem, welches in diesem Falle sogar linear ist, mit dem modifizierten Newton-Kantorovich-Verfahren

$$ y^{\nu+1}_{n+1} = y^\nu_{n+1} - W^{-1}Fy_{n+1}^\nu, \qquad\hbox{mit}\quad Fy_{n+1} = (1-H)y_{n+1} - y_n, $$

so gelten die Ableitungen, wie untenstehend. Man beachte, daß natürlich das unmodifizierte Newton-Kantorovich-Verfahren die Lösung des Nullstellenproblem $Fy_{n+1}=(1-H)y_{n+1}-y_n=0$, in einem Iterationsschritt, bis auf Rundungsfehler sofort exakt liefern würde. Die Iterationsmatrix $W$ ist hier $W=1-\tilde H$. Man berechnet nun für das modifizierte Newton-Kantorovich Iterationsverfahren

$$ y_{n+1}^1 = y_{n+1}^0 - {(1-H)y_{n+1}^0 - y_n\over1-\tilde H} % = \biggl(1 - {1-H\over1-\tilde H}\biggr)y_{n+1}^0 + {1\over1-\tilde H}y_n = {H-\tilde H\over1-\tilde H}y_{n+1}^0 + {1\over1-\tilde H}y_n. $$

Dies zeigt schon im wesentlichen folgendes: Der Hindmarsh-Test ist sinnvoll, da er testet in wie weit sich $H$ und $\tilde H$, also $h\gamma$ und $\tilde h\tilde\gamma$ entfernt haben, seit der letzten Refaktorisierung der Iterationsmatrix $W$, hier also der Reziprokenbildung. In den Programmen GEAR, EPISODE und LSODE, in denen keine Jacobimatrix $J$ separat gespeichert wird, kann man also aus $H-\tilde H$, den Faktor $J$ herausfaktorisieren. In dem Programm TENDLER ist dieses Ausklammern, aufgrund der gesonderten Speicherung von $J$, nicht so einfach möglich. Hier sind dann weitere Überprüfungen nötig, z.B. durch Betrachten der Konvergenzgeschwindigkeit, was auch tatsächlich getan wird.

Interessanter ist nun, daß es sich bei der obigen Gleichung um eine eingliedrige Rekurrenzgleichung handelt der Form $a_{n+1}=pa_n+q$, mit der allgemeinen Lösung

$$ a_{n+1} = p^{n+1}a_0 + \sum_{i=0}^n p^{n-i}q. $$

Man ersetze 1 durch $\nu+1$ und 0 durch $\nu$. Der Wert

$$ p = {H-\tilde H\over1-\tilde H}, $$

ist für genügend kleinen Zähler stets kleiner als eins zu kriegen und es gilt dann näherungsweise

$$ \sum_{i=0}^n p^i = {1-p^{n+1}\over1-p} \approx 1, \qquad\hbox{falls}\quad p\approx 0. $$

Aufgrund dieser Herleitung wird nun verständlich den Wert

$$ q = {1\over1-\tilde H}y_n $$

grob als Näherung für die Lösung der impliziten Gleichung mit dem Newton-Kantorovich Iterationsverfahren anzusehen, unter den oben erwähnten, stark vereinfachenden Voraussetzungen.

Der Fall $W=I$, also Picard Iteration, führt direkt auf

$$ y^{\nu+1}_{n+1} = Hy^\nu_{n+1}+y_n $$

und damit auf

$$ y_{n+1}^{\nu+1} = H^{\nu+1}y_{n+1}^0 + \biggl(\sum_{i=0}^\nu H^{\nu-i}\biggr)y_n, \qquad\hbox{also}\qquad y_{n+1}^{\nu+1}\approx H^{\nu+1}y_{n+1}^0 + {y_n\over1-H}. $$

Man ersieht natürlich sofort, daß $\mathopen|H\mathclose|<1$ eine notwendige und hinreichende Bedingung für Konvergenz ist. Beachtenswert ist hier, daß der erste Iterationsschritt genauer geworden wäre, bei Verwendung eines nicht konsistenten Prädiktors. Dies gilt natürlich ersteinmal nur für diese Testgleichung.

Die Verallgemeinerung auf den $n$-dimensionalen Fall mit $\gamma\ne0$, also für ein echt implizites Verfahren, ist nun einfach. Der Vollständigkeit halber sei dieses Ergebnis hier kurz aufgeschrieben.

Hat man vorgelegt $y_{n+1}=\gamma hf(y_{n+1})+\psi$ und betrachtet man die Testgleichung $f=Jy$, mit einer Diagonalmatrix $J$, so geht man analog wie oben vor. Man setzt $H=hJ$, $\tilde H=\tilde h\tilde J$ und die Iterationsmatrix $W$ ist $W=I-\tilde\gamma\tilde H$. Das zu lösende Gleichungssystem ist weiter $Fy_{n+1}=(I-\gamma H)y_{n+1}-\psi$ und man erhält dann

$$ y^{\nu+1}_{n+1} = W^{-1}(\gamma H-\tilde\gamma\tilde H)y^\nu_{n+1} + W^{-1}\psi. $$

Hier wäre $p=W^{-1}(\gamma H-\tilde\gamma H)$ und $q=W^{-1}\psi$.

Für die Picard Iteration, also $\tilde H=\bf 0$ und somit $W=I$, ergibt sich

$$ y^{\nu+1}_{n+1} = \gamma Hy_{n+1}^\nu + \psi $$

und daher

$$ y^{\nu+1}_{n+1} = (\gamma H)^{\nu+1}y_{n+1}^0 + \biggl(\sum_{i=0}^\nu (\gamma H)^i\biggr)\psi \overset{(\nu\to\infty)}{\rightarrow} (I-\gamma H)^{-1}\psi. $$

Konvergenzbedingung für die Neumannsche Reihe ist natürlich $\rho(\gamma H)<1$, welches für genügend kleine Schrittweite $h$ immer zu erreichen ist.

Biographisch:

  1. Neumann, Carl Gottfried (1832–1925)

Um den 40%/100%-Test in dem Programm TENDLER besser zu verstehen, seien hier zur Illustriation genauere numerische Werte geliefert. Für die BDF sieht man in der nachstehenden Tabelle, wie sich aufeinander folgende Werte $\gamma_p$, für $p=1,\ldots,6$ verhalten.

$p$ $\gamma_p$ $\gamma_p $ $\gamma_{p+1}/\gamma_p$ $\gamma_p/\gamma_{p+1}$ $1.1\gamma_{p+1}/\gamma_p$ $1.1\gamma_p/\gamma_{p+1}$
1 1 1 0.67 1.50 0.73 1.65
2 0.67 2/3 0.82 1.22 0.90 1.34
3 0.55 6/11 0.88 1.15 0.97 1.27
4 0.48 12/25 0.91 1.10 1.00 1.21
5 0.44 60/137 0.93 1.07 1.03 1.18
6 0.41 60/147 n/a n/a n/a n/a

Es wird jetzt deutlich, wie in gewisser Hinsicht günstig der Wert 30% bei dem Hindmarsh-Test im Programm GEAR gewählt ist. Man beachte, daß wenn die Schrittweite von der Schrittweiten- und Ordnungssteuerung geändert wird, sich diese um mindestens 10% geändert hat. Dieser Mindestanstieg der Schrittweite wird auch in den meisten anderen Programmen, basierend auf den BDF, wie DIFSUB, EPISODE, GEAR und LSODE, gefordert. Für die zyklischen Formeln von J.M. Tendler erhält man leicht anders aussehende Verhältnisse. U.a. ist die Mindestschrittweitenvergrößerung, die in der Schrittweiten- und Ordnungssteuerung gefordert wird, leicht größer.

7. Der Konvergenztest im Programm STINT

Nachdem nun dargestellt wurde, wie in dem Programm TENDLER der Konvergenztest durchgeführt wird, ist es von Interesse, wie in anderen Programmen dieser Test durchgeführt wird. Zudem wurde auch schon deutlich gemacht, daß es Kritikpunkte an bisherigen Konvergenztests gab, hier ist wieder Shampine (1980) aufzuführen. Auch in dem Program STINT ist der Konvergenztest ein einflußreiches Segment, welches dazu dient festzustellen, wann die Korrektoriteration abgebrochen werden kann. Da jeder Iterationsschritt eine Funktionsauswertung und eine Normauswertung benötigt, stellt er einen rechenintensiven Part des gesamten Programmes dar.

Bei diesem Test wird nun geprüft, ob die gewichtete Korrektur, also das Pseudoresiduum für die $z$-Iterierten, genügend klein ist. Sei nun $\delta_{m\ell+i}^\nu$ die gewichtete Norm dieses Pseudoresiduums im $\nu$-ten Iterationsschritt, also

$$ \delta_{ml+i}^\nu = \mathopen\|z^{\nu+1}_{m\ell+i} - z^\nu_{m\ell+i}\mathclose\|_w = \left\| W^{-1} \left[ hf(t_{ml+i}, {\mskip 3mu}y_{ml+i}^\nu) - z_{ml+i}^\nu \right] \right\|_w, $$

wobei $y_{ml+i}^\nu$ und $z_{ml+i}^\nu$ entsprechend die $\nu$-ten Iterierten von $y_{ml+i}$ bzw. $h\dot y_{ml+i}$ sind und $|\cdot |_w$ die gewichtete RMS-Norm bezeichnet. In dem Programm STINT wird die gewichtete RMS-Norm benutzt zu

$$ \mathopen\|y\mathclose\|_w = \sqrt{{1\over n}\sum_{i=1}^n \left({y_i\over w_i}\right)^2}. $$

$z^\nu$ wird hierbei durch Erfülltsein der Formel $y=\gamma z+\psi$ berechnet, also $z=(y-\psi)/\gamma$. Polynomiale Extrapolation wird für die Bestimmung der nullten Iterierten verwendet. Diese Aufgabe wurde schon vom Prädiktor übernommen. Sei nun weiter

$$ c_\nu = \max\left\{ {9 \over 10}c_{\nu-1},{\mskip 3mu} \left(\delta_{ml+i}^\nu\over\delta_{ml+i}^{\nu-1}\right)^2 \right\}, \qquad\hbox{mit}\quad\nu \gt 1 $$

die modifizierte Konvergenzrate zwischen der $(\nu-1)$-ten und $\nu$-ten Iterierten. Es wird der folgende empirisch bestimmte Faktor

$$ \tau = {1 \over 2} \cdot {1 \over c_{1,p+1}/\alpha_{11}} \cdot {\varepsilon \over p+2}, $$

im Konvergenztest verwendet, wobei $\varepsilon$ die vom Anwender vorgegebene Fehlertoleranz pro Schritt ist und $p$ die Ordnung des benutzten Verfahrens ist. Es ist

$$ c_{i,p+1} = {1\over (p+1)!} C_{p+1,k}\pmatrix{\alpha\cr \beta\cr}, $$

der Fehlerfaktor der ersten Stufe. Die Vektoren $\alpha$ und $\beta$ kennzeichnen wieder die Koeffizienten des linearen Mehrschrittverfahrens (pro Stufe), und $C_{p+1,k}\in\mathbb{Z}^{(p+2)\times(2k+2)}$ ist die Konsistenzmatrix.

Es wird nun angenommen, daß der Korrektor konvergierte, falls

$$ \delta_{ml+i}^\nu{\mskip 5mu}\cdot{\mskip 5mu}\min\{ 1,{\mskip 5mu}\sqrt{2c_\nu} \} \le \tau. \tag{Kgz} $$

Man beachte, daß hierbei $c_0$ die Konvergenzrate der vorherigen Stufe des Zykluses ist. Wenn jetzt die obige Bedingung (Kgz) erfüllt ist, so geht man davon aus, daß Korrektorkonvergenz erreicht wurde oder doch zumindestens im nächsten Schritt erreicht wird. Auf jeden Fall wird dann das modifizierte Newton-Kantorovich Verfahren abgebrochen.

Bei der Überprüfung der Norm des Pseudoresiduums $\delta_{m\ell+i}^\nu$ geht man davon aus, daß lineare Konvergenz vorliegt.

Der Faktor $\sqrt2$ im Konvergenztests (Kgz) wird nur deswegen angewendet, um sicherzustellen, daß im nächsten Iterationsschritt die Ungleichung $\delta_{ml+i}^{\nu+1}\le\tau$ auf alle Fälle erfüllt ist. Dadurch, daß der nächste Iterationsschritt nicht mehr durchgeführt wird, erspart man sich weiteres Auswerten der Funktion $f$.

Weiter wird in der Ungleichung $\tau$ und nicht $\varepsilon$ verwendet. Dies dient dazu die Gleichung etwas genauer zu lösen als gefordert. Die Wahl von $\tau$ geht auf C.W. Gear zurück und wurde beibehalten, da sie zufriedenstellend arbeitet.

Der Nenner des zweiten Faktors bei $\tau$ wurde deswegen eingefügt, weil in (Phase 2) der Fehlerkontrolle dieser Term ebenfalls durchgelassen wird. Daher sollte der Korrektor diese Fehler auch tolerieren. In der Phase 1 der Fehlerkontrolle wird jedoch nicht der Quotient ${c_{1,p+1}/\alpha_{11}}$ verwendet, sondern $\delta_{\ell-1,p+1}$. Da die Kehrwerte $\alpha_{11}/c_{1,p+1}$ stets größer sind als die Werte $\delta_{\ell-1,p+1}$ und zwar für alle Ordnungen $p=1,\ldots,7$, fällt damit der Faktor $\tau$ korrekt aus. Ausserdem gilt

$$ {\alpha_{11}\over c_{1,p+1}} \gt {\alpha_{ii}\over c_{i,p+1}}, \qquad i=2,3,4, \quad p=1,\ldots,7, $$

außer für die Ordnung $p=6$. Hier gilt

$$ {\alpha_{11}\over c_{1,p+1}} \lt {\alpha_{33}\over c_{1,p+1}},\qquad p=6. $$

Jedoch sind in diesem Ausnahmefall die beiden Werte in den Vorkommastellen gleich und unterscheiden sich lediglich geringfügig in der ersten Nachkommastelle.

Im Programm STINT kommen u.a. die folgenden Strategien zur Anwendung:

  1. Es werden grundsätzlich niemals mehr als drei Iterationen ausgeführt. Nach diesen drei Iterationen wird weiterverfahren, so, als ob der Korrektor nicht konvergierte. An dieser Stelle wird dann die Iterationsmatrix neuberechnet, falls sie noch von vorherigen Zyklen stammen sollte. Anschliessend wird erneut maximal 3-mal iteriert.
  2. Wurde das Newton-Kantorovich-Verfahren schon wiederholt, so wird die Schrittweite um den Faktor $3/10$ vermindert. Dies natürlich nur, falls die minimal mögliche Schrittweite dabei nicht unterschritten wird. Die Iterationsmatrix $W$ wird neu berechnet. Dabei braucht die Jacobimatrix $J$ nicht mehr erneut ausgewertet werden, da sie in der Matrix rj gespeichert wird. Dies kostet zwar nicht unerheblich Speicherplatz, beschleunigt aber die Integration merkbar.
  3. Sollte nun tatsächlich die Schrittweite die minimal mögliche nicht überschreiten, so wird die Ordnung auf 2 gesetzt, falls sie vorher größer als 2 war und auf 1 gesetzt, falls sie vorher 2 gewesen ist.

In den Fällen (2) und (3) wird der gesamte Zyklus wiederholt. Im Gegensatz zu einstufigen Verfahren, wo lediglich die nicht akzeptierte Stufe erneut durchlaufen wird, wird hier der komplette Zyklus neu berechnet. Schrittzurückweisungen sind jedoch die Ausnahme, sodaß dieser Punkt nicht in dem starken Maße Bedeutung hat, wie man möglicherweise befürchten könnte. Programmiermässig wäre es möglich nur die nicht akzeptierte Stufe zu wiederholen, jedoch wird dies aus theoretischen Überlegungen heraus nicht genutzt. Jedoch sprechen auch Gründe wie Effizienz und Speicherplatzbedarf gegen eine stufenweise Wiederholung. Dies gilt auch für das Programm TENDLER, nicht jedoch für das Programm ODIOUS. Sollte sich der Zyklus nicht wiederholen lassen, d.h. die Ordnung und die Schrittweite sind so klein als möglich, so wird das Unterprogramm verlassen und eine entsprechende Markierung in den Rückgabeparametern gesetzt. Die obigen Überprüfungen werden bei den Marken 780 und 800 durchgeführt.

Die Jacobimatrix $J=f_y$ der Differentialgleichung wird nur an drei Stellen ausgewertet und zwar einmal natürlich ganz zu Anfang, also beim allerersten Aufruf des Unterprogrammes STINT. Dann jedoch nur, falls $h{\beta_{ii}/\alpha_{ii}}$ sich um mehr als 80% geändert hat, oder aber falls wie in (1), nach drei Iterationen der Konvergenztest nicht zufriedenstellend verlaufen war und gleichzeitig dabei die Iterationsmatrix noch aus vergangenen Zyklen stammte.

8. Der Konvergenztest im Programm DIFSUB

Voreingestellt ist die Rechnung in doppelter Genauigkeit, jedoch wird die Invertierung der Iterationsmatrix $W=I-h\gamma J$ in einfacher Genauigkeit vorgenommen. Dies hat den Vorteil, daß weniger Speicherplatz in Anspruch genommen werden muß. Die Multiplikation mit der Inversen erfordert jedoch stets die Umwandlung von einfacher Genauigkeit in doppelte Genauigkeit. Das ist nicht unaufwendig, wenn es wiederholt durchgeführt werden muß. Diese Rückmultiplikation wird natürlich wesentlich häufiger durchgeführt, als die Invertierung. C.W. Gear wählt absichtlich den Weg der Invertierung der Matrix anstatt einer $LU$-Zerlegung, zur Vermeidung des Aufrufes eines Unterprogrammes, hier der Faktorisierung der Matrix $W$. Die Berechnung des Matrix-Vektor Produktes erfolgt inline.

Diese Programmiertaktik ist nicht nachahmenswert. Die Rechenarbeit steckt in Gleitkommaoperationen und nicht in der Verwaltung des Kellerspeichers und der Übergabe von Parametern.

Die Anzahl der Gleitkommaoperationen ist bei einer Invertierung für jede Dimension $n\ge2$ größer als bei einer $LU$-Zerlegung. Jedoch verfolgt Tendler (1973) in seinem Programm DIFJMT den gleichen eingeschlagenen Weg, wie in dem Programm DIFSUB. In dem Programm STINT wird dann allerdings, wie in dem Programm GEAR, die Invertierung ersetzt durch eine $LU$-Zerlegung nach dem Programm von Moler (1972). Die Programme EPISODE, GEAR, LSODE, $\ldots$ wählen alle den Weg der $LU$-Zerlegung.

In dem Programm DIFSUB wird die Iteration abgebrochen, falls jede Komponente des Pseudoresiduums kleiner ausfällt als

$$ \tau = {\varepsilon\over(2p+4)n} \mathtt{YMAX(I)}. $$

In DIFSUB wird also die Iteration abgebrochen, falls gilt

$$ \left\| y_{n+1}^{\nu+1} - y_{n+1}^\nu \right\|_w {\mskip 5mu}\le{\mskip 5mu} {\varepsilon\over 2n(p+2)}. $$

Dabei ist $w$ ein Gewicht und es wird die gewichtete Maximumnorm benutzt, also

$$ \mathopen\|y\mathclose\|_w = \max_{i=1}^n \left|y_i\over w_i\right|. $$

Diese gewichtete Maximumnorm ist insoweit erstaunlich, als für die Fehlerkontrolle und die Schrittweiten- und Ordnungssteuerung eine andere Norm benutzt wird, nämlich

$$ \mathopen\|y\mathclose\|_w = \sqrt{\sum_{i=0}^n\left(y_i\over w_i\right)^2} $$

Der programmiermässige Vorteil liegt nun darin, daß wenn auch nur einmal festgestellt wird, daß keine Konvergenz in allen Komponenten vorliegen kann, der Konvergenztest sofort abgebrochen werden kann. Da der Test natürlich komponentenweise von unteren Indices zu größeren Indices vollzogen wird, genügt schon das Nichterfülltsein des obigen Konvergenzkriteriums für eine einzige Komponente. Diese Möglichkeit der geringen Effizienzsteigerung wird im Programm DIFSUB nicht genutzt. Jedoch ist der Fall der nicht-akzeptierten Korrektoriteration so selten, was natürlich auch sehr erwünscht ist, daß dieser geringe Rechenvorteil selbstverständlich nicht ins Gewicht fällt und daher von neueren Programmen auch nicht mehr so durchgeführt wird. Dies erleichtert auch sehr den Austausch von Normen. Auch auf Vektorrechnern hätte dieses Vorgehen u.U. Nachteile.

Konvergenzraten werden in dem Programm DIFSUB nicht gebildet und daher auch nicht benutzt. Sehr wohl werden Konvergenzraten in den Programmen EPISODE, GEAR, LSODE, $\ldots$ benutzt. Man vgl. hier auch Shampine (1978) für weitere Bemerkungen bzgl. der Evolution von Programmen, maßgeblich an der dortigen Stelle Adams-Verfahren.

Es kommen nun folgende Heuristiken zur Anwendung:

  1. Wird selbst nach drei Iterationen mit Hilfe des modifizierten Newton-Kantorovich-Verfahrens die oben angegebene Bedingung nicht erfüllt, so wird beim ersten Konvergenzversagen die Jacobimatrix $J$ der Differentialgleichung und insbesondere die Iterationsmatrix $W=I-h\gamma J$ neu berechnet.
  2. Bei zwei Konvergenzversagen hintereinander wird die Schrittweite $h$ geviertelt.
  3. Ist die Schrittweite minimal und wurde beim modifizierten Newton-Kantorovich-Verfahren schon die Iterationsmatrix $W$ neu berechnet, so wird das Unterprogramm DIFSUB mit einer entsprechenden Fehlermarkierung in einer Statusvariablen verlassen.

Beim Picardschen Iterationsverfahren, bei der die Jacobimatrix $J$ nicht benutzt wird, entfällt die Möglichkeit der Auffrischung der Iterationsmatrix und es wird daher sofort die Schrittweite $h$ um den Faktor vier vermindert. Da in dem Programm DIFSUB lediglich Speicherplatz für eine einzige Matrix bereitsteht, können auch die beiden Matrizen $W$ und $J$ nicht getrennt behandelt werden, wie dies in dem Programm TENDLER geschieht.

Wird also in den Programmen DIFSUB, EPISODE, GEAR, LSODE, $\ldots$ eine Auffrischung auch nur einer der Matrizen $J$ oder $W$ nötig, so müssen beide Matrizen neu berechnet werden. In diesem Zusammenhang sind die Bemerkungen von Enright (1978) von Interesse, aber auch Shampine äußert sich zu der getrennten oder gemeinsamen Speicherung der Matrizen in Shampine (1980), Shampine (1981a) und Shampine (1982b).

9. Der Konvergenztest in DSTIFF, GEAR und EPISODE

Es sei nocheinmal darauf hingewiesen, daß der Faktor $\tau$ sehr stark direkt auf die Effizienz des gesamten Programmes einwirkt. Zum einen bestimmt ja der Korrektor ausschließlich die Anzahl der Funktionsauswertungen. Muß nun aufgrund eines sehr strengen Konvergenztests mehr als nur einmal iteriert werden, so muß auch entsprechend häufig die Differentialgleichung, also die Funktion ausgewertet werden. Zum anderen müssen für den Konvergenztest selber mindestens eine Norm ausgerechnet werden. Dies ist nicht selten aufwendiger, als die Auswertung der Funktion selber. Für das Schalten sogar müssen insgesamt zwei Normen berechnet werden, wenn der Benutzer den Schaltmodus angewählt hat. Der aufwendigste Teil jedoch ist die Rücksubstitution. Hier müssen quadratisch in der Dimension, also $n^2$ viele arithmetische Operationen durchgeführt werden. Die Aufmerksamkeit dem man deswegen diesem kritischen Parameter schenkt, ist aus diesen Gründen daher sehr verständlich.

1. Der eigentliche, zentrale Test lautet in dem Programm GEAR:

$$ \delta^\nu\cdot\min\left(1,\sqrt2{\mskip 3mu}c_\nu\right) \le{\varepsilon\over2(p+2)}\xi, $$

wobei $\xi$ eine verfahrensabhängige Konstante ist, die also abhängt von der Ordnung $p$ und den Koeffizienten der Formeln, hier also entweder Adams-Formeln oder BDF. $c_\nu$ ist hierbei eine modifizierte Konvergenzrate, welche ganz zu Anfang der Integration auf 1 gesetzt wird.

Programmiert wird dieser Test in der Form

$$ \hat c_\nu \gets \max\left[{9\over10}c_{\nu-1}, \left(\hat\delta^\nu\over\hat\delta^{\nu-1}\right)^2\right],\qquad \hbox{falls}\quad\nu\gt 1, $$

mit

$$ \hat\delta^\nu \gets \mathopen\|y^\nu-y^{\nu-1}\mathclose\|_w^2, \qquad \mathopen\|y\mathclose\|_w = \sqrt{\sum_{i=1}^n\left(y_i\over w_i\right)^2}. $$

Der Konvergenztest gilt als bestanden, wenn

$$ \hat\delta^\nu\cdot\min\left(1,2\hat c_\nu\right) \le n\left(\mathtt{tq[4]}{\mskip 3mu}\varepsilon\right)^2, $$

mit

$$ \mathtt{tq[4]} \gets {\mathtt{tq[2]}\over2(p+2)},\qquad\mathtt{tq[2]}\gets\mathtt{pertst[}p,\mathtt{meth},2\mathtt{]} $$

und der Konstantentryade $\mathtt{pertst[][][]}\in\mathbb{C}^{12\times2\times3}$. Dabei wird $2c^\nu\delta^\nu$ als eine Schätzung der Korrektur für die nächste Iteration angesehen. Auch hier ist natürlich $c^\nu$ am Anfang der Iteration nicht bekannt. Deswegen wird der letzte Wert $c^\nu$ vom vorherigen Schritt genommen.

2. Konvergenztest. Der Konvergenztest im Programm EPISODE ist gegenüber dem Test im Programm GEAR leicht verändert. Dies liegt zum einen daran, daß das Programm EPISODE Formeln verwendet, bei denen die Koeffizienten von den Schrittweitenverhältnissen abhängig sind und desweiteren, weil es Intention bei der Entwicklung des Programmes EPISODE war, möglichst große Schrittweitenänderungen zuzulassen.

Sei $c_\nu$ die modifizierte Konvergenzrate (und nicht das Quadrat der modifizierten Konvergenzrate). Aktualisiert wird diese Größe durch die Rechenvorschrift

$$ c_\nu\gets\max\left({1\over\sqrt{10}}, {\mskip 3mu}{\delta^\nu\over\delta^{\nu-1}}\right),\qquad\nu\ne0 $$

und die gewichtete Norm des Pseudoresiduums $\delta^\nu$ im $\nu$-ten Iterationsschritt werde berechnet durch

$$ \delta^\nu \gets \mathopen\|y^\nu-y^{\nu-1}\mathclose\|_w, \qquad \mathopen\|y\mathclose\|_w = \sqrt{{1\over n}\sum_{i=1}^n\left(y_i\over w_i\right)^2}. $$

Der eigentliche Konvergenztest gilt als erfolgreich bestanden, wenn

$$ \delta^\nu\cdot\min\left(1,c_\nu\right)\le{\varepsilon\over10}\xi. $$

Der Wert $\xi$ ist wiederum eine vom Verfahren und von den Schrittweitenverhältnissen abhängige Größe. Erneut wird $c_\nu$ ganz zu Anfang der Integration auf 1 gesetzt. Auffälligste Veränderung ist hierbei der Faktor $1/\sqrt{10}$, welcher in der Tat dafür sorgt, daß sich die modifizierte Konvergenzrate wesentlich schneller an die wahren Konvergenzverhältnisse anpassen kann. Bei Probeläufen mit dem Programm TENDLER wurde gefunden, daß dies von großem Vorteil ist. Das Programm LSODAR verwendet hier entsprechend den Faktor $1/5=0.2$ und nicht mehr $3/\sqrt{10}\approx0.95$. Diese Beobachtung scheint daher universellere Gültigkeit zu haben.

Das Programm DSTIFF verwendet genau wie das Programm GEAR diesen Test. Das Programm DSTIFF ist dem Programm GEAR überhaupt sehr ähnlich. Vereinfacht ausgedrückt ist das Programm DSTIFF gleich dem Programm GEAR mit veränderten Formeln und gewissen veränderten Abstimmungsparametern.

10. Der Konvergenztest im Programm LSODE

In dem Programm LSODE nun wurde dieser Test noch weiter abgeändert. Der Wert $\delta^\nu$ wird wie oben berechnet, aber es wird jetzt die gewichtete RMS-Norm statt der diskreten ${\cal L}_2$-Norm verwendet. Auf Konvergenz wird jetzt wie folgt getestet:

$$ {3\over 2} c^\nu\delta^\nu p! {\mskip 3mu} k_{m+1} \ell_m {\mskip 5mu} \le {\mskip 5mu} {\varepsilon \over 2(p+2)} {\mskip 5mu}\: . $$

Hierbei ist $k_{m+1}$ der lokale Diskretisierungsfehler des benutzten Verfahrens und $\ell_m$ ist die letzte Komponente des Korrektor-Vektors $\ell$ (Nordsieck-Darstellung). Sollte nun die Konvergenzrate $c^\nu$ in der zweiten oder dritten Iteration größer als 2 sein, so geht das Programm LSODE davon aus, daß die Iteration divergiert.

Die Beschreibung der Konvergenztests in den Programmen LSODA und LSODAR wird später, bei der Beschreibung dieser beiden Programme, gegeben. Sie sind jedoch dem obigen Test sehr ähnlich. Dennoch wurden geringfügige Modifikationen durchgeführt. Insbesondere betrifft dies den Aspekt der Rechnung nahe der Maschinengenauigkeit und natürlich bzgl. des Schaltens. Bei beiden Programmen ist es möglich, daß der Konvergenztest eine Entscheidung trifft, ob und wie geschaltet wird.

Es gibt gute Gründe, warum gerade nur maximal dreimal iteriert wird. Umfangreiche Testläufe, u.a. auch von verschiedensten Autoren, ergaben das folgende Bild:

  1. Falls vier oder mehr Iterationen erlaubt sind, ist die Anzahl der Fälle, in denen keine Konvergenz erzielt wurde nicht wesentlich geringer und
  2. falls nur zwei Iterationen erlaubt sind, ist die Anzahl der Fälle, wo keine Konvergenz erreicht wurde jedoch bedeutend größer.

Zu diesen Ergebnissen beachte man auch den Startfehler-Satz von Liniger. Dies deckt sich mit Erfahrungen wie sie auch C.W. Gear gemacht hat, siehe Gear (1971a), Gear (1971b), Gear (1971c).

Bibliographisch:

  1. Gear, Charles William: “Experience and Problems with the Software for the Automatic Solution of Ordinary Differential Equations”, in “Mathematical Software”, ACM Monograph Series, Editor John Rice, Academic Press, New York London, 1971, pp.211–227
  2. Gear, Charles William: “The Automatic Integration of Ordinary Differential Equations”, CACM, Vol 14, No 3, March 1971, pp.176–179
  3. Gear, Charles William: “Algorithm 407. DIFSUB for Solution of Ordinary Differential Equations [D2]”, CACM, Vol 14, No 3, March 1971, pp.185–190
  4. Gear, Charles William: “Numerical Initial Value Problems in Ordinary Differential Equations”, Prentice-Hall, Englewood Cliffs (New Jersey), 1971, xvii+253 S., Kopie

11. Die Korrektoriteration im Programm DDASSL

1. Die grundlegende Idee beim Programm DDASSL zur Lösung von differential-algebraischen Gleichungen der Form

$$ F(t,y,\dot y)=0,\qquad y(t_0)=y_0,\qquad \dot y(t_0)=\dot y_0, $$

besteht darin, $\dot y$ zu ersetzen durch eine Näherung für die Ableitung unter Zuhilfenahme zurückliegender Werte, also Werten, welche schon berechnet wurden. Für Diskretisierungen solcher Gleichungen vgl. man auch Liniger (1979). Beispielsweise löst man bei der Ordnung 1, in jedem Zeitschritt das nicht-lineare Gleichungssystem

$$ F\left(t_n,y_n,{y_n-y_{n-1}\over h_n}\right) = 0, $$

nach $y_n$ auf. Für die Approximationen der Ableitung $\dot y$ werden in dem Programm DDASSL die BDF$i$, mit $i=1,\ldots,5$, verwendet. Nach einem gelungenen Schritt wird überprüft, ob die Ordnung oder die Schrittweite verändert werden soll, in Abhängigkeit der berechneten Näherungslösung an zurückliegenden Punkten.

Die Darstellung der Lösung erfolgt in Form von modifizierten rückwärtsgenommenen Differenzen, wobei zusätzlich darauf geachtet wird, daß der führende Koeffizient bei der BDF, bei wechselnden Schrittweitenverhältnissen, immer gleich bleibt.

Durch Diskretisierung der Ableitung $\dot y$ mittels der BDF, entsteht ganz allgemein das Problem, das i.d.R. nicht-lineare Gleichungssystem

$$ F(t,y,\hat\alpha y+\beta) = 0 $$

in jedem Zeitschritt zu lösen. $\hat\alpha$ ist eine skalare Konstante, welche sich verändert kann, falls sich die Ordnung oder die Schrittweite ändert. $\beta$ ist ein Vektor, welcher Werte an zurückliegenden Zeitpunkten zusammenfasst. Die obige nicht-lineare Gleichung wird durch ein gedämpftes, modifiziertes Newton-Kantorovich Iterationsverfahren, in jedem Zeitschritt gelöst. Das Iterationsverfahren hat mithin die Form

$$ y^{\nu+1} = y^\nu-\omega W^{-1}F(t,y^\nu,\hat\alpha y^\nu+\beta). $$

Die Iterationsmatrix $W$ ist gegeben durch

$$ W = {\partial\over\partial y}F(t,y,\hat\alpha y+\beta) = \left(F_t,{\mskip 3mu}F_y,{\mskip 3mu}F_{\dot y}\right)\pmatrix{0\cr I\cr \hat\alpha I\cr} = F_y+\hat\alpha F_{\dot y}. $$

$\omega$ ist ein skalarer Dämpfungsparameter und definiert durch

$$ \omega = \cases{ {2/(1+\hat\alpha/\alpha)}, & falls $\hat\alpha\ne\alpha$, \cr 1, & falls $\hat\alpha=\alpha$. \cr } $$

2. Die Modifizierung des gedämpften Newton-Kantorovich Iterationsverfahrens besteht darin, daß die Iterationsmatrix $W$ auch noch in nachfolgenden, und häufiger auch noch in zahlreichen weiteren Zeitschritten, benutzt wird. Das Newton-Kantorovich Iterationsverfahren benutzt ein und dieselbe Iterationsmatrix $W$, über alle Iterationsinidizes $\nu$, im Gegensatz zum Newton-Raphson Iterationsverfahren, welches auch bei jedem neuen Iterationsschritt $\nu$, die Iteratiosmatrix $W=W(y^\nu)$ neu berechnen und faktorisieren würde.

Die Iterationsmatrix $W$ wird im Programm DDASSL genau dann neu berechnet und faktorisiert, falls

  1. die Iteration nicht konvergiert und dabei die Iterationsmatrix “alt” war,
  2. wenn sich die führenden Koeffizienten der benutzten Formeln, welche von der Schrittweite $h$ abhängen, “stark” geändert haben. Dies enstpricht dem üblichen Hindmarsh-Test.

Die Iterationsmatrix $W$ wird dabei in der Form $(1/\hat\alpha)F_y+F_{\dot y}$ gespeichert. Zur Lösung der linearen Gleichungssysteme werden LINPACK-Routinen benutzt, sowohl für vollbesetzte, als auch für schwachbesetzte Systeme mit Bandstrukturen.

Als Vektornorm benutzt das Programm DDASSL die RMS-Norm, welche berechnet wird in der Form

$$ \mathopen\|x\mathclose\| = M\sqrt{{1\over n}\sum_{i=1}^n\left(x_i\over Mw_i\right)^2} = \sqrt{{1\over n}\sum_{i=0}^n\left(x_i\over w_i\right)^2} $$

und $\mathopen|x\mathclose|=0$, falls $M=0$. Das letzte ist natürlich äquivalent mit $x_1=x_2=\ldots=x_n=0$. Der Gewichstvektor $w$ bestimmt sich zu

$$ w_i = \varepsilon_r\mathopen|y_i\mathclose|+\varepsilon_a, $$

wobei $\varepsilon_a$ die geforderte absolute Genauigkeit ist und $\varepsilon_r$ die vom Benutzer angeforderte relative Genauigkeit bezeichnet. Der Gewichstvektor hängt also von der Genauigkeitsforderung des Benutzers ab. Das Programm LSODE benutzt die gleiche Norm wie das Programm DDASSL, während hingegen die beiden Programme LSODA und LSODAR eine gewichtete Maximumnorm verwenden und dies um eine damit passende Vektornorm zu einer entsprechend gewählten Matrixnorm zu erhalten.

Wird mehr als einmal iteriert, so wird eine Konvergenzrate $\rho$ berechnet, mit

$$ \rho = \left(\|y^{\nu+1}-y^\nu\|\over\|y^1-y^0\|\right)^{1/\nu} $$

Es sei ferner

$$ s = \cases{ {\rho/(1-\rho)}, &falls $\rho\le9/10$ , \cr 100, &sonst . \cr } $$

Auf Konvergenz, oder genauer auf Abbruch der Iteration, wird wie folgt getestet. Zuerst wird überprüft, ob

$$ \mathopen\|y^{\nu+1}-y^\nu\mathclose\| \le 100 u \mathopen\|y^0\mathclose\|. $$

Hierbei ist $u$, wie üblich, die kleinste positive Zahl, die auf der Rechenanlage darstellbar ist, sodaß $1+u>1$ gilt. Anschliessend wird getestet, ob $\rho>9/10$ ausfällt. Ist dies der Fall, so wird die Iteration abgebrochen, und es wird also dann angenommen, daß keine Konvergenz statt gefunden hat. Zum Schluß des eigentlichen Tests Schließlich folgt ggf. der zentrale Test, ob

$$ s \mathopen\|y^{\nu+1}-y^\nu\mathclose\| \le {33\over100}. $$

Iteriert wird maximal 4 Male. Nach der vierten Iteration wird ebenfalls Divergenz angenommen.

Bibliographisch:

  1. Liniger, Werner: “Multistep and One-Leg Methods for Implicit Mixed Differential Algebraic Systems”, IEEE Transactions on Circuits and Systems, Vol CAS-26, No 9, September 1979, pp.755–762