, 7 min read

Das Schalten im Programm NTI/SIMPLE

Denn der Mensch denkt nicht gern logisch, dazu ist er nicht geschaffen. Von Natur aus denkt er assoziativ. Logisch zu denken muß er erst lernen, und es fällt ihm schwer. Mein Vater, der einmal in der Mathematischen Gesellschaft der Berliner Universität einen Vortrag hielt über ein numerisches Verfahren zur Berechnung höherer Differentialgleichungen, hatte sich zu dem Einwand zu äußern, das Verfahren sei doch maßlos umständlich. Er erwiderte, das sei wohl zugegeben, aber das komme daher, daß der liebe Gott dem Menschen den Verstand nicht zum Berechnen von Differentialgleichungen gegeben habe, sondern zur Futtersuche.

W.T. Runge (1966)

Konzeptionell wird bei dem Programm NTI/SIMPLE der gleiche Weg gewählt, der auch bei dem Programm TENDLER gewählt wurde: Man schaltet nicht zwischen verschiedenen Integrationsverfahren hin und her, sondern nur zwischen zwei verschiedenen Iterationsarten. Hier wird also zwischen modifiziertem Newton-Kantorovich-Verfahren und Picard Iteration hin und her gewechselt. Die nötige Information zum Schalten wird während der Iteration gewonnen, falls mit dem Newton-Kantorovich Verfahren gearbeitet wird und es wird keinerlei Information herangezogen beim Wechsel von nicht-steif nach steif.

Im letzteren Falle wird einfach bei Divergenz der Fixpunktiteration zur Newton-Kantorovich Iteration übergewechselt. Das Programm TENDLER benutzt während der Picard Iteration sehr wohl weitere Informationen, die während dieser Iterationsart anfallen.

1. Die Schaltentscheidung im Programm NTI/SIMPLE

Das Programm NTI/SIMPLE von Nørsett/Thomsen (1986) verwendet, wie das Programm TENDLER, beim Wechsel vom steifen Modus in den nicht-steifen Modus, den Wert QRP einer Newton-Kantorovich-artigen Iteration zur Auflösung impliziter Differenzengleichungen. Der Wert QRP wird im Programm NTI/SIMPLE allerdings “stiffness-ratio” genannt wird, im Gegensatz zur üblichen Definition von

$$ \mu = {\max_i\mathopen|\lambda_i\mathclose| \over \min_k\mathopen|\lambda_k\mathclose|}, \qquad \min \mathopen|\lambda_k\mathclose| \ne 0. $$

Mit diesem Wert des Steifheitsmaßes steht das QRP im Falle steifer Differentialgleichung zwar sinngemäß im Zusammenhang, jedoch besteht keine einfache, direkte funktionale Abhängigkeit.

Die Strategie für die Schaltung in dem Programm NTI/SIMPLE ist:

  1. Wechsele von der Newton-Kantorovich zur Picard Iteration über, falls das QRP kleiner als 2 ist.
  2. Wechsele von der Fixpunktiteration zum Newton-Kantorovich Iterationsverfahren über, falls Divergenz bei einfacher Iteration gemeldet wurde.
  3. Warte mindestens 10 Schritte bevor überhaupt an ein erneutes Schalten gedacht werden kann.

Die Rechtfertigung für den ersten Punkt ist der folgende Sachverhalt, der in gewisser Hinsicht eine notwendige Voraussetzung für Konvergenz der Fixpunktiteration darstellt: Beachtet man, daß $\hgJ < 1$ eine hinreichende Bedingung für Konvergenz der Picard Iteration ist, so erhält man mit der Forderung, daß das QRP kleiner ist als 2, eine notwendige Bedingung dafür, daß die hinreichende Konvergenzbedingung erfüllt ist. Das Programm NTI/SIMPLE basiert auf SDIRK Formeln. Der endgültige im Programm NTI/SIMPLE benutzte Wert QRP$^*$ wird als Maximum aller QRP genommen, gebildet über alle Stufen des Runge-Kutta Verfahrens.

In der Iterationsmatrix $W$ ist $h\gamma J$ enthalten. Wenn nun $\hgJ$ groß ist, so muß $\mathopen|W\mathclose|$ “nahe” bei $\hgJ$ liegen, und wenn $\hgJ$ nahe oder unterhalb von 1 ist, so muß $\mathopen|W\mathclose|$ in etwa in der Nähe von 2 sein. Dies ist eine grobe und einfache Motivation für die Benutzung des QRP einer Newton-Kantorovich-artigen Iteration im Programm NTI/SIMPLE. Allerdings benutzt das Programm NTI/SIMPLE, wie das Programm TENDLER, die “verbotene” Rückrichtung. Während im Programm TENDLER hierfür versucht wird teilweise eine Kompensation zu erhalten durch weitere Informationen aus dem Schrittweiten- und Ordnungssteuerungssegment, so wird dies in dem Programm NTI/SIMPLE nicht getan. Allerdings sollte hinzugefügt werden, daß z.B. Minimierungsverfahren für Funktionen $F:\mathbb{R}^n\to\mathbb{R}$ häufig ebenfalls nur Punkte ansteuern mit $F'(x)=0$ und somit wird auch dort letztlich das Bestehen einer “verbotenen” Rückimplikation angenommen.

Divergenz während der Picard Iteration wird auch schon dann angenommen, wenn die Konvergenzrate über 0.8 liegt. Allerdings darf die Konvergenzrate ein einziges Mal doch über 0.8 liegen, da dann von Nørsett/Thomsen (1986) angenommen wird, daß dies durch Rundungsfehler entstanden sein könnte. Dies deutet darauf hin, daß in dem Programm NTI/SIMPLE dem Rechnen nahe an der Grenze der Maschinengenauigkeit nicht die Aufmerksamkeit zuteil wurde, wie dies beispielsweise bei dem Programm LSODA der Fall ist. Die Konvergenzrate während der Picard Iteration ist eine untere Schranke für die Lipschitkonstante $L$ der Funktion $f$:

$$ {1\over \mathopen|h\gamma\mathclose|} {\mathopen\|y^{\nu+1}-y^\nu\mathclose\| \over \mathopen\|y^\nu-y^{\nu-1}\mathclose\|} \le L, $$

wie durch Subtraktion von

$$ y^{\nu+1} = h\gamma f(y^\nu)+\psi, \qquad y^\nu = h\gamma f(y^{\nu-1})+\psi $$

sofort folgt.

2. Bewertung des Programmes NTI/SIMPLE

1. Das Programm NTI/SIMPLE wurde an den Testgleichungen C1, C2, C3, C4, C5, D2, D5 und NA4, NA5, NB5, NE4, ND1, ND2, ND5 aus DETEST untersucht und mit der nicht schaltfähigen Version des Programmes verglichen. Gründe für das Auslassen von Testgleichungen aus DETEST und die benutzte Rechenanlage werden in Nørsett/Thomsen (1986) nicht genannt.

Die Ergebnisse sind nun im einzelnen wie folgend. Hierbei werden die ermittelten Werte für die Genauigkeit $\varepsilon=10^{-2}$ nicht mit angegeben, da bei dieser niedrigen Genauigkeitsanforderung der globale Fehler nicht gut unter Kontrolle blieb. Man findet diese Werte bei Nørsett/Thomsen (1986). Diese Erfahrungen mit den geringen Genauigkeitsansprüchen wurden auch mit dem Programm TENDLER gemacht.

Zur Vermeidung solcher Effekte, wurde in dem Programm TENDLER dann soweit gegangen, daß diese geringe Genauigkeitsanforderung von einem Treiber als Fehleingabe abgefangen wird. Die Tabelle gilt für die Genauigkeitsanforderung $\varepsilon=10^{-4}$.

Dgl. #fix nfe nje #$LU$ $Ux=c$ CPU code
C1 41 460 3 16 195 0.85 neu
  0 421 1 25 416 1.03 alt
C2 35 427 1 18 210 0.85 neu
  0 407 1 26 402 1.01 alt
C3 38 412 1 18 193 0.83 neu
  0 389 1 26 384 0.97 alt
C4 229 3506 11 58 1916 6.66 neu
  0 2992 4 63 2972 7.20 alt
D2 14 608 2 14 508 1.02 neu
  0 583 3 21 571 1.05 alt
D5 0 114 4 20 99 0.16 neu
  0 111 4 20 99 0.15 alt

Deutlich wird der, wenn auch geringe, Rechenzeitgewinn ggü. dem nicht schaltfähigen Programm.

Bei beiden Programmen war der beobachtete globale Fehler am Endpunkt größenordnungsmässig gleich, sodaß die Werte gut vergleichbar sind.

Stellenweise auffällig ist die höhere Anzahl an Jacobimatrixauswertungen des schaltfähigen Programmes. Ebenso sind die benötigten Funktionsauswertungen i.d.R. höher, jedoch sind dafür die Rücksubstitutionen (in der Spalte $Ux=c$) drastisch zurückgegangen.

Die Anzahl der $LU$-Zerlegungen ist ebenfalls vermindert worden. Wenn man beachtet, daß gerade die Kosten für die lineare Algebra in Programmen zur Lösung steifer Differentialgleichungen einen wichtigen Beitrag leisten, so ist dieser Gesichtspunkt hier von besonderem Interesse.

In den Rechenzeiten schlägt sich dieser Gewinn nicht in der erwarteten Weise nieder, da es sich bei den ausgewählten Testdifferentialgleichungen um im wesentlichen kleindimensionale Probleme handelt: alle Dimensionen liegen bei, oder unterhalb von 4.

Bei den Gleichungen NA4, NA5, NB5, NE4, ND1, ND2, ND5 wird kein Vergleich in Nørsett/Thomsen (1986) angegeben, wodurch man die Kosten der Schaltung und der Schaltentscheidung erkennen könnte. Deutlich wird allerdings, daß bei diesen nicht-steifen Differentialgleichungen tatsächlich keine einzige Jacobimatrixauswertung und keine $LU$-Zerlegung benötigt wird, mit Ausnahme bei der Gleichung NA4. Diese Gleichungen werden also korrekt als nicht-steife Gleichungen erkannt. Diese Ergebnisse gelten nicht für die Genauigkeitsanforderung von $10^{-2}$, für die schon mehrfach darauf hingewiesen wurde, daß man einer solchen Genauigkeitsanforderung des Benutzers zumindestens vorsichtig gegenüber stehen sollte. Bei dieser niedrigen Genauigkeit werden häufig $LU$-Zerlegungen und Jacobimatrixauswertungen benötigt, weiterhin ist der globale Fehler inakzeptabel.

2. Zusätzlich zu den obigen Testgleichungen wurden Vergleichsmessungen an der van der Polschen, zweidimensionalen Differentialgleichung durchgeführt. Die hier benutzte van der Polsche Gleichung lautet

$$ \eqalign{ \dot y_1 &= y_2,\cr \dot y_2 &= 100(1-y_1^2)y_2 - y_1,\cr } \qquad\hbox{mit}\quad y(0)=\pmatrix{2\cr 0\cr}\qquad\hbox{für}\quad t\in[0,100]. $$

Anhand dieser Gleichung ergab sich nun für eine Genauigkeitsanforderung von $\varepsilon=10^{-4}$

code #steps #fix $\Vert y_N-y(100)\Vert$ nfe nje #$LU$ $Ux=c$ CPU
neu 187 150 $2.5e-3$ 1172 4 36 261 1.15
alt 185 0 $4.9e-4$ 1172 4 71 1160 1.49

Man erkennt den leichten Rechenzeitgewinn bei zumindestens vergleichbaren globalen Fehlern am Endzeitpunkt. Die Erklärung hierfür ist, daß bei den Spitzen der Lösungen der van der Polschen Differentialgleichung, der wesentlich effizientere Modus der Fixpunktiteration gewählt wurde.

Es ist Zufall, daß die Anzahl der Funktionsauswertungen gleich geblieben ist. Die hohe Anzahl an Fixpunktiterationen deutet darauf hin, daß erfolgreich hin und her geschaltet wurde.

Dies ist ein Hinweis dafür, daß mit Hilfe des QRP eine Möglichkeit gegeben ist im Falle von Runge-Kutta Formeln automatisch zu schalten.

3. Bibliographie

  1. Nørsett, Syvert Paul (1944–2025)
  2. Nørsett, Syvert Paul und Thomsen, Per Grove: “Switching between Modified Newton and Fix-Point Iteration for Implicit ODE-Solvers”, BIT, Vol 26, 1986, pp.339–348
  3. Thomsen, Per Grove