14th April 2026, 29 min read

Das Schalten im Programm LSODA

• 1. Kurzbeschreibung des Programmes LSODA
• 2. LSODAR und Programmtransparenz
• 3. Bewertung und Vergleich des Programmes LSODA
• 4. LSODA im Vergleich mit anderen Lösern aus ODEPACK
• 5. Literatur

Das Programm LSODE und seine Modifikationen zählt zu den leistungsfähigsten und erfolgreichsten Programmpaketen überhaupt, zur Lösung steifer und nicht-steifer Anfangswertaufgaben bei gewöhnlichen Differentialgleichungen, z.B. siehe Seifert (1987). Allerdings muß vor der eigentlichen Integration der Charakter der Differentialgleichung, also steif oder nicht-steif, bekannt sein. Dies gilt insbesondere für die nicht-steifen Löser. Das Programm LSODA, siehe Petzold (1983a), und seine Modifikationen nimmt diese zu treffende Vorentscheidung und Vorabklassifikation dem Benutzer ab. Geschaltet wird in dem Programm LSODA zwischen den Adams-Formeln und den BDF. Die Schaltentscheidung wird maßgeblich vom Verhalten der Korrektoriteration abhängig gemacht, wobei allerdings auch das Verhalten während der Schrittweitensteuerung mit beachtet wird. Interessant ist, daß, obwohl das Programm LSODE sehr durchdacht und schon vergleichsweise verwickelt ist, es trotzdem möglich ist, eine weitere im gewissen Rahmen mögliche Leistungssteigerung hinzuzufügen. LSODA und LSODAR sind wie LSODE Bestandteil von ODEPACK, siehe Hindmarsh (1983). Aber auch NAG bietet einen schaltfähige Löser.

1. Kurzbeschreibung des Programmes LSODA

1. Die Programme LSODA und LSODAR benutzen als Vektornorm

wobei $w=(w_1,\ldots,w_2)$ ein Vektor von Gewichten ist. Die Division bei der Normbildung wird durch Multiplikation mit dem Kehrwert umgangen. Dies geschieht der Effizienz wegen. Man denke z.B. an zurückgewiesene Schritte. Dies ist jedoch im weiteren unerheblich. Wie in dem Programm TENDLER werden die Gewichte auf Positivität getestet. Bei $w_i=0$ drucken die Programme LSODA und LSODAR eine Fehlermeldung aus, die nicht unterdrückt werden kann, es sei denn durch Modifikation der beiden Programme. Die Kehrwerte werden vor jedem Schritt aktualisiert und zwar durch die Formel

wobei $\varepsilon_r$ die Toleranzanforderung des Benutzers an die relative Genauigkeit ist und $\varepsilon_a$ entsprechend die Toleranzanforderung an die absolute Genauigkeit ist.

Das Gewicht hängt also von den Genauigkeitswünschen des Benutzers ab. $\varepsilon_r$ und $\varepsilon_a$ können selber Vektoren sein (der Länge $n$). Dadurch ist es möglich für jeweils verschiedene Komponenten auch verschiedene Genauigkeitsanforderungen zu stellen, oder m.a.W., innerhalb der Normroutinen gewisse Komponenten verschieden stark zu gewichten. Dies setzt allerdings i.d.R. eine gewisse a priori Kenntnis der Lösung voraus.

Die benutzte Matrixnorm in dem Programmen LSODA und LSODAR ist

Die Normen wurden ggü. dem Programm LSODE geändert. Das Programm LSODE benutzt in der Routine vnorm die Rechenvorschrift

also die gewichtete RMS-Norm, mit allerdings dem gleichen Gewicht w.o. angegeben. Insbesondere hängt das Gewicht wieder von der Genauigkeitsanforderung $\varepsilon$ ab.

Von L.R. Petzold stammt auch das Programm DDASSL zur Lösung von differential-algebraischen Gleichungen, also Gleichungen, die neben den Differentialgleichungen noch zusätzlich algebraische Restriktionen enthalten können, u.U. nichtlinear gekoppelt. Damit sind auch implizite Differentialgleichungen erfasst. Bei der Lösung von differential-algebraischen Gleichungen handelt es sich um die Behandlung von Problemen der Form

Hierbei sind $F$, $t_0$, $y_0$ und $\dot y_0$ fest vorgegeben. Die Funktion $y$ ist gesucht. In dem Programm DDASSL wird als Vektornorm die RMS-Norm verwendet. Programmiert wird die Auswertung dieser Norm in der Form

Ist $M=0$, so wird $\left\Vert v\right\Vert$ zu Null gesetzt. Die Vorabberechnung der Maximumnorm des gewichteten Vektors geschieht aus Genauigkeits- und Sicherheitsgründen. Aus Sicherheitsgründen, da ohne die implizite Berechnung der Maximumnorm $\mathopen|v\mathclose|_{w,\infty}$, es möglicherweise zu einem Überlauf oder Unterlauf kommen könnte, und dies obwohl das Endergebnis auf der Rechenanlage darstellbar ist.

Viele Programme, so u.a. auch DISFUB und LSODE, schränken unnötigerweise den Darstellungsbereich stark ein, aufgrund der Tatsache, daß diese Programme der Normberechnung nicht die Aufmerksamkeit widmen, die die Normberechnung verdient. Allerdings kostet die Vorabberechnung der Maximumnorm zusätzliche Rechenleistung. Die Berechnung der Norm verbraucht einen nicht zu unterschätzenden Teil der gesamten Rechenzeit, bei der numerischen Lösung von Differentialgleichungen.

Bei linearen Differentialgleichungen mit konstanten Koeffizienten und konstanter Inhomogenität überwiegt die Rechenzeit für die Normberechnung bei weitem die Rechenzeit für die Funktionsauswertung.

In dem Programm DDASSL wird der Vektor mit den Komponenten $v_i/w_i$ nicht gespeichert, sondern muß zweimal berechnet werden. Die Speicherung hätte eine Erhöhung des Speicherplatzbedarfes erfordert. Deutlich wird, daß die Berechnung der RMS-Norm (korrekt programmiert) aufwendiger ist als die Berechnung der Maximumnorm.

2. In den beiden Programmen LSODA und LSODAR wird vollständig zwischen zwei verschiedenen Integrationsarten hin und her gewechselt.

Bei beiden Programmen werden die verwendeten linearen Mehrschrittverfahren (Adams-Formeln und BDF) stets im $P(EC)^i$-Modus ($2\le i\le3$) betrieben. Die Anzahl der Iterationsschritte wird abhängig gestaltet vom beobachteten Konvergenzverhalten. Insbesondere kann sich die Anzahl der Iterationen in jedem Schritt verändern. Die maximale Anzahl der Iterationen beträgt hier $2\le i\le3$, vom Sonderfall des Rechnens nahe der Grenzgenauigkeit abgesehen, wo auch $1\le i\le3$ gelten darf. Der Benutzer kann, anders als in dem Programm TENDLER, nicht im voraus angeben, wie oft iteriert werden soll, d.h. die Anzahl der Iterationen wird stets von dem Programm LSODA entschieden.

Das Programm LSODAR ist eine Erweiterung des Programmes LSODA; dieses Programm baut also auf dem Programm LSODA auf. Das Programm LSODAR bietet die Möglichkeit der Lösung von Differentialgleichungsproblemen mit impliziter Zeitvorgabe ($g$-stop-facility, oder auch implicit output genannt), d.h., diejenige Stelle $t$, zu der hinintegriert werden soll, ist nicht explizit bekannt, sondern ergibt sich als Lösung eines Gleichungssystemes, in denen die berechneten Lösungswerte $y(t)$ eingehen.

Dies tritt z.B. bei der Untersuchung von Grenzzyklen auf, wo man wissen möchte, wann die Lösung die Poincaré-Ebene das erste Mal trifft. Mit dem Programm LSODAR kann man also Probleme lösen, der Form

Man nennt das auch g-stop.

Biographisch:

1. Poincaré, Jules Henri (1854–1912)

3. Das Schalten zwischen den impliziten Adams-Verfahren und den BDF beruht maßgeblich auf den folgenden Überlegungen. Ausgangspunkt bei der Lösung der impliziten Differenzengleichung ist

Im steifen Modus werden die BDF bis zur Ordnung 5 benutzt, und es wird die Norm der Jacobimatrix $J$, also $K=\mathopen|J\mathclose|$, direkt berechnet. Grundannahme dieser Maßnahme ist, daß die Faktorisierung der Iterationsmatrix $W$ wesentlich aufwendiger ist, als die Berechnung der Matrixnorm, die der obigen Vektornorm zugeordnet ist. Dies ist für Differentialgleichungen mit $n>1$ der Fall. Da die Jacobimatrix $J$ über mehrere Schritte konstant gelassen wird, also nicht stets neu berechnet wird, muß man mit Fehlern rechnen.

Eine Neuberechnung der Jacobimatrix $J$ wird angeordnet

1. durch Entscheidungen während des Hindmarsh-Testes,
2. falls Konvergenzversagen auftrat und
3. mindestens auf jeden Fall alle 20 Schritte.

Da nun in den Programmen LSODA und LSODAR drei Maßnahmen vorhanden sind, die Jacobimatrix $J$ neu zu berechnen, falls sich diese Matrix “stark” ändert, nimmt man in beiden Programmen an, daß die Fehler, die man durch Gebrauch einer zu alten Jacobimatrix begeht, nicht ins Gewicht fallen.

Zudem wird durch Sicherheitsfaktoren eine Überbewertung vermieden, bzw. abgemildert.

4. Im nicht-steifen Falle werden die impliziten Adams-Verfahren mit den Adams-Bashforth Formeln als Prädiktor, von der Ordnung 1 bis zur Ordnung 12 benutzt. Im nicht-steifen Teil verschafft man sich durch Betrachten der Konvergenzrate während der Picard Iteration eine untere Näherung für den Spektralradius der Jacobimatrix $J$. Es ist nämlich

was direkt folgt durch Subtraktion der beiden Gleichungen

Wird mehr als zweimal iteriert, so nimmt man das Maximum der so gebildeten Konvergenzraten. Die so gewonnene untere Schranke für $L$ heiße $K$. Weitere generelle Informationen über Konvergenzraten im Zusammenhang mit der Korrektoriteration findet man bei Shampine (1980)2. Es wird mindestens zweimal pro Schritt mit dem Korrektor iteriert, dreimal jedoch höchstens nur, falls dies nötig ist.

Beim vierten Male wird Divergenz angenommen. Der Zwang mindestens zweimal zu iterieren, also mindestens im $P(EC)(EC)$ Modus zu arbeiten, verteuert die gesamte Integration merkbar, jedoch nicht wesentlich. Dies ist der (geringe) Preis, den man für das automatische Schalten in dem Programm LSODA bzw. LSODAR zu zahlen bereit sein muß. Die obige Rechenvorschrift zur Bestimmung einer unteren Schranke für $L$ kann in speziellen Fällen als von Misessches Iterationsverfahren für den dominaten Eigenwert aufgefasst werden. Ist nämlich

so folgt

Besitzt $M$ nicht-reelle dominante Eigenwerte, so kann bekanntlich die von Misessche Iteration unrichtige Ergebnisse liefern.

Biographisch:

1. Mises, Richard Edler von (1883–1953)
2. Lipschitz, Rudolf (1832–1903)

5. Für den nicht-steifen Teil wird nun verlangt, daß die Schrittweite $h_N$ für diesen Teil also, wie folgt beschränkt ist:

um genügend schnelle Konvergenz während der Fixpunktiteration zu erhalten; hier also mit einer Konvergenzrate von \frac1/2. $K$ ist hierbei die untere Schranke für die Lipschitzkonstante, w.o. beschrieben. Weiterhin wird verlangt, daß

Man fordert hierdurch, daß die Schrittweite $h_N$ im Stabilitätsgebiet, angegeben durch den Radius $r$, bleibt. Der Faktor 2 im Nenner ist ein Sicherheitfaktor. Diese letzte Bedingung wird nicht nur während des Schaltens gefordert, sondern überhaupt auch während des gesamten Rechnens im nicht-steifen Teil. Es gibt hier verschiedene Radien, die von Bedeutung sind. Zum einen hängt der Radius $r$ von der Ordnung $p$ ab und zum anderen von der Anzahl der Iterationen während der Picard Iteration. L.R. Petzold wählt hier den Radius für den $P(EC)E$-Modus in Abhängigkeit der Ordnung $p$, also $r=r_p$. Es sei nocheinmal darauf hingewiesen, daß im $P(EC)E$-Modus niemals gearbeitet wird, sondern stets nur im $P(EC)^i$-Modus. Der Wert dieses Radiuses $r$ sei keine kritische Größe in dem Progamm, solange er von der richtigen Größenordnung sei, schreibt Petzold (1983a).

6. Die Entscheidung des Wechselns wird nun im Schrittweiten- und Ordnungssteuerungssegment gefällt. Man schaltet von den Adams-Formeln auf die BDF, wenn

1. die Schrittweite $h_N$ für den nächsten Schritt im nicht-steifen Falle günstiger ist als die Schrittweite $h_S$ für den nächsten Schritt mit den BDF,
2. die obigen beiden Restriktionen für die Schrittweite im nicht-steifen Falle erfüllt sind, und
3. die Ordnung nicht höher als 5 ist. Ist umgekehrt die Schrittweite für die BDF fünfmal günstiger als die Schrittweite für die Adams-Formeln, so schaltet man von den BDF auf die Adams-Formeln.

Man beachte, daß das Programm LSODA tatsächlich beim Schalten verschiedene Integrationsmethoden auswählt. Deshalb ergibt es Sinn, den Schalttest im Schrittweiten- und Ordnungssteuerungs-Segment auszuführen. Bei einem Wechsel werden stets die Koeffizienten des Verfahrens neu berechnet, so, wie dies das Programm LSODE auch bei einem gänzlich neuen Integrationsbeginn tut. Man beachte, daß z.B. der erste Schritt mit einer BDF, nach Benutzung des Adams-Verfahrens, kein “echter” Schritt mit einer BDF ist. Da jedoch Maßnahmen in dem Programm LSODA vorhanden sind, um zu häufiges Wechseln zu vermeiden, ist dieser Umstand nicht von übergeordneter Bedeutung. Er zeigt jedoch eine grundsätzliche Problematik auf.

Es gibt eine Reihe von Vorsichtsmaßnahmen, die weiterhin getroffen wurden. Um möglichen Instabilitäten bedingt durch häufiges Schalten zumindestens ansatzweise entgegenzutreten, wird verlangt, daß höchstens alle 20 Schritte ein Schalten in Betracht gezogen wird.

Weiterhin wird der Fall, daß an der Grenzgenauigkeit der Rechenanlage gearbeitet wird, viel Aufmerksamkeit gewidmet, und es werden hier besondere Vorkehrungen getroffen, die stellenweise auch das Schalten berühren.

7. Das Programm LSODA (und natürlich auch LSODAR) benutzen einen anderen Konvergenztest, als dies das Programm TENDLER tut. Der Test, der zwischen den Marken 400 und 410 in dem Unterprogramm stoda durchgeführt wird, ist ggü. dem des Programmes LSODE unwesentlich verändert. Die Konvergenzrate darf jede Iteration, aber erst natürlich nur im zweiten Iterationsschritt, höchstens um den Faktor 1/5 abnehmen.

Der eigentliche Konvergenztest ist nun

wobei $r$ die Konvergenzrate ist und ganz zu Anfang der Integration auf $0.7$ gesetzt wird. Die Konstante $\hbox{const}$ hängt von der Genauigkeit $\varepsilon$ und den Diskretisierungsfehlern der Verfahren ab. Der Konvergenztest ist insoweit modifiziert worden, als daß eine zusätzliche Normberechnung eingeführt wurde um weitere Iterationen nahe der Grenzgenauigkeit zu vermeiden.

2. LSODAR und Programmtransparenz

1. Es wurde nun untersucht, in wie weit Bemühungen unternommen worden sind, das Programm LSODAR lesbar und modular zu gestalten. Als ein Indiz, von tatsächlich sehr vielen Indizien, kann man die Länge der Parameterlisten werten. Zum Begriff der Programmlesbarkeit sollte man hinzufügen, daß dies teilweise eine Frage des persönlichen Geschmackes ist. Dies betrifft insbesondere Benennung von Programmvariablen, Einrückungstiefe, Absatztiefe, etc. Allerdings ist und bleibt ein einziges Konglomerat von Unterprogrammen, einhergehend mit einem schlecht durchdachten Konzept, ein Programm, welches dem Leser Schwierigkeiten bereiten wird, beim Verständnis und insbesondere bei der Modifikation.

Auch Tischer/Gupta/Maeder (1986) betrachten Differentialgleichungslöser nicht nur unter einem rein mathematischen Gesichtspunkt. Es gibt eine Reihe von Programm-Entwurfsverfahren, u.a. das Modul-Schnittstellen-Konzept. Hierzu schreibt Richter (1977):

Die Aufgaben $\ldots$ werden zu Beginn der Entwurfsphase so detailliert wie möglich in Funktionen zerlegt, und jede Funktion wird durch ein Modul beschrieben, dessen Schnittstellen zu angrenzenden Modulen definiert werden. Hauptvorteil dieses Verfahrens ist der geringe initiale Planungsaufwand, der beträchtlich unter dem der vorangehenden Verfahren liegt.

Weiter äußert sich nun Richter zu den Nachteilen:

Die Zeitvorgaben bei der Implementierung eines solchen Systemes können schwer eingehalten werden. Häufig wird die Komplexität der Module auf Kosten der Komplexität der Schnittstellen minimiert. Nach diesem Verfahren realisierte Systeme verhalten sich gegen Änderungen außerordentlich sensitiv, da bereits relativ geringfügige Modifikationen unter Umständen schon ziemlich viele Module betreffen können.

Man beachte, daß das Programm LSODAR eher nach dem Prinzip der Nukleus-Erweiterung entworfen und entwickelt wurde, mit LSODE als Nukleus. Die hauptsächlichen Entwurfsentscheidungen waren also damit schon getroffen worden.

2. Weiterhin hängt die Parameterliste aber auch ab von der Aufgabenstellung, von der Problemschwierigkeit und von der Programmiersprache. Es kann weiterhin vorkommen, daß ein Programm zwar zahlreiche Unterprogramme mit sehr langen Parameterlisten enthält (statische Sichtweise), diese jedoch bei der Programmausführung keine Rolle spielen (dynamische Sichtweise), oder aber in etwa gleich viele Unterprogramme mit langen und kurzen Parameterlisten vorhanden sind, jedoch maßgeblich nur die Unterprogramme mit den kurzen Parameterlisten für das Programm wichtig sind, die anderen Unterprogramme also lediglich z.B. Spezialfälle berücksichtigen (problemorientierte Sichtweise). Der Effekt des Einflusses der Programmiersprache auf das Gesamtprogramm wird deutlich anhand des Programmes Dhrystone, siehe Weicker (1984), welches einmal in der Programmiersprache Ada und einmal in der Programmiersprache Pascal codiert wurde. Es existiert auch eine C Version.

Der Rest verteilt sich auf statische Variablen, formale Parameter (in, out, in out), etc.

Bei neuzeitlichen Systemprogrammen trifft man es häufig an, daß die Parameterlisten i.d.R. kurz sind, man vgl. hierzu z.B. Tanenbaum (2005). Bei der Unterprogrammsammlung NAG hingegen ist es üblich, daß die Parameterlisten eher vergleichsweise lang sind.

Tanenbaum (2005) gibt die folgenden Werte für ausgewählte Systemprogramme an:

3. Stellenweise unter diesem Aspekt sollen daher die Schnittstellen des Programmes LSODAR vorgestellt werden. Untersucht wird hier die Version vom 30.März.1987. Enthalten in dem Programm LSODAR sind insgesamt 4 benannte common-Blöcke, mit

Die wichtigsten Unterprogramme im Gesamtprogramm LSODAR haben nun die folgende Anzahl an Parametern:

Auf die LINPACK und BLAS Routinen dgefa, dgesl, dgbfa, daxpy, dscal, idamax, ddot, dcopy und d1mach wird hier nicht eingegangen, da sie ausserhalb der Entscheidungen standen, beim Entwurf des Programmes LSODAR. Die durchschnittliche Parameterlänge für die einzelnen Unterprogramme in dem Programm LSODAR beträgt nun knapp 8 Paramter pro Unterprogramm bzw. Funktion, mit einer Standardabweichung von rund 5.

Auffällig sind hier die maßgeblichen Routinen lsodar, stoda und roots, die die hauptsächlichen algorithmischen Leistungen am gesamten Programm durchführen.

3. Bewertung und Vergleich des Programmes LSODA

1. Nicht-steife Gleichungen.

Für die nicht-steifen Testdifferentialgleichungen aus DETEST ergaben sich für das Programm LSODA die folgenden Ergebnisse. Unter den jeweiligen Genauigkeitsanforderungen $10^{-3}$, $10^{-6}$ und $10^{-9}$ sind sämtliche Laufzeiten zusammengefasst. Die Vergleiche wurden auf einer Rechenanlage vom Typ CDC 6600, in wie üblich einfacher Genauigkeit, durchgeführt. Alle Differentialgleichungen wurden mit einer Startschrittweite von $h_0=10^{-12}$ integriert. Beide Programme erzielten vergleichbare globale Fehler.

Für das Programm LSODE mit der Einstellung mf=10, wurden die Laufzeiten wie angegeben ermittelt. Die Einstellung mf=10 benutzt die Adams-Verfahren im $P(EC)^i$-Modus, mit $i\in\{1,2,3\}$.

Für das Programm LSODE mit der Einstellung mf=22, ergaben sich nun die folgenden Werte. Es wurde also mit BDF, gekoppelt mit der Newton Iteration und einer Jacobimatrix, die durch numerische Differentiation berechnet wurde, gearbeitet.

Vergleicht man nur die Rechenzeiten, also diejenige Größe, die bei sehr großen Differentialgleichungen von entscheidender Bedeutung wird, so beobachtet man, daß das Programm LSODE mit mf=10 leicht besser abschneidet, als das Programm LSODA. Dennoch muß der Benutzer um diese Einstellung im voraus wissen. Bei der Einstellung des für den Benutzer am bequemsten Wertes mf=22, erkennt man, daß sich die Rechenzeit mehr als verdoppelt hat.

Der Zwang mindestens zweimal im nicht-steifen Falle mit Hilfe der Picard Iteration zu iterieren, bemerkt man an der Anzahl der Funktionsauswertungen des Programmes LSODA. Diese sind ggü. dem Programm LSODE um über 40% gestiegen. Die Gesamtrechenzeit ist jedoch nicht um diesen Anteil gestiegen. Dies liegt u.a. daran, daß das Programm LSODA offensichtlich die größeren Schrittweiten wählen kann, wie man an der Spalte der Anzahl der Schritte erkennen kann, allerdings auch daran, daß die Jacobimatrix durch numerische Differentiation bestimmt wurde. Hier zeigt sich, daß das Programm LSODA allen anderen Einstellungen eindeutig überlegen ist. Diese Aussage gilt natürlich wohlgemerkt nur für die betrachteten Testdifferentialgleichungen, die jedoch einen recht repräsentativen Querschnitt bilden. Dies heißt somit nicht, daß grundsätzlich, für jede beliebige Differentialgleichung LSODA dem Programm LSODE überlegen ist, sondern zeigt lediglich, daß für eine recht große Klasse von Problemen, das Schalten sehr sinnvoll sein kann.

Weiter fällt auf, daß sich bei sehr scharfer Genauigkeitsanforderung, die Rechenzeitunterschiede zwischen den beiden Programmen LSODA und LSODE mit mf=10 verwischen. Mit diesen Zahlen im Hinterkopf muß man das vom Programm LSODA gewählte Schaltschema als durchaus sehr erforlgreich bezeichnen. Bei allen zahlenmässigen Vor- und Nachteilen, beachte man den Bequemlichkeitsgewinn, denn man durch Einsatz des Programmes LSODA erhält. Die a priori Klassifikation als steif oder nicht-steif wird automatisch vom Programm erledigt und braucht nun nicht mehr vom Benutzer vorgenommen zu werden. Man beachte auch den Sicherheitsgewinn. Es kann einem Benutzer nun nicht mehr so einfach passieren, daß er einen Löser, konzipiert für nicht-steife Differentialgleichungen, ansetzt auf steife Gleichungen. Es ist dieser letzte Fall, der zu besonders untragbaren Rechenzeiten führt.

Nebenbei sei erwähnt, daß es bei nicht-steifen Differentialgleichungen zumindestens akzeptabel ist, diese Gleichungen mit einem Programm, konzipiert für steife Differentialgleichungen, zu lösen.

Eine Verdoppelung oder ggf. mehr, der Rechenzeiten ist bei durchweg geringen Laufzeiten, im Bereich weniger Sekunden, durchaus annehmbar. Die längsten Wartezeiten fallen dann sowieso in den Bereich der Ein- und Ausgabe. Bei großdimensionalen nicht-steifen Differentialgleichungen ist ein steifer Löser ggü. einem nicht-steifem Löser schon eher benachteiligt. Bei nicht-steifen Differentialgleichungen verändert sich die Jacobimatrix häufig signifikant, zudem sind u.U. größere Schrittweitenvariationen vonnöten, als bei steifen Gleichungen. Dies wirkt sich gravierend auf die Faktorisierungsarbeit aus, die der steife Löser zu leisten hat.

2. Steife Gleichungen.

Das Programm LSODA wurde auch an den steifen Testdifferentialgleichungen von DETEST erprobt und mit dem Programm LSODE (mf=22) verglichen. Petzold (1983a) schreibt hier, daß bei den Differentialgleichungen B5 und E4 geringfügige Schwierigkeiten beobachtet wurden. Insbesondere das Problem B5 hat Eigenwerte der (konstanten) Jacobimatrix mit $-10\pm100i$. Die Probleme werden damit erklärt, daß die untere Schranke für den Spektralradius flukturierte.

Bei echt komplexen (nicht reellen) Eigenwerten liefert bekanntermaßen das von Misessche Iterationsverfahren für den dominanten Eigenwert weniger brauchbare Ergebnisse. Für das Programm LSODA ergab sich nun hierbei:

Für das Programm LSODE mit der Einstellung mf=22 ergaben sich die folgenden Resultate.

Mit dem Programm LSODE mit der Parametereinstellung mf=10 wurde natürlich nicht verglichen, da dies zu rechenzeitaufwendig gewesen wäre. Ein Programm, konzipiert für nicht-steife Differentialgleichungen, welches man auf steife Probleme ansetzt, braucht exzessiv viel Rechenzeit und liefert dann dennoch nicht hochgenaue Ergebnisse.

Man erkennt an den beiden letzten Tabellen, daß die Rechenzeit für das schaltfähige Programm ggü. dem nicht schaltfähigen Programm geringer ist.

Dies ist zwar nur gering, lediglich rund 20%, man muß jedoch behalten, daß die Testdifferentialgleichungen aus STIFF-DETEST wirklich sämtlich steif sind, abgesehen von Ausnahmen mit nur “geringer” Steifheit, wie z.B. das Problem B3 und B4. Die vollen Rechenzeitvorteile eines schaltfähigen Programmes kommen daher potentiell nicht zur gänzlichen Geltung. Die Ersparnisse hängen daher mehr oder weniger stark von der Startphase ab.

Erneut sieht man, daß die Anzahl der Funktionsauswertungen des Programmes LSODA ggü. LSODE (mf=22), gerade bei den höheren Genauigkeitsanforderungen, leicht ungünstiger ausfällt. Die Anzahl der Jacobimatrixauswertungen des Programmes LSODA ist immer geringer, als die entsprechende Zahl bei dem Programm LSODE. Ebenso verhält es sich mit den Schritten, die zur Vollendung der Integration benötigt werden. Hier kann das Programm LSODA, über das gesamte Integrationsintervall verteilt gesehen, größere Schrittweiten wählen, als dies das Programm LSODE tun kann.

3. Um auch die Unterschiede zu beleuchten, welche die beiden Programme liefern würden, wenn sie mit einem Differentialgleichungsproblem konfrontiert werden, welches den Charakter, also steif oder nicht-steif, während der Integration häufiger ändert, wurden Ergebnisse gemessen für die

van der Polsche Differentialgleichung der Form

Hier ergaben sich für das Programm LSODA die Werte:

und für das unmodifizierte Programm LSODE fand Petzold:

Deutlich wird die geringere Anzahl an Schritten. Die Anzahl an Jacobimatrixauswertungen ist besonders gesunken ggü. dem Programm LSODE. Die Anzahl der Funktionsauswertungen sinkt nicht im gleichen Maße, wie die Anzahl der Schritte, was verständlich ist, wenn man beachtet, daß das Programm LSODA im nicht-steifen Falle stets zwei Iterationsschritte durchführt, um eine Näherung der Lipschitzkonstanten zu erhalten. Das Programm LSODA arbeitet also im nicht-steifen Falle immer mindestens im $P(EC)(EC)$-Modus. Allerdings sollten die obigen letzten Ergebnisse nicht übergewichtet werden, da nicht auszuschliessen ist, daß das Programm LSODA auf diese Differentialgleichung hin getrimmt wurde. Bei anderen Differentialgleichungen mit wechselnder Steifheit lassen sich solche eindrucksvollen Ergebnisse nicht erzielen, wie die nachfolgenden Ausführungen zeigen.

4. Nicht immer ist das Programm LSODA der gewöhnlichen Nutzung des Programmes LSODE (mf=21) überlegen. Sogar bei einem Differentialgleichungsproblem, welches seinen Charakter ändert, und zwar einmalig von steif nach nicht-steif, kann das Programm LSODE ggü. dem Programm LSODA dennoch (leichte) Vorteile aufweisen. Byrne/Hindmarsh (1987) geben für das Problem P5 folgende Rechenzeiten auf einer Rechenanlage Cray-1S an, wobei

Es ergab sich nun:

mf=21 bei LSODE bedeutet Benutzung der BDF im Newton-Kantorovich Iterationsmodus, wobei die Jacobimatrix vom Benutzer bereitgestellt wird, und mf=10 bedeutet wie üblich Adams-Formeln mit Picard Iteration. jt=1 bei LSODA zeigt an, daß die Jacobimatrix vom Benutzer bereitgestellt wurde.

mf=21$\to$10 gibt an, daß vom Benutzer bei $t=4.9\cdot10^5$ ein entsprechender Wechsel der Variablen mf vorgenommen wurde. Dies setzt natürlich voraus, daß man das Verhalten der Lösung der Differentialglewichung P5 schon im voraus kennt, oder zumindestens einen gewissen Überblick besitzt.

Bei den obigen Angaben bzgl. der Differentialgleichungen von DETEST, bzw. NSDTST und STDTST, wurden die Jacobimatrizen stets durch numerische Differentiation gewonnen und damit nivellierten sich die Unterschiede zwischen den beiden Programmen LSODA und LSODE im großen und ganzen aus, aufgrund der Startvorteile eines schaltfähigen Programmes basierend auf linearen Mehrschrittverfahren, welches mit sehr kleiner Schrittweite startet. Bei dem zweidimensionalen Problem P5 ist die Auswertung der Jacobimatrix sehr einfach, sodaß der Zwang von LSODA stets im $P(EC)(EC)$-Modus zu arbeiten, hier bei diesen Gegebenheiten durchschlägt.

Würde man das Problem P5 ein zweites Mal lösen, wobei die Jacobimatrix diesmal durch numerische Differentiation gewonnen würde, so würde man erwarten, daß sich erneut die Laufzeitunterschiede verkleinern. In diesem Falle jedoch liegt ein “Versagen” des Programmes LSODA vor, da dieses Programm nicht erkennt, daß man den hochoszillatorischen Verlauf der Lösung nach ca. $t=4.5\cdot10^5$ am günstigsten mit den Adams-Verfahren mit Picard Iteration erhält. “Versagen” allerdings hier nur in dem Sinne, daß ein Vorteil nicht genutzt wurde.

Die Wahl des steifen Modus führt ja durchaus zu einer korrekten Lösung. Die Rechenzeiten von LSODE (mf=21) zeigen allerdings, daß man auch mit den BDF im Newton-Kantorovich Iterationsmodus günstigere Laufzeiten erhalten kann. Jedoch sind solche Überlegungen bei Rechenzeiten im Subsekundenbereich von eher geringerem Interesse, wobei man aber die zugrunde liegende Rechenanlage zu berücksichtigen hat. Auf einer kleineren Rechenanlage stellt sich das Problem P5 als durchaus aufwendig dar, was auch an dem vergleichsweise langen Integrationsintervall liegt.

Da bei der Simulation großer elektrischer Schaltkreise ein solches sporadisch auftauchendes, oszillatorisches Verhalten möglich und realistisch ist, deutet das Differentialgleichungsproblem P5 eine gewisse Schwachstelle der zugrundeliegenden Schaltentscheidung an. Dennoch zeigt sich bei einem Vergleich mehrerer schaltfähiger Programme, daß von diesen das Programm LSODA zu den Wirkungsvollsten gehört, siehe Bruder/Strehmel/Weiner (1988).

4. LSODA im Vergleich mit anderen Lösern aus ODEPACK

Für die Rechenanlagen CDC-7600 und Cray-1 wurden mit den vier Programmen LSODE, LSODA, LSODES und GEARBI das Differentialgleichungsproblem P8 gelöst. Alle vier Programme verwenden die völlig gleichen Formeln, nämlich die BDF$i$, mit $i=1,2,3,4,5$, wobei LSODA zusätzlich die Möglichkeit hat die Adams-Verfahren zu benutzen, falls das Programm LSODA glaubt, daß die Differentialgleichung im Augenblick nicht-steif ist. Das Problem P8 wurde einmal mit einer Ortsdiskretisierung von $10\times10$ und einmal mit einer von $20\times20$ gelöst. Zu den nachfolgenden Tabellen und Erläuterungen vgl. man Hindmarsh (1983).

Bei den nachstehenden Tabellen bedeutet 'CPU' die verbrauchte CPU-Zeit in Sekunden auf dem angegebenen Rechner. nstep gibt die Anzahl der benötigten Schritte, nfe gibt die Anzahl der Funktionsauswertungen an, nje entsprechend die Anzahl der Jacobimatrix Auswertungen. Die Spalte 'LU' weist die Anzahl der $LU$-Faktorisierungen aus und schließlich 'mem' informiert über den benötigten Speicherplatzbedarf, welcher gerade bei großdimensionalen Problemen, wie hier P8, von entscheidender Bedeutung ist. Die Abkürzung 'usj' (user supplied jacobian) gibt an, daß die Jacobimatrix vom Benutzer programmiert wurde und 'jdq' (jacobian via difference quotient) gibt an, daß die Jacobimatrix intern durch numerische Differentiation, mit Hilfe eines Differenzenquotienten, berechnet wurde.

Als erstes die Ergebnisse für P8 mit einem $10\times10$-Gitter auf der Rechenanlage CDC-7600.

Anschliessend wurde die gleiche Differentialgleichung, also P8 mit $10\times10$-Ortsgitter, auf der Rechenanlage ^{Cray-1} gelöst.

Zuletzt die Ergebnisse für die Differentialgleichung P8 mit einem $20\times20$-Ortsgitter auf der Rechenanlage Cray-1.

Anhand dieser und weiterer Daten gelangt Hindmarsh (1983) zu den nachstehenden Überlegungen.

1. Für jedes der beiden Probleme variiert die Anzahl der Schritte nicht stark unter den einzelnen Lösern. Die Anzahl der Schritte wird maßgeblich von der Genauigkeitsanforderung bestimmt. Durch Vergleich der Ergebnisse für das $20\times20$-Gitter mit den Ergebnissen für das $10\times10$-Gitter erkennt man, daß letztere in etwa einen Fehler aufgrund der Ortsdiskretisierung von 2% aufweisen.

2. Für jedes Problem sind die Leistungsdaten für die beiden Programme LSODE und LSODA in etwa gleich, wobei LSODA häufig günstigere Werte aufzuweisen hat. Den größten Gewinn erzielt das Programm LSODA während der Startphase, da diese mit den wesentlich weniger aufwendigen Adams-Verfahren im Picard Iterationsmodus gelöst werden. Da jedoch das Programm LSODA für die Schätzung einer Näherung für die Lipschitzkonstante bei jedem Schritt mindestens zwei Korrektoriterationen und damit Funktionsauswertungen einschließlich Normberechnungen durchführen muß, verwischen sich die Unterschiede insgesamt.

3. Für die beiden Programme LSODA und LSODE bedeutet die Benutzung der Jacobimatrix, generiert durch numerische Differentiation, einen deutlichen Mehraufwand. Auf der Rechenanlage CDC-7600 beträgt dieser Mehraufwand nicht mehr als 25%, während hingegen auf der Rechenanlage Cray-1 liegt der Mehraufwand zwischen 65% und 118%. Dies liegt u.a. maßgeblich an der Möglichkeit zur Vektorisierung. Z.B. ist auf der Cray-1 die Bandelimination 10-mal schneller als auf der CDC-7600. Die LINPACK-Routinen sind u.a. speziell auf die Vektorrechenfähigkeiten der Cray-1 angepaßt. Gleichzeitig sind die Funktionsauswertungen auf der Cray-1 nur doppelt so schnell wie auf der CDC-7600. Daher nehmen die Funktionsauswertungen und damit die Jacobimatrixauswertungen durch numerische Differentiation einen größeren prozentualen Teil der Gesamtrechenzeit auf sich. Z.B. beträgt für LSODE (usj) auf dem $10\times10$-Gitter der Anteil der Rechenzeit für die Funktionsauswertungen an der Gesamtrechenzeit etwa 4% auf der CDC-7600, hingegen jedoch 19% auf der Cray-1.

4. Auf der Rechenanlage CDC-7600 lässt sich durch die Benutzung des Programmes LSODES ein Geschwindigkeitsvorteil von $1.6$ bis $2.1$ erzielen. Auf der Cray-1 ergibt sich kein, oder ein nur vernachlässigbarer Vorteil. Im Gegensatz zu der Bandelimination, welche auf der Cray-1 einen großen Geschwindigkeitszuwachs verzeichnet, liegt der Vorteil bei der Benutzung von speziellen Algorithmen für dünn-besetzte Matrizen auf der Cray-1 bei lediglich 2 zu 1. Die Geschwindigkeitsvorteile des Programmes LSODES auf der CDC-7600 liegen maßgeblich an der Art und Weise der Speicherung der Iterationsmatrix $W=I-h\gamma J$. Dadurch lassen sich zahlreiche Jacobimatrixauswertungen einsparen, welche auf der CDC-7600 einen maßgeblichen Rechenzeitanteil haben. Bei dem Programm LSODES wird die Jacobimatrix $J$ zwischen 26 und 49 Schritte hintereinander mehrmals benutzt. Hingegen bei den Programmen LSODA und LSODE wird die Jacobimatrix $J$ nur zwischen 3 und 6 Schritte hintereinander mehrmals benutzt. Man beachte, daß bei den beiden Programmen LSODA und LSODE eine Refaktorisierung der Iterationsmatrix $W$ automatisch auch zu einer Neuauswertung der Jacobimatrix $J$ führt. Für Probleme, die ähnlich sind wie das Problem P8 und bei denen zusätzlich die Funktionsauswertungen noch aufwendiger sind, empfiehlt sich daher eher die Benutzung von LSODES.

5. Auf beiden Rechenanlagen und für beide Ortsgitter beträgt der Mehraufwand für die Generierung der Jacobimatrix $J$ durch numerische Differentiation für das Programm LSODES weniger als 10%. Dies liegt zum einen daran, daß das Programm LSODES wesentlich weniger Jacobimatrixauswertungen benötigt und zum anderen daran, daß die Berechnung der Jacobimatrix durch numerische Differentiation nur 8 Funktionsauswertungen kosten und zwar unabhängig von der Gitterweite der Ortsdiskretisierung.

6. Der Speicherplatzbedarf für das Programm LSODES ist insbesondere für feine Orstdiskretisierungen geringer als bei den beiden Programmen LSODA und LSODE. Der Speicherplatzminderbedarf beträgt auf dem $10\times10$-Gitter etwa 12% und 41% auf dem $20\times20$-Gitter. Für sehr grobe Gitter hat das Programm LSODES keinerlei Speicherplatzvorteil, ja sogar einen Speicherplatzmehrbedarf, aufgrund der benötigten Zusatzinformationen für die Besetztheitsstrukturen und aufgrund der Tatsache, daß sowohl $W$ als auch die $LU$-Zerlegung von $W$ getrennt gespeichert werden. Wegen der Dünnbesetztheit bedeutet das Speichern von zwei quadratischen Matrizen keinen unvertretbaren Speichermehraufwand. Dies weist darauf hin, daß es zweckmässig ist, die Speicherung von zwei Matrizen in dem Programm TENDLER beizubehalten (unter Berücksichtigung spezieller Eigenheiten der Besetztheit).

7. Bei den hier vorgestellten Ergebnissen verzeichnet das älteste Programm GEARBI die besten Ergebnisse. Aufgrund der regelmässigen Blockgestalt der Jacobimatrix $J$, kann die Block-SOR Iteration bei dem Programm GEARBI hier großen Nutzen ableiten. Die $LU$-Zerlegungen sind hier lediglich die Zerlegungen der $2\times2$-Blöcke. Die Gesamtanzahl der Block-SOR Iterationen auf dem $10\times10$-Gitter betrug 607, also im Mittel weniger als 2 Iterationen pro Schritt. Für das $20\times20$-Ortsgitter stieg die Gesamtanzahl der Block-SOR Iterationen auf 927, also im Mittel $2.7$ Iterationen pro Schritt. Aufgrund des geringen Vektorisierungsgrades in dem Programm GEARBI, sinkt die Rechenzeit bei einem Wechsel der Rechenanlage von der CDC-7600 auf die Cray-1 nur um den Faktor $2.1$. Insbesondere hat das Programm GEARBI ggü. dem Programm LSODA nur noch geringere Rechenzeitvorteile, wenn überhaupt. Der Speicherplatzvorteil des Programmes GEARBI ist jedoch in jedem Falle enorm. Auf dem $10\times10$-Gitter beträgt der Gewinn $4.7$, und auf dem feineren Ortsgitter von $20\times20$ steigt der Vorteil auf $8.7$ an. Für ähnliche Differentialgleichungen wie P8, welche SOR Iterationsverfahren nahelegen, ist es angebracht, diese auch tatsächlich zu benutzen.

5. Literatur

1. Bruder, Jürgen und Strehmel, Karl und Weiner, Rüdiger: “Partitioned Adaptive Runge-Kutta Methods for the Solution of Nonstiff and Stiff Systems”, Numerische Mathematik, Vol 52, Fasc 6, 1988, pp.621–638
2. Byrne, George Dennis
3. Byrne, George D. und Hindmarsh, Alan C.: “Stiff ODE Solvers: A Review of Current and Coming Attractions”, Journal of Computational Physics, Vol 70, No 1, May 1987, pp.5—62
4. Hindmarsh, Alan C. (*1942)
5. 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
6. Petzold, Linda Ruth (*1954)
7. Petzold, Linda Ruth: “Automatic Selection of Methods for Solving Stiff and Non-Stiff Systems of Ordinary Differential Equations”, SIAM Journal on Scientific and Statistical Computing, Vol 4, No 1, March 1983, pp.136–148
8. Richter, Lutz (1936–2017)
9. Richter, Lutz: “Betriebssysteme”, B.G. Teubner Verlag, Stuttgart, 1977, 152 S.
10. Seifert, Peter (*1941)
11. Seifert, Peter: “Computational Experiments with Algorithms for Stiff ODEs”, Computing, Vol 38, 1987, pp.163–176
12. Shampine, Lawrence Fred, Math Genealogy
13. Shampine, Lawrence Fred: “Implementation of Implicit Formulas for the Solution of ODEs”, SIAM Journal on Scientific and Statistical Computing, Vol 1, No 1, March 1980, pp.103–118
14. Tanenbaum, Andrew Stuart (*1944)
15. Tanenbaum, Andrew Stuart: Operating Systems: Design and Implementation
16. Tanenbaum, Andrew Stuart: “Operating Systems: Design and Implementation”, Prentice-Hall, Englewood Cliffs (New Jersey), 1987, xvi+719 S.
17. Tischer, Peter E. und Gupta, Gopal K. und Maeder, A.J.: “A Software Engineering Approach for ODE Solvers”, in “Computational Techniques and Applications”, CTAC-85, Editor John Noye and Robert May, North Holland, Amsterdam New York Oxford Tokyo, 1986, pp.299–308
18. Weicker, Reinhold P.: “Dhrystone: A Synthetic Systems Programming Benchmark”, CACM, Vol 27, No 10, October 1984, pp.1013–1030