, 21 min read
Erfahrungsmaterial mit bisherigen Lösern
- 1. Strategien und Rechnergenauigkeit
- 2. STINT vers. GEAR und EPISODE
- 3. STINT versus LSODE
- 4. Rationalisierte Schrittweiten und dessen Einfluß
Es gibt eine Fülle von Differentialgleichungslösern, sowohl für steife, als auch für nicht-steife Gleichungen. Allerdings sind nicht alle Löser auch gleich gut für jeden Zweck geeignet. Es gibt auch Löser für Gleichungen mit speziellen Eigenschaften. Teilweise konnte man mit den Programmen zur Lösung von gewöhnlichen Differentialgleichungen besondere Auffälligkeiten feststellen, welche nicht nur löserspezifisch sind, sondern allgemeiner beobachtet werden konnte. Hier wird u.a. über Testergebnisse berichtet, die mit dem Programm STINT durchgeführt wurden. Diese Ergebnisse sind auch für das Programm TENDLER von Interesse, da ja beide Programme, wie das Programm DIFJMT, die gleichen zyklischen Formeln verwenden. Zudem kamen zahlreiche Ergebnisse und Erfahrungsberichte, die mit dem Programm STINT gewonnen wurden, dem Programm TENDLER zugute.
Es zeigt sich, daß mit dem Programm STINT ein Löser vorliegt, der sich mit anderen bekannten Lösern in seiner Leistungsfähigkeit durchaus messen lassen kann. Dennoch offenbart das Programm STINT einige besonders hervorhebenswerte Andersartigkeiten, denen hier Raum gegeben wird. Zuerst wird die Unstetigkeit der Schrittweiten- und Ordnungssteuerung anhand bekannter Testergebnissen dargestellt. Anschliessend werden die Testdaten von Tendler/Bickart/Picel (1978) in knapper Form angegeben und interpretiert. Schließlich wird auf die Testläufe von Gaffney (1984) eingegangen. Zum Abschluß wird ein wenig bekannter Effekt von rationalisierter Schrittweite anhand eines Beispieles dargestellt.
1. Strategien und Rechnergenauigkeit
Als Testdifferentialgleichung sei gewählt das abgewandelte Kontrollproblem, zurückgehend auf E.J. Davison,
mit
wobei hier $a={3/2}$, $b=1/10$ und $c=1/100$ ist. Die $80\times80$ Matrix $A$ enthält also auf der Diagonalen $a_{ii}=-a^{80-i}$, $i=1,\ldots80$, auf der Subdiagonalen $a_{i+1,i}=b$ und der Rest ist gleichmässig mit $c$ gefüllt. Die Steuerung des Problems geschieht mit
Die Funktion $u({}\cdot{})$ ist hierbei eine Näherung für eine Rechteckkurve mit der Amplitude 1 und der Periode $2\pi$. Es gibt bessere Approximationen für diese Funktion, die jedoch weiterhin die Unstetigkeit vermeiden.
Enright (1978) und Enright/Kamel (1979) geben nun die folgenden Meßresultate an für das Programm EPISODE und MOD-EPISODE.
Letzteres benutzt Enright's update technique.
| Problem | $\varepsilon$ | Zeit A |
Zeit B |
#Fkt A |
#Fkt B |
#Jac A |
#Jac B |
|---|---|---|---|---|---|---|---|
| EPISODE | $10^{-2}$ | 5.72 | 10.27 | 44 | 43 | 11 | 13 |
| " | $10^{-4}$ | 20.19 | 39.15 | 322 | 350 | 31 | 44 |
| " | $10^{-6}$ | 65.36 | 77.48 | 817 | 743 | 107 | 84 |
| MOD-EPISODE | $10^{-2}$ | 2.82 | 2.68 | 42 | 41 | 1 | 1 |
| " | $10^{-4}$ | 11.42 | 13.40 | 250 | 272 | 1 | 1 |
| " | $10^{-6}$ | 31.81 | 32.53 | 670 | 597 | 2 | 4 |
Der Buchstabe A kennzeichnet die Rechenanlage IBM 360/165 und der Buchstabe B bezeichnet die Anlage IBM 3033. Benutzte Rechengenauigkeit, verwendeter Übersetzer und globale Fehler werden in Enright (1978) und Enright/Kamel (1979) nicht genannt. Man ersieht deutlich, daß bei einem Wechsel der Rechenanlage die Ergebnisse sich verändern und dies sogar bei Rechenanlagen vom gleichen Hersteller. Sind also die beiden Steuerungen einmal aus dem Takt, so finden sie nicht mehr zusammen. Das gleiche gilt auch für die Korrektoriteration und die Fehlerkontrolle, wo ebenfalls unstetige Strategien zur Anwendung gelangen.
Ähnliche Ergebnisse wurden auch mit dem Programm STINT gemacht und zwar maßgeblich von Gaffney (1984). Die Unterschiede sind jedoch bei dem Programm STINT erheblich ausgeprägter, als dies hier lediglich ansatzweise zutage tritt.
Noch bemerkenswerter sind die folgenden Resultate, die mit dem Programm STINT gemacht wurden auf einer Rechenanlage Siemens 7.890-E.
Ein und dasselbe Programm, hier STINT, wurde viermal laufen gelassen und zwar zweimal mit dem Siemens-Fortran-77 Übersetzer und zwei andere Male mit dem Waterloo-Fortran-77 Übersetzer. Beide Übersetzer laufen unter dem gleichen Betriebssystem und natürlich auf dem gleichen Rechner. Dennoch zeigte sich, daß die Ergebnisse nicht immer die gleichen waren. Es ergab sich vielmehr für das Problem E4:
| $\varepsilon=10^{-7}$ | exakte J Siemens |
num. J Siemens |
exakte J WATFOR |
num. J WATFOR |
|---|---|---|---|---|
| #Schritte | 587 | 592 | 585 | 610 |
| #Fktn | 810 | 920 | 848 | 986 |
| #Jac | 16 | 17 | 16 | 19 |
| #LU | 60 | 54 | 63 | 58 |
Diese Werte wurden ermittelt bei fast identischen globalen Fehlern und und zwar in allen Komponenten. Die Unterschiede sind für sich betrachtet nicht bemerkenswert groß, jedoch die Tatsache, daß alle Ergebnisse auf ein und derselben Rechenanlage mit ein und demselben Programm ermittelt wurden, mit lediglich verschiedenen Übersetzern, verdeutlicht den hohen Grad an Unstetigkeit der genannten Segmente und den nicht vernachlässigbaren Einfluß von Geringfügigkeiten. Zugleich wird sichtbar, daß ein Vergleich verschiedener Programme und eine sich daran anschliessende Bewertung, diesen Gegebenheiten Rechnung tragen muß. Die Frage nach einem “besten” Programm wird daher weiter noch verkompliziert, da eine Garantie der Reproduzierbarkeit von Ergebnissen nicht mehr gewährleistet werden kann. Ähnliche Schlüsse und Ergebnisse berichten auch Parlett/Wang (1975) bei Programmen für Berechnungen in der linearen Algebra.
Bei der Lösung elliptischer, partieller Differentialgleichungen mit der Programmsammlung bzw. Methodenbank ELLPACK, welches dem Benutzer zahlreiche Diskretisierungsarten und verschiedene Lösungsverfahren für lineare und u.U. nichtlineare Gleichungssysteme bereitstellt, siehe Rice/Boisvert (1985), tritt ebenfalls das Problem auf, wie sich diese Sammlung bei einem Wechsel der Rechenanlage und bei einem Wechsel der Übersetzer verhält. Beobachtet werden konnte in Rice (1983) eine Streuung von bis zu 50%.
Bibliographisch:
- Beresford Neill Parlett
- Y Wang
- John Richard Rice, lived 1934–2024, complete bibliography
- Rice, John Richard: “Machine and Compiler Effects on the Performance of Elliptic PDE Software”, 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.97–101
- Ronald F. Boisvert
2. STINT vers. GEAR und EPISODE
However, the step size chosen by STINT is almost constant for long intervals of the integration period. The main reason for this is that the STINT code has failed to recognize that ^{methods with inadequate stability regions were used} over the interval $0.01\le t\le 25.0$. This is a consequence of the fact that the step size and order selection in STINT is based on error control and not on stability analysis. However, since none of the codes monitors stability, it is a little surprising that STINT exhibits this difficulty.
Hier soll anhand der Test- und Meßresultate von STIFF DETEST gezeigt werden, daß ein Programm basierend auf den zyklischen Formeln von J.M. Tendler durchaus erfolgreich mit anderen bewährten Lösern konkurrieren kann und in speziellen Fällen sogar günstigere Ergebnisse liefern kann. Über Schwächen des Programmes STINT wird u.a. auch noch berichtet. Da das Programm TENDLER bei bisher allen getesteten Differentialgleichungen effizienter arbeitet, unter günstigen Ümständen sogar stellenweise bis zu dreimal weniger Funktionsauswertungen benötigt, als das Programm STINT, werden ähnliche Ergebnisse auch für das Programm TENDLER erwartet. Zusätzlich liefert das Programm TENDLER mehr Informationen zurück, als dies das Programm STINT tut. Die erweiterten Informationen beziehen sich auf die generierten Statistiken und die rückwärtsgenommenen Differenzen. Die weiter unten stehenden Tabellen wurden dem Aufsatz von Tendler/Bickart/Picel (1978) entnommen. Dort findet man vielfach ausführlicher, weitere Tabellen, die zahlreichen Aufschluß über die Wirkungsweise des Programmes STINT liefern.
Biographisch:
- Joel Marvin Tendler (1943–2022)
- Zdenek Picel
- Theodora A. Bickart
1. Während der Entwicklung des Programmes STINT wurden eine Reihe von Test-Problemen entworfen und das Programm STINT wurde derart abgestimmt, daß es bei dieser Klasse von Problemen sehr effizientes Verhalten zeitigte. Obwohl diese Probleme nur einen kleinen Querschnitt durch alle überhaupt auftretenden Differentialgleichungen darstellen, wurde jedoch festgestellt, daß das Programm STINT durchaus leistungsfähig ist und die Löser GEAR und EPISODE in diesem Falle deutlich überbietet. U.a. wurde das Programm STINT anhand des folgenden Beispieles getestet:
mit der exakten Lösung
und den Eigenwerten von $A$ zu $\lambda_{1/2} = \pm i$ und $\lambda_{3/4} = -10\pm 90i$. Für Integrationen über längere Intervalle wäre also ein Widlund-Winkel von $\alpha>83.6^\circ$ nötig, also maximal die BDF3 oder der Zyklus dritter Ordnung. Das Steifheitsmaß beträgt hier ungefähr $90$, wobei das Steifheitsmaß festgelegt ist zu
Integriert man von $t_0=0$ bis $t=15$ und einer Genauigkeitsanforderung von $\varepsilon=10^{-5}$ des lokalen Fehlers, so konnten die nachstehenden Leistungsdaten gemessen werden:
| Programm | Schritte | #Fktn | #Jac | #LU |
|---|---|---|---|---|
| STINT | 947 | 1016 | 6 | 23 |
| GEAR | 1590 | 1856 | 97 | 97 |
| EPISODE | 1622 | 2281 | 111 | 111 |
Reduziert man die Genauigkeit auf $\varepsilon=10^{-3}$ und begrenzt man, wegen der großen Steifheit die höchstzulässige Ordnung auf 4, so erhält man das folgende Bild für das Programm STINT:
| Schritte | #Fktn | #Jac | #LU |
|---|---|---|---|
| 552 | 622 | 5 | 14 |
mit einem maximalen absoluten und maximalen relativen Fehler von $0.0144$. Alle Rechnungen wurden auf einer IBM 370/165 und einem nicht näher genannten Übersetzer durchgeführt.
2. In den nun folgenden Tabellen gibt die Spalte Genauigkeit an, mit welcher Genauigkeitsanforderung $\varepsilon$ das Programm gestartet wurde. Die Spalte `Overhead' kennzeichnet dabei die Differenz der Gesamtrechenzeit abzüglich der Zeit für
- die Auswertung der Funktion,
- die Berechnung der Jacobimatrix und
- die Bestimmung der $LU$-Zerlegung.
Hierbei beachte man, daß zahlreiche Gleichungen in STDTST sehr einfach auszuwertende Funktionen und Jacobimatrizen besitzen und Rücksubstitutionen, Normberechnungen, etc., dem Overhead zugeschlagen werden.
Die Testresultate wurden unabhängig voneinander einmal für relative Genauigkeit und beim anderen Male für die absolute Genauigkeit ermittelt. Das Programm STINT ist konzeptionell besonders für relative Genauigkeit angelegt und auch so programmiert. Vergleicht man die unten folgenden Tabellen, so stellt man fest, daß das Programm STINT für die Anzahl der Jacobimatrixauswertungen besonders günstige Werte aufweisen kann. Dies liegt natürlich an der besonderen Behandlung der Jacobimatrix, nämlich deren Speicherung.
3. Absolute Fehler. Für das Programm STINT erhielt man die folgende Tabelle. Einmal war hierbei der globale Fehler am Endpunkt größer als 1000. Zwangsabbrüche aufgrund von Konvergenzversagen oder Fehlertestversagen traten nicht auf.
| $\varepsilon$ | CPU-Zeit | Overhead | #Fktn | #Jac | #LU | Zyklen |
|---|---|---|---|---|---|---|
| $10^{-2}$ | 4.236 | 3.613 | 5585 | 334 | 708 | 1027 |
| $10^{-4}$ | 8.613 | 7.765 | 10325 | 375 | 945 | 1969 |
| $10^{-6}$ | 16.426 | 15.040 | 18670 | 483 | 1516 | 3410 |
| $10^{-8}$ | 32.757 | 30.448 | 35181 | 626 | 2423 | 6723 |
| $\Sigma$ | 62.033 | 56.886 | 69761 | 1818 | 5592 | 13129 |
Für das Programm GEAR ergab sich das folgende. Hier war einmal der globale Fehler am Endpunkt größer als 1000. Abbruch der Integration aufgrund von Konvergenzschwierigkeiten bzw. aufgrund von Fehlertestschwierigkeiten trat nicht auf.
| $\varepsilon$ | CPU-Zeit | Overhead | #Fktn | #Jac | #LU | Schritte |
|---|---|---|---|---|---|---|
| $10^{-2}$ | 1.997 | 1.925 | 3618 | 454 | 454 | 2343 |
| $10^{-4}$ | 5.995 | 5.830 | 8946 | 751 | 751 | 6768 |
| $10^{-6}$ | 9.031 | 8.783 | 13430 | 1006 | 1006 | 10581 |
| $10^{-8}$ | 17.671 | 17.271 | 21718 | 1398 | 1398 | 18409 |
| $\Sigma$ | 34.694 | 33.809 | 47712 | 3609 | 3609 | 38101 |
Für das Programm EPISODE ergab sich nun die folgende Tabelle. Einmal versagte die Fehlerkontrolle und führte zu einem Abbruch der gesamten Integration, und einmal versagte das Programm gleich beim Start weg vollständig.
| $\varepsilon$ | CPU-Zeit | Overhead | #Fktn | #Jac | #LU | Schritte |
|---|---|---|---|---|---|---|
| $10^{-2}$ | 5.684 | 5.570 | 5884 | 821 | 821 | 4075 |
| $10^{-4}$ | 8.390 | 8.204 | 9624 | 1182 | 1182 | 6585 |
| $10^{-6}$ | 12.389 | 12.093 | 15350 | 1567 | 1567 | 9834 |
| $10^{-8}$ | 22.950 | 22.488 | 24764 | 2098 | 2098 | 17445 |
| $\Sigma$ | 49.413 | 48.356 | 41849 | 5668 | 5668 | 37939 |
Es folgen nun die Resultate für die relativen Fehler. Hier wurden also die gleichen Testdifferentialgleichungen wie oben, also A1 bis F5 gelöst und die dabei sich ergebenden Werte tabelliert.
4. Relative Fehler. Für das Programm STINT ergaben sich die folgenden Ergebnisse, die in der untenstehenden Tabelle aufgegliedert sind. Sechsmal war der globale Fehler am Endpunkt größer als 1000, und einmal brach die Integration ab, weil der Korrektor nicht konvergierte.
| $\varepsilon$ | CPU-Zeit | Overhead | #Fktn | #Jac | #LU | Zyklen |
|---|---|---|---|---|---|---|
| $10^{-2}$ | 2.201 | 1.713 | 3117 | 283 | 507 | 598 |
| $10^{-4}$ | 4.980 | 4.273 | 6363 | 356 | 711 | 1251 |
| $10^{-6}$ | 11.411 | 10.315 | 13365 | 412 | 1019 | 2727 |
| $10^{-8}$ | 21.108 | 19.269 | 23789 | 498 | 1616 | 4594 |
| $\Sigma$ | 39.701 | 35.570 | 46634 | 1549 | 3853 | 9170 |
Das Programm GEAR hatte fünfmal am Ende des Integrationsintervalls einen Fehler produziert, der größer als 1000 war. Auch hier gab es keine Korrektor-Konvergenzprobleme, die den Löser zum Halten zwangen.
| $\varepsilon$ | CPU-Zeit | Overhead | #Fktn | #Jac | #LU | Schritte |
|---|---|---|---|---|---|---|
| $10^{-2}$ | 1.359 | 1.311 | 2375 | 361 | 361 | 1584 |
| $10^{-4}$ | 4.842 | 4.713 | 6975 | 652 | 652 | 5346 |
| $10^{-6}$ | 6.580 | 6.390 | 9992 | 822 | 822 | 8057 |
| $10^{-8}$ | 11.766 | 11.470 | 15721 | 1089 | 1089 | 13155 |
| $\Sigma$ | 24.546 | 23.884 | 35063 | 2924 | 2924 | 28142 |
Das Programm EPISODE startete einmal überhaupt nicht, einmal konvergierte die Newton-Kantorovich Iteration für den Korrektor nicht und resultierte somit in einem Programmabbruch, und schließlich fünfmal war der Fehler zum Schluß des Intervalls höher als 1000.
| $\varepsilon$ | CPU-Zeit | Overhead | #Fktn | #Jac | #LU | Schritte |
|---|---|---|---|---|---|---|
| $10^{-2}$ | 4.568 | 4.473 | 4832 | 744 | 744 | 3740 |
| $10^{-4}$ | 5.930 | 5.782 | 7534 | 1011 | 1011 | 5177 |
| $10^{-6}$ | 8.746 | 8.515 | 11685 | 1355 | 1355 | 7981 |
| $10^{-8}$ | 14.940 | 14.592 | 17798 | 1558 | 1558 | 12756 |
| $\Sigma$ | 34.184 | 33.361 | 41849 | 4668 | 4668 | 29654 |
5. Den oben abgedruckten Zahlen, aber auch den ausführlicheren Tabellen, die in dem Aufsatz von Tendler/Bickart/Picel (1978) zu finden sind, entnimmt man die folgenden Interpretationen der Werte. Das Verhältnis von Zyklusanzahl zu Jacobimatrixauswertungen beträgt für das Programm STINT $2:1$ bis $9:1$. Bei höheren Genauigkeiten ist das Verhältnis höher. Dies erklärt sich leicht daraus, daß bei höheren Genauigkeitsanforderungen höhere Ordnungen benutzt werden, die ihrerseits in dem Programm STINT dazu führen, daß die Schrittweite weniger häufig gewechselt wird. Andererseits sind bei scharfen Genauigkeitsanforderungen die Schrittweiten generell kleiner und die Änderungen ebenfalls. Dies berifft dann die Iterationsmatrix $W=I-h\gamma J$, welche sich dann weniger häufig ändert und weniger häufig einer Refaktorisierung bedarf. Damit erhält man das Ergebnis, daß der Hindmarsh-Test in dem Programm STINT dazu führt, daß in grobem Zahlen ausgedrückt, alle 2 bis 9 Zyklen die Jacobimatrix neu berechnet wird.
Das Verhältnis von $LU$-Zerlegungen zu Jacobimatrixauswertungen beträgt ungefähr $2:1$ bis $4:1$. Wiederum bei höheren Genauigkeiten vergrößert sich dieses Verhältnis. Mit anderen Worten: Doppelt so häufig bis viermal so häufig wird die Iterationsmatrix $W=I-h\gamma J$ faktorisiert, bevor die Jacobimatrix $J$ neu berechnet wird. Das Verhältnis von Zyklen zu $LU$-Zerlegungen beträgt ungefähr $1:1$ bis $3:1$. In Worten: Entweder jeden neuen Zyklus, bis zu jedem dritten Zyklus wird die Iterationsmatrix $W$ faktorisiert. Dies verdeutlicht, daß der Hindmarsh-Test vergleichsweise früh anzieht.
In den Programmen GEAR und EPISODE findet man mit zum erstenmal den Hindmarsh-Test. Das Programm DIFSUB, siehe Gear (1971), kennt keinen Hindmarsh-Test und faktorisiert die Iterationsmatrix genau dann, wenn die Ordnung gewechselt wird oder Konvergenzschwierigkeiten auftreten.
A.C. Hindmarsh verwendete ursprünglich einen 30%-Test. In dem Programm STINT wird ein 40%/80%-Test verwendet und schließlich das Programm TENDLER benutzt einen erweiterten 40%/100%-Test. Die Erweiterungen liegen hierbei, daß weitere Umstände, wie beispielsweise angesammelte Minuspunkte, die Faktorisierung der Iterationsmatrix $W$ beeinflussen.
Bibliographisch:
- Alan C. Hindmarsh (*1942)
- Charles William Gear (1935–2022)
- C.W. Gear: “The automatic integration of ordinary differential equations” Communications of the ACM, Volume 14, Issue 3, pp. 176 - 179, https://doi.org/10.1145/362566.362571
3. STINT versus LSODE
In particular, we note the sensitivity of STINT and BLEND to changes in computer precision. From time to time we shall use LSODE as a yardstick by which to evaluate the above codes.
Nach Erscheinen des Aufsatzes von Tendler/Bickart/Picel (1978) sind eine Reihe von weiteren Programmen und Verfahren entwickelt worden, und zudem ist das Programm GEAR weiter verfeinert und verbessert worden zu dem Programm LSODE. Da das Programm TENDLER dieselben Formeln benutzt wie das Programm STINT, ist es aufschlußreich die Ergebnisse und Daten, die durch Vergleich des Programmes STINT mit anderen Lösern ermittelt wurden, zu sammeln, auszuwerten und schließlich mit in die Entwicklung des Programmes TENDLER hineinfliessen zu lassen. Getestet und beurteilt wurde das Programm STINT erneut von Gaffney (1984). Leitgedanke dieser Messungen war nicht allein schwerpunktmässig ein Vergleich. Jedoch liefern die gewonnenen Daten für eine Gegenüberstellung ebenfalls im begrenzten Umfange Anhaltspunkte.
P.W. Gaffney testete die Programme:
- STRIDE von Burrage, implizite RK-Verfahren
- BLEND von Skeel und Kong (blended linear multistep-methods)
- STINT
- DIRK von Alexander, implizite RK-Verfahren
- LSODE von Hindmarsh
Gaffney untersuchte diese Programme nicht wie Enright/Hull (1975) an 19 verschiedenen Testdifferentialgleichungen oder Fox (1972) an 6 verschiedenen Gleichungen, sondern nur an einer einzigen, die jedoch durch Parameterveränderung stark variiert werde konnte. Hiermit ergab sich eine einparametrige Schar von Testgleichungen. P.W. Gaffney betrachtete das folgende Anfangswertproblem, analog den Problemen B2 bis B5 $$ \dot y = \pmatrix{ -10& \omega & 0 & 0 & 0 & 0\cr -\omega&-10 & 0 & 0 & 0 & 0\cr 0 & 0 & -4 & 0 & 0 & 0\cr 0 & 0 & 0 & -1 & 0 & 0\cr 0 & 0 & 0 & 0 & -1/2 & 0\cr 0 & 0 & 0 & 0 & 0 & -1/10\cr } y, \quad y(0) = \pmatrix{1\cr 1\cr 1\cr 1\cr 1\cr 1}, $$ wobei $\omega$ nicht nur zwischen $3$ und $100$ lief, sondern den Bereich zwischen $1$ und $1000$ umfaßte. Die beiden nicht-reellen Eigenwerte $\lambda_{1/2}=-10\pm i\omega$ der Matrix überprüfen, ob die verwendeten Programme auch noch hinreichend effizient arbeiten, wenn die Lösungen der Differentialgleichung gedämpfte, aber stark oszillierende Schwingungen beschreiben. Dies erfordert von den verwendeten Formeln u.a. entweder große Widlund-Winkel, oder die Ordnungssteuerungen der Programme müssen schnell genug in der Lage sein die gewählten Ordnungen zu verringern. Man beachte, daß auch bei den 19 Testdifferentialgleichungen von Enright/Hull (1975), einige verschiedene Gleichungen nur durch Variation von Parametern entstehen.
Nicht allein entscheidend war der Vergleich der Programme, sondern die Auswirkung der Wahl der Rechenanlage und die Wahl der Maschinengenauigkeit auf die Ergebnisse, die dann die obigen Programme für eine partielle Differentialgleichung lieferten. Diese zeitabhängige partielle Differentialgleichung wurde in eine gewöhnliche Differentialgleichung “umgewandelt”, die nun ähnliche Eigenschaften besitzt, wie die oben angegebene Testgleichung, insbesondere ist sie steif und erfordert sehr große Widlund-Winkel. Die semidiskretisierte partielle Differentialgleichung ähnelt von ihrer Schwierigkeit her den Problemen B2 bis B5, und dies ist gerade der Grund, warum die obige Testgleichung gewählt wurde.
Gaffney erstellte nun u.a. die folgenden Vergleichswerte, die hier leicht gekürzt wiedergegeben werden. In der unten gedruckten Tabelle sind zur Gegenüberstellung u.a. die beiden Programme STINT und LSODE angegeben. Zu erwähnen ist, daß das Programm STINT bei der Beschränkung auf die Ordnung $p=2$ genau die gleichen Formeln benutzt wie das Programm LSODE, nämlich die BDF1 und BDF2. Hier zeigt sich deutlich der Einfluß der verschieden gewählten Darstellungsformen der Lösung, die verschiedenen zur Anwendung gelangenden Strategien. Das Programm STINT speichert zusätzlich zu der Iterationsmatrix $W=I-h\gamma J$ auch die Jacobimatrix $J$, während hingegen das Programm LSODE nur alleine die Iterationsmatrix $W$ speichert. Das Programm LSODE besitzt daher nicht die Möglichkeit gezielt nur die Iterationsmatrix zu refaktorisieren, ohne dabei gleichzeitig die Jacobimatrix $J$ mitzuberechnen. Die Refaktorisierung der Iterationsmatrix $W$ hat aber häufig ihre alleinige Ursache in der dominanten Veränderung der Werte $h\gamma$, nicht jedoch in einer Veränderung der Matrix $J$. Dennoch muß dann aus Speicherplatzgründen die Jacobimatrix $J$ mitberechnet werden. Die Weiterentwicklung des Programmes LSODE zu dem Programm VODE trägt diesen Ergebnissen Rechnung, als daß nun auch in VODE eine Option vorhanden ist, die es erlaubt zu bestimmen, ob die Jacobimatrix $J$ getrennt gespeichert werden soll oder nicht. Die Vergleiche, die Hindmarsh/Brown/Byrne (1989) veröffentlichen, zeigen erneut, daß eine separate Speicherung von $J$ sehr nützlich ist.
| $\omega$ | Routine | $p$ | $\overline p$ | #Fktn | Schritte | #Jac | $\varepsilon_{\max}$ | $\varepsilon_{\rm end}$ | CPU |
|---|---|---|---|---|---|---|---|---|---|
| 500 | STINT | 3 | 2.97 | 2057 | 1923 | 9 | 0.048 | 0.59E-3 | 0.2830 |
| " | LSODE | 2 | 2.00 | 2869 | 2053 | 116 | 0.056 | 0.23E-3 | 0.4597 |
| " | BLEND | 4 | 3.98 | 2834 | 1685 | 30 | 0.0014 | 0.11E-5 | 0.3203 |
| " | STRIDE | 15 | 5.00 | 3473 | 264 | 31 | 2.3E-5 | 0.92E-7 | 1.0872 |
| " | DIRK | 5 | 5.00 | 3931 | 256 | 28 | 5.4E-6 | 0.88E-8 | 0.2629 |
| 800 | STINT | 2 | 2.00 | 3265 | 3198 | 5 | 0.081 | 0.14E-3 | 0.4244 |
| " | LSODE | 2 | 2.00 | 4322 | 3120 | 171 | 0.089 | 0.23E-3 | 0.6969 |
| " | BLEND | 4 | 3.98 | 3281 | 2590 | 26 | 0.0065 | 0.19E-5 | 0.4103 |
| " | STRIDE | 15 | 4.90 | 5103 | 411 | 36 | 3.9E-5 | 0.11E-6 | 1.6158 |
| " | DIRK | 5 | 5.00 | 5987 | 393 | 33 | 8.4E-6 | 0.14E-7 | 0.4000 |
| 1000 | STINT | 2 | 2.00 | 3956 | 3882 | 6 | 0.101 | 0.21E-3 | 0.5140 |
| " | LSODE | 2 | 2.00 | 5237 | 3790 | 206 | 0.110 | 0.23E-3 | 0.8458 |
| " | BLEND | 4 | 3.97 | 4144 | 3327 | 26 | 0.084 | 0.76E-6 | 0.5231 |
| " | STRIDE | 15 | 5.10 | 6443 | 494 | 41 | 4.9E-5 | 0.22E-8 | 2.0431 |
| " | DIRK | 5 | 5.00 | 7337 | 483 | 37 | 1.0E-4 | 0.23E-7 | 0.4887 |
Hierbei ist $p$ die maximal zulässige Ordnung, und $\overline p$ bezeichnet die mittlere Ordnung, die vom Programm verwendet wurde. $\varepsilon_{\max}$ ist der maximale Fehler, der überhaupt gemacht wurde, und $\varepsilon_{\rm end}$ ist der globale Fehler am Integrationsende, in diesem Falle bei $t=64$. Die angegebenen Werte wurden auf einer Rechenanlage Cray-1 ermittelt.
Vergleicht man nun die obigen Resultate des Programmes STINT mit den anderen Programmen, so ergibt sich das folgende Bild für STINT.
- Von der CPU-Zeit her ist das Programm STINT zuerst zweitschnellstes Programm, danach dritt- und anschließend wieder zweit-schnell-stes. Das Programm DIRK ist hier stets am effizientesten. Hier gilt aber zu beachten, daß die Funktionsauswertungen und die Jacobimatrixauswertungen ausgesprochen einfach sind.
- Bei den Jacobimatrixauswertungen bietet das Programm STINT deutliche Vorteile ggü. allen anderen Programmen. Dieses geht aber zu Lasten des Speichers, was bei großen Problemen zu berücksichtigen ist. Allerdings sind bei großdimensionalen Problemen häufig Bandstrukturen vorhanden, sodaß sich der zusätzliche Speicherbedarf nicht im entsprechenden Maße auswirkt.
- Bei den Funktionsauswertungen zeigt das Programm STINT generell deutlich Vorteile, verglichen mit sämtlichen anderen gegenübergestellten Unterprogrammen. Das CPU-Zeit-günstigste Programm DIRK schneidet bei dieser Betrachtung hier jedesmal am wenigsten günstig ab. Dies unterstreicht nocheinmal, daß bei dem gewählten Testbeispiel die Funktion sehr einfach auszuwerten ist, es sind nämlich nur 8 Gleitkommamultiplikationen und nur 2 Gleitkommaadditionen nötig. Die Jacobimatrix ist sogar konstant. Bei den in der Einleitung erwähnten Homotopieverfahren können die Funktionsauswertungen u.U. relativ aufwendig sein, sodaß das Runge-Kutta-Verfahren DIRK in diesen Fällen nicht mehr CPU-Zeit-mässig brauchbare Werte liefern könnte.
Man beachte die sehr verschiedenen globalen Fehler, die von den untersuchten Programmen geliefert wurden. Dies relativiert die obigen Gegenüberstellungen, jedoch sind die beiden Programme STINT und LSODE gut miteinander vergleichbar. Besondere Aufmerksamkeit gebe man den sehr wenigen Schritten von den beiden Programmen STRIDE und BLEND, d.h. es wurden sehr große Schrittweiten gewählt. Auffällig ist, daß das Programm LSODE deutlich am meisten Jacobimatrixauswertungen benötigte und zwar sogar mehr, als alle anderen Programme zusammen.
Ein wichtiges Anliegen der Vergleiche, die Gaffney durchführte, war der Einfluß der Rechenanlage auf die Ergebnisse der benutzen Programme. Hierbei zeigte das Programm STINT besondere Auffälligkeiten. Bei der nun folgenden Tabelle beobachte man das---abweichend vom Üblichen---% Verhalten des Programmes STINT, einmal bei Wechsel der Genauigkeit und zum anderen beim Wechsel der Rechenanlage. Es ist hier $\omega = 100$. ^{PDP10/SP} beziehungsweise ^{PDP10/DP} kennzeichnet die gewählte Genauigkeit, also single precision oder double precision auf einer Rechenanlage vom Typ PDP10.
| $\varepsilon$ | Rechenanlage | PDP10/SP | Cray-1 | CDC-7600 | PDP10/DP |
|---|---|---|---|---|---|
| $10^{-4}$ | Schritte | 2795 | 977 | 977 | 977 |
| " | #Fktn | 2946 | 1114 | 1114 | 1114 |
| " | #Jac | 8 | 8 | 8 | 8 |
| " | CPU-Zeit | 6.680 | 0.161 | 0.576 | 3.506 |
| $10^{-5}$ | Schritte | 815 | 2917 | 1186 | 2278 |
| " | #Fktn | 1069 | 3125 | 1329 | 2489 |
| " | #Jac | 17 | 11 | 12 | 12 |
| " | CPU-Zeit | 2.138 | 0.470 | 0.774 | 8.670 |
| $10^{-6}$ | Schritte | 703 | 3690 | 3690 | 3690 |
| " | #Fktn | 976 | 3867 | 3867 | 3867 |
| " | #Jac | 12 | 12 | 12 | 12 |
| " | CPU-Zeit | 2.047 | 0.572 | 2.060 | 13.495 |
Zum Vergleich seien hier die entsprechenden Werte für das Programm BLEND angegeben. Die Werte für die anderen Programme findet man in dem schon erwähnten Aufsatz von Gaffney (1984). Es zeigt sich dabei, daß lediglich das Programm STINT besonders auffällige Anzeichen aufweist. Wie bei der obigen Tabelle ist hier $\omega = 100$.
| $\varepsilon$ | Rechenanlage | PDP10/SP | Cray-1 | CDC-7600 | PDP10/DP |
|---|---|---|---|---|---|
| $10^{-4}$ | Schritte | 3345 | 366 | 366 | 366 |
| " | #Fktn | 740 | 8834 | 8834 | 8834 |
| " | #Jac | 32 | 28 | 28 | 28 |
| " | CPU-Zeit | 1.513 | 0.100 | 0.493 | 2.757 |
| $10^{-5}$ | Schritte | 389 | 389 | 389 | 389 |
| " | #Fktn | 871 | 902 | 902 | 902 |
| " | #Jac | 27 | 27 | 27 | 27 |
| " | CPU-Zeit | 1.947 | 0.106 | 0.530 | 3.046 |
| $10^{-6}$ | Schritte | 547 | 728 | 728 | 936 |
| " | #Fktn | 1122 | 1251 | 1251 | 1532 |
| " | #Jac | 26 | 32 | 32 | 34 |
| " | CPU-Zeit | 2.734 | 0.182 | 0.976 | 7.157 |
Biographisch:
4. Rationalisierte Schrittweiten und dessen Einfluß
In einfacher Genauigkeit sind die IBM 360 Systeme Maschinen mit Basis 16 und abgehackter Arithmetik, die etwa sieben Dezimalstellen mitführen. Die Wortlängen und Rundungsfehlereigenschaften sind so schlecht, daß einfache Genauigkeit vollständig nutzlos für die ernsthafte Lösung von Differentialgleichungen ist (im Unterschied zu Lösungen, die als Beispiel im Unterricht vorgeführt werden). Andererseits ist die doppelte Genauigkeit von etwa 16 Dezimalstellen hardwaremässig eingebaut und sie kostet nur etwa 25% mehr als einfache Genauigkeit. Aus diesen Gründen sollte man gewöhnlich auf dieser Serie von Maschinen vollständig mit doppelter Genauigkeit rechnen. (Manche wundern sich, wie billig doppelte Genauigkeit ist, andere, warum eine so schlechte einfache Genauigkeit vorgesehen wurde.)
Lawrence F. Shampine, M.K. Gordon, "Computer Solutions" (1984)
Filippi/Kraska (1973) machten auf den folgenden Sachverhalt aufmerksam. Wählt man die Startschrittweite derart, daß sie im ersten Viertel des Rechnerwortes exakt darstellbar ist und ist ebenfalls der Anfangspunkt $t_0$ so gewählt, daß er exakt im Rechnerwort darstellbar ist und wird nun weiterhin die Schrittweite lediglich halbiert oder verdoppelt, so ergibt sich eine sehr wünschenswerte Genauigkeitsanhebung. Die Einschränkung auf ein Viertel der Wortlänge dient nur zur Illustriation und um genügende Reserven bei der Schrittweitenverdopplung oder Halbierung zu haben.
Die nachstehende Tabelle enthält die Werte für die Differentialgleichung $$ \dot y = \pmatrix{0&15\cr -0.6&0\cr}y, \qquad y(0)=\pmatrix{0\cr 1\cr}, \qquad t\in[0,50], $$ mit der exakten Lösung $$ y(t) = \pmatrix{5\sin3t\cr \cos3t\cr}, \qquad t\in[0,50]. $$ Diese Testdifferentialgleichung wurde einmal mit dual rationalisierter Schrittweite und einmal ohne diesen Effekt integriert und zwar mit von S. Filippi und E. Kraska entworfenen Hybrid-Verfahren. Da bei der obigen Differentialgleichung das Integrationsintervall vergleichsweise lang ist, tritt der erwünschte Effekt besonders deutlich zutage. Die Ergebnisse wurden auf einer Rechenanlage vom Typ CDC 6400 mit 48 bit Mantissenlänge gewonnen.
| $\varepsilon$ | $h_0$ | Fehler in $y_1(50)$ | Fehler in $y_2(50)$ |
|---|---|---|---|
| $10^{-12}$ | einfach | $7.659\cdot10^{-13}$ | $1.657\cdot10^{-13}$ |
| " | angepaßt | $1.411\cdot10^{-13}$ | $3.886\cdot10^{-14}$ |
| $10^{-13}$ | einfach | $6.371\cdot10^{-13}$ | $1.298\cdot10^{-13}$ |
| " | angepaßt | $1.184\cdot10^{-14}$ | $2.851\cdot10^{-15}$ |
| $10^{-14}$ | einfach | $6.262\cdot10^{-13}$ | $1.272\cdot10^{-13}$ |
| " | angepaßt | $9.773\cdot10^{-16}$ | $2.100\cdot10^{-16}$ |
Wollte man diese Rationalisierung der Schrittweite in das Programm TENDLER einbauen, so wären
- die Schrittweiten- und Ordnungssteuerung,
- die Fehlerkontrolle und
- der Konvergenztest davon betroffen.
Man sollte vermerken, daß dieser Effekt mehr für sehr hohe Genauigkeitsanforderungen konzipiert ist. Wichtig ist, daß die Startzeit $t_0$ eine günstige Interndarstellung besitzt. Ist dies nicht der Fall, so muß die erste Schrittweite so bemessen werden, daß man durch $t_0+h_0$ auf eine günstige Interndarstellung kommt. Dies zu erreichen ist einfach, wenn man Gleitkommazahlen maskieren könnte. Weiß man um die Interndarstellung der Gleitkommazahlen im Hauptspeicher, im Cache und in den Registern, so genügt u.U. ein einfaches §union§.
Biographisch.