15th April 2026, 12 min read

Das Schalten im Programm DEASY

• 1. Kurzbeschreibung des Programmes DEASY
• 2. Die nichtlineare von Misessche Iteration
• 3. Bewertung und Vergleich des Programmes DEASY
• 4. Literatur

Das Programm DEASY, siehe Shampine (1984)2, ist eine Zusammenfassung zweier sonst völlig voneinander unabhängiger Differentialgleichungslöser und zwar einem Programm basierend auf linearen Mehrschrittverfahren und einem zweiten Programm basierend auf einem Runge-Kutta Verfahren. Von den hier vorgestellten schaltfähigen Verfahren ist es das einzige, welches zwischen so verschiedenartigen Methoden hin- und herwechselt. Geschaltet wird auch vollständig, so wie dies auch das Programm LSODA tut. Es erfolgt keine unterschiedliche Behandlung der einzelnen Komponenten der Differentialgleichungen, wie dies beispielsweise bei Suleiman/Hall (1985) geschieht, bei dem für das nicht-steife Teilsystem die Adams-Formeln und für das steife Teilsystem die BDF verwendet werden. Von seiner ursprünglichen Konzeption her, ist das Programm DEASY maßgeblich für den unerfahrenen Benutzer erdacht, welcher schnell und unkompliziert seine Gleichungen lösen möchte und sich hierbei nicht um irgendwelche Besonderheiten kümmern will, wie z.B. Dünnbesetztheiten der Jacobimatrix, oder eben Steifheit oder Nicht-Steifheit. Das Programm DEASY ist, wie DGERK, Bestandteil der Unterprogrammsammlung DEPAC. DEPAC ist als Bestandteil von SLATEC darin aufgegangen. Die Geschichte von SLATEC findet man hier.

1. Kurzbeschreibung des Programmes DEASY

1. Das Programm DEASY löst gewöhnliche Differentialgleichungen der Form

Der Fall $t\in[b,a]$ ist zugelassen. Es sei aber für die folgenden Beschreibungen $a<b$ vorausgesetzt.

Als Matrixnorm von $\mathopen|J\mathclose|$ wird benutzt

Diese Norm wird deswegen gebraucht, weil zum einem $\left\Vert J\right\Vert_2 \le \left\Vert J\right\Vert_E$ gilt. Zum anderen weil die Frobenius-Norm mit der euklidischen Norm verträglich ist und im Programm DEASY die euklidische Norm als Vektornorm verwendet wird.

Bei der euklidischer Norm (auch 2-Norm) und bei der Froebenius Norm kann man unitäre Matrizen invariant aus der Norm herausfaktorisieren. Es gilt nämlich

2. Lemma: Für jede Matrix $A\in\mathbb{C}^{n\times n}$ und für beliebige unitäre $X,Y\in\mathbb{C}^{n\times n}$ gilt

Beweis: Für die euklidische Norm gilt

Für die Frobenius-Norm rechnet man

wegen $\mathop{\rm tr} UV = \mathop{\rm tr} VU$, $\forall U\in\mathbb{C}^{m\times n}$, $\forall V\in\mathbb{C}^{n\times m}$.     ☐

3. Lemma: Für eine beliebige Matrix $A\in\mathbb{C}^{n\times n}$ gilt $\left\Vert A\right\Vert_2 \le \left\Vert A\right\Vert_E \le \sqrt n \left\Vert A\right\Vert_2$.

Beweis: Sei $A=USV$ eine Singulärwertzerlegung von $A$ mit unitären $U,V\in\mathbb{C}^{n\times n}$ und $S=\mathop{\rm diag}(\sigma_1,\ldots,\sigma_n)$, wobei o.B.d.A. $\sigma_1\ge\ldots\ge\sigma_n\ge0$ angenommen werden kann. Dann gilt

und

    ☐

4. In dem Programm DEASY wird stets die vollständige Jacobimatrix durch numerische Differentiation bestimmt. Behandlung von Dünnbesetztheiten in der Jacobimatrix findet nicht statt, da das Programm DEASY vorwiegend für den ungeübten Benutzer gedacht ist und so entsprechend konzipiert wurde.

Geschaltet wird in Abhängigkeit einer Näherung für die Lipschitzkonstante $L$. Es wird vollständig zwischen zwei sonst völlig unabhängigen Programmen geschaltet, nämlich DERKF, einem (wenn man so will, dem) Runge-Kutta-Fehlberg Verfahren der Ordnung 4(5) und DEBDF, einer Modifikation des Programmes LSODE.

Beim Wechsel startet das neu ausgewählte Programm, also insbesondere DEBDF, völlig neu, so, als ob das Programm gänzlich neu aufgerufen wurde. Dies heißt also für DEBDF mit der Ordnung $1$ und einer vom Programm berechneten Anfangsschrittweite $h_0$. Gestartet wird stets mit DERKF.

Dem Programm DEASY werden beim Aufruf die beiden Intervallgrenzen $a$ und $b$ mit angegeben. Diese sind somit bekannt und können mit zur Schaltentscheidung herangezogen werden.

In dem Programm TENDLER ist dies nicht der Fall. Das Programm TENDLER benutzt für das Schalten keinerlei Informationen bzgl. der Intervallgrenzen.

Es sei $a<t<b$ vorausgesetzt. Der Wert $L$ bezeichne hier eine vom Programm DEASY berechnete Näherung für die Lipschitzkonstante der Funktion $f$. Der Schaltteil sieht nun wie folgt aus:

1. Falls $(b-t)L\ge300$ und $t+15h\ge b$, dann schalte von DERKF nach DEBDF.
2. Falls $(b-t)L\le100$ und $t+15h\ge b$, dann schalte von DEBDF nach DERKF.

Die Bedingung $t+15h\ge b$ in (1) und (2) erzwingt, daß nicht geschaltet wird, falls wahrscheinlich nur noch $15$ Schritte bis zum Integrationsende durchgeführt werden müssen. Die Bedingungen für $L$ testen, ob die Lipschitzkonstante “groß” oder “klein” ist. Dabei geht bei diesem Test der noch zurückzulegende Teil bis zum Integrationsende in die Überlegungen ein.

Die Näherung für die Lipschitzkonstante $L$ wird in DEBDF, also im steifen Falle, “direkt” berechnet über die Norm $\left\Vert J\right\Vert_E$, während hingegen sie im nicht-steifen Falle, also DERKF, über die von Misessche Iteration gewonnen wird. Im steifen Modus wird $\left\Vert J\right\Vert_E$ nur bei einer Refaktorisierung der Iterationsmatrix berechnet. Die Überprüfung, ob nun geschaltet werden soll oder nicht, findet beim Runge-Kutta-Fehlberg Verfahren alle 15 Schritte statt. Dies dient dazu einerseits die Integration nicht durch übermäßig vieles Schalten zu destabilisieren und andererseits den Aufwand zur Erkennung von Steifheit niedrig zu halten.

Der Shampine-Test findet in dem Programm DEASY keine Verwendung. Der Shampine-Test, der sich in das Programm TENDLER vergleichsweise einfach einbauen ließ, zeigte als Entscheidungshilfe zum Schalten einen eher gering einzuschätzenden Einfluß. Dies zeigen die generierten Statistiken. Dennoch eignet sich der Shampine-Test u.a. für die Beurteilung ob die Differentialgleichung insgesamt als steif oder als nicht-steif anzusehen ist. Dieses Ergebniss steht jedoch erst nach der vollständigen Integration fest, nicht jedoch während der Integration, also an der Stelle, an der es für eine Schaltentscheidung wünschenswert wäre.

2. Die nichtlineare von Misessche Iteration

1. Die nachstehende Überprüfung, die maximal alle 15 Schritte durchgeführt wird, erfolgt im nicht-steifen Falle, also Benutzung von DERKF. Man möchte eine Näherung für die Lipschitzkonstante $L$ erhalten. W.o. angegeben, wird ja aufgrund einer Näherung einer Lipschitzkonstanten geschaltet. Die Funktion $f$ wird dabei höchstens 5-mal ausgewertet. Man erhält eine untere Schranke für die Lipschitzkonstante $L$ der Funktion $f$. Shampine (1984)2 geht nun wie folgt vor: Die erste Iterierte $y^1$ wird definiert durch

wobei $y^0=y_k$ und $f^0=f(t_k,y_k)$ ist, mit $a\le t_k\le b$. Die skalare Größe $\rho_0$ wird derart gewählt, daß die Differenz $y^1-y^0$ “klein” ist. Hinreichende Bedingungen für Konvergenz der von Misesschen Iteration gegen einen betragsgrößten Eigenwert sind, daß ein eindeutiger dominanter Eigenwert vorliegt und man nicht im Komplementärraum zum dominanten Eigenvektor startet. Diese Bedingungen führen zu gewissen Restriktionen an die Jacbimatrix von $f$. Shampine (1984)2 äußert sich nicht weiter zu der Wahl dieses Wertes $\rho_0$. Für $m=1,2,3,4,5$ wird nun weiter definiert

und schließlich wird $y^{m+1}$ gewonnen durch

Nach obiger Festlegung der $\rho_m$ gilt stets, daß $\rho_m\le L$, also ist $\rho_m$ eine untere Schranke für die Lipschitzkonstante $L$ für die Funktion $f$. Weiter ist nach Definition der $\rho_m$

2. Um den Zusammenhang mit der von Misesschen Iteration, zur Bestimmung des betragsgrößten Eigenwertes einer Matrix, besser zu erkennen, sei auf die Gleichung $(*)$ in jeder Komponente der Mittelwertsatz für $\mathbb{R}^n\to\mathbb{R}$ Funktionen angewandt. Dann ergibt sich

wobei $J_m$ die “Jacobimatrix” der Funktion $f$ ist, ausgewertet bei $t_k$. Die Zeilen der Jacobimatrix $J$ werden an jeweils unbekannten Zwischenstellen auf der Strecke von $y^0$ nach $y^m$ ausgewertet. Hält man nun $y^m$ nahe bei $y^0$, im Sinne einer Norm, so ist es naheliegend $J\approx J_m$ anzunehmen, für alle $m$. Man hat also nun

Jetzt wird deutlich, daß man tatsächlich den dominanten Eigenwert der Matrix $J$ mit dem obigen Verfahren auszurechnen versucht. Natürlich gilt dies nur bei gewissen Voraussetzungen an die Matrix $J$ und damit letztlich an die Funktion $f$.

3. Alternativ kann man den Grenzwert obiger Iteration, falls vorhanden, mit der Lipschitzkonstanten wie folgt in Zusammenhang bringen. Man erkennt, daß $\rho_m^2$ der Rayleigh-Quotient für den größten Eigenwert von $J_m^\top J_m = \left\Vert J_m\right\Vert_2^2$ ist. Dies zeigt erneut die Verbindung zu einer Lipschitzkonstanten.

3. Bewertung und Vergleich des Programmes DEASY

Die Ergebnisse, die mit dem Programm DEASY gewonnen wurden, zeigen, daß für den Zweck für den das Programm DEASY konzipiert wurde, das Programm auch geeignet ist.

• Im Falle nicht-steifer Differentialgleichungen (nonstiff DETEST problem set) nimmt die Rechenzeit nur um 5.2% zu. Dies ist ein sehr günstiger Wert für den ungünstigen Falle, nämlich vollständig nicht-steife Differentialgleichungen.

• Im Falle steifer Differentialgleichungen (restricted stiff DETEST problem set) nimmt die Rechenzeit lediglich um 2.2% zu.

Diese beiden Zahlen verdeutlichen, daß selbst ein recht einfacher Weg des Schaltens dem Benutzer hilft, das Erkennen und Vorhersagen von Steifheit und Nichtsteifheit im großen und ganzen völlig zu ignorieren, falls der konzipierte Verwendungszweck des Programmes DEASY beachtet wird.

Vergleicht man die beiden Programme DERKF und DEASY an den nicht-steifen Differentialgleichungen von DETEST, so ergibt sich das untenstehende Bild.

Die Rechenzeiten gelten für eine CDC6600. Es wurde stets mit absoluter Fehlertoleranz gearbeitet.

Beim Vergleich an den steifen Differentialgleichungen von DETEST ergab sich das folgende Bild. Die Probleme B2, B3, E2 und E4 wurden aus der Gruppe der Testdifferentialgleichungen herausgenommen. Die ersten drei Gleichungen sind nicht wirklich steif. Sie liefern daher daher keine zusätzlichen Informationen. Das Problem E4 reagiert sehr empfindlich auf schon kleinere Änderungen der Anfangswerte, ist also in diesem Sinne mathematisch nicht korrekt gestellt. Dies sind die Gründe sie nicht in den Vergleich mit aufzunehmen.

Shampine (1981b) äußert sich in größerer Ausführlichkeit zu den Testdifferentialgleichungen in DETEST.

Es ergab sich nun:

Selbstverständlich gibt es Differentialgleichungen —und Shampine (1984)2 gibt diese auch an— für die sich die Rechenzeiten verdoppeln können, überhaupt also die Benutzung des optimalen Lösers erheblich effizienter ist als die Benutzung des schaltfähigen Programmes. Dies geschieht gerade bei denjenigen Differentialgleichungen, die nicht nur eine “große” Lipschitzkonstante besitzen, sondern hinzu noch instabil sind, wie z.B. das Differentialgleichungssystem ND5, welches gerade das Zweikörperproblem in elliptischer Bewegung darstellt. Die Exzentrizität beträgt dort gerade $\varepsilon=0.9$. Diese Differentialgleichung hat am Anfang der Integration die folgenden Eigenwerte der Jacobimatrix

Die Norm der Jacobimatrix beträgt hier $\left\Vert J\right\Vert_2 = 2000$. Zudem ist diese Gleichung instabil an den Stellen nahe $t\approx6.5$ und nahe $t\approx12.7$; Shampine bemerkt hierzu:

$\ldots$ which code is best in a region of instability is not clear. In this instance using the BDF code was acceptable but not optimal.

Speziell für das Problem ND5 gibt Shampine (1984)2 die Werte in der nachfolgenden Tabelle an.

Man erkennt deutlich, daß in diesem Falle nicht die günstigste Wahl des Integrators getroffen wurde. Jedoch sind bei so unbedeutenden Rechenzeiten, diese Überlegungen von untergeordneter Bedeutung. Auf Kleinrechnern und bei nicht-trivialen Differentialgleichungen wäre diese Aussage allerdings wieder abzuändern. Auffallend ist, daß die Anzahl der Funktionsauswertungen in beiden Fällen annähernd gleich ist. Bei sehr aufwendig auszuwertenden Differentialgleichungen, wie sie z.B. bei Homotopieverfahren auftreten, würden sich die Rechenzeitunterschiede weiter verwischen.

Viel wichtiger ist der Bequemlichkeits- und Sicherheitsgewinn den der Benutzer erhält.

1. Bequemlichkeitsgewinn deswegen, weil der Benutzer sich zum einen nur mit einer Programmdokumentation vertraut machen muß und zum anderen
2. Sicherheitsgewinn deswegen, weil der Benutzer nicht mit einem Programm konzipiert für nicht-steife Differentialgleichungen auf eine steife Differentialgleichung losgehen kann.

Es ist gerade der letzte Fall, der rechenzeitmässig besonders ineffizient ist. Hier liegen die Geschwindigkeitsfaktoren bei 100 bis 1000 und noch höher.

Der Vorteil des Kennenlernens nur einer Programmdokumentation ist nicht zu unterschätzen. Eine Differentialgleichung, die mit dem Programm DE/STEP gelöst wurde, kann nicht mit dem gleichen Aufrufprogramm auch mit dem Programm GEAR gelöst werden. Dies obwohl beide Löser die Differentialgleichung in gleicher Standardform, also $\dot y=f(t,y)$, akzeptieren und sogar letztlich mit den völlig gleichen Formeln lösen. Es ist sogar weitaus mehr als nur die Argumentenliste zu verändern. Die zahlreichen “Flags”, Fehlermarkierungen, Rücksprünge und Speicherverwaltungen unterscheiden sich nicht unerheblich.

Im weiteren werden schaltfähige Programme vorgestellt, die maßgeblich auf Runge-Kutta Verfahren basieren. Stellenweise wurden hierzu originär neue Integrationsformeln entwickelt, was bei den Programmen LSODA, DEASY und TENDLER nicht der Fall war. Damit wird allerdings nicht nur eine neue Idee des Schaltens verglichen, sondern auch die neuentwickelte Formel samt ihrer programmiermässigen Realisation. Hier treten dann wieder übliche Probleme der Bewertung auf. Nämlich, ist die Idee des Schaltens in der verwendeten Form und in der verwendeten Programmrealisierung erfolgreich, oder liegt das gute oder schlechte Funktionieren nur an den Formeln? Wenn das Programm wenig wirkungsvoll arbeitet, liegt dies an der aufgezwungenen Schaltung?

4. Literatur

1. Mises, Richard Edler von (1883–1953)
2. Rayleigh, Lord John William Strutt (1842–1919)
3. Shampine, Lawrence Fred, Math Genealogy
4. Shampine, Lawrence F.: “Evaluation of a Test Set for Stiff ODE Solvers”, ACM TOMS, Vol 7, No 4, December 1981, pp.409–420
5. Shampine, Lawrence Fred: “Stiffness and the Automatic Selection of ODE Codes”, Journal of Computational Physics, Vol 54, No 1, April 1984, pp.74–86
6. Suleiman, M.B. und Hall, George: “A Single Code for the Solution of Stiff and Nonstiff ODE's”, SIAM Journal on Scientific and Statistical Computing, Vol 6, No 3, July 1985, pp.684–697