, 13 min read

TENDLER: 2. Benutzung des Programmes

Das Programm TENDLER dient zur näherungsweisen Lösung von gewöhnlichen Differentialgleichungssystemen erster Ordnung und zwar von Anfangswertproblemen. Differentialgleichungen höherer Ordnung lassen sich leicht in Gleichungen erster Ordnung überführen durch Anhebung der Dimension. I.a. werden die vom Programm TENDLER gelieferten Näherungswerte “etwas ungenauer” sein, als eingebaute Funktionen, wie z.B. $\sin$, $\cos$, etc. I.a. werden dafür auch Lösungen schneller berechnet. Z.B. die Differentialgleichung, welche $\sin$ und $\cos$ als Lösung besitzt, kann schneller gelöst werden, als die mehrfache Auswertung der eingebauten Sinus- und Kosinusfunktion. Dies ist beispielsweise interessant für graphische Anwendungen. Diese Ausführungen sollten allerdings nur als grober Hinweis angesehen werden.

Herkömmliche Löser, wie z.B. LSODE, DE/STEP, u.s.w., offerieren dem Benutzer letztlich ein einziges Unterprogramm, welches wiederholt aufgerufen wird. Dieses eine Unterprogramm wird dann durch eine Vielzahl von Parametern gesteuert. Bei dem Programm TENDLER ist dies anders. In CVODE ist dies auch anders. Definition eines Problems und Lösung desselben werden beim Programm TENDLER auch als solches unterschiedlich behandelt. Also für jeden Bereich existiert auch ein separates Unterprogramm, welches demzufolge weitaus weniger Parameter enthält. Dies hat weitreichende Folgen. Zum einen ist diese Vorgehensweise etwas effizienter (weniger Rechenzeit) und zum anderen ist diese Vorgehensweise für den Anwender auch wesentlich bequemer.

Das Programm TENDLER liefert lediglich Näherungswerte für die exakte Lösung des Anfangswertproblems. Eine Garantie, daß die Lösung bis auf Maschinengenauigkeit exakt berechnet wird, kann nicht gegeben werden. Die Programme LSODE, CVODE, DE/STEP, etc., versprechen dies übrigens auch nicht.

Sämtliche TENDLER Unterprogramme fangen mit dem Präfix tdl an.

1. Definition der Differentialgleichung

1. Zweck. Das Programm TENDLER dient in der augenblicklichen Fassung zur Lösung von

$$ \dot y = f(t,y), \qquad y(a) = y_0 \in \mathbb{R}^n, \qquad (a,t\in\mathbb{R}). $$

Dabei ist zu fordern, daß $f$ mindestens acht Mal stetig differenzierbar ist. Bei unstetigen rechten Seiten kann es zu Schwierigkeiten kommen. Für unstetige $f$, für Randwertaufgaben, für partielle Differentialgleichungen, für stochastische Differentialgleichungen ist das Programm TENDLER nicht direkt einsetzbar.

  1. Bei Unstetigkeiten zweiter Art (Sprüngen) kann man sich, wenn man weiß wo die Sprungstellen liegen, u.U. dadurch behelfen, daß man den Löser bei den Sprungstellen neu startet.
  2. Für Randwertaufgaben kann man das Programm TENDLER als Unterprogramm in einem Mehrfachschießverfahren einsetzen.
  3. Gewisse partielle Differentialgleichungen lassen sich durch Diskretisierung umformulieren zu einer gewöhnlichen Differentialgleichung, entweder als Rand- oder aber als Anfangswertproblem.
  4. Stochastische Differentialgleichung lassen sich u.U. durch wiederholtes Lösen angehen.

Die Erweiterung des Programmes TENDLER für die Behandlung von $g$-stop Problemen der Form

$$ g\colon \mathbb{R}\times\mathbb{R}^n\times\mathbb{R}^n \to \mathbb{R}^m, \qquad g(t,y,\dot y) = 0, $$

ist berücksichtigt, aber nicht eingebaut. Differentialgleichungsprobleme mit impliziter Zeitvorgabe werden gelegentlich auch $g$-stop Probleme genannt.

Das Programm TENDLER besteht aus mehr als 100 Unterprogrammen. Davon sind jedoch nur ein Teil für den Anwender i.d.R. direkt von Interesse. Dies sind maßgeblich tdldefsys(), tdldefivp() und tdlundefsys() zur An- und Abmeldung eines Anfangswertproblems. Die C-Prototypen hierzu befinden sich in "ode.h".

2. dd = tdldefsys (n, func, jac, m, g) :

  1. n ist die Dimension der Differentialgleichung.
  2. func ist der Name des Unterprogrammes zur Berechnung der Funktion $f$.
  3. jac ist der Name des Unterprogrammes zur Berechnung der Jacobimatrix $f_y$.
  4. m ist $m$ bei $g$, und g ist das Unterprogramm für die Berechnung von $g$.

tdldefsys() dient also zur Definition der Differentialgleichung ohne die Anfangswerte. Intern wird Speicherplatz via calloc() reserviert.

Der zurückgelieferte Wert dd dient im weiteren zur Identifikation des Problems, differential equation descriptor. Jetzt wird auch ersichtlich, warum es so einfach ist, beliebig viele Differentialgleichungen gleichzeitig zu lösen.

Man würde wie folgt vorgehen:

$$ \eqalign{ \mathtt{dd1} &= \mathtt{tdldefsys} (n_1, f_1, J_1, m_1, g_1);\cr \mathtt{dd2} &= \mathtt{tdldefsys} (n_2, f_2, J_2, m_2, g_2);\cr } $$

Wenn die Jacobimatrix nicht berechnet werden soll oder kann (man vgl. z.B. NC5, nonstiff C5), so übergibt man nojac, also

dd = tdldefsys(n, f, nojac, m, g);

Da z.Z. keine $g$-stop Probleme gelöst werden können, muß $m$ immer Null und $g$ auf jeden Fall nogstop sein, also

dd = tdldefsys(n, f, J, 0, nogstop);

Ist $m\ne0$ und/oder $g\ne\mathtt{nogstop}$, so erzeugt die augenblickliche Fassung des Programmes TENDLER eine Fehlermeldung. Durch eine Fehlspezifikation werden also niemals falsche Werte erzeugt.

$n$ und $f$ müssen immer angegeben werden. Inkorrekt wäre also

dd = tdldefsys(n, NULL, J, 0, nogstop);    /* wrong */

Gute Übersetzer würden dies auch schon anhand der Prototypen erkennen.

3. tdldefivp (dd, a, b, y0) :

  1. dd ist die obige Identifzierung um welche Differentialgleichung es sich handelt.
  2. a und b sind Integrationsgrenzen und y0 der Anfangswertvektor.
  3. a ist $a$.
  4. b ist an für sich unerheblich.

Es dient lediglich dazu, den Größenordnungsbereich und die Integrationsrichtung (links oder rechts) anzugeben. Es darf über b hinaus integriert werden. Umgekehrt muß b auch nicht unbedingt erreicht werden. Intern wird b zur Berechnung einer Anfangs- und Maximalschrittweite benutzt. Zum Schalten wird b nicht benutzt. Es ist möglich und auch effizient, bei einer Differentialgleichung, aber mehreren verschiedenen Anfangswerten, das Unterprogramm tdldefivp() mehrfach aufzurufen, z.B. im Rahmen eines Schießverfahrens. tdldefsys() wird hierzu nicht benötigt.

4. tdlundefsys (dd) : Dies ist das Gegenstück zu tdldefsys(). Es meldet die Kennung dd ab, insbesondere wird der belegte Speicherplatz wieder frei gegeben. Die Funktion tdlundefsys() ruft hierzu free() auf.

Das Aufrufen von tdlundefsys() kann entfallen, wenn nur eine Differentialgleichung gelöst wird und anschließend im Hauptprogramm nichts weiter berechnet wird. Dies ist eine Eigenschaft von calloc() und free() und wird nicht durch das Programm TENDLER erzwungen.

tdlundefsys() ist eigentlich nur bei großdimensionalen Problemen wirklich interessant. tdldefsys() und tdlundefsys() entsprechen in ihrer Funktionalität in etwa fopen() und fclose().

2. Lösung der Differentialgleichung

Bisher wurde lediglich das Problem definiert, also noch kein einziger Wert berechnet. Es könnte so aussehen, als hätte das Programm TENDLER bisher nichts Sinnvolles geleistet. In Wirklichkeit wurden die Daten auf Fehlerfreiheit und innere Konsistenz untersucht und intern eine Reihe von Voreinstellungen wirksam. Nun sollen tatsächlich Näherungen berechnet werden. Es gibt hier zwei Möglichkeiten.

  1. Die Punkte $t_i$, an denen Näherungen ausgegeben werden sollen, sind vom Benutzer vorgegeben. Hierzu dient tdlsolve().
  2. Die Gitterpunkte $t_i$ werden vom Programm TENDLER gewählt. Hierzu dienen tdlmeagersol() (“Magerstufe”) und tdlautosol().

1. y = tdlsolve (dd, t) :

  1. dd ist w.o. die Kennung (differential equation descriptor), und
  2. t ist ein beliebiger Wert, an dem die Lösung gesucht wird.

Zurückgeliefert wird ein Vektor $y\in\mathbb{R}^n$. Wünscht man an der Stelle $t=\pi$ einen Näherungswert, so schriebe man

y = tdlsolve(dd,3.141592653589793238462643383279502884197);

Bei zwei voneinander unabhängigen Gleichungen $y_1$ und $y_2$ schriebe man

$$ \eqalign{ y_1 &= \mathtt{tdlsolve} (\mathtt{dd1}, t);\cr y_2 &= \mathtt{tdlsolve} (\mathtt{dd2}, t);\cr } $$

Der Wert $t=a$ ist zulässig, ebenso $t>b>a$ (analog: $t<b<a$). tdlsolve() entspricht in seiner Funktionalität in etwa fgetc().

2. y = tdlmeagersol (dd, &t) :

  1. Genau wie tdlsolve(), bloß, daß in t ein vom Programm TENDLER erzeugter Gitterpunkt übergeben wird.
  2. In y befindet sich wiederum der Lösungsvektor $\in\mathbb{R}^n$. Die vom Programm TENDLER gewählten Gitterpunkte brauchen nicht äquidistant zu sein.

Die Funktion tdlmeagersol() ist nützlich, wenn man über die Lösung, außer ihrer Existenz, nichts weiter weiß. Beispielsweise ist für die Sinuskurve auf $\left[0,100\pi\right]$ ein Gitter $0,\pi,2\pi,3\pi,\ldots$ i.a. nicht so interessant. Dies weiß man häufig allerdings erst nach der Lösung einer Differentialgleichung. Der erste von tdlmeagersol() zurückgegebene Wert ist nicht $y_0$ und $t\ne a$.

3. y = tdlautosol (dd, &t) : Genau wie tdlmeagersol(). Das Gitter liegt nur enger, ist i.a. auch wieder nicht äquidistant. tdlautosol() ist gut geeignet für Plotten. tdlmeagersol() wiederum ist besser geeignet für zahlenmässiges Drucken. Die von den Funktionen tdlsolve(), tdlmeagersol() und tdlautosol() zurückgegebenen Vektoren können zwischen Aufrufen dieser Funktionen beliebig überschrieben werden. Will man jedoch einen Vektor über mehrere Aufrufe hinweg speichern, so muß man den Vektor puffern, da das Programm TENDLER diesen Vektor u.U. überschreiben kann.

4. Mit den bisher beschriebenen Routinen tdldefsys(), tdldefivp() und tdlsolve(), lässt sich ein Großteil von Differentialgleichungen lösen. Bei diesen Routinen wurde nur das an Angaben verlangt, was das mathematische Problem beschreibt. Interessant ist hierbei, daß dies unabhängig von der verwendeten Lösungsmethode ist, hier zyklischen Mehrschrittverfahren. Ein Programm basierend auf Runge-Kutta Formeln könnte man mit völlig gleichlautenden Unterprogrammen genauso aufrufen. Von Dingen wie Anfangsschrittweite, Steifheit, Speicherplatzbedarf, etc., wurde der Benutzer entlastet. Es ist möglich, daß die Genauigkeit, die das Programm TENDLER liefert, zu hoch ist, also man letztlich daran interessiert ist, die Rechengeschwindigkeit zu erhöhen. Umgekehrt kann die gelieferte Genauigkeit zu gering sein. Eine Anhebung der Genauigkeit ist prinzipiell möglich. Es gibt weitere Unterprogramme im Programm TENDLER, mit denen sich dies bewerkstelligen lässt. Zuvor jedoch ein vollständiges Beispiel in ANSI C.

5. Beispiel: Zu lösen sei $\dot y=-y$, $y(0)=1$. Die Funktion $f(t,y)=-y$ ist

int f (double t, double y[], double ydot[])  {
    ydot[0] = -y[0];
}

Gute Übersetzer, oder lint, werden eine Warnung ausgeben, daß in f der Parameter t nicht benutzt wird und kein Rückgabewert explizit übergeben wird. Diese Warnungen können in diesem Falle natürlich ignoriert werden. Gesucht sei $y(2)$.

#include <stdio.h>      /* für printf() */
#include "ode.h"        /* für tdldefsys(), tdldefivp(), tdlsolve(), nojac, nogstop */

int main (void) {
    double *y, yi[]={1.0};
    void *dd = tdldefsys (1, f, nojac, 0, nogstop);

    tdldefivp (dd, 0.0, 3.0, yi);
    y = tdlsolve (dd, 2.0);
    printf ("%e\n", y[0]);
    return 0;
}

Der Wert 3 bei tdldefivp() ist im großen und ganzen willkürlich. Der Wert 0.5 wäre hier ebenfalls akzeptabel gewesen.

Ganz allgemein haben die Funktionen func und jac in ANSI C das Aussehen

int func (double t, double y[], double ydot[]);

bzw.

int jac (double t, double y[], double J[]);

Tatsächlich ist J[] und nicht J[][] anzugeben!

Die Elemente von $J$, also $J_{ik}$, sind spaltenweise in J[] abzulegen, d.h.

$$ \mathtt{J[i+n*k]} = J_{ik}, \qquad\hbox{falls } 0\le i,k\lt n, $$

bzw.

$$ \mathtt{J[(i-1)+n*(k-1)]} = J_{ik}, \qquad\hbox{falls } 1\le i,k\le n. $$

Diese Form der Speicherung ist nicht die standardmäßige Speicherungsform für C Matrizen! Vielmehr entspricht sie der Speicherung wie sie in Fortran üblich ist.

6. Beispiel: Für die lineare Differentialgleichung

$$ \dot y = \pmatrix{ -10 & \alpha\cr -\alpha & -10\cr } y %{-10,\: \alpha\choose -\alpha,{\mskip 3mu}-10}y $$

würde man jac wie folgt programmieren:

int jac (double t, double y[], double J[]) {
    J[0] = J[3] = -10,  J[1] = -α, J[2] = α;
}

Auch hier wäre eine Warnung bzgl. t, y[] und des Rückgabewertes ignorierbar. Ggf. würde man vor den Zuweisungen noch einfügen:

static int loaded = 0;
if (loaded) return;
loaded = 1;

3. Extras

1. Der erfahrene Benutzer wird ggf. wünschen die Genauigkeitsanforderung für den lokalen Fehler, die Maximalordnung, die Anfangsschrittweite, die Maximalschrittweite, etc., selber bereitzustellen. Dies ist möglich. Hierzu dienen die Unterprogramme tdlsetaccuracy(), tdlsetstep() und tdlsetmodes(). Die Parameterlisten erklären hinreichend, was gesetzt wird:

  1. tdlsetaccuracy (dd, tol, ordmax, nfemax),
  2. tdlsetstep (dd, minh, maxh, nexth),
  3. tdlsetmodes (dd, jacm, iterm, printmode, errcm).

I.a. wird lediglich tdlsetaccuracy() aufgerufen und zwar in der Form

tdlsetaccuracy (dd, newtol, AUTO_SELECT, AUTO_SELECT);

Die anderen Optionen, die durch tdlsetstep() und tdlsetmodes() verändert werden können, erhalten durch das Programm TENDLER Vorbesetzungen, die es sich i.a. nicht lohnt zu ändern.

Ist $f$ lediglich $k$-mal stetig differenzierbar, so ist eine Begrenzung der Maximalordnung auf $k-1$ anzuraten, wobei hier

$$ k\in\{2,\ldots,8\}. $$

Eine Begrenzung der Maximalordnung ordmax empfiehlt sich u.U. auch dann, wenn man Kenntnis vom maximal zulässigen Widlund-Winkel besitzt, was bei linearen, zeitinvarianten Differentialgleichungen nicht gänzlich ungewöhnlich wäre.

Diese Unterprogramme sollten i.a. nach tdldefivp() verwendet werden. Sie können zwischen Aufrufen von tdlsolve() oder tdlmeagersol()/tdlautosol() beliebig oft, in beliebiger Reihenfolge aufgerufen werden. I.a. ist dies nicht ratsam, es sei denn, man hat bereits einen sehr detaillierten Kenntnisstand über den Lösungsverlauf.

Übliche Werte

  1. für die Toleranz sind $10^{-4}$ bis $10^{-12}$,
  2. für die Maximalordnung 3 bis 5,
  3. für die maximale Anzahl der Funktionsauswertungen 1000–30000.

Minimale und maximale Schrittweite hängen von $t$ ab. Bei einem Integrationsintervall von $10^{30}$ bis $10^{60}$ ist minh=$10^3$ als sehr klein zu betrachten.

Überhaupt sind diese Einstellungen abhängig von der zu lösenden Differentialgleichung.

Für tdlsetaccuracy() ist $\mathtt{tol}>0$ und $1\le\mathtt{ordmax}\le7$ einzuhalten.

Für tdlsetstep() muß gelten

$$ 0 \lt \mathtt{minh} \le \mathtt{nexth} \le \mathtt{maxh} . $$

tdlsetmodes() ist unterteilt in Exklusivoptionen und Inklusivoptionen. Die Exklusivoptionen sind aus den Wertemengen

$$ \eqalign{ \mathtt{jacm} &\in \{\mathtt{EXACTJ}, \mathtt{APPROXJ}\},\cr \mathtt{iterm} &\in \{\mathtt{MNI}, \mathtt{MNIFLOAT}, \mathtt{FIX}, \mathtt{FIXFLOAT}, \mathtt{FIXPECE}, \mathtt{SWITCHING}\},\cr \mathtt{errcm} &\in \{\mathtt{ECEVSTAGE}, \mathtt{ECEOC}\},\cr } $$

und für die Inklusivoptionen muß gelten

$$ \mathtt{printmode} \in \{\mathtt{PMCVGCE}, \mathtt{PMSSOC}, \mathtt{PMSTAT}, \mathtt{PMWARN}, \mathtt{PMERREX}, \mathtt{PMDUMP}\} . $$

Exklusivoptionen sind Schaltervariablen, die nur einen der angegebenen Werte annehmen können. iterm=MNI und gleichzeitig iterm=FIX ist unzulässig.

Eine suggestive Bedeutungsbeschreibung für Inklusivoptionen lautet: Inklusivoptionen sind Schaltervariablen, die mehrere Werte gleichzeitig inne haben können.

Beispielsweise ist printmode=PMSSOC und gleichzeitig printmode=PMERREX zulässig und sinnvoll.

Genauer: Inklusivoptionen sind Teilmengen aus der Potenzmenge der o.a. angegebenen Menge, also mengenwertige Schaltervariablen.

2. Mit dem Unterprogramm tdlsetaccuracy() lässt sich leicht eine Schätzung des globalen Fehlers realisieren, durch Vergleich von Lösungen berechnet mit verschiedenen Genauigkeiten.

Man definiert sich zwei Kennungen für eine einzige Differentialgleichung: $g$ (genau) und $u$ (ungenau). Man programmiert

$$ \eqalign{ &g = \mathtt{tdldefsys} (n, f, \mathtt{nojac}, 0, \mathtt{nogstop});\cr &u = \mathtt{tdldefsys} (n, f, \mathtt{nojac}, 0, \mathtt{nogstop});\cr &\mathtt{tdldefivp} (g, a, b, y_0);\cr &\mathtt{tdldefivp} (u, a, b, y_0);\cr &\mathtt{tdlsetaccuracy} (g, \varepsilon_g, \mathtt{AUTO\_SELECT}, \mathtt{AUTO\_SELECT});\cr &\mathtt{tdlsetaccuracy} (u, \varepsilon_u, \mathtt{AUTO\_SELECT}, \mathtt{AUTO\_SELECT});\cr } $$

Es sei hierbei $\varepsilon_g < \varepsilon_u$.

Man löst somit die Differentialgleichung $\dot y=f(t,y)$, $y(a)=y_0$ zweimal:

  1. einmal mit der Genauigkeit $\varepsilon_g$ und
  2. einmal mit $\varepsilon_u$.

Der Vorteil ein und dieselbe Differentialgleichung zweimal gleichzeitig zu lösen, liegt darin, daß man nun sehr leicht die Abweichung berechnen kann, nämlich durch Berechnen der Normdifferenz der beiden Vektoren

$$ \eqalign{ y_g &= \mathtt{tdlsolve} (g, t);\cr y_u &= \mathtt{tdlsolve} (u, t);\cr } $$

Würde man die Differentialgleichung zweimal hintereinander lösen, so fielen Pufferungsprobleme an. I.a. wird die Lösung korrespondierend zu $\varepsilon_g$ genauer sein. Allerdings insbesondere bei sehr hohen Genauigkeitsanforderungen (Grenzgenauigkeit) bringt eine weitere Verkleinerung von $\varepsilon_g$ grundsätzlich keine weitere Steigerung der Genauigkeit. Diese Grenzgenauigkeit hängt von der Differentialgleichung ab und wird vom Programm TENDLER ständig berechnet um zu überprüfen, ob sie auch nicht unterschritten wird.

3. iterm dient dazu die Iterationsart fest vorzugeben, außer natürlich bei iterm=SWITCHING, wo es gerade erlaubt ist zu schalten. Dieses Fixieren auf eine Iterationsart kann nützlich sein, wenn man z.B. weiß, daß die Differentialgleichung nur schwach steif ist und es daher ratsam ist, eine der Iterationsarten FIX, FIXFLOAT, FIXPECE von vorneherein vorzugeben.

Es ist erlaubt, daß in tdldefsys() die Jacobimatrix angegeben wird, aber mit tdlsetmodes() numerische Differentiation vorgeschrieben wird. In gelegentlichen Fällen ist dies nützlich. Z.B. wenn die Jacobiamtrix, die in tdldefsys() übergeben, nicht die eigentliche Jacobimatrix $f_y$ ist, sondern eine geschickt gewählte Näherung hiervon, die darauf bedacht sein könnte, Bandstrukturen, Rechenvorteile, o.ä. auszunutzen. Mit der numerischen Differentiation kann man auf einem kurzen Intervall schnell überprüfen, ob die Näherungen mit exakter Jacobimatrix mit den Näherungen der abgekürzten Jacobimatrix hinreichend übereinstimmen. Wie man erkennt, ist dies ein sehr spezieller Fall.

4. Statistiken lassen sich erzeugen durch den Aufruf von tdlprtstat(dd). Die Ausgabe der Statistiken ist in mehrere Abschnitte eingeteilt.

Ein Dump wird durch tdldddump(dd) erzeugt.

Für den sehr an Details interessierten Benutzer ist folgende Information möglicherweise interessant:

  1. Es ist zulässig und sinnvoll ptdlrtstat() bzw. dtdlddump() aufzurufen, aber in tdlprintmode die Schalterstellungen PMSTAT bzw. PMDUMP zu deaktivieren. Es werden dann keine Statistiken bzw. es wird kein Dump gedruckt.
  2. Sobald die Kennung dd mit tdldefsys() erlangt wurde, kann zu jeder Zeit tdlprtstat() aufgerufen werden. Interessant sind diese Daten i.a. nur am Ende einer Integration.

Die Statistiken sind für jede Kennung verschieden und völlig unabhängig voneinander. Dies gilt auch für die Schaltervariablen, ob Inklusiv- oder Exklusivoptionen. Bei dem obigen Beispiel der simplen Schätzung des globalen Fehlers mit Hilfe von $g$ und $u$, werden, falls gewünscht, zwei verschiedene Statistiken ausgegeben. Man erkennt dann, daß i.a. die Berechnung von $g$ aufwendiger ist, d.h. es werden i.d.R. mehr Funktionsauswertungen, Jacobimatrixauswertungen, $LU$-Zerlegungen, etc., durchgeführt.