, 19 min read

Divergenz der Korrektoriteration: Theorie und Experimente

Original post is here eklausmeier.goip.de/blog/2024/06-17-divergenz-und-korrektoriteration-theorie-und-experimente.


Einer der ganz zentralen Bestandteile eines Programmes, basierend auf Verfahren mit impliziten Stufen, ist die Auflösung der Implizitheit durch ein Iterationsverfahren zur Lösung von Gleichungssystemen. Diese Gleichungssysteme sind i.d.R. nichtlinear. Bei steifen Differentialgleichungen und bei linearen Mehrschrittformeln, als auch bei impliziten Runge-Kutta-Verfahren, ist der Aufwand, der hierfür getrieben wird, enorm und stellt den maßgeblichen Anteil an der Gesamtrechenzeit dar. Man ist jedoch nicht so sehr an der exakten Lösung der zahlreich auftauchenden nichtlinearen Gleichungssysteme interessiert, sondern an der zügigen Integration der Differentialgleichung. Die Gleichungssysteme sind hier nur ein Weg dorthin. Es stellt sich nun heraus, daß das zur Anwendung gelangende Newton-Kantorovich Iterationsverfahren unter gewissen Umständen nicht immer konvergiert. Diese Divergenz bleibt häufig unbemerkt und äußert sich lediglich in größeren globalen Fehlern. Da die Schaltentscheidung im Programm TENDLER, wie auch in anderen Programmen, während der Korrektoriteration gefällt wird, ist diesem Punkte in schaltfähigen Programmen besondere Aufmerksamkeit zu widmen.

Zuerst wird das Newton Iterationsverfahren in Beziehung gesetzt zu bekannten Iterationsverfahren für lineare Gleichungssysteme und umgekehrt. Anschliessend erlauben die Ergebnisse von Hughes Hallett (1984) eine einfache Antwort darauf, ob überhaupt gewisse Modifikationen des Newton Iterationsverfahren konvergieren können. Die Erfahrungen von Byrne/Hindmarsh/Jackson/Brown (1977) zeigen anhand praktischer Beispiele, daß die Korrektoriteration tatsächlich wesentlichen Einfluß auch auf die Genauigkeit des Verfahrens haben und nicht nur auf die Effizienz.

1. Das modifizierte Newton-Verfahren und Spezialisierungen

Hier seien kurz die wichtigsten Iterationsarten zur Lösung linearer Gleichungsysteme angeschrieben. Später, bei einem Erklärungsversuch für auffallende Verhalten der Korrektoriterationen, werden die Untersuchungen wichtig werden. Bei der Erklärung des sehr großen globalen Fehlers bei den beiden Programmen GEAR und EPISODE, bei der Wahl einer Diagonalapproximation der Jacobimatrix, sind die Divergenzsätze von Hughes Hallett (1984) ein möglicher Ansatzpunkt. Eine Charakteristik der Erklärung ist, daß zwischen den beiden möglichen Interpretationen, modifizierten Newton-Verfahren und Iterationsverfahren für lineare Gleichungssysteme, hin und her gewechselt wird.

Das gedämpfte modifizierte Newton-Verfahren lautet

$$ x^{\nu+1} = x^\nu - \omega W^{-1}\varphi(x^\nu,x^{(\nu+1)}), \qquad\nu=1,2,\ldots $$

zur Lösung des Nullstellenproblems $f(x)=0$, $x={\mskip 5mu}?$ Hierbei ist die Matrix $W$ in der Regel eine Näherung für die Jacobimatrix $J=f'(x^\nu)$. Insbesondere hängt die Matrix $W$ häufig auch von den vorhergehenden Iterierten $x^1, x^2, \ldots{\mskip 3mu}$ ab, also $W=W(x^1,\ldots,x^\nu)$. Dies muß jedoch nicht unbedingt so sein. Um die Schreibweise handlich und übersichtlich zu halten, seien im weiteren gewisse Abhängigkeiten in der Schreibweise nicht gesondert erwähnt. Den Wert $f(x^\nu)$ bezeichnet man als das Residuum und die Differenz zweier aufeinander folgender Iterierten $x^{\nu+1}-x^\nu$ wird Pseudoresiduum genannt.

Die beiden Begriffe Residuum und Pseudoresiduum werden auch für die Beträge, bzw. Normen dieser Größen benutzt, also $|f(x^\nu)|$ und $|x^{\nu+1}-x^\nu|$. Wird also von genügend kleinem Residuum gesprochen, so ist damit natürlich die entsprechende Norm gemeint.

Die beiden wichtigsten Spezialfälle dieser recht allgemeinen Iterationsklasse sollen jetzt kurz angegeben werden. Die Auffassung als Spezialisierung muß nicht immer vorgenommen werden. In anderem Zusammenhang kann es durchaus günstiger sein, anders vorzugehen.

Es seien allgemein Zerlegungen der invertierbar vorausgesetzten Matrix $A$ betrachtet: $A=L+D+R$, mit Subdiagonalmatrix $L$, Diagonalmatrix $D$ und Superdiagonalmatrix $R$. Voraussetzung sei weiter, daß keine der Diagonalelemente der Diagonalmatrix $D$ verschwindet, also $\det D\ne0$. Eine andere Zerlegung, allgemeinere Zerlegung der Matrix $A$, sei $A=M-N$, mit invertierbarer Matrix $M$. Gegeben sei das Nullstellenproblem $Ax-b=0$, $x={\mskip 5mu}?$

Zwei typische Vertreter für Iterationsarten zur Lösung obiger Gleichung sind nun die Jacobi-Überrelaxation (JOR-Verfahren) und die sukzessive Überrelaxation (SOR-Verfahren).

Die erstere der beiden Iterationen findet sich auch unter anderen Benennungen.

Die Iterationsvorschrift für das JOR-Verfahren lautet nun

$$ x^{\nu+1} = x^\nu - \omega D^{-1}(Ax^\nu - b), \qquad\nu=1,2,\ldots $$

mit der Iterationsmatrix $J_\omega = D^{-1}(D-\omega A) = I - \omega D^{-1}A$. Hierbei ist $J_\omega$ die Matrix, die bei der Schreibweise auftritt

$$ x^{\nu+1} = J_\omega x^\nu + \omega D^{-1}b. $$

Die weitere wichtige Spezialisierung für den Dämpfungs- bzw.
Verstärkungsfaktor $\omega$, nämlich $\omega=1$, liefert die gewöhnliche [Jacobi Iteration]. Jacobi, Carl Gustav (1804--1851).

Die Iterationsvorschrift für das SOR-Verfahren lautet

$$ x^{\nu+1} = x^\nu - \omega D^{-1}\left(Lx^{\nu+1} + (D+R)x^\nu - b\right), \qquad\nu=1,2,\ldots $$

mit der Iterationsmatrix

$$ R_\omega = (D+\omega L)^{-1}\bigl((1-\omega)D-\omega R\bigr) = (I+\omega D^{-1}L)^{-1}\bigl((1-\omega)I-\omega D^{-1}R\bigr). $$

Die weitere Spezialisierung $\omega=1$, liefert das Gauß-Seidel Iterationsverfahren. Gauß, Carl Friedrich (1777--1855), von Seidel, Philipp Ludwig (1821--1896).

Nebenläufig sei erwähnt, daß die obige Form auch diejenige Gestalt hat, die sich für die Programmierung am besten eignet, wobei man die Vorkonditionierung mit der Inversen der Diagonalmatrix $D$ natürlich am Anfang der Rechnung und nur einmal durchführt. Wichtig ist zu vermerken, daß man bei der obigen Iterationsvorschrift für die einzelnen Komponenten von unteren Indices zu höheren Indices durchläuft, also $x_i^\nu, i=1,\ldots,n$. Die umgekehrte Reihenfolge für den Durchlaufsinn der Indices ist natürlich ebenso möglich; hier vertauschen sich dann lediglich die Rollen der Matrizen $L$ und $R$.

Auch noch auf eine andere Art und Weise läßt sich das Gauß-Seidel Iterationsverfahren als modifiziertes Newton-Kantorovich Verfahren interpretieren. Mit der Zerlegung

$$ (L+D)x^{\nu+1} = b - Rx^\nu $$

erhält man

$$ \eqalign{ x^{\nu+1} &= (L+D)^{-1}(b-Rx^\nu) = x^\nu - (L+D)^{-1}\bigl((L+D)x^\nu + Rx^\nu - b\bigr)\cr &= x^\nu - (L+D)^{-1}(Ax^\nu - b).\cr } $$

Die Iterationsmatrix $W=L+D$ läßt hier jetzt die Interpretationen zu als Näherung für die Jacobimatrix $J$ der Funktion $f$, mit $f(x)=Ax-b$. Hier ersieht man die Asymmetrie die zutage tritt. Der Subdiagonalgehalt der Matrix $J$ wird hier vollständig gewertet, während hingegen der Superdiagonalgehalt, also hier $R$, gar nicht auftritt. Dies ist nicht verwunderlich an sich, soll jedoch vermerkt werden. Bei anderer Reihenfolge des Durchlaufens der Komponenten $x_i^{\nu+1}$ des Vektors $x^{\nu+1}$, hat man natürlich hier $W=D+R$, wie schon oben erwähnt.

Die letzte Bemerkung gilt selbstverständlich ganz allgemein. Mit der oben angegebenen Zerlegung $A=M-N$ im Kopfe, ist

$$ x = x - M^{-1}(Ax-b) = x - M^{-1}\bigl((M-N)x - b\bigr) = M^{-1}Nx + M^{-1}b, $$

also $Mx=Nx+b$. Für die Zerlegungsmatrizen $M$ und $N$ der wichtigsten Iterationsarten erhält man

Verfahren $M$-Matrix $N$-Matrix
JOR $M={1\over\omega}D$ $N={1\over\omega}D-A$
SOR $M=L+{1\over\omega}D$ $N={1\over\omega}\bigl(-\omega R+(1-\omega)D\bigr)$

Sowohl das JOR-Verfahren als auch das SOR-Verfahren lassen sich als konvergenzbeschleunigte Verfahren auffassen. In diesen beiden Fällen wird nur die einfachst denkbare Konvergenzbeschleunigung aus zwei Iterierten gewählt, man erhält also lediglich das gewichtete Zweiermitel der beiden letzten Iterierten, somit

$$ x^{\nu+1}\gets x^\nu + \omega(x^{\nu+1} - x^\nu) = \omega x^{\nu+1} + (1-\omega)x^\nu. $$

Dies ist natürlich nichts anderes als das Verfahren von Euler-Knopp. Euler, Leonhard (1707--1783), Knopp, Konrad Hermann Theodor (1882--1957).

Allgemein arbeitet man mit der unendlichen Block-Dreiecksmatrix

$$ \mathbb{P} = \pmatrix{ B_{11} & & & 0\cr B_{21} & B_{22}\cr B_{31} & B_{32} & B_{33}\cr \vdots & \vdots & \vdots & \ddots\cr }, $$

wobei die Matrizen $B_{ij}\in\mathbb{R}^{n\times n}$ sind und $B_{ij}=\bf 0$, falls $j>i$. Blockimplizite Verfahren werden hier nicht weiter betrachtet. Natürlich summieren sich die Zeilensummen von $\mathbb{P}$ stets auf eins auf, man hat also das Erfülltsein der Konsistenzbedingung

$$ \sum_{j=1}^i B_{ij} = I, \qquad\hbox{für alle } i=1,2,\ldots $$

Man rechnet dann $\mathbb{P}\bf x$, mit $x=(x^0, x^1,\ldots)$.

Diese Art der Konvergenzbeschleunigung hat nun eine Reihe von Bezeichnungen. Im folgenden werde der Bezeichnungsweise von Albrecht/Klein (1984) gefolgt. Autoren Prof. Dr. Peter Albrecht und Prof. Dr. Marcelo Pinheiro Klein. Die zuletzt angeführte allgemeine Art der linearen Konvergenzbeschleunigung heiße $\mathbb{P}$-Extrapolation. Ist $G$ die grundlegende Iterationsmatrix, also $x^{\nu+1}=Gx^\nu+c$, so heiße das beschleunigte Verfahren

$$ x^{\nu+1}\gets Bx^{\nu+1} + (I-B)x^\nu = \bigl(BG+(I-B)\bigr)x^\nu + Bc $$

entsprechend $B$-Extrapolation. Reduziert sich die Matrix $B\in\mathbb{R}^{n\times n}$ zu einem Skalar $\gamma\in\mathbb{R}$, so werde einfach von $\gamma$-Extrapolation gesprochen.

Wenn hier von beschleunigt gesprochen wird, so deutet dies lediglich die Veränderung der Iterationsvorschrift an, nicht jedoch ist damit automatisch gesagt, daß das neue Verfahren auch tatsächlich bei jeder Wahl der Extrapolationsparameter wirklich schneller konvergiert. Dies ist gerade eine der Aufgaben, nicht die einzige, durch geeignete Wahl der Parameter dies zu erzielen. Es ist durchaus zweckmässig und auch sinnvoll auf schnellere Konvergenz zu verzichten, dafür aber ein größeres Konvergenzgebiet zu erzielen und damit für eine wesentlich größere Klasse von Problemen die Garantie zu erhalten, daß das verwendete Iterationsverfahren überhaupt konvergiert. Man vgl. zu diesen Beschreibungen auch den Aufsatz von Varga/Eiermann/Niethammer (1987) und natürlich die entsprechenden Monographien, beispielsweise von Varga (1962) oder Birkhoff/Lynch (1984)](https://www.amazon.com/Numerical-Solution-Elliptic-Problems-Mathematics-6/dp/0898711975). Zu den Autoren: Richard Steven Varga (1928--2022), Michael Eiermann, Wilhelm Niethammer (1933--2023), Garrett Birkhoff (1911--1996), Robert E. Lynch (*1932).

2. Die Divergenzsätze von Hughes Hallett

Es ist nun von Wichtigkeit in Erfahrung zu bringen, unter welchen Umständen man durch $B$-Extrapolation überhaupt Konvergenz erzielen kann. Aussagen hierzu liefert u.a. Hughes Hallett (1984), Andrew Jonathan Hughes Hallett (1947--2019). Eine einfache Strukturaussage macht zuerst

1. Lemma: $A$-Extrapolation der $B$-Extrapolation (direkt hintereinander) liefert eine $AB$-Extrapolation. Insbesondere für kommutierende Matrizen $A$ und $B$ sind $AB$- und $BA$-Extrapolation dasselbe.

Von generellem Interesse ist das folgende Lemma, welches aber auch bei einem späteren Satz die Schlüsselrolle spielen wird.

2. Lemma: (1) Die $\gamma$-Extrapolation konvergiert für gewisse $\gamma>0$ genau dann, wenn sämtliche Realteile der Eigenwerte der Iterationsmatrix $G$ kleiner als eins sind, es muß also gelten $\mathop{\rm Re}\nolimits \lambda_i<1$, für alle $i=1,\ldots,n$ und den $\lambda_i$ aus dem Spektrum $\sigma(G)$. {\parskip=0pt plus 2pt\parindent=20pt\par} (2) Analog konvergiert die $\gamma$-Extrapolation für $\gamma<0$ genau dann, wenn $\mathop{\rm Re}\nolimits \lambda_i>1$, für alle $i=1,\ldots,n$.

Eines der wesentlichen Ergebnisse von Hughes Hallett ist nun der folgende Satz.

3. Proposition: Für eine beliebige Iterationsmatrix $G$ läßt sich nicht immer eine reelle Diagonalmatrix $R$ finden, so, daß die $R$-Extrapolation konvergent ist.

Der Beweis benutzt das oben angegebene Lemma. Die Beweisführung konzentriert sich nun darauf zu zeigen, daß die Eigenwerte von $R(I-G)$ nicht stets in der linken Halbebene $\mathbb{C}^{{}-{}}$ liegen können. Dies wiederum legt es nahe, bekannte Sätze aus der Stabilitätstheorie zu benutzen, wie z.B. die Sätze von Liénard-Chipart und Hurwitz. Hurwitz, Adolf (1859--1919), Liénard, Alfred-Marie (1869--1958), Chipart, A.H..

Beweis: (zum Diagonalextrapolationssatz von Hughes Hallett) Sei $G$ Grunditerationsmatrix. Die $R$-Extrapolation hat als Iterationsmatrix $RG+I-G = I+R(G-I)$, also hat die $\gamma$-Extrapolation der $R$-Extrapolation die Iterationsmatrix $I+\gamma R(G-I)$. Nach dem Lemma muß also gelten: Die Realteile der Eigenwerte von $R(G-I)$ müssen in der linken komplexen Halbebene $\mathop{\rm Re}\nolimits z<0$ liegen ($R\gets\gamma R$). Das charakteristische Polynom

$$ \left|R(G-I)-\lambda I\right| = a_0\lambda^n + a_1\lambda^{n-1} + \cdots + a_{n-1}\lambda + a_n $$

hat die Koeffizienten

$$ \eqalign{ a_1 &= - \sum_{1\le i\le n} \left(R(G-I)\right)_i^i = - \sum_{1\le i\le n} r_i (G-I)_i^i, \cr a_2 &= + \sum_{i_1\lt i_2} \left(R(G-I)\right)_{i_1i_2}^{i_1i_2} = + \sum_{i_1\lt i_2} r_{i_1} r_{i_2} (G-I)_{i_1i_2}^{i_1i_2}, \cr a_3 &= - \sum_{i_1\lt i_2\lt i_3} r_{i_1} r_{i_2} r_{i_3} (G-I)_{i_1i_2i_3}^{i_1i_2i_3}, \cr \vdots & \qquad\qquad\vdots\cr a_n &= (-1)^n r_1 r_2 \ldots r_n \left|G-I\right|. } $$

Nach dem Stabilitätskriterium von Liénard/Chipart ist aber abwechselnde Positivität geeigneter $a_i$ vonnöten. Dies lässt sich durchbrechen, falls $(G-I)_i^i=0$ $\forall i$, oder $(G-I)_{i_1i_2}^{i_1i_2}=0$ $\forall i_1<i_2$, oder $(G-I)_{i_1i_2i_3}^{i_1i_2i_3}=0$ $\forall i_1<i_2<i_3$, oder $\ldots$ $\left|G-I\right|=0$ ist.     ☐

Es ist nun leicht möglich auf den nichtlinearen Fall zu verallgemeinern. Benutzt man das Verfahren $x^{\nu+1}=g(x^\nu)$ so ist ganz entsprechend wie im linearen Falle, die $R$-Extrapolation

$$ x^{\nu+1} = Rg(x^\nu) + (I-R)x^\nu. $$

Hughes Hallett (1984) zeigt nun ganz analog:

4. Satz: (1) Für die oben angegebene $R$-Extrapolation ist es nicht immer möglich eine reelle Diagonalmatrix $R$ zu finden, der Gestalt, daß das resultierende Verfahren dann konvergent ist.

(2) Dennoch ist es grundsätzlich möglich eine Nichtdiagonalmatrix zu finden, so, daß das Verfahren für einen Startwert, der nahe genug an der gewünschten Lösung liegt, konvergiert.

(3) Weiterhin ist es möglich eine Folge von Matrizen zu bestimmen, derart, daß das Verfahren für beliebige Funktionen $g({}\cdot{},{}\cdot{})$ konvergiert.

Der Beweis wird durch eine Linearisierung auf den vorhergehenden Satz zurückgespielt.

Es sind nun diese Ergebnisse mit den Vorbereitungen von oben, die ein zusätzliches Licht auf das Konvergenzversagen des modifizierten Newton-Verfahrens werfen. Dies komplementiert auch die Beobachtungen von Shampine (1980), welcher nicht generelle Konvergenzunmöglichkeiten in Erwägung zog, sondern seine Überlegungen auf den Konvergenztest fokussierte. Lawrence F. Shampine.

Shampine (1980) zitiert als experimentelle Stütze den Aufsatz von Byrne/Hindmarsh/Jackson/Brown (1977), in dem deutlich wird, daß das modifizierte Newton-Verfahren in den beiden Programmen GEAR und EPISODE, bei Verwendung der Richtungsableitung der Funktion $f$ als Näherung für die Diagonalelemente der Jacobimatrix $J$, nicht mehr zufriedenstellend arbeitet. Autoren: Byrne, George D., Hindmarsh, Alan C. (*1942), Jackson, Kenneth R., Brown, H. Gordon. Shampine (1980)](https://epubs.siam.org/doi/epdf/10.1137/0901005 schreibt:

Regardless of the Jacobian approximation, if the convergence test is reliable, the codes should deliver a good solution to the problem. Of course the efficiency is affected, but the accuracy of the results should not be. … it appears that the convergence test is unreliable, and that the potential unreliability can sometimes be exhibited as the result of a very poor approximate Jacobian.

Hierbei ist zu beachten, daß die Programme meistens maximal dreimal iterieren und damit lediglich drei Tests zur Erkennung von Divergenz durchgeführt werden. In dem Programm SYMPOL, zur Lösung von polynomialen nichtlinearen Gleichungssytemen, beispielsweise, wird bis zu 25-mal iteriert. In dem Programm BRENTM zur Lösung allgemeiner nichtlinearer Gleichungssysteme wird maximal 50-mal iteriert. Das Programm COLNEW iteriert pro Gitter nicht mehr als 40-mal.

3. Die Experimente von Byrne/Hindmarsh/Jackson/Brown

Das an anderer Stelle angegebene zweidimensionale Differentialgleichungsproblem P3 von Byrne/Hindmarsh/Jackson/Brown (1977), welches seinen Ursprung in der chemischen Kinetik hat, wurde nun von Byrne, Hindmarsh, Jackson und Brown mit den verschiedensten Parametereinstellungen der beiden Programme GEAR und EPISODE ausgetestet. Die dabei gewonnenen Erfahrungen und Ergebnisse sind für die weitere Analyse sehr wertvoll. Daher soll dieses Beispiel kurz näher untersucht und interpretiert werden. Die folgenden Daten wurden gemessen bei dem Problem P3.

$\varepsilon$ code T $\Vert y-Y\Vert$ nst ${\rm nfe\over nst}$ nfe nje ${\rm nst\over nje}$ $J$-CPU $%T$ $f$-CPU %T
$10^{-3}$ E21 2.33 1.89 4428 1.8 7979 859 5.2 0.0904 4 0.389 17
" E22 2.37 3.70 4418 1.8 7884 845 5.2 0.159 7 0.384 16
" E23 1.92 1120 4337 1.7 7350 893 4.9 0 0 0.358 19
" G13 2.74 2520 11681 1.1 13087 1362 8.6 0 0 0.638 23
" G21 2.25 1.67 5619 1.6 9145 751 7.5 0.0790 4 0.446 20
" G22 2.30 1.44 5573 1.7 9390 754 7.4 0.142 6 0.458 20
" G23 1.45 11400 4532 1.5 6578 662 6.8 0 0 0.321 22
$10^{-6}$ E21 5.18 9.22 10507 1.6 16594 1263 8.3 0.133 3 0.809 16
" E22 5.33 22.2 10610 1.6 16903 1337 7.9 0.251 5 0.824 15
" E23 5.67 2190 13075 1.6 21099 2715 4.8 0 0 1.03 18
" G21 5.75 5.32 14992 1.5 22919 1579 9.5 0.166 3 1.12 19
" G22 5.95 4.68 14984 1.6 23229 1571 9.5 0.295 5 1.13 19
" G23 4.77 3360 15187 1.5 22304 2129 7.1 0 0 1.09 23
$10^{-9}$ E21 16.9 74.2 31282 1.3 41622 2465 12.7 0.259 2 2.03 12
" E22 17.1 50.3 31413 1.3 41794 2532 12.4 0.476 3 2.04 12
" E23 18.7 4310 39078 1.4 55623 6304 6.2 0 0 2.71 14
" G21 16.7 24.4 41058 1.4 56423 3581 11.5 0.377 2 2.75 16
" G22 17.4 26.8 40963 1.4 56513 3613 11.3 0.679 4 2.76 16
" G23 14.4 3840 42199 1.4 56990 4358 9.7 0 0 2.78 19

Dabei bezeichnet $T$ die Gesamtrechenzeit auf einer Rechenanlage vom Typ ^{IBM 370/195} in doppelter Genauigkeit. Der Eintrag code bezeichnet welches Programm, mit welcher Parametereinstellung benutzt wurde. Hierbei deuten die beiden letzten Ziffern den Wert von mf an. Beispielsweise wurde bei E21 das Programm EPISODE mit der Parametereinstellung mf=21 benutzt. Der durch beiderseitige Einrahmung hervorgehobene Eintrag $\|y-Y\|$ kennzeichnet den globalen Fehler, der begangen wurde. nstep gibt die Anzahl der Schritte an zur Vollendung der Integration. nfe gibt die Anzahl der Funktionsauswertungen an, und der Bruch $\hbox{nfe}\over\hbox{nst}$ gibt an, wieviele Funktionsauswertungen pro Schritt getätigt wurden. Schließlich wird durch nje die Anzahl der Jacobimatrixauswertungen angezeigt. Der Bruch $\hbox{nst}\over\hbox{nje}$ zeigt an, über wieviele Schritte im Mittel die Jacobimatrix konstant gelassen wurde. $J$-CPU bedeutet, wieviel Rechenzeit verbraucht wurde zur Ausführung des Unterprogrammes pset und aller Routinen, die dieses Unterprogramm selbst aufruft. Die sich hieran direkt anschließende Spalte $%T$ gibt den prozentualen Anteil an der Gesamtrechenzeit wieder. "$f$-CPU" offenbart die Gesamtrechenzeit, die zur Funktionsauswertung benötigt wurde, allerdings ohne die Zeit zur numerischen Approximation der Jacobimatrix. Die sich ebenfalls direkt rechts daneben anschließende letzte Spalte gibt den prozentualen Anteil am Gesamtrechenzeitbedarf an.

Der Eingabeparameter mf besteht stets aus zwei Ziffern. Die erste Ziffer ist immer aus $\{1,2\}$, wobei 1 das Adams-Verfahren anwählt und 2 die BDF anzeigt. Die zweite Ziffer ist immer aus $\{0,1,2,3,4,5\}$. Für die obige Tabelle sind nur wichtig die Ziffern 1, 2 und 3. Zur Vollständigkeit seien jedoch alle Ziffern kurz erläutert. Es sei erwähnt, daß die beiden letzten Ziffern 4 und 5 noch nicht in den beiden Programmen GEAR und EPISODE zur Verfügung standen, sondern erst in der späteren Version LSODE und den entsprechenden Modifikationen dieses Programmes. Dennoch, da an anderer Stelle, nämlich bei der Beschreibung von LSODA und LSODAR, kurz nocheinmal diese unterschiedlichen Paramtereinstellungen zur Sprache kommen, seien hier sämtliche Einstellungen aufgelistet.

0: Fixpunktiteration. Es wird keine Jacobimatrix ausgewertet und vom Benutzer muß auch keine Routine hierfür bereitgestellt werden.

1: Modifiziertes Newton Verfahren mit einer vollen Jacobimatrix. Diese Matrix muß vom Benutzer durch ein separat programmiertes Unterprogramm dem Programm zur Verfügung gestellt werden. Ob diese Matrix tatsächlich die passende Jacobimatrix $J$ zur Funktion $f$ ist, wird nicht überprüft.

2: Wie 1, jedoch wird die Jacobimatrix intern durch numerische Differentiation berechnet. Diese Einstellung ist für den Benutzer wesentlich bequemer und einfacher. Bei großen Differentialgleichungen ist diese Einstellung jedoch i.a. nicht so effizient, wie diejenige bei 1.

3: Diagonalapproximation der Jacobimatrix. Diese Paramter-Einstellung ist sehr speicherplatzökonomisch. Bei der Diagonalapproximation handelt es sich nicht wirklich um eine Diagonal_approximation_, sondern es wird lediglich eine gewichtete Richtungsableitung der Funktion $f$ zur Newton Iteration herangezogen.

4: Modifiziertes Newton-Verfahren mit einer vom Benutzer bereitgestellten Bandmatrix. Diese Matrix muß wie bei der Parameter-Einstellung 1 durch ein getrenntes Unterprogramm bereitgestellt werden, und wie bei 1 wird diese Matrix nicht überprüft.

5: Wie 4, nur wird hier die Bandmatrix durch numerische Differentiation ermittelt. Dies ist erneut für den Benutzer die bequemste und einfachste Art die Jacobimatrix zu spezifizieren.

Die Interpretation der oben angegebenen Tabelle ist nun wie folgt. Durchschnittlich werden von GEAR über alle Parametereinstellungen betrachtet, durchschnittlich 1.5 Funktionsauswertungen pro Schritt durchgeführt (Standardabweichung unter 0.2), und für EPISODE ergibt sich hier der Wert 1.6 (Standardabweichung ebenfalls unter 0.2). Dies heißt, daß beide Programme zu einem großem Teil im $PEC$-Modus arbeiten. Ferner bedeutet dies, daß eine Konvergenzrate sehr häufig nicht gebildet werden kann und somit die Konvergenzrate aus einem vergangenen Schritt genommen wird, bei dem Konvergenztest in beiden Programmen.

Der Rechenzeitanteil der Funktionsauswertungen an der Gesamtrechenzeit liegt bei dem Programm GEAR bei durchschnittlich 20% (Standardabweichung unter 3%) und für EPISODE bei durchschnittlich 15% (Standardabweichung ebenso unter 3%). Diese Zahlen sind für ein Mehrschrittverfahren mit variabler Schrittweite und variabler Ordnung nicht ungewöhnlich. Für das Programm TENDLER ergeben sich Zahlen ähnlicher Größenordnung.

Für die Strategien bei der Neuauswertung der Jacobimatrix $J$ ergibt sich das folgende Bild. Das Programm GEAR wartet im Mittel ungefähr 9 Schritte, bevor die Jacobimatrix neu ausgewertet wird (Standardabweichung unter 2 Schritten), während hingegen das Programm EPISODE im Mittel rund 8 Schritte wartet (Standardabweichung: ca. 3 Schritte). Erwartungsgemäß wird bei schärferer Genauigkeitsanforderung $\varepsilon$ länger gezögert, bis eine Neuauswertung vorgenommen wird, da ja die gewählten Schrittweiten naturgemäß kleiner sind und damit die Iterationsmatrix $W=I-h\gamma J$ sich nicht so stark ändert. Die Neuberechnung der Jacobimatrix $J$ wird ja entscheidend im Hindmarsh-Test gesteuert und beeinfußt, weniger durch Konvergenzversagen.

Auffällig an den angegebenen Daten in der Tabelle sind nun die eingerahmten globalen Fehler $\|y-Y\|$. Bei der Wahl des Eingabesteuerungsparameters mf=23 für beide Programme, also GEAR und sowohl auch EPISODE, ergeben sich deutlich auffallende, ungewöhnlich große Fehler. Für die Wahl mf=13 für das Programm GEAR gelten die Bemerkungen analog. Hierzu bemerken Byrne, Hindmarsh, Jackson und Brown:

mf=23 does not control the error well for this problem for either code.

Im Lichte der oben vorbereiteten Bemerkungen ist dies nicht mehr so überraschend. Es sind diese Gründe, die bewogen ein Iterationsverfahren mit Diagonalapproximation der Jacobimatrix nicht in das Programm TENDLER mit aufzunehmen. Insbesondere wird beim Schalten nicht zwischen modifizierter Newton-Kantorovich Iteration und Newton-Jacobi Iteration hin und her gewechselt, sondern lediglich zwischen modifiziertem Newton-Kantorovich Iteration und Picard Iteration.

Denkbar wäre, die Diagonalelemente von $J$ weiterzubenutzen, bei einem Wechsel der Iterationsart, also Verwendung des modifizierten Newton-Jacobi Verfahrens, siehe Ortega (1972). Lawrence Shampine (1982) schlägt dies ebenso vor. Er bemerkt hierzu:

We propose using simple iteration until the first Jacobian is formed and thereafter using Jacobi-iteration.

Autoren: Lawrence F. Shampine, James M. Ortega (*1932).

Dennoch, die $n$ Divisionen, $n$ Multiplikationen und $n$ Additionen und die zusätzlichen Speicherstellen sind zu bedenken. Das Programm LSODA bietet dem Benutzer nicht mehr an, entgegen den Offerierungen beim Programm LSODE, die Jacobimatrix durch eine Diagonalapproximation zu ersetzen. Bei einfach auszuwertenden Funktionen ist es vorteilhafter Picard Iteration zu verwenden, als ältere Informationen aus den Diagonalelementen der Jacobimatrix, ja sogar letztlich nur ältere Informationen aus mehr oder minder guten Näherungen für diese Diagonalelemente. In dem Programm von {Suleiman, M.B.}{Hall, George}Suleiman/Hall (1985) wird insoweit besonderer Gebrauch von den Diagonalelementen der Jacobimatrix $J$ gemacht, als daß die Spur der Jacobimatrix mit zum Schalten herangezogen wird. Das Programm TENDLER benutzt diese Informationen nicht. Bei Diagonaldominanz ergibt sich offensichtlich sofort ein Vorteil, bei der Benutzung des modifizierten Newton-Jacobi Verfahrens.

Es ist der zusätzliche Speicherbedarf und die Unsicherheit des Vorteils, die bewogen keinen Gebrauch der modifizierten Newton-Jacobi Iteration zu machen. Vom Speicherplatz sind zyklische lineare Verfahren einstufigen linearen Verfahren in linearen Termen der Dimension der Differentialgleichung _unterlegen:, da ja für jede Stufe eine geeignete Differenzentabelle bereit gehalten werden muß. Beispielsweise beträgt der lineare Anteil des Speichers für die drei Programme LSODA, STINT und TENDLER:

TENDLER STINT LSODA
$37n$ bzw. $42n$ $41n$ $11n$

Angeführt ist hier jeweils der Speicherplatzbedarf, falls die Höchstordnung auf 5 begrenzt wird. Dies liegt daran, daß die Programme LSODE, LSODA und LSODAR im steifen Modus maximal die BDF5 verwenden können. Die BDF6 wird aufgrund ihres kleinen Widlund-Winkels nicht in diesen Programmen benutzt. Die BDF$i$ ($i>6$) sind nicht $D$-stabil. Die zyklischen Verfahren von J.M. Tendler sind hingegen bis zur Ordnung 7 einsetzbar. Das Programm TENDLER lässt sich leicht dahingehend modifizieren, daß nur noch $37n$ Speicherzellen benutzt werden. Allerdings sind dann unter gewissen Umständen, während der Schrittweiten- und Ordnungssteuerung, Doppelberechnungen durchzuführen. Diese Doppelberechnungen müssen im Programm STINT auf jeden Fall durchgeführt werden.