, 25 min read

TENDLER: 1. Grobaufbau und prinzipielle Überlegungen


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”.

D.E. Knuth (1974)

C procides the infinitely abusable goto statement, and labels to branch to. Formally, the goto is never necessary, and in practise it is almost easy to write code without it. We have not used goto in this book.

B.W. Kernighan, D.M. Ritchie (1978)

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.

  1. Donald Ervin Knuth (*1938)
  2. Brian Wilson Kernighan (*1942)
  3. Dennis Mike Ritchie (1941–2011)
  4. Bjarne Stroustrup (*1950)

1. Die Sprache C

C was chosen as the base language for C++ because it

  1. is versatile, terse, and relatively low-level;
  2. is adequate for most system programming tasks;
  3. runs everywhere and on everything; and
  4. 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.

Jack F. Clemons (1984)

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:

  1. Daten: Jacobimatrix $J$, Iterationsmatrix $W$, Startiterationswert $y^0$
  2. Operation: Lieferung einer Näherung $y^\nu$ für $y^*$, mit $y^*=h\gamma f(y^*)+\psi$
  3. 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:

  1. 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})$
  2. Operationen: Aktualisierung der Differenzentabelle, Bereitstellen von Schätzwerten für die Ableitungen der Ordnung $p-1$, $p$ und $p+1$, Interpolation
  3. 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:

  1. sie ist i.a. wesentlich genauer als die numerische Differentiation,
  2. i.d.R. leichter zu programmieren und
  3. 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:

  1. Max E. Jerrel, ResearchGate

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:

  1. 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 im common Block einige Konstanten gesetzt werden und gleich das gesamte Programm GEAR mitübersetzt werden.
  2. Veränderte heuristische Parameter beim Korrektor-Konvergenztest (Hindmarsh-Test) und bei der Fehlerkontrolle.
  3. 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 dieses save Bereiches 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.)
  4. $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.

  1. Erweiterte Fehlerabprüfungen und Fehlerbehandlungen. Allgemein erweitertes “user-interface”. Numerisch wird jedoch genau das gleiche Verfahren verwendet.
  2. Geringfügig veränderte heuristische Parameter, wie Anfangskonvergenzrate, Strategiefaktoren bei der Fehlerkontrolle und bei der Schrittweiten- und Ordnungssteuerung.
  3. Benutzung von LINPACK-Routinen anstatt optimierter, “selbstgemachter” $LU$-Zerlegungen.
  4. 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.

  1. Das ganze Programm ist in der Programmiersprache C programmiert, anstatt in Fortran.
  1. Starke Einkapselung von Aufgaben und entsprechende Kommentierung. Data hiding. Profiling und statistics.
  2. 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.
  3. Der Benutzer erhält wesentlich mehr Informationen zurück.
  4. 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.
  5. Portabel: Übersetzer-spezifische Details (wie NULL, size_t, etc.) sind in einer Header-Datei zentral zusammengefaßt.
  6. Über 78 statistische Zähler inklusive einer Auswertung der hierduch gewonnenen Daten.
  7. Es können nun beliebig viele Differentialgleichungen gleichzeitig gelöst werden.
  8. Das Programm TENDLER ist type-insensitive.
  9. Vollständige Vektorisierung.
  10. 5 verschiedene Korrektorarten in einem Programm. Hinzufügung weiterer Iterationsarten einfach: Funktionszeiger.
  11. 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:

  1. Joel Marvin Tendler: “A Stiffly Stable Integration Process Using Cyclic Composite Methods”, Ph.D. Diss., Syracuse University, Syracuse, New York, 26-Feb-1973
  2. 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
  3. 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.

flowchart TD DIFSUB(DIFSUB, 1971) DIFBDF(DIFBDF) DIFJMT(DIFJMT, 1973) STINT(STINT, 1976) TENDLER(TENDLER V1, 1990) TENDLER2(TENDLER V2, 2026) GEAR(GEAR, GEARB, GEARS, GEARBI, GEARBIL, GEARIB, GEARV, GEARST, 1976) LSODE(LSODE, LSODA, LSODAR, LSODES, 1980) VODE(VODE, 1989) CVODE(CVODE, 1995) IDA(IDA, 1999) %% Dependencies DIFSUB --> TendlerPhD TendlerPhD --> STINT STINT --> TENDLER TENDLER --> TENDLER2 subgraph TendlerPhD direction LR DIFBDF --> DIFJMT end DIFSUB --> GEAR GEAR --> LSODE --> VODE --> CVODE --> IDA

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

  1. Prädiktor,
  2. Korrektor,
  3. Fehlerkontrolle und
  4. Schrittweiten- und Ordnungssteuerung.

Jeder dieser vier Segmente wird im weiteren ausführlicher beschrieben.

Der dynamische Aufbau ist jetzt für ein vierstufiges Verfahren

  1. $P(EC)^{i_1}$,
  2. $P(EC)^{i_2}$, Fehlerkontrolle,
  3. $P(EC)^{i_3}$, Fehlerkontrolle,
  4. $P(EC)^{i_4}$, Fehlerkontrolle und Schrittweiten- und Ordnungssteuerung,

mit

$$ 1 \le i_k \le \mathtt{MAX}\alpha\mathtt{ITERAT}, \quad \alpha\in\{\mathtt{mni}, \mathtt{fix}\}, \qquad k=1,\ldots,4. $$

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:

  1. Peter Jochen Deuflhard (1944 – 2019)
  2. Gustafsson, Kjell
  3. Lundh, Michael
  4. Söderling, Gustaf

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

$$ \varphi_{p+1} \gets \mathtt{yp}-\varphi_{p+2},\quad \varphi_{p+2} \gets \varphi_{p+1}-\varphi_{p+2},\qquad \varphi_i \mathrel{\mathtt{+=}} \varphi_{p+1},\quad i=1,\ldots,p. $$

DIFSUB: Hier sind $p-2$ VOPs durchzuführen und zwar

$$ y_i\gets y_i+a_i{\mskip 3mu}\mathtt{error}_i, \qquad i=3,\ldots,p+1. $$

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

$$ \left(y_{p+j-j_1} \mathrel{\mathtt{+=}} y_{p+j-j_1+1},\quad j_1=j,\ldots,p+1\right) \qquad j=2,\ldots,p+1. $$

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

$$ \mathtt{save}_9 \gets y_2-h{\mskip 3mu} \mathtt{save}_{10},\quad y_1 \mathrel{\mathtt{+=}} a_1{\mskip 3mu} \mathtt{save}_9,\quad y_2\mathrel{\mathtt{-=}} \mathtt{save}_9,\quad \mathtt{error} \mathrel{\mathtt{+=}} \mathtt{save}_9 $$

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

$$ (p-1)(p+1)+{(p+1)p\over2}+1 = {3\over2}p^2+{1\over2}p. $$

DE/STEP: Sei $s$ die Anzahl der Schritte mit konstanter Schrittweite $h$. Falls $p>s$, dann $p-(s+1)$ VOPs, wegen

$$ \bigl(\varphi_{s+1},{\mskip 3mu}\ldots,{\mskip 3mu}\varphi_p\bigr)\cdot \mathop{\rm diag}\left(\beta_{s+1},\ldots,\beta_p\right). $$

Im Grenzfalle hat man also

$$ (p-2)+(p-3)+\cdots+2+1={(p-2)(p-1)\over2}={\cal O}\left(p^2\over2\right). $$

DIFSUB: Die Nordsieck-Darstellung zeigt hier deutliche Vorteile. Stets sind $p$ VOPs nötig, wegen

$$ \bigl(y_2,{\mskip 3mu}\ldots,{\mskip 3mu}y_{p+1}\bigr)\cdot \mathop{\rm diag}\left(\omega,\omega^2,\ldots,\omega^p\right) $$

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

$$ \left(y_{p+j_1-j_2} \mathrel{\mathtt{-=}}y_{p+j_1-j_2+1}, \quad j_2=j_1,\ldots,p\right) \quad j_1=1,\ldots,p. $$

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

$$ \ell=3+\left\lfloor p\over5\right\rfloor=\cases{ 3,&falls $p\le4$,\cr 4,&falls $p\ge5$,\cr }\qquad p\le p_\max=7. $$

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.