, 13 min read

Das Schalten im Programm RKF4RW

Das Programm RKF4RW, siehe Rentrop (1985), ist hier nun eines der Programme, welche ausschliesslich Runge-Kutta-Fehlberg und Rosenbrock-Wolfbrandt-Verfahren benutzen und zwischen diesen beiden Verfahren hin und her schalten. Neu ist hierbei ebenfalls, daß nicht nur lediglich zwischen bewährten klassischen Verfahren geschaltet wird, sondern das zusätzlich speziell neue Formeln entwickelt und gesucht wurden. Das neugefundene ROW-Verfahren enthält das klassische Runge-Kutta-Fehlberg Verfahren in sich eingebettet. Damit ist schon eine gewisse Strategie der Fehlerkontrolle und der Schrittweitensteuerung nahe gelegt.

Entgegen den beiden bisher vorgestellten linearen Verfahren, bei denen durch die Schaltfähigkeit eine weitere Effizienzsteigerung, wenn auch sehr geringe, erreicht werden konnte, ist dies bei dem Programm RKF4RW nicht mehr in solchem Maße der Fall. Sogar das Gegenteil tritt ein.

Von seiner Grundkonzeption her ist das Programm RKF4RW eher nur für geringe bis mittlere Genauigkeitsanforderungen gedacht, also Genauigkeiten, die sich i.d.R. im Bereich von $\varepsilon=10^{-3}$ bis $\varepsilon=10^{-6}$ befinden. Für Runge-Kutta Verfahren ist es unüblich die Ordnung zu wechseln. Dies muß der Benutzer selber erledigen, indem er ein anderes Programm aufruft. Bei den üblichen Programmen basierend auf linearen Mehrschrittverfahren, wie z.B. LSODA, DE/STEP, aber auch TENDLER, geschieht die Wahl einer Ordnung vollautomatisch. Vom Benutzer werden keinerlei Sondermaßnahmen gefordert. Extrapolationsverfahren, wie z.B. das Programm METAN1, siehe Deuflhard (1983), siehe Deuflhard (1985), passen die Ordnung ebenfalls den Genauigkeitsforderungen an, neben der Schrittweite, die von allen bisher erwähnten Programmen variiert wird.

Shampine (1983) erweiterte das Programm METAN1 zu einem Programm METAN3, welches ggü. dem Programm METAN1 eine deutliche Verbesserung darstellt. Zum einem ist das Programm METAN3 schaltfähig, basierend auf der direkten Berechnung von $\left\Vert J\right\Vert_1$, und zum anderen wurden gewisse Ineffizienzen vermieden.

Runge-Kutta Verfahren und Extrapolationsverfahren arbeiten jedoch nicht recht wirkungsvoll, wenn man viele Ausgabepunkte wünscht oder benötigt, die vom Benutzer vorgegeben werden. In diesem Falle sind Programme basierend auf linearen Mehrschrittverfahren i.d.R. vorzuziehen.

1. Kurzbeschreibung des Programmes RKF4RW

Das Verfahren von Runge-Kutta. Es handelt sich hier um ein für die numerische Rechnung sehr brauchbares Näherungsverfahren, das letzten Endes auf eine Verallgemeinerung der Simpsonschen Formel hinausläuft. Ich will mich damit begnügen, Ihnen die Formel anzuschreiben und die Durchführung des Verfahrens an einem Beispiel zu illustrieren; auf die Begründung, die einen erheblichen und nicht sehr instruktiven Rechenaufwand erfordert, verzichte ich.

Adalbert Duschek (1953)

1. Das Programm RKF4RW löst Differentialgleichungen der Form

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

Das Programm RKF4RW schaltet zwischen dem klassischen Runge-Kutta-Fehlberg Verfahren der Ordnung 4(5) mit 6 Stufen und einem eingebettetem Rosenbrock-Wolfbrandt Verfahren (ROW-Verfahren) der Ordnung 3(4).

Geschaltet wird, wenn überhaupt, nur nach $ \def\hgJ{\left\Vert h\gamma J\right\Vert} k = \max (n, 16) $ Schritten. Zu den Werten $n$ und 16 vermerkt Rentrop (1985), daß beide Programme, sowohl PRK4 als auch RKF4RW, nicht sehr empfindlich auf Änderungen dieses Faktors reagieren. Werte zwischen 8 und 20 seien durchaus effizient.

2. Vom nicht-steifen Teil wird in den steifen Teil gewechselt, falls die Schrittweite, die das ROW-Verfahren liefern würde, um den Faktor hfak = $1+(n/6)$ besser ist, als der Schrittweitenvorschlag durch das Runge-Kutta-Fehlberg Verfahren. Dieser Faktor hfak beeinflußt maßgeblich die Effizienz der beiden Programme. Dennoch ist die “Einstellung” dieses Faktors nicht einfach. Ein kleinerer Wert würde dazu führen, daß die Differentialgleichung stets als steif betrachtet würde, also somit nicht geschaltet würde und ein vergrößerter Faktor hätte die Eigenschaft, daß viel zu spät in den steifen Modus zurückgeschaltet würde.

Die Schrittweitensteuerung begrenzt Verkleinerungen und Vergrößerungen der Schrittweite auf Faktoren liegend im Intervall $[0.5, 1.5]$. Gewählter Sicherheitsfaktor für die Gewichtung des geschätzten lokalen Fehlers ist $9/10$.

3. Für den Wechsel von steif nach nicht-steif verlangt man mehr, als nur einen Schrittweitenvorteil. Rentrop (1985) gibt an, daß er während der $LU$-Zerlegung der Iterationsmatrix $W=I-h\gamma J$, gleichzeitig auch die Maximumnorm $\left\Vert W\right\Vert_\infty$ während der Pivotisierung erhält. Diese Norm wird nun dazu benutzt, um festzustellen, ob $\left\Vert h\gamma\right\Vert \rho(J) < 2.4$ ist, wobei $2.4$ der Radius des Stabilitätsbereiches des Runge-Kutta-Fehlberg Verfahrens ist. Es wird nun verlangt, daß sogar die leicht gröbere Forderung $\hgJ_\infty < 2.4$ erfüllt wird. Als Matrixnorm wird die Maximumnorm $|\cdot|_\infty$ verwendet. Mit $\gamma=1/2$ erhält man dann als notwendige, nicht hinreichende, Bedingung für Arbeiten im Stabilitätsbereich, daß

$$ \left\Vert W\right\Vert_\infty \le 1+\hgJ_\infty = 1+0.5\cdot 2.4 = 2.2 $$

Dieser Test, zusammen mit einem Schrittweitenvorteil, führt zum Schalten. Der Schrittweitenvorteilsfaktor ist, w.o. auf hfak = $1+n/6$ gesetzt.

4. Das verwendete (klassische) Runge-Kutta-Fehlberg Verfahren der Ordnung 4(5) lautet, wie in der Tabelle sichtbar.

$$ \begin{array}{c|cccccc} 0 & & & & & &\\[6pt] \displaystyle\frac{1}{4} & \displaystyle\frac{1}{4} & & & & &\\[6pt] \displaystyle\frac{3}{8} & \displaystyle\frac{3}{32} & \displaystyle\frac{9}{32} & & & &\\[6pt] \displaystyle\frac{12}{13} & \displaystyle\frac{1932}{2197} & \displaystyle\frac{7296}{2197} & 0 & & &\\[6pt] 1 & \displaystyle\frac{439}{216} & -8 & \displaystyle\frac{3680}{513} & \displaystyle-\frac{845}{4104} & &\\[6pt] \displaystyle\frac{1}{2} &\displaystyle-\frac{8}{27} & 2 & \displaystyle-\frac{3544}{2565} & \displaystyle\frac{1859}{4104} & \displaystyle-\frac{11}{40} &\\[6pt] \hline y_1 & \displaystyle\frac{25}{216} & 0 & \displaystyle\frac{1408}{2565} & \displaystyle\frac{2197}{4104} & \displaystyle-\frac{1}{5} & 0\\[6pt] \hat y_1 & \displaystyle\frac{16}{135} & 0 & \displaystyle\frac{6656}{12825} & \displaystyle\frac{28561}{56430} & \displaystyle-\frac{9}{50} & \displaystyle\frac{2}{55} \\[6pt] \end{array} $$

Hierbei werden die Werte

$$ 0, \quad\frac{1}{4}, \quad\frac{3}{8}, \quad\frac{12}{13}, \quad 1 \quad\text{und}\quad \frac{1}{2} $$

nicht benutzt, da vom Benutzer des Programmes RKF4RW verlangt wird, die Differentialgleichung in autonomer Form bereitzustellen, also die Gleichung die Form hat $\dot y=f(y)$.

Shampine (1982a), “Implementation of Rosenbrock Methods”, entschliesst sich in seinem schaltfähigen Differentialgleichungslöser DGERK, welcher Bestandteil des Paketes DEPAC ist, die nicht-autonome Form, wie üblich, beizubehalten. Eine Umwandlung einer nicht-autonomen Gleichung in eine autonome Gleichung hat folgende Nachteile:

  1. Das übliche software-interface wird zerstört, und es ist ungewohnt. Dies ist wichtig für größere Unterprogrammbibliotheken, wie z.B. DEPAC.
  2. Bandstruktur kann verloren gehen und jede zusätzliche Gleichung lässt den Aufwand quadratisch ansteigen.
  3. Ein lineares Problem kann nicht-linear werden, man vergl. beispielsweise das Problem D1.
  4. Neben $f_y$ muß zusätzlich $f_t$ berechnet werden, wenn der Benutzer die Jacobimatrix direkt bereitstellen will.
  5. Ein Fehlerkriterium für die unabhängige Variable $t$ ist nicht leicht angebbar.
  6. Die Norm $\left\Vert J\right\Vert_1 = \max(\left\Vert f_y\right\Vert_1, \left\Vert f_t\right\Vert_1)$ kann sich ändern. Dies ist für ein schaltfähiges Programm, welches die Schaltentscheidung auf der Norm der Jacobimatrix $J$ beruhen lässt, ein sehr unerwünschter Effekt.

Rentrop gelang es in das obige Runge-Kutta-Fehlberg Verfahren der Ordnung 4(5) mit 6 Stufen, ein ROW-Verfahren der Ordnung 3(4) einzubetten. Das Rosenbrock-Verfahren ist dabei sogar $A$-stabil, jedoch nicht $L$-stabil. Das Runge-Kutta-Verfahren 4(5) wird wie üblich ohne lokale Extrapolation benutzt, also es wird mit der Ordnung 4 weitergerechnet, und das Verfahren 5.ter Ordnung wird nur zur Schrittweiten-Steuerung benutzt.

5. Das ROW-Verfahren, welches in dem Programm RKF4RW benutzt wird, lässt sich beschreiben durch die Formeln

$$ \eqalign{ y_{n+1} &= y_n+(c_1k_1+c_2k_2+c_3k_4+c_4k_4+c_5k_5+c_6k_6),\cr \hat y_{n+1} &= y_n+(\hat c_1k_1+\hat c_2k_2+\hat c_3k_3+\hat c_4k_4+ \hat c_5k_5+\hat c_6k_6),\cr } $$

mit

$$ (I-h\gamma J)k_i = hf(y_n+\sum_{j=1}^{i-1}\alpha_{ij}k_j) + hJ\cdot\sum_{j=1}^{i-1}\gamma_{ij}k_j, \qquad\hbox{für}\quad i=1,2,3,4,5,6 $$

und

$$ J=\cases{\hbox{Jacobimatrix von $f$}, &falls steif;\cr \hbox{Nullmatrix} , &falls nicht steif.\cr} $$

Die Werte $\alpha_{ij}$, $c_i$ und $\hat c_i$ sind w.o., und für die Werte

$$ \tilde\gamma_{ij} := {\gamma_{ij}\over\gamma},\qquad\gamma={1\over2}, $$

erhält man mit dem Parameterschema

$$ \begin{array}{|ccccc} \tilde\gamma_{21}& & & & & \\ \tilde\gamma_{31}&\tilde\gamma_{32}& & & & \\ \tilde\gamma_{41}&\tilde\gamma_{42}&\tilde\gamma_{43}& & & \\ \tilde\gamma_{51}&\tilde\gamma_{52}&\tilde\gamma_{53}&\tilde\gamma_{54}& & \\ \tilde\gamma_{61}&\tilde\gamma_{62}&\tilde\gamma_{63}&\tilde\gamma_{64}&\tilde\gamma_{65}& \\ \hline \end{array} $$

die von Rentrop berechneten Werte

$$ \begin{array}{|ccccc} -3.7 & & & & & \\ -0.05395888512 &-0.4960411149 & & & & \\ 3.592176012 &6.782987834 &-16.22024878 & & & \\ 8.404898124 &15.19965650 &-37.17118671 &1.411793372 & & \\ 3.022508944 &-8.091206674& 4.763352827 &0.6940545809 &1.15 \\ \hline \end{array} $$

Mit den gegebenen Konstanten $\alpha_{ij}$, $c_i$ und $\hat c_i$ durch das Runge-Kutta-Fehlberg Verfahren 4(5), mußte Rentrop 9 zusätzliche nichtlineare Gleichungen lösen:

Er hatte jedoch 15 freie Parameter $\gamma_{ij}$, und $\gamma$ ist ebenso noch frei verfügbar.

Die nichtlinearen Gleichungen wurden nun durch einen gemischten Ansatz von analytischer Lösungstechnik und numerischer Lösung, mit Hilfe eines Newton-Verfahrens angegangen. Rentrop (1985) entwickelte auch ein Verfahren, welches komponentenweise Steifheit behandelt. Hier sind dann andere Ableitungen erforderlich. Man vgl. hier auch das Buch von Hairer/Wanner/Nørsett (1987), S.276–285.

2. Bewertung des Programmes RKF4RW

Nur erwähnt werden soll hier das Verfahren von Runge und Kutta, über das man in zahlreichen Lehrbüchern nachlesen kann, auf dessen Wiedergabe aber hier aus Raumgründen verzichtet werden muß. Diese Methode der numerischen Integration von Differentialgleichungen verschieden hoher Ordnung hat den Differenzenmethoden gegenüber einige bemerkenswerte Vorteile, vor allem den, daß ein Übergang auf andere (insbesondere kleinere) Schrittweiten während des Integrationsprozesses ohne Komplikationen gelingt, während dies bei den Differenzenverfahren immer umständlich und lästig ist. Dagegen besteht der Nachteil, daß der Ablauf des Runge-Kuttaschen Verfahrens undurchsichtig bleibt und daß die Rechnung der vielfachen Kontrollmöglichkeiten ermangelt, die das Differenzenverfahren auszeichnet. Schließlich ist bei Runge-Kutta Verfahren auch die Genauigkeit der Resultate bei gleichem Aufwand merklich geringer, wenn auch Varianten dieser Methode (z.B. Runge-Kutta-Fehlberg) diesem Übelstand abzuhelfen bemüht sind.

Karl Stumpff (1965)

Rentrop (1985) vergleicht seine beiden Programme RKF4RW und PRK4 mit den Programmen RKF4 und GRK4. Zusammenfassend stellt er fest, daß das Programm RKF4RW besser arbeitet als das Programm PRK4 mit der komponentenweisen Steifheitsbehandlung, obwohl der Unterschied nicht allzu groß ist. Die Testdifferentialgleichungen entnimmt er DETEST. Im einzelnen wählt Rentrop (1985) die Testgleichungen NA2, NB3, NC5, ND1 und NE3 als nicht-steife Gleichungen und A2, B5, C4, D3 und E3 als steife Testsysteme. Schließlich testet er seine Programme auch noch an der van der Polschen Gleichung

$$ \ddot y = p(1-y^2)\dot y - y, \qquad \left\{\eqalign{y(0)&=2,\cr \dot y(0)&=0,\cr}\right. \qquad t\in[0,50] $$

für die Werte $p\in\{5, 10, 1000\}$ in Rentrop (1984) und äußert sich zu der Benutzung des Programmes auf Kleinstrechnern. Es werden nicht sämtliche Testdifferentialgleichungen aus DETEST benutzt. Warum nicht sämtliche Testgleichungen aus DETEST durchgegangen werden, wird nicht weiter begründet.

Getestet wird nur mit der Genauigkeit von $\varepsilon=10^{-4}$. Als Anfangsschrittweite wird stets $h_0=10^{-3}$ gewählt. Aussagen zum globalen Fehler werden in Rentrop (1984) und in Rentrop (1985) nicht gemacht.

Bei der Besprechung des Programmes PAI4 von Bruder/Strehmel/Weiner (1988), wird nocheinmal kurz auf Vergleichsergebnisse mit dem Programm RKF4RW eingegangen. Hierbei stellt sich dann heraus, daß das Programm RKF4RW fast durchweg weniger befriedigendere Ergebnisse liefert als das Programm PAI4. Ebenso zeigt sich, daß LSODA vergleichsweise bessere Ergebnisse liefert und zudem wesentlich mehr Einstellungsmöglichkeiten im Rahmen der linearen Algebra Routinen zu bieten hat. Dennoch werden in dem Programm RKF4RW mit zum ersten Male Schaltfähigkeit und ROW-Verfahren miteinander verbunden und ihre Leistungsfähigkeit untersucht und getestet.

Im einzelnen ermittelt Rentrop (1985) nun im Vergleich der beiden Programme RKF4RW und PRK4, die folgenden Werte. In dem Programm PRK4, welches versucht komponentenweise Steifheit zu erkennen und zu behandeln, wird ein Verfahren dritter Ordnung in ein $A$-stabiles und sogar $L$-stabiles Verfahren vierter Ordnung eingebettet. Die bei dem Programm RKF4RW benutzten Verfahren sind $A$-, aber nicht $L$-stabil, also nicht $A_\infty^0$- bzw. $S_\infty^0$-stabil. Die Ergebnisse wurden auf einer CYBER 175 (48 bit Mantissenlänge${}\approx14$ Dezimalstellen) gewonnen.

  Dgl. A2 B5 C4 D3 E3 NA2 NB3 NC5 ND1 NE3
RKF4RW: CPU 0.24 0.44 0.23 0.12 0.23 0.01 0.04 0.29 0.18 0.23
  nfe 1165 2424 1370 780 1532 169 302 174 1270 1638
PRK4: CPU 0.20 0.22 0.44 0.11 0.17 0.01 0.04 0.60 0.26 0.25
  nfe 840 1296 2839 671 1189 181 333 329 1684 1753

Man ersieht, daß der Unterschied der beiden Programme nicht beträchtlich ist. Dennoch, durch die eingeschränkte Wahl der Testprobleme aus DETEST sind größere Unterschiede auch nicht zu erwarten, da die größer dimensionalen Testdifferentialgleichungen, wie z.B. NC1, NC2, NC3 und NC4, fehlen. Bei der größten Differentialgleichung, hier NC5, dem 5 Körperproblem beschrieben durch insgesamt 30 Gleichungen, wird das Profil deutlicher.

Summiert man Funktionsauswertungen und Rechenzeiten auf, so erhält man die Tabelle wie nachfolgend.

  steif nicht-steif Summe
RKF4RW: CPU 1.26 0.75 2.01
  nfe 7271 3553 10824
PRK4: CPU 1.14 1.16 2.30
  nfe 6835 4280 11115

Von Interesse ist jedoch nicht nur der Vergleich zweier schaltfähiger Programme, sondern auch der Vergleich von Programmen, die speziell nur für einen einzigen Aufgabentyp konzipiert sind. Hier sind die Ergebnisse abweichend von demjenigen Verhalten, wie man es bei den Programmen DEASY, siehe Shampine (1984)2, und LSODA, siehe Petzold (1983a), her kannte.

Dieses abweichende Verhalten fällt umso mehr auf, da der zwei Jahre früher erschienene Artikel von Petzold (1983a) deutlich machte, daß nicht nur lediglich ein Bequemlichkeitsgewinn erzielt werden kann, sondern u.U. (abhängig von der Jacobimatrixauswertung) sogar ein bei bestimmten Testdifferentialgleichungen (wenn auch nur geringer) Rechenzeitvorteil möglich ist. Auch die ein Jahr vorher gewonnenen Ergebnisse von Shampine (1984)2, deuteten in diese Richtung.

Verglichen wurde nun von Rentrop (1985):

  1. das Programm RKF4RW, welches nur im steifen Modus betrieben wurde, mit dem Programm GRK4A. Das Programm GRK4 von Kaps und Rentrop benutzt ebenfalls Rosenbrock-Verfahren. Auch wurde das Programm RKF4RW mit dem Programm PRK4 verglichen, welches ebenso nur im steifen Modus betrieben wurde.
  2. Für die nicht-steifen Differentialgleichungen, aus dem w.o. angegebenen eingeschränkten Vorrat an Testgleichungen von DETEST, wurde das Programm RKF4 mit dem Programm PRK4 verglichen. Das Programm PRK4 wurde beim obigen Vergleich ausschließlich im nicht-steifen Modus betrieben.

Das Programm RKF4 verwendet die gleichen Formeln wie das Programm RKF4RW im nicht-steifen Modus, nämlich das Runge-Kutta-Fehlberg Verfahren der Ordnung 4(5). I.d.S. kann man RKF4 auffassen als “RKF4RW-nicht-steif”. Es ergab sich hierbei:

Dgl. Programm CPU nfe
steif GRK4A 0.43 3694
  PRK4-steif 0.88 5093
  RKF4RW-steif 0.97 6881
nicht-steif RKF4 0.28 1395
  PRK4-nicht-steif 1.06 5686

Bedenkt man, daß die Laufzeiten von “RKF4-steif” günstiger ausfallen, als die Laufzeiten von RKF4RW, so erkennt man, daß im steifen Falle mit Laufzeitverlängerungen um ca. den Faktor 2 gerechnet werden muß. Für nicht-steife Gleichungen ist mit Geschwindigkeits_einbußen_ mit bis zum Faktor 5 zu rechnen. Diese Resultate sind im Lichte mit den anderen schaltfähigen Programmen, die hier vorgestellt werden, als eher atypisch zu bezeichnen. Rentrop bemerkt hierzu:

The increase of computing time up to a factor 3 seems to be tolerable.

3. Bibliographie

  1. Bruder, Jürgen und Strehmel, Karl und Weiner, Rüdiger: “Partitioned Adaptive Runge-Kutta Methods for the Solution of Nonstiff and Stiff Systems”, Numerische Mathematik, Vol 52, Fasc 6, 1988, pp.621–638
  2. Deuflhard, Peter Jochen (1944 – 2019)
  3. Deuflhard, P.: “Order and Stepsize Control in Extrapolation Methods”, Numerische Mathematik, Vol 41, 1983, pp.399–422
  4. Deuflhard, P.: “Recent Progress in Extrapolation Methods for Ordinary Differential Equations”, SIAM Reviews, Vol 27, No 4, December 1985, pp.505–535
  5. Duschek, Adalbert Ludwig (1895–1957)
  6. Hairer, Ernst (*1949)
  7. Hairer, E., Wanner, G., Nørsett, S.: Solving Ordinary Differential Equations I – Nonstiff Problems. Springer Berlin, Heidelberg (2008), https://doi.org/10.1007/978-3-540-78862-1
  8. Nørsett, Syvert Paul (1944–2025)
  9. Petzold, Linda Ruth (*1954)
  10. Petzold, Linda Ruth: “Automatic Selection of Methods for Solving Stiff and Non-Stiff Systems of Ordinary Differential Equations”, SIAM Journal on Scientific and Statistical Computing, Vol 4, No 1, March 1983, pp.136–148
  11. Rentrop, Peter (*1948)
  12. Rentrop, Peter: “The Performance of a Type-Insensitive Runge-Kutta Code on a Personal Computer”, TUM-I8410, Juli 1984, München, Institut für Informatik, ii+13 S.
  13. Rentrop, Peter: “Partitioned Runge-Kutta Methods with Stiffness Detection and Stepsize Control”, Numerische Mathematik, Vol 47, 1985, pp.545–564
  14. Shampine, Lawrence Fred
  15. Shampine, Lawrence Fred: “Implementation of Rosenbrock Methods”, ACM TOMS, Vol 8, No 2, June 1982, pp.93–113
  16. Shampine, Lawrence Fred: “Type-Insensitive ODE Codes Based on Extrapolation Methods”, SIAM Journal on Scientific and Statistical Computing, Vol 4, No 4, December 1983, pp.635–644
  17. Shampine, Lawrence Fred: “Stiffness and the Automatic Selection of ODE Codes”, Journal of Computational Physics, Vol 54, No 1, April 1984, pp.74–86
  18. Stumpff, Karl (1895–1970)
  19. Stumpff, Karl: “Himmelsmechanik Band II — Das Dreikörperproblem”, VEB Deutscher Verlag der Wissenschaften, Berlin, 1965, 682 S.
  20. Wanner, Gerhard (*1942)