, 25 min read
TENDLER: 1. Grobaufbau und prinzipielle Überlegungen
- 1. Die Sprache C
- 2. Objektorientierte Programmierung
- 3. Evolution von Differentialgleichungslösern
- 4. Statische und dynamische Organisation des Programmes TENDLER
- 5. Programmrealisierungen von Mehrschrittformeln
- 6. Speicherplatzaufwand bei Mehrschrittformeln
In the late 1960's we witnessed a “software crisis”, which many people thought was paradoxical because programming was supposed to be so easy. As a result, people are now beginning to renounce every feature of programming that has become guilty by association with difficulties; not only goto statements are being questioned, but we also hear complaints about floating-point calculations, global variables, semaphores, pointer variables, even assignment statements. Soon we might be restricted to only a dozen or so programs that are sufficiently simple as to be allowable; then we will be almost certain that these programs cannot lead us into trouble, but of course we won't be able to solve many problems. $\ldots$ and I do consider assignment statements and pointer variables to be among computer science's “most valuable treasures”.
C procides the infinitely abusable goto statement, and labels to branch to. Formally, the
gotois never necessary, and in practise it is almost easy to write code without it. We have not usedgotoin this book.
Hier werden zuerst einige Eigenschaften der Sprache C in tabellarischer Form aufgelistet um das Verständnis weiterer Abschnitte zu erleichtern. Auch ein sehr kurzer Überlick über den Entstehungsverlauf moderner Löser für gewöhnliche Differentialgleichungen wird angegeben. Anschliessend wird der grundsätzliche Aufbau des Programmes TENDLER beschrieben, welches auf zyklischen Formeln basiert. Der hier gegebene Überblick ist lediglich kursorisch, da erst in den folgenden Abschnitten genauer auf die Einzelheiten der verschiedenen Segmente eingegangen wird. Schließlich werden Analogien zu Regelsystemen aufgezeigt.
Differentialgleichungslöser benötigen nicht unerheblichen Speicherplatz zur Erfüllung der ihnen aufgetragenen Aufgaben. An dieser Stelle soll genauer aufgeschlüsselt werden wie sich der Speicherplatzbedarf ändert in Abhängigkeit von der Höchstordnung und der Anzahl der Stufen innerhalb eines Zykluses. Zugleich wirken die verschiedenen Datenorganisationen zurück auf die Effizienz mit der gewisse Operationen, die notgedrungen bei einer vollautomatischen Lösung anfallen, wie z.B. Interpolation, durchgeführt werden können. Zahlreiche Differentialgleichungslöser haben sich zu umfangreichen Programmpaketen entwickelt. Diese Entwicklung geschah nicht auf einen Schlag, sondern die Programme entwickelten sich immer stückchenweise weiter. Es ist aus diesem Grunde von Interesse zu beleuchten wie sich verschiedene Versionen eines Programmes unterscheiden.
Biographpisch.
- Donald Ervin Knuth (*1938)
- Brian Wilson Kernighan (*1942)
- Dennis Mike Ritchie (1941–2011)
- Bjarne Stroustrup (*1950)
1. Die Sprache C
C was chosen as the base language for C++ because it
- is versatile, terse, and relatively low-level;
- is adequate for most system programming tasks;
- runs everywhere and on everything; and
- fits into the UNIX programming environment.
B. Stroustrup (1986)
There are several morals to this story. First don't pick the most complex language to add nonstandard arithmetic to. Secondly, don't depend on compiler writers. Thirdly, develop the package so it can be adapted to other systems as easily as possible.
C.B. Dunham (1971)
Entwicklungssprache des Programmes TENDLER ist C.
Konzeptionell wurden jedoch zahlreiche Funktionen von der Sprache C++ beeinflußt.
Für die Sprache ist das Buch “The C Programming Language” von Kernighan/Ritchie (1988) maßgeblich. Die Sprache C++ wird beschrieben in dem Buch von Stroustrup (2013), erste Auflage von 1986.
Eine vergleichsweise ausführliche Bibliographie über UNIX und die Sprache C findet man in dem Übersichtsaufsatz von A.N. Maslov (1987).
Die Liste der reservierten Schlüsselwörter in der Sprache C und C++ lautet:
asm auto break case char
class(2) const(1) continue default delete(2)
do double else enum(1) extern
float for friend(2) goto if
inline(2) int long new(2) operator(2)
overload(2) public(2) register return short
sizeof static struct switch this(2)
typedef union unsigned virtual(2) void(1)
while
(1) nicht in oldstyle C
(2) nur in C++
Die Schlüsselwörter signed und volatile sind in C++ für spätere
Versionen reserviert, sind aber augenblicklich undefiniert.
Das Schlüsselwort entry wurde von oldstyle C für spätere Versionen
reserviert, ist aber nicht benutzt.
fortran ist bei einigen Übersetzern ebenfalls reserviert.
asm ist in C++ reserviert, in C allerdings lediglich potentiell reserviert.
In C ist u.U. noalias ebenfalls als Schlüsselwort reserviert.
volatile ist in C reserviert, nicht jedoch in oldstyle C.
Gewisse nicht-graphische Zeichen können durch die folgenden Ersatzsymbolfolgen dargestellt werden.
| Zeichen | Abk. | C code |
|---|---|---|
| Zeilenvorschub | NL (LF) | \n |
| Horizontaler Tabulator | HT | \t |
| Vertikaler Tabulator | VT | \v |
| Zeichen zurück | BS | \b |
| Zeilenvorschub | CR | \r |
| Seitenvorschub | FF | \f |
| akustische Meldung | BEL | \a |
| umgekehrter Schrägstrich | \ |
\\ |
| Fragezeichen | ? |
\? |
| Hochkomma | ' | \' |
| Anführungsstriche | " |
\" |
| Bitmuster (oktal) | 0ddd |
\ddd |
| Bitmuster (sedezimal) | 0xddd |
\xddd |
Hiervon wird in sehr eingeschränktem Maße bei den Statistiken Gebrauch gemacht.
Wichtig ist ferner die Operatorenreihenfolge.
Hierbei gilt: Unäre Operatoren und Zuweisungsoperatoren sind
rechts-assoziativ, alle anderen sind links-assoziativ.
D.h. a=b=c bedeutet a=(b=c), a+b+c bedeutet (a+b)+c,
und *p++ bedeutet *(p++), keinesfalls (*p)++.
Die Vorrangregeln und die Reihenfolgen des Auswertens lauten wie nachstehend.
| Operator | Assoziativität |
|---|---|
() [] -> . |
links nach rechts |
! ++ -- - (Typ) * & sizeof |
rechts nach links |
* / % |
links nach rechts |
+ - |
links nach rechts |
<< >> |
links nach rechts |
< <= > >= |
links nach rechts |
== != |
links nach rechts |
& |
links nach rechts |
^ |
links nach rechts |
| |
links nach rechts |
| |
links nach rechts |
|| |
links nach rechts |
?: |
rechts nach links |
= += -= etc. |
rechts nach links |
, |
links nach rechts |
2. Objektorientierte Programmierung
The crew was about to go into this dump sequence when all four of our flight computer machines locked up and went “catatonic.” Had this been the real thing, the Shuttle would probably have had difficulty landing. $\ldots$ but there it was: Our machines all stopped. Our greatest fear had materialized---a generic software problem.
1. Nicht nur bei interaktiven Programmen, häufig unter Umgebungen mit graphischer Oberfläche, sondern auch bei anderen komplizierteren Software-Systemen ist das Prinzip der objektorientierten Programmierung sehr nützlich. Als Beispiel sei das Programm TENDLER angeführt. Hier gibt es eine Fülle von Anwendungsmöglichkeiten, von denen hier einige genannt werden sollen. Von Funktionszeigern wird schon jetzt sehr stark Gebrauch gemacht, nämlich bei der Auswahl der Iterationsarten, bei der Steuerung des Hauptprogrammes und natürlich bei den Aufrufen der Funktion und der Jacobimatrix.
Die objektorientierte Sichtweise für die Korrektoriteration:
- Daten: Jacobimatrix $J$, Iterationsmatrix $W$, Startiterationswert $y^0$
- Operation: Lieferung einer Näherung $y^\nu$ für $y^*$, mit $y^*=h\gamma f(y^*)+\psi$
- Nebenindikatoren: Versagen, fehlerhafte Eingaben, Schaltinformation
Ausführung in Pseudo-C++.
Drei Punkte “...” heißt in diesem Falle nicht variable Anzahl an
Argumenten, sondern deutet hier lediglich nicht weiter interessierende
Programmteile an.
class correct {
real *rw, // Iterationsmatrix $W=I-h\gamma J$
*rj; // Jacobimatrix $J$
index* intpivot;
friend ecfunc(...); // ...
// ...
public:
iwtype setconv(...,itr);
virtual itfunc(...,real tau);
// ...
correct(...); // Konstruktor
~correct(...); // Destruktor
};
Zu beachten ist, daß ein anderer Programmierstil die Probleme nicht selbstständig löst, mag der Programmierstil auch noch so neumodisch klingen. Ein neues Programmierparadigma ermöglicht jedoch u.U. eine bessere und leichtverständlichere Gliederung des Programmes. Bei der objektorientierten Programmierung geht es maßgeblich um eine Gliederung und das Verdeutlichen einer übergreifenden Ordnung. Nicht selten jedoch werden die Programme auch erheblich kürzer.
Objektorientierte Sichweise für die Fehlerkontrolle und die Schrittweiten- und Ordnungssteuerung:
- Daten: Differenzentabelle (Differenzenmatrix) $(\nabla^0 y,{\mskip 3mu}\ldots,\nabla^{p+2}y)$ in unmodifizierter Darstellung oder modifizierte Darstellung in der Form $(\phi_1,{\mskip 3mu}\ldots,\phi_{p+2})$
- Operationen: Aktualisierung der Differenzentabelle, Bereitstellen von Schätzwerten für die Ableitungen der Ordnung $p-1$, $p$ und $p+1$, Interpolation
- Nebenindikatoren: Versagen, fehlerhafte Eingabe
Problem: Berücksichtigung von Statistiken.
Zum einen könnte man sie als globale Variablen auffassen, oder aber besser
man deklariert sie als friend aller wichtigen Klassen.
Ein alternativer Realisierungsvorschlag für die Korrektoriteration wäre beispielsweise:
class correct {
friend ecfunc(...); // ...
// ...
public:
iwtype setconv(..., int itr);
// ...
virtual itfunc(...);
correct(...);
correct(...);
};
Die abgeleitete Klasse kontrolliert den Zugriff auf $J$ und $W$, beispielsweise
class mni: public correct {
real *rw, *rj;
index* intpivot;
// ...
public:
itfunc(..., real tau);
// ...
mni(...); // Konstruktor für mni
~mni(...); // Destruktor für mni
}
2. Siehe Bjarne Stroustrup (2013).
Worin besteht nun der Unterschied zwischen einer Klasse und einem Modul?
Das Unsichtbarmachen von Daten und Operationen, und eine gegliederte
Benutzerschnittstelle erhält man auch durch Module (in C und C++: files as
module technique; in Modula: MODULE).
Eine Klasse ist ein Datentyp, im abstrakten Sinne. Um eine Klasse zu benutzen muß man eine Inkarnation der Klasse schaffen, also Objekte dieser Klasse erzeugen, sogenannte Instanzen. Man kann hierbei beliebig viele Objekte dieser Klasse erzeugen.
Im Gegensatz hierzu ist ein Modul selbst ein Objekt. Um ein Modul zu benutzen ist es lediglich zu initialisieren, allerdings gibt es dann auch nur ein einziges Objekt.
Durch Klassen wird eine Vielzweckverwendbarkeit erreicht, wie man sie mit Modulen nicht erreichen kann.
Z.B. in einem Übersetzer, in dem man sich eine Klasse deklariert hat, welche eine Symboltabelle implementiert, kann man sich danach beliebig viele Symboltabellen durch Definition erzeugen: eine für globale Variablen, eine für Schlüsselwörter, eine für lokale Variablen, u.s.f. Die ringzyklische Liste für $(\nabla^iy)_i$ und für $z=hf$ ließen sich gut als abgeleitete Klasse einer generischen, zyklischen Ringliste darstellen. Dies als primitives Beispiel.
3. Ein weiteres Beispiel für Klassen, abgeleitete Klassen und virtual
functions bei interaktiven graphischen Arbeiten folgt im weiteren,
man vgl. auch Stroustrup (2013).
Es soll ein Programm geschrieben werden, welches geometrische Figuren
auf dem Bildschirm darstellt.
Die gemeinsamen Eigenschaften aller Figuren seien repräsentiert durch
eine Klasse shape, und die den einzelnen Figuren inhärenten
Eigenschaften sind in abgeleiteten Klassen abgelegt.
class shape {
point center;
color col;
// ...
public:
void move(point to) { center=to; draw(); }
point where() { return center; }
virtual void draw();
virtual void rotate(int);
// ...
};
Funktionen, welche keinen besonderen Bezug auf die speziellen Eigenschaften
der Figuren nehmen, können jetzt wie üblich deklariert werden.
Der Rest wird als virtual deklariert, d.h. die eigentliche
Definition erfolgt in einer abgeleiteten Klasse.
Beispielsweise für einen Kreis
class circle: public shape {
int radius;
public:
void draw();
void rotate(int i); { /* do nothing */ }
// ...
};
Wenn nun shape_vec[] ein Vektor von Figuren ist, so kann man wie folgt
vorgehen, um sämtliche Figuren auf dem Schirm um $45^\circ$ zu drehen.
for (int i=0; i<no_of_shapes; ++i)
shape_vec[i].rotate(45);
Die typische Operation, die ein Benutzer vollzieht, ist, daß er auf ein Objekt zeigt und fragt was es ist und dieses Objekt sich entprechend verhalten soll, ohne daß man jedes Detail genauer spezifizert und ohne daß man Typusinformationen beisteuert. Das Programm muß diese selber ausfindig machen.
4. Jerrell (1990) hat eine sehr elegante
Möglichkeit vorgeschlagen, in der Sprache C++ automatisch zu differenzieren,
ohne Rückgriff auf numerische Verfahren oder Formelmanipulationssprachen.
Zwei Eigenschaften machen diese Realisierung in C++ ggü. anderen Sprachen
besonders vorteilhaft: einmal das Operator overloading und das
class-Konzept.
Die grundsätzliche Idee ist hierbei einfach: Funktionen wie $\sin$, $\cos$,
etc., und Operatoren $\pm$, u.s.w., werden einfach überladen.
Genauer:
overload log, exp, sin, cos;
class ad { // class for automatic differentiation
double value; // $f(x)$
double gradient[NVAR]; // $\nabla f$
double hessian[NVAR][NVAR]; // $\mathop{\rm Hess}f=(\partial^2f/\partial x_i\partial y_j)_{ij}$
public
ad (double, int); // Initialisierung der Variablen
...
friend ad operator + (ad, ad);
friend ad operator + (double, ad);
friend ad operator + (ad, double);
...
double log (double);
friend ad log (ad);
...
}
Für die automatische Differentiation von $f(x_1,x_2)=x_1+x_2+\ln x_1x_2$ an der Stelle $(x_1,x_2)=(1.0,2.0)$ schriebe man einfach
#define NVAR 2
#include "ad.hpp"
int main (void) {
ad x1(1.0, 1), x2(2.0, 2);
ad z = x1 + x2 + log(x1 * x2);
...
}
Für $f(x_1,x_2)=\exp(x_1^2+x_2^2)$ schriebe man ganz entsprechend
ad z = exp (x1^2 + x2^2);
Diese automatische Erzeugung der Ableitungen hat viele Vorteile:
- sie ist i.a. wesentlich genauer als die numerische Differentiation,
- i.d.R. leichter zu programmieren und
- fehlerunanfälliger für den Anwender.
Wird $f$ korrekt ausgewertet, so wird gleichzeitig auch $f'$ und $f''$ korrekt ausgewertet.
Für das Programm TENDLER wären offensichtliche Modifikationen nötig, um zum einem nicht ständig bei jeder Funktionsauswertung eine Jacobimatrixauswertung zu veranlassen und zum anderen werden die zweiten Ableitungen bei den zyklischen Formeln von J.M. Tendler grundsätzlich nicht benötigt. Diese Modifikationen sind völlig problemlos zu realisieren, genauso einfach wie die notwendige variable Dimensionierung.
Die grundsätzlichen Ideen lassen sich auch für die Programmiersprachen C, Pascal und Fortran übernehmen, allerdings müssen gewisse Umwege in Kauf genommen werden, siehe Jerrell (1989).
Die grundsätzliche Idee des obigen Ansatzes ist schon älter und in der
Literatur mehrfach wiederentdeckt worden.
Bei obiger Klasse ad ist hessian[][] eine quadratische Matrix.
Wenn $f\colon\mathbb{R}^n\to\mathbb{R}$ zweimal stetig differenzierbar ist (es genügen
aber auch leicht schwächere Voraussetzungen), so ist $\mathop{\rm Hess} f$
eine symmetrische Matrix.
Es werden daher nur, wenn es darauf ankommt, nur $(n+1)n/2$ Speicherstellen
benötigt.
Biographisch:
3. Evolution von Differentialgleichungslösern
Although there are a few limitations of computers such as the amount of memory and processor speed, these limitations do not matter for a large collection of programs. Instead, the limitations encountered in programming are most often related to our own, very human capabilities: The principal limitations on the programs we write are often imposed by our inability to comprehend the design — not by physical capabilities of computers.
W.A. Wulf, M. Shaw, P.N. Hilfinger (1981)
Der Übergang vom Programm DIFSUB zum Programm LSODE geschah nicht ruckartig und plötzlich, sondern vollzog sich in einzelnen Zwischenschritten, wobei das Programm GEAR einen Zwischenschritt markiert. Das Programm GEAR ist allerdings nicht in irgend einer Form lediglich ein Zwischenprodukt, sondern selbst universell einsetzbar.
Der Übergang von DIFSUB zu GEAR ist maßgeblich gekennzeichnet durch die folgenden Stichworte:
- Die variable Dimensionierung im Programm DIFSUB wurde aufgegeben
zugunsten einer statischen Dimensionierung in
common-Blöcken. Eine Differentialgleichungsdimension von größer als 20 erfordert eine völlige Neuübersetzung des Programmes GEAR. Dies ist nicht so wesentlich, wie man vermuten könnte, da zum einen das Programm GEAR vergleichsweise kurz ist und zum anderen, da eine neue Funktion ohnehin eine Übersetzung dieser (einzelnen) Funktion benötigt. Im letzteren Falle können dann auch imcommonBlock einige Konstanten gesetzt werden und gleich das gesamte Programm GEAR mitübersetzt werden. - Veränderte heuristische Parameter beim Korrektor-Konvergenztest (Hindmarsh-Test) und bei der Fehlerkontrolle.
- Kein anfängliches Abspeichern in
save. Dies wurde so im Programm DIFSUB gemacht, was natürlich relativ schnell geht (Datentransport, keine Rechnung). Bei einem Konvergenztestversagen oder einem Fehlertestversagen ist dann lediglich eine Kopie diesessaveBereiches zurückzuspeichern. Im Programm GEAR werden die alten Werte in der Nordsieck-Darstellung zurückberechnet. Der Vorteil liegt hierbei darin, daß weniger Speicher benötigt wird. Der Nachteil ist darin zu suchen, daß durch zusätzliche Rechnung zum einen die Werte im Nordsieck-Vektor weiter verfälscht werden und zum anderen damit ein Versagen etwas aufwendiger wird. Doch ist es gerade die Aufgabe der Schrittweiten- und Ordnungssteuerung dies zu vermeiden. Auf Vektorrechnern hätte das Programm DIFSUB aufgrund seiner gewählten Datenorganisation ggü. dem Programm LSODE stellenweise prinzipielle, leichte Vorteile. Natürlich wird das Programm LSODE auf einer Cray erheblich schneller ablaufen als das Programm DIFSUB. Zum einen weil das Programm LSODE Compiler-Direktiven für die Cray enthält und zum anderen weil die Unterprogrammsammlung LINPACK vektorisiert auf der Cray zur Verfügung steht. (LSODE verwendet LINPACK.) - $LU$-Zerlegung anstatt Matrixinvertierung bei steifen Verfahren.
Der Übergang von GEAR zu LSODE. Gelegentlich wird von Benutzern dieser beiden Programme berichtet, daß sich die Ergebnisse beider Löser fast überhaupt nicht unterscheiden. Hierzu Skip Thompson (1986):
There is no real reason to switch solvers in programs using the older codes, however.
Innerlich unterscheiden sich diese beiden Programme jedoch enorm. U.a. ist das Programm LSODE erheblich umfangreicher als das Programm GEAR.
- Erweiterte Fehlerabprüfungen und Fehlerbehandlungen. Allgemein erweitertes “user-interface”. Numerisch wird jedoch genau das gleiche Verfahren verwendet.
- Geringfügig veränderte heuristische Parameter, wie Anfangskonvergenzrate, Strategiefaktoren bei der Fehlerkontrolle und bei der Schrittweiten- und Ordnungssteuerung.
- Benutzung von LINPACK-Routinen anstatt optimierter, “selbstgemachter” $LU$-Zerlegungen.
- Stärkere Untergliederung in Unterprogramme und Funktionsunterprogramme.
Der Übergang von dem Programm STINT zu dem Programm TENDLER ist in gewisser Hinsicht wesentlich markanter. Hauptentwicklungsziel war neben einem schaltfähigen Programm, ein leicht zu veränderndes Programm. Dies ist u.a. interessant bei dem Vergleich neuerer Versuche, Taktiken, etc., bei der Programmierung zyklischer linearer Mehrschrittformeln.
- Das ganze Programm ist in der Programmiersprache C programmiert, anstatt in Fortran.
- Starke Einkapselung von Aufgaben und entsprechende Kommentierung. Data hiding. Profiling und statistics.
- Statt einer Triade (Hypermatrix) wird eine ringzyklische Liste als zentrale Datenstruktur verwendet. Dadurch wird ein Fehlschlag der Korrektoriteration oder der Fehlerkontrolle etwas einfacher zu behandeln.
- Der Benutzer erhält wesentlich mehr Informationen zurück.
- Bequemere Benutzung durch eine Vielzahl von Treibern. Die Parameterlisten zum Aufruf von Treibern sind kurz. Irgendeine Einschränkung bei der Dimension der Differentialgleichung besteht allerdings keineswegs, wie etwa bei dem Programm GEAR. Rigorose Kontrolle von Eingabedaten.
- Portabel: Übersetzer-spezifische Details (wie
NULL,size_t, etc.) sind in einer Header-Datei zentral zusammengefaßt. - Über 78 statistische Zähler inklusive einer Auswertung der hierduch gewonnenen Daten.
- Es können nun beliebig viele Differentialgleichungen gleichzeitig gelöst werden.
- Das Programm TENDLER ist type-insensitive.
- Vollständige Vektorisierung.
- 5 verschiedene Korrektorarten in einem Programm. Hinzufügung weiterer Iterationsarten einfach: Funktionszeiger.
- 1000 Zeilen $\longrightarrow$ 10000 Zeilen.
Wie das Programm LSODE und die entsprechenden Modifikationen dieses Programmes, z.B. LSODAR, so hat auch das Programm TENDLER von den beiden Vorgängerprogrammen DIFJMT, siehe Tendler (1973), und STINT, siehe Tendler/Bickart/Picel (1976), siehe Tendler/Bickart/Picel (1978), als auch von weiteren Programmen, Ideen und Konzepte übernommen, erweitert und verfeinert.
Bibliographisch:
- Joel Marvin Tendler: “A Stiffly Stable Integration Process Using Cyclic Composite Methods”, Ph.D. Diss., Syracuse University, Syracuse, New York, 26-Feb-1973
- Tendler, Joel Marvin und Bickart, Theodore A. und Picel, Zdenek: “STINT: STiff ordinary differential equations INTegrator. Program and User Manual”, Electrical and Computer Engineering Department, Syracuse University, Syracuse, New York 13210, USA, and Department of Electrical Engineering, Delft University of Technology, Delft, The Netherlands, Technical Report TR-76-12, December 1976, ii+85 pages
- Tendler/Bickart/Picel (1978): A Stiffly Stable Integration Process Using Cyclic Composite Methods
Der Entstehungsbaum für die beiden Programme LSODE und TENDLER ist in untenstehendem Diagramm angegeben.
Die Jahresangaben in Klammern geben nicht notwendigerweise das Entstehungsjahr der Erstversion an.
4. Statische und dynamische Organisation des Programmes TENDLER
The program consists of a number of small, essentially unconnected modules for code generation, a simple lexical analyzer, a canned parser which we did not have to write, and some miscellany associated with input files and the macro facility. There are only about 15 global variables. The program is currently about 1000 lines of C. About 20 percent of these lines are “print” statements, generating the output code. To our everlasting shame, there are two GOTOs in the program.
B.W. Kernighan, L.L. Cherry (1976)
Der statische, programmmässige Aufbau ist
- Prädiktor,
- Korrektor,
- Fehlerkontrolle und
- Schrittweiten- und Ordnungssteuerung.
Jeder dieser vier Segmente wird im weiteren ausführlicher beschrieben.
Der dynamische Aufbau ist jetzt für ein vierstufiges Verfahren
- $P(EC)^{i_1}$,
- $P(EC)^{i_2}$, Fehlerkontrolle,
- $P(EC)^{i_3}$, Fehlerkontrolle,
- $P(EC)^{i_4}$, Fehlerkontrolle und Schrittweiten- und Ordnungssteuerung,
mit
In welcher Stufe eine Fehlerkontrolle durchgeführt werden soll, kann ebenfalls gezielt gesteuert werden. In diesem Falle wird deutlich, daß die erste Stufe anders behandelt wird als die restlichen. Der $P(EC)^{i_k}$-Modus ist lediglich illustrativ. Ebenso gut hätte hier PECE stehen können. Der Benutzer des Programmes TENDLER kann entscheiden, welcher Modus verwendet werden soll, oder aber kann diese Entscheidung dem Programm selber überlassen.
Der Start spielt hierbei eine Sonderrolle. Die Statistiken überwachen das gesamte Geschehen, bleiben aber selbst stets inaktiv, werden also nicht zu Entscheidungen mit hinzugezogen. Der Schaltteil erkennt Steifheit, reagiert auf Steifheit, antwortet auf Veränderungen und zwar in potentiell jedem Schritt.
Der Vergleich mit einem Regelkreis liegt sehr nahe. Dieser Gedanke wird sogar auch fruchtbar angewendet und zwar bei den Arbeiten von Deuflhard (1983) und bei Gustafsson/Lundh/Söderling (1988). Siehe auch das Buch “Scientific Computing with Ordinary Differential Equations” von Deuflhard/Bornemann (2002).
| Terminologie Regelkreis | Löser |
|---|---|
| Regelgröße/Zustandsgröße | vom Programm berechnete Näherungswerte |
| Führungsgröße | exakte Lösung |
| Stellgrößen | Schrittweite, Ordnung, Iterationsarten |
| Latenzzeit/Totzeit | 1–8 Stufen |
| Hilfsstellgrößen/Hilfsregelgrößen | Fehlerkontrolle |
Wie üblich ist die exakte Lösung i.a. nicht verfügbar. Es sei bemerkt, daß sich Schweinezyklen bzgl. des Schaltens im Programm TENDLER gut beobachten lassen.
Biographisch:
5. Programmrealisierungen von Mehrschrittformeln
Es gilt das Gesetz von der Erhaltung der Mühsal.
Prof. Dr. Joachim Treusch (1984)
Neben den verschiedenen Formeln, welche lineare Mehrschrittverfahren in zahlreiche Lager trennt, erfolgt eine weitere Untergliederung durch die Art und Weise der Programmrealisierung. Hier sei vorwiegend an die Darstellung der Lösung gedacht. Maßgeblich sind hier die Darstellung in Form klassischer rückwärtsgenommener Differenzen, modifizierten rückwärtsgenommenen Differenzen und die Nordsieck-Darstellung. Jeder dieser Darstellungen entfaltet in gewissen Situationen seine ganz besonderen Vorzüge. Jedoch hat auch jede Darstellung seine Nachteile. Anhand ganz konkreter Programme, wie TENDLER, DIFSUB und DE/STEP, sollen die benötigte Anzahl gewisser signifikanter Rechenoperationen verglichen werden.
VOP heiße im folgenden Vektoroperation.
1. Aktualisierungsaufwand pro Schritt/Stufe. Zur Pflege und Erhaltung benötigt jede Darstellung der Lösung einen gewissen Verwaltungsaufwand. Dieser berücksichtigt das Fortschreiten der Lösung in der Zeit. Im einzelnen fallen hier jeweils an:
TENDLER: $\nabla^{i+1}y_1=\nabla^iy_1-\nabla^iy_0$, für $i=1,\ldots,p$. Insgesamt also $p$ VOPs.
DE/STEP: Hier sind stets $p+2$ VOPs durchzuführen und zwar
DIFSUB: Hier sind $p-2$ VOPs durchzuführen und zwar
2. Prädiktorbildung. Zur Berechnung eines Startwertes für die Picard- oder Newton-Kantorovich Iteration wird ein explizites Mehrschrittverfahren benutzt, genannt Prädiktor.
TENDLER: Hier sind $2p+4$ VOPs nötig.
DE/STEP: Ebenfalls $2p+4$ VOPs.
DIFSUB: Es sind hier am meisten Rechenoperationen nötig, nämlich $(p+1)p/2$ VOPs. Berechnet wird
3. Korrektoriteration. Die eigentliche Iteration erfordert ebenfalls gewissen Rechenaufwand, neben der Auswertung der Funktion. Zur Vereinfachung sei zuerst hier lediglich der Fall der $P(EC)E$-Picard Iteration betrachtet.
TENDLER: 3 VOPs, für $y^{\nu+1}\gets y^\nu+\gamma(hf(y^\nu)-z^\nu)$
DE/STEP: 3 VOPs, für $x\gets y+\lambda(z-u)$.
DIFSUB: ca. 4 VOPs, für
4. Neue, veränderte Schrittweite nach erfolgreichem Schritt. Nach erfolgreichem Durchlaufen eines Schrittes oder eines gesamten Zykluses, kann die Schrittweiten- und Ordnungssteuerung vorschlagen, daß die Schrittweite und die Ordnung geändert werden sollen. Hierfür sind dann jedoch Startwerte auf einem anderen Gitter nötig, oder die Formel muß sich ändern. Erneut geschieht diese Änderung nicht kostenfrei, sondern erfordert gewisse Rechenschritte.
TENDLER:
Das Programm TENDLER berechnet die neuen benötigten Startwerte durch
Interpolation.
Es ist hierbei $p-1$ mal zu interpolieren, wobei $p$ die neue Ordnung ist.
Nachdem die neuen zurückliegenden Werte durch Interpolation bestimmt
wurden, ist die Differenzentabelle neu aufzubauen.
$y_{m\ell}$ ist schon bekannt, darum muß man nur $(p-1)$-mal interpolieren,
anstatt $p$-mal.
Hierzu wird $p-1$ mal die Funktion newtform() aufgerufen, welche
jeweils $p+1$ VOPs benötigt.
Die Neugenerierung der Differenzentabelle in der Funktion formbackdiff()
benötigt $(p+1)p/2$ VOPs.
Schließlich ist für den nächsten Prädiktorschritt noch $z=hf$ auf die
neue Schrittweite anzupassen: Die Funktion adapthz() benötigt genau eine
VOP.
Insgesamt also
DE/STEP: Sei $s$ die Anzahl der Schritte mit konstanter Schrittweite $h$. Falls $p>s$, dann $p-(s+1)$ VOPs, wegen
Im Grenzfalle hat man also
DIFSUB: Die Nordsieck-Darstellung zeigt hier deutliche Vorteile. Stets sind $p$ VOPs nötig, wegen
5. Zurückgewiesener Schritt. Wird ein Schritt von der Fehlerkontrolle, oder im Falle von steifen Differentialgleichungen von der Newton-Kantorovich-Korrektor Iteration, zurückgewiesen, so ist eine neue Schrittweite und ggf. eine neue Ordnung zu wählen. Ähnlich wie bei einem akzeptierten Schritt mit einer von der Schrittweiten- und Ordnungssteuerung vorgeschlagenen Schrittweite und Ordnung, sind jetzt ebenfalls neue Startwerte zu besorgen. In diesem Falle ist nicht nur die theoretische Darstellung von Interesse, sondern zusätzlich die Art der Speicherung, also die Speicherorganisation.
TENDLER: Gleiche Anzahl wie bei neuer, veränderter Schrittweiter nach erfolgreichem Schritt, also insgesamt ${\cal O}(3p^2/2)$ VOPs.
DE/STEP: Genau wie bei dem Programm TENDLER, also ${\cal O}(3p^2/2)$, im ungünstigsten (!) Falle.
DIFSUB: Hier sind lediglich $(p+1)$ VOPs nötig und zwar lediglich eine Umspeicherung einer Matrix.
GEAR: Obwohl das Programm GEAR die gleichen Formeln und die gleiche Darstellungsform wählt, sind hier aufgrund eines kompakteren Speicherschematas, nun in diesem Falle mehr VOPs nötig. Zu berechnen ist jetzt
Dies wird auch genau so in dem Programm LSODE und seinen Varianten vollzogen. Der Vorteil ist hierbei, daß weniger Speicher benötigt wird, als in dem Program DIFSUB. Der Nachteil ist, daß durch diese zusätzlichen Rechnungen Rundungsfehler eingeschleppt werden und es natürlich mehr Rechenarbeit kostet. Jedoch sind zurückgewiesene Schritte eher selten. Wären sie auffallend häufig, so wäre dies ein deutliches Indiz dafür, daß die Schrittweiten- und Ordnungssteuerung nicht zur vollen Zufriedenheit arbeitet.
6. Speicherplatzaufwand bei Mehrschrittformeln
Jede verschiedene Programmausführung eines linearen Mehrschrittverfahrens, oder eine zyklische Kombination eines solchen Verfahrens, benötigt zur Verwaltung der interenen Abläufe, wie Fehlerkontrolle und Korrektoriteration beispielsweise, einen gewissen Speicherraum. Dieser benötigte Speicherbedarf ist sogar bei Programmen, die die gleichen Formeln benutzen, unterschiedlich. Maßgeblich wird der Speicherplatzaufwand von der Darstellungsart der Lösung beeinflußt. Bei der numerischen Lösung steifer gewöhnlicher Differentialgleichungen überwiegt natürlich der Speicherplatzbedarf zur Speicherung der Iterationsmatrizen und ggf. noch der Jacobimatrix.
Es ist $n$ stets die Dimension der Differentialgleichung, $p$ bezeichnet immer die Ordnung, $p_\max$ ist die maximal mögliche Ordnung, $\ell$ ist die Anzahl der Stufen innerhalb eines Zykluses und $\ell_\max$ ist die Höchstanzahl der Stufen innerhalb eines Zykluses über alle Ordnungen betrachtet, die zugelassen sind. Bei einem Programm basierend auf den BDF ist die maximal mögliche Ordnung immer durch 6 nach oben beschränkt, da oberhalb der Ordnung 6, die BDF nicht mehr nullstabil sind, also erst recht nicht mehr steif-stabil sind. Die meisten Programme, die auf den BDF basieren, begrenzen sogar die Höchstordnung auf 5, da die BDF6 einen sehr kleinen Widlund-Winkel besitzt.
1. Der Speicherplatzbedarf in Abhängigkeit von $p$ und $\ell$. Durch Wahl entsprechender Eingabeargumente kann der Benutzer die Höchstordnung $p_\max$ und damit einhergehend ggf. die Höchstanzahl der Stufen $\ell_\max$ begrenzen. Der Speicherplatzbedarf hängt direkt von diesen beiden Größen ab.
TENDLER: Maximal werden $2n^2+52n$ Speicherzellen benötigt, falls $p=p_\max=7$ und $\ell=\ell_\max=4$. Hierbei ist
Im einzelnen schlüsselt sich der Speicherplatzbedarf wie folgt auf:
| Variable | Speicherbedarf |
|---|---|
rw |
$n^2$ |
rj |
$n^2$ |
yr |
$(p+2)(\ell+1)n$ |
hfr |
$(\ell+1)n$ |
intpivot |
$n$ |
ymax |
$n$ |
Dies macht zusammen dann also $2n^2+\bigl[(\ell+1)(p+3)+2\bigr]n$.
DIFSUB: Im Höchstfalle werden $n^2+23n$ Speicherzellen benötigt. Dies schlüsselt sich auf in
| Variable | Speicherbedarf |
|---|---|
pw |
$n^2$ |
error |
$n$ |
y |
$(p+1)n$ |
save |
$[(p+1)+4]n$ |
ipiv |
$n$ |
ymax |
$n$ |
Insgesamt ergibt dies $n^2+\bigl[2(p+1)+7\bigr]n$.
In dem ursprünglichen Programm DIFSUB wurde aus vermeintlichen
Effizienzgründen keine $LU$-Zerlegung durchgeführt, sondern explizit
die Inverse der Iterationsmatrix bestimmt und dann das
Matrix-Vektor-Produkt inline berechnet.
Der besseren Vergleichbarkeit wegen wurde deshalb ein Indizesvektor ipiv
hinzugefügt, welcher die Permutationsindices der $LU$-Zerlegung
enthält.
STINT: Im Höchstfalle werden $2n^2+51n$ Speicherzellen benötigt. Es ist $p_\max=7$ und $\ell_\max=4$.
| Variable | Speicherbedarf |
|---|---|
rw |
$n^2$ |
rj |
$n^2$ |
y |
$(p+1)\ell n$ |
dy |
$\ell n$ |
save |
$[(p+1)+3+2]n$ |
ipivot |
$n$ |
ymax |
$n$ |
Dies summiert sich zu $2n^2+\bigl[(\ell+1)(p+1)+\ell+7\bigr]n$.
DE/STEP: Es ist bei diesem Programm nicht vorgesehen, daß der Benutzer den Speicherplatzbedarf der Höchstordnung anpasst. Benötigt werden $22n$ Speicherzellen. Es ist $p_\max=12$. Die Aufschlüsselung des Speicherbedarfs ist
| Variable | Speicherbedarf |
|---|---|
y, yy, wt, p, yp, ypout |
$n$ |
phi |
$(p+4)n$ |
Also insgesamt $(p+10)n$.
Nur de() ruft intrp() auf.
LSODAR: Im Höchstfalle werden $4n_g+18n+n^2$ Speicherzellen benötigt.
| Variable | Speicherbedarf |
|---|---|
y, iwork |
$n$ |
jroot |
$n_g$ |
rwork |
$\displaystyle{ \max(16,n+9)n+3n_g=\cases{ (p+4)n+3n_g,&falls Adams\cr (p+4)n+n^2+3n_g,&falls BDF\cr } }$ |
Der Speicherplatz kann sich um $2n$ Speicherzellen erhöhen, wenn man für
jede Komponente der Lösung getrennte relative (rtol) und absolute (atol)
Genauigkeitsbedingungen stellt.
DDASSL: Im Höchstfalle werden $n^2+12n$ Speicherzellen benötigt. Es ist $p_\max=5$.
| Variable | Speicherbedarf |
|---|---|
y |
$n$ |
yprime |
$n$ |
rwork |
$(p+4)n+n^2$ |
iwork |
$n$ |
Dies summiert sich auf zu $n^2+(p+7)n$. Auch hier ggf. ein Mehr von $2n$ Zellen bei besonderen Genauigkeitsforderungen.
2. Der temporäre Speicherplatzbedarf. Zwischen den Aufrufen wird ein Teil des angeforderten Speichers wieder frei. Dieser temporäre Speicherraum eignet sich daher maßgeblich nur für Zwischenrechnungen, welche zwischen den Aufrufen des Lösers anfallen. Insbesondere beim Programm TENDLER liegen im vermeintlich temporären Speicher Nutzinformationen vor. Jedoch kann der Benutzer diese Nutzinformationen überschreiben, ohne die Fortsetzbarkeit der Integration zu behindern.
TENDLER: Insgesamt $13n$ Speicherzellen werden frei.
Von yr werden $(p+2)n$ und von hfr werden $\ell n$ Speicherzellen frei.
Insgesamt also $(\ell+p+2)n$ Speicherzellen.
DIFSUB: Hier werden $12n$ Zellen frei, wegen save, mit seinen
$[(p+1)+4]n$ Zellen.
STINT: Hier werden $16n$ Zellen frei, wegen dy mit $(\ell-1)n$ Zellen und
save mit $(p+6)n$ Zellen.
DE/STEP: Insgesamt $3n$ Speicherzellen, wegen y, yp, $\varphi_{15}$ und
$\varphi_{16}$ mit jeweils $n$ Zellen.
LSODAR: $n+n_g$ Zellen werden frei, wegen y mit $n$ und jroot mit
$n_g$ Zellen.