18th April 2026, 27 min read

Das Schalten im Programm RKF4RW

Das Programm RKF4RW, siehe Rentrop (1985), Rentrop, Peter, 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, P., Deuflhard (1983)2, siehe Deuflhard, P., Deuflhard (1985), passen die Ordnung ebenfalls den Genauigkeitsforderungen an, neben der Schrittweite, die von allen bisher erwähnten Programmen variiert wird. _{Shampine, Lawrence F.}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 $\mathopen|J\mathclose|_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

\ifabzwick\else \beginepigram_{Duschek, Adalbert}% 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. \author A. Duschek (1953) \endepigram

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 $ k = \max (n, 16) $ Schritten. Zu den Werten $n$ und 16 vermerkt Rentrop (1985)_{Rentrop, Peter}, 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 $\mathopen|W\mathclose|\infty$ während der Pivotisierung erhält. Diese Norm wird nun dazu benutzt, um festzustellen, ob $\mathopen|h\gamma\mathclose| \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ß

$$ \mathopen\|W\mathclose\|_\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. %

$$ \vbox{\offinterlineskip \halign {\fracstrut $\displaystyle{#}$\quad& \vrule#& \quad$\displaystyle{#}$\quad& \quad$\displaystyle{#}$\quad& \quad$\displaystyle{#}$\quad& \quad$\displaystyle{#}$\quad& \quad$\displaystyle{#}$\quad& \quad$\displaystyle{#}$\quad\cr 0 && & & & & &\cr 1\over4 && 1\over4 & & & & &\cr 3\over8 && 3\over32 & 9\over32 & & & &\cr 12\over13 && 1932\over2197 & 7296\over2197 & 0 & & &\cr 1 && 439\over216 & -8 & 3680\over513 & -{845\over4104} & &\cr 1\over2 &&-{8\over27} & 2 & -{3544\over2565} & 1859\over4104 & -{11\over40} &\cr \noalign{\hrule} y_1 && 25\over216 & 0 & 1408\over2565 & 2197\over4104 & -{1\over5} & 0\cr \hat y_1 && 16\over135 & 0 & 6656\over12825& 28561\over56430 & -{9\over50} & 2\over55\cr } % } $$

Hierbei werden die Werte 0, \frac1/4, \frac 3/8, \frac12/{13}, 1 und \frac1/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, Lawrence F.}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:

Das übliche software-interface wird zerstört, und es ist ungewohnt. Dies ist wichtig für größere Unterprogrammbibliotheken, wie z.B. DEPAC.
Bandstruktur kann verloren gehen und jede zusätzliche Gleichung lässt den Aufwand quadratisch ansteigen.
Ein lineares Problem kann nicht-linear werden, man vergl. beispielsweise das Problem D1. %tatsächlich stiff D1!
Neben $f_y$ muß zusätzlich $f_t$ berechnet werden, wenn der Benutzer die Jacobimatrix direkt bereitstellen will.
Ein Fehlerkriterium für die unabhängige Variable $t$ ist nicht leicht angebbar.
Die Norm $\mathopen|J\mathclose|_1 = \max(\mathopen|f_y\mathclose|_1, \mathopen|f_t\mathclose|_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

$$ \vbox{\offinterlineskip \halign{\strut #&\vrule#& \quad$#$\quad& \quad$#$\quad& \quad$#$\quad& \quad$#$\quad& \quad$#$\quad& \quad$#$\quad\cr &&\tilde\gamma_{21}& & & & &\cr &&\tilde\gamma_{31}&\tilde\gamma_{32}& & & &\cr &&\tilde\gamma_{41}&\tilde\gamma_{42}&\tilde\gamma_{43}& & &\cr &&\tilde\gamma_{51}&\tilde\gamma_{52}&\tilde\gamma_{53}&\tilde\gamma_{54}& &\cr &&\tilde\gamma_{61}&\tilde\gamma_{62}&\tilde\gamma_{63}&\tilde\gamma_{64}&\tilde\gamma_{65}&\cr \noalign{\hrule} }} $$

die von Rentrop berechneten Werte

$$ \vbox{\offinterlineskip \halign{\strut #&\vrule#& \quad$#$\quad& \quad$#$\quad& \quad$#$\quad& \quad$#$\quad& \quad$#$\quad& \quad$#$\quad\cr &&-3.7 & & & & &\cr &&-0.05395888512 &-0.4960411149 & & & &\cr &&3.592176012 &6.782987834 &-16.22024878 & & &\cr &&8.404898124 &15.19965650 &-37.17118671 &1.411793372 & &\cr &&3.022508944 &-8.091206674& 4.763352827 &0.6940545809 &1.15\cr \noalign{\hrule} }} $$

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: 5 Gleichungen zur Erreichung der Ordnung 4, 2 Gleichungen für das Erreichen der Ordnung 3 und schließlich 2 Gleichungen zur Sicherung der $A$-Stabilität. 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){Rentrop, Peter} entwickelte auch ein Verfahren, welches komponentenweise Steifheit behandelt. Hier sind dann andere Ableitungen erforderlich. Man vgl. hier auch das Buch von Hairer/Wanner/N\o rsett (1987), {Hairer, Ernst}{Wanner, Gerd}{N\o rsett, Syvert Paul}% S.276--285.

2. Bewertung des Programmes RKF4RW

\ifabzwick\else \beginepigram 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. \author K. Stumpff (1965)\silentreference{Stumpff (1965)} \endepigram

% Rentrop (1985)_{Rentrop, Peter} 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\null. 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)_{Rentrop, Peter} 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, Peter}Rentrop (1985) nicht gemacht.

Bei der Besprechung des Programmes PAI4 von {Strehmel, Karl}{Bruder, Jürgen}_{Weiner, Rüdiger}% 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){Rentrop,Peter} 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 \ifabzwick $L$-stabil. \else $L$-stabil, also nicht $A\infty^0$- bzw. $S_\infty^0$-stabil. \fi Die Ergebnisse wurden auf einer CYBER 175 (48 bit Mantissenlänge${}\approx14$ Dezimalstellen) gewonnen.

$$ \vbox{\offinterlineskip\halign { \strut# & \quad# & \vrule# & \quad# & \quad# & \quad# & \quad# & \quad# & \vrule# & \quad# & \quad# & \quad# & \quad# & \quad#\cr & Dgl. && A2 & B5 & C4 & D3 & E3 && NA2 & NB3 & NC5 & ND1 & NE3\cr \noalign{\hrule} RKF4RW: & CPU && 0.24 & 0.44 & 0.23 & 0.12 & 0.23 && 0.01 & 0.04 & 0.29 & 0.18 & 0.23\cr & nfe && 1165 & 2424 & 1370 & 780 & 1532 && 169 & 302 & 174 & 1270 & 1638\cr \noalign{\hrule} PRK4: & CPU && 0.20 & 0.22 & 0.44 & 0.11 & 0.17 && 0.01 & 0.04 & 0.60 & 0.26 & 0.25\cr & nfe && 840 & 1296 & 2839 & 671 & 1189 && 181 & 333 & 329 & 1684 & 1753\cr }} $$

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.

$$ \vbox{\offinterlineskip\halign{ \strut# & \quad# & \vrule# & \quad# & \quad# & \vrule# & \quad#\cr & && steif & nicht-steif && Summe\cr \noalign{\hrule} RKF4RW: & CPU && 1.26 & 0.75 && 2.01\cr & nfe && 7271 & 3553 && 10824\cr \noalign{\hrule} PRK4: & CPU && 1.14 & 1.16 && 2.30\cr & nfe && 6835 & 4280 && 11115\cr }} $$

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, Shampine, Lawrence F., und LSODA, siehe Petzold (1983a), Petzold, Linda Ruth, her kannte. Dieses abweichende Verhalten fällt umso mehr auf, da der zwei Jahre früher erschienene Artikel von _{Petzold, Linda Ruth}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, Lawrence F.}Shampine (1984)2, deuteten in diese Richtung.

Verglichen wurde nun von _{Rentrop, Peter}Rentrop (1985): (1) das Programm RKF4RW, welches nur im steifen Modus betrieben wurde, mit dem Programm GRK4A\null. 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:

$$ \vbox{\offinterlineskip\halign{ \strut# & \quad# & \quad# & \quad#\cr Dgl. & Programm & CPU & nfe\cr \noalign{\hrule} steif & GRK4A & 0.43 & 3694\cr & PRK4-steif & 0.88 & 5093\cr & RKF4RW-steif & 0.97 & 6881\cr \noalign{\hrule} nicht-steif & RKF4 & 0.28 & 1395\cr & PRK4-nicht-steif & 1.06 & 5686\cr }} $$

Bedenkt man, daß die Laufzeiten von “RKF4-steif\thinspace” 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 Geschwindigkeitseinbuß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: \beginquote \llap{“}The increase of computing time up to a factor 3 seems to be tolerable.” \endquote

date: "2026-04-18 21:00:00" title: "Das Schalten im Programm NTI/SIMPLE" MathJax: true TeXSource: "/home/klm/Dipl/t4grk.tex" lang: "de" categories: ["mathematics"] tags: ["Matrixpolynome", "Jordanketten"] index: true author: "Elmar Klausmeier"

$ \def\multisub#1#2{{\textstyle\mskip-3mu{\scriptstyle1\atop\scriptstyle#2_1}{\scriptstyle2\atop\scriptstyle#2_2}{\scriptstyle\ldots\atop\scriptstyle\ldots}{\scriptstyle#1\atop\scriptstyle#2_#1}}} \def\multisup#1#2{{\textstyle\mskip-3mu{\scriptstyle#2_1\atop\scriptstyle1}{\scriptstyle#2_2\atop\scriptstyle2}{\scriptstyle\ldots\atop\scriptstyle\ldots}{\scriptstyle#2_{#1}\atop\scriptstyle#1}}} \def\multisubsup#1#2#3{{\textstyle\mskip-3mu{\scriptstyle#3_1\atop\scriptstyle#2_1}{\scriptstyle#3_2\atop\scriptstyle#2_2}{\scriptstyle\ldots\atop\scriptstyle\ldots}{\scriptstyle#3_{#1}\atop\scriptstyle#2_{#1}}}} \def\diag{\mathop{\rm diag}} \def\tridiag{\mathop{\rm tridiag}} \def\col{\mathop{\rm col}} \def\row{\mathop{\rm row}} \def\dcol{\mathop{\rm col\vphantom {dg}}} \def\drow{\mathop{\rm row\vphantom {dg}}} \def\rank{\mathop{\rm rank}} \def\grad{\mathop{\rm grad}} \def\adj#1{#1^*} \def\iadj#1{#1^*} \def\tr{\mathop{\rm tr}} \def\mapright#1{\mathop{\longrightarrow}\limits^{#1}} \def\fracstrut{} $

\ifabzwick\else \beginepigram_{Runge, W.T.}% % Denn der Mensch denkt nicht gern logisch, dazu ist er nicht geschaffen. Von Natur aus denkt er assoziativ. Logisch zu denken muß er erst lernen, und es fällt ihm schwer. Mein Vater, der einmal in der Mathematischen Gesellschaft der ^{Berliner Universität} einen Vortrag hielt über ein numerisches Verfahren zur Berechnung höherer Differentialgleichungen, hatte sich zu dem Einwand zu äußern, das Verfahren sei doch maßlos umständlich. Er erwiderte, das sei wohl zugegeben, aber das komme daher, daß der ^{liebe Gott} dem Menschen den Verstand nicht zum Berechnen von Differentialgleichungen gegeben habe, sondern zur ^{Futtersuche}. \author W.T. Runge (1966) \endepigram

% Konzeptionell wird bei dem Programm NTI/SIMPLE_{NTI/SIMPLE} der gleiche Weg gewählt, der auch bei dem Programm TENDLER gewählt wurde: Man schaltet nicht zwischen verschiedenen _Integrations_verfahren hin und her, sondern nur zwischen zwei verschiedenen _Iterations_arten. Hier wird also zwischen modifiziertem Newton-Kantorovich-Verfahren und Picard Iteration hin und her gewechselt. Die nötige Information zum Schalten wird während der Iteration gewonnen, falls mit dem Newton-Kantorovich Verfahren gearbeitet wird und es wird keinerlei Information herangezogen beim Wechsel von nicht-steif nach steif. % % Im letzteren Falle wird einfach bei Divergenz der Fixpunktiteration zur Newton-Kantorovich Iteration übergewechselt. Das Programm TENDLER benutzt während der Picard Iteration sehr wohl weitere Informationen, die während dieser Iterationsart anfallen.

3. Die Schaltentscheidung im Programm NTI/SIMPLE

Das Programm NTI/SIMPLE von N\o rsett/Thomsen (1986) {N\o rsett, Syvert Paul}{Thomsen, Per Grove}% verwendet, wie das Programm TENDLER, beim Wechsel vom steifen Modus in den nicht-steifen Modus, den Wert $QRP$ einer Newton-Kantorovich-artigen Iteration zur Auflösung impliziter Differenzengleichungen. Der Wert $QRP$ wird im Programm NTI/SIMPLE allerdings “^{stiffness-ratio}” genannt wird, im Gegensatz zur üblichen Definition von

$$ \mu = {\max_i\mathopen|\lambda_i\mathclose| \over \min_k\mathopen|\lambda_k\mathclose|}, \qquad \min \mathopen|\lambda_k\mathclose| \ne 0. $$

Mit diesem Wert des Steifheitsmaßes steht das $\it QRP$ im Falle steifer Differentialgleichung zwar sinngemäß im Zusammenhang, jedoch besteht keine einfache, direkte funktionale Abhängigkeit.

Die Strategie für die Schaltung in dem Programm NTI/SIMPLE ist:

Wechsele von der Newton-Kantorovich zur Picard Iteration über, falls das $QRP$ kleiner als 2 ist.
Wechsele von der Fixpunktiteration zum Newton-Kantorovich Iterationsverfahren über, falls Divergenz bei einfacher Iteration gemeldet wurde.
Warte mindestens 10 Schritte bevor überhaupt an ein erneutes Schalten gedacht werden kann.

Die Rechtfertigung für den ersten Punkt ist der folgende Sachverhalt, der in gewisser Hinsicht eine notwendige Voraussetzung für Konvergenz der Fixpunktiteration darstellt: Beachtet man, daß $\hgJ < 1$ eine hinreichende Bedingung für Konvergenz der Picard Iteration ist, so erhält man mit der Forderung, daß das $QRP$ kleiner ist als 2, eine notwendige Bedingung dafür, daß die hinreichende Konvergenzbedingung erfüllt ist. Das Programm NTI/SIMPLE basiert auf SDIRK Formeln. Der endgültige im Programm NTI/SIMPLE benutzte Wert $QRP^*$ wird als Maximum aller $QRP$ genommen, gebildet über alle Stufen des Runge-Kutta Verfahrens.

In der Iterationsmatrix $W$ ist $h\gamma J$ enthalten. Wenn nun $\hgJ$ groß ist, so muß $\mathopen|W\mathclose|$ “nahe” bei $\hgJ$ liegen, und wenn $\hgJ$ nahe oder unterhalb von 1 ist, so muß $\mathopen|W\mathclose|$ in etwa in der Nähe von 2 sein. Dies ist eine grobe und einfache Motivation für die Benutzung des $QRP$ einer Newton-Kantorovich-artigen Iteration im Programm NTI/SIMPLE\null. Allerdings benutzt das Programm NTI/SIMPLE, wie das Programm TENDLER, die “verbotene” Rückrichtung. Während im Programm TENDLER hierfür versucht wird teilweise eine Kompensation zu erhalten durch weitere Informationen aus dem Schrittweiten- und Ordnungssteuerungssegment, so wird dies in dem Programm NTI/SIMPLE nicht getan. Allerdings sollte hinzugefügt werden, daß z.B. Minimierungsverfahren für Funktionen $F:\mathbb{R}^n\to\mathbb{R}$ häufig ebenfalls nur Punkte ansteuern mit $F'(x)=0$ und somit wird auch dort letztlich das Bestehen einer “verbotenen” Rückimplikation angenommen.

Divergenz während der Picard Iteration wird auch schon dann angenommen, wenn die Konvergenzrate über 0.8 liegt. Allerdings darf die Konvergenzrate ein einziges Mal doch über 0.8 liegen, da dann von N\o rsett/Thomsen (1986) angenommen wird, daß dies durch Rundungsfehler entstanden sein könnte. Dies deutet darauf hin, daß in dem Programm NTI/SIMPLE dem Rechnen nahe an der Grenze der Maschinengenauigkeit nicht die Aufmerksamkeit zuteil wurde, wie dies beispielsweise bei dem Programm LSODA der Fall ist. Die Konvergenzrate während der Picard Iteration ist eine untere Schranke für die Lipschitkonstante $L$ der Funktion $f$:

$$ {1\over \mathopen|h\gamma\mathclose|} {\mathopen\|y^{\nu+1}-y^\nu\mathclose\| \over \mathopen\|y^\nu-y^{\nu-1}\mathclose\|} \le L, $$

wie durch Subtraktion von

$$ y^{\nu+1} = h\gamma f(y^\nu)+\psi, \qquad y^\nu = h\gamma f(y^{\nu-1})+\psi $$

sofort folgt.

4. Bewertung des Programmes NTI/SIMPLE

1. Das Programm NTI/SIMPLE wurde an den Testgleichungen C1, C2, C3, C4, C5, D2, D5 und NA4, NA5, NB5, NE4, ND1, ND2, ND5 aus DETEST untersucht und mit der nicht schaltfähigen Version des Programmes verglichen. Gründe für das Auslassen von Testgleichungen aus DETEST und die benutzte Rechenanlage werden in N\o rsett/Thomsen (1986) {N\o rsett, Syvert Paul}{Thomsen, Per Grove}% nicht genannt.

Die Ergebnisse sind nun im einzelnen wie folgend. Hierbei werden die ermittelten Werte für die Genauigkeit $\varepsilon=10^{-2}$ nicht mit angegeben, da bei dieser niedrigen Genauigkeitsanforderung der globale Fehler nicht gut unter Kontrolle blieb. Man findet diese Werte bei N\o rsett/Thomsen (1986). {N\o rsett, Syvert Paul}{Thomsen, Per Grove}% Diese Erfahrungen mit den geringen Genauigkeitsansprüchen wurden auch mit dem Programm TENDLER gemacht. % % Zur Vermeidung solcher Effekte, wurde in dem Programm TENDLER dann soweit gegangen, daß diese geringe Genauigkeitsanforderung von einem Treiber als Fehleingabe abgefangen wird. Die Tabelle gilt für die Genauigkeitsanforderung $\varepsilon=10^{-4}$.

$$ \vbox{\offinterlineskip\halign{ \strut# && \quad#\cr Dgl. & \#fix & nfe & nje & \#$LU$ & $Ux=c$ & CPU & code\cr \noalign{\hrule} C1 & 41 & 460 & 3 & 16 & 195 & 0.85 & neu\cr & 0 & 421 & 1 & 25 & 416 & 1.03 & alt\cr \xx C2 & 35 & 427 & 1 & 18 & 210 & 0.85 & neu\cr & 0 & 407 & 1 & 26 & 402 & 1.01 & alt\cr \xx C3 & 38 & 412 & 1 & 18 & 193 & 0.83 & neu\cr & 0 & 389 & 1 & 26 & 384 & 0.97 & alt\cr \xx C4 & 229 & 3506 & 11 & 58 & 1916 & 6.66 & neu\cr & 0 & 2992 & 4 & 63 & 2972 & 7.20 & alt\cr \xx D2 & 14 & 608 & 2 & 14 & 508 & 1.02 & neu\cr & 0 & 583 & 3 & 21 & 571 & 1.05 & alt\cr \xx D5 & 0 & 114 & 4 & 20 & 99 & 0.16 & neu\cr & 0 & 111 & 4 & 20 & 99 & 0.15 & alt\cr }} $$

Deutlich wird der, wenn auch geringe, Rechenzeitgewinn ggü. dem nicht schaltfähigen Programm. Bei beiden Programmen war der beobachtete globale Fehler am Endpunkt größenordnungsmässig gleich, sodaß die Werte gut vergleichbar sind. Stellenweise auffällig ist die höhere Anzahl an Jacobimatrixauswertungen des schaltfähigen Programmes. Ebenso sind die benötigten Funktionsauswertungen i.d.R. höher, jedoch sind dafür die Rücksubstitutionen (in der Spalte $Ux=c$) drastisch zurückgegangen. Die Anzahl der $LU$-Zerlegungen ist ebenfalls vermindert worden. Wenn man beachtet, daß gerade die Kosten für die lineare Algebra in Programmen zur Lösung steifer Differentialgleichungen einen wichtigen Beitrag leisten, so ist dieser Gesichtspunkt hier von besonderem Interesse. In den Rechenzeiten schlägt sich dieser Gewinn nicht in der erwarteten Weise nieder, da es sich bei den ausgewählten Testdifferentialgleichungen um im wesentlichen kleindimensionale Probleme handelt: alle Dimensionen liegen bei, oder unterhalb von 4.

Bei den Gleichungen NA4, NA5, NB5, NE4, ND1, ND2, ND5 wird kein Vergleich in N\o rsett/Thomsen (1986) angegeben, wodurch man die Kosten der Schaltung und der Schaltentscheidung erkennen könnte. Deutlich wird allerdings, daß bei diesen nicht-steifen Differentialgleichungen tatsächlich keine einzige Jacobimatrixauswertung und keine $LU$-Zerlegung benötigt wird, mit Ausnahme bei der Gleichung NA4. Diese Gleichungen werden also korrekt als nicht-steife Gleichungen erkannt. Diese Ergebnisse gelten nicht für die Genauigkeitsanforderung von $10^{-2}$, für die schon mehrfach darauf hingewiesen wurde, daß man einer solchen Genauigkeitsanforderung des Benutzers zumindestens vorsichtig gegenüber stehen sollte. Bei dieser niedrigen Genauigkeit werden häufig $LU$-Zerlegungen und Jacobimatrixauswertungen benötigt, weiterhin ist der globale Fehler inakzeptabel.

2. Zusätzlich zu den obigen Testgleichungen wurden Vergleichsmessungen an der van der Polschen, zweidimensionalen Differentialgleichung durchgeführt. Die hier benutzte van der Polsche Gleichung lautet_{Pol, van der}_{van der Pol}

$$ \eqalign{ \dot y_1 &= y_2,\cr \dot y_2 &= 100(1-y_1^2)y_2 - y_1,\cr } \qquad\hbox{mit}\quad y(0)=\pmatrix{2\cr 0\cr}\qquad\hbox{für}\quad t\in[0,100]. $$

Anhand dieser Gleichung ergab sich nun für eine Genauigkeitsanforderung von $\varepsilon=10^{-4}$

$$ \vbox{\offinterlineskip\halign{ \strut# & \vrule# && \quad#\cr code && \#steps & \#fix & $\|y_N-y(100)\|$ & nfe & nje & \#$LU$ & $Ux=c$ & CPU\cr \noalign{\hrule} neu && 187 & 150 & $2.5e-3$ & 1172 & 4 & 36 & 261 & 1.15\cr alt && 185 & 0 & $4.9e-4$ & 1172 & 4 & 71 & 1160 & 1.49\cr }} $$

Man erkennt den leichten Rechenzeitgewinn bei zumindestens vergleichbaren globalen Fehlern am Endzeitpunkt. Die Erklärung hierfür ist, daß bei den Spitzen der Lösungen der van der Polschen Differentialgleichung, der wesentlich effizientere Modus der Fixpunktiteration gewählt wurde. Es ist Zufall, daß die Anzahl der Funktionsauswertungen gleich geblieben ist. Die hohe Anzahl an Fixpunktiterationen deutet darauf hin, daß erfolgreich hin und her geschaltet wurde. Dies ist ein Hinweis dafür, daß mit Hilfe des $QRP$ eine Möglichkeit gegeben ist im Falle von Runge-Kutta Formeln automatisch zu schalten.

date: "2026-04-18 21:00:00" title: "Das Schalten im Programm PAI4" MathJax: true TeXSource: "/home/klm/Dipl/t4grk.tex" lang: "de" categories: ["mathematics"] tags: ["Matrixpolynome", "Jordanketten"] index: true author: "Elmar Klausmeier"

Das Programm PAI4, welches man beschrieben findet in dem Aufsatz von Bruder/Strehmel/Weiner (1988), {Bruder, Jürgen}{Strehmel, Karl}_{Weiner, Rüdiger}% verwendet ebenso wie das Programm RKF4RW, neu entwickelte Formeln. Getestet wird also nicht nur die Form einer programmmässigen Realisierung eines Verfahrens, sondern zugleich wird auch eine neue Formel auf ihre praktische Eignung hin überprüft. Beurteilt man daher das Programm PAI4, so beurteilt man nicht notwendig als solches die neuen Formeln, sondern die Formeln im Zusammenspiel mit dem Rest des Programmes. Die Art der Schaltung birgt in sich nun keine allzu neuen Ideen, welche nicht schon vorher beschrieben wurden. Wiederum zieht man zur Schaltentscheidung die Norm der Jacobimatrix maßgeblich heran. Dies hat auch seinen Grund darin, daß die Jacobimatrix ohnehin für das Verfahren benötigt wird. Zum Schluß werden die drei Programme PAI4, LSODA und RKF4RW einander gegenübergestellt. Erneut erkennt man, daß das zuerst vorgestellte Programm LSODA sehr erfolgreich ist.

Von seiner Grundkonzeption her ist das Programm PAI4 nur für sehr geringe Genauigkeitsanforderungen gedacht, also Anforderungen im Bereich $\varepsilon=10^{-2}$ bis $\varepsilon=10^{-5}$. Für Anwendungen im Bereich der Simulation ist dies durchaus ausreichend: Die Differentialgleichung selber beschreibt die Wirklichkeit nicht genauer, und die zahlreich einfliessenden Parameter sind ebenfalls nicht hochgenau. Für ein schaltfähiges Programm kann jedoch dieser Genauigkeitsbereich auch vielfach zu einschränkend sein. Nicht-steife Differentialgleichungen verlangen nicht selten genaue Rechnung. Der Grund für die Benutzung eines schaltfähigen Programmes liegt nun u.a. maßgeblich darin, daß man sich vor plötzlich auftauchender Steifheit wappnen will. Es zeigt sich, daß das Programm PAI4 bei den Testgleichungen aus STDTST und NSDTST dem Programm LSODA, bei sehr geringen Genauigkeitsanforderungen ebenbürtig ist, falls auch die Dimensionen der Differentialgleichungen klein sind. Bei Gleichungen mit stark variierender Jacobimatrix kann das Programm PAI4 versagen.

5. Kurzbeschreibung des Programmes PAI4

1. Im Rahmen seiner Dissertation entwickelte J. Bruder_{Bruder, Jürgen} das schaltfähige Programm PAI4. % Dieses Programm benutzt für den nicht-steifen Teil das Runge-Kutta Verfahren 3.ter Ordnung mit dem Parameterschema

$$ \vbox{\offinterlineskip \halign { \fracstrut$\displaystyle{#}$\quad & \vrule# && \quad $\displaystyle{#}$ \quad\cr 1\over2 && 1\over2\cr 1 && -1 & 2\cr \noalign{\hrule} 1 && 1\over6 & 4\over6 & 1\over6\cr }} $$

und für den steifen Teil ist vorgesehen das Verfahren dritter Ordnung zu

$$ \vbox{\offinterlineskip \halign { \fracstrut$\displaystyle{#}$\quad & \vrule# && \quad $\displaystyle{#}$ \quad\cr 1\over2 && {1\over2}R_1^2\cr 1 && -40R_1^3+78R_2^3 & 41R_1^3-78R_2^3\cr \noalign{\hrule} && {1\over2}R_1^4-R_2^4+R_3^4 & {4\over3}R_1^4-2R_3^4 & -{2\over3}R_1^4+R_2^4+R_3^4\cr }} $$

mit

$$ \eqalign{ R_0^i(z) &= {1+(1-2\gamma_i)z+({1\over2}-2\gamma_i+\gamma_i^2)z^2 \over (1-\gamma_iz)^2}, \qquad\hbox{für}\quad i\in\{2,3\},\cr R_0^4(z) &= {1+(1-3\gamma)z+({1\over2}-3\gamma+\gamma^2)z^2 \over (1-\gamma z)^3}, \qquad\hbox{mit $\gamma = 0.435 866 521 508 4592$} \cr } $$

und

$$ \gamma_i = {\gamma\over c_i}, \qquad\hbox{für}\quad i\in\{2,3\}. $$

Die Parameter $c_i$ treten in den (nichtlinearen) Ordnungsbedingungen auf. Die obige Wahl der $\gamma_i$ bietet Ersparnisse bei den Rechenoperationen bei der anschließenden Implementierung dieses Verfahrens. % Dabei ist das angegebene Parameterschema die Charakterisierung für das ARK-Verfahren\ifabzwick\else und darf nicht mit dem gewöhnlichen Parameterschema verwechselt werden\fi, man siehe hierzu Bruder/Strehmel/Weiner (1988), {Bruder, Jürgen}{Strehmel, Karl}_{Weiner, Rüdiger}% für das auch auf die Herleitung und Begründung dieser Verfahren verwiesen sei. Die Ableitung von ARK-Verfahren höherer Ordnung, wird von den Autoren als recht prohibitiv bezeichnet. Das Verfahren für den steifen Teil ist $L$-stabil.

Wie schon in dem Programm LSODA wird vollständig zwischen zwei _Integrations_arten hin und her gewechselt. Es wird also nicht nur zwischen den _Iterations_arten geschaltet, wie dies in dem Programm TENDLER geschieht.

2. % Der eigentliche Schalttest besteht nun darin zu überprüfen, ob gilt

$$ \|hf_y\| \le c\cdot r . $$

Die Konstante $c$ dient dazu, die Ungenauigkeiten in der Benutzung von $|f_y|$ anstatt $\rho(f_y)$ wegzudämpfen. Diese Konstante $c$ wurde als $c=1.1>1$ gewählt. Der zweite Faktor $r$ ist der Radius des Stabilitätsgebietes des verwendeten Verfahrens.

Im weiteren wird dann nach den folgenden Regeln verfahren:

Im nicht-steifen Falle, also bei Benutzung des expliziten Runge-Kutta Verfahrens dritter Ordnung, wird die Jacobimatrix neuberechnet, falls sich die Schrittweite insgesamt seit der letzten Berechnung der Jacobimatrix entweder gefünftelt oder aber verachtfacht hat. Auf jeden Fall wird wird die Jacobimatrix nach 16 Schritten erneut ausgewertet.
Im steifen Falle, also Benutzung des ARK-Verfahrens, wird die Jacobimatrix aufgefrischt, falls sich die Schrittweite verkleinert hat, oder in jedem Falle alle 8 Schritte.
Nach jedem Schritt wird geprüft, ob ein Wechsel der Integrationsart möglich ist.
Als Matrixnorm, für die Berechnung der für das Schalten so wichtigen Norm der Jacobimatrix, wird sowohl die 1-Norm $|{}\cdot{}|1$ als auch die Maximumnorm $|{}\cdot{}|\infty$ benutzt, und es wird dann das Minimum dieser beiden Werte genommen.
Die Schrittweite wird in dem Programm PAI4 über Richardson-Extrapolation gesteuert. Die Schrittweite darf sich höchstens verfünffachen oder fünfteln. %Ein Schritt wird dann mit halber Schrittweite durchgeführt.

6. Bewertung und Vergleich des Programmes PAI4

Man beachte, daß die Verfahren im Programm PAI4 nur für geringe Genauigkeitsanforderungen konzipiert wurden. Die Rechenzeitersparnis ggü. dem Programm LSODA von Petzold ist vernachlässigbar: lediglich 6%. D.h., aufgrund einer solchen, möglichen Rechenzeitersparnis, wird ein Benutzer des Programmes LSODA wohl nicht auf das Programm PAI4 überwechseln, da die nötigen Änderungen im aufrufenden Hauptprogramm und das Austesten, diese Ersparnisse i.d.R. nicht aufwiegen. Diese Daten und die nachfolgende Tabelle beziehen sich auf eine angeforderte Genauigkeit von $\varepsilon=10^{-4}$, d.h. für eine vergleichsweise geringe Genauigkeitsanforderung. Gegenüber dem Programm RKF4RW von Rentrop (1985)_{Rentrop, Peter} beträgt die Rechenzeitersparnis jedoch schon 43%.

$$ \vbox{\offinterlineskip \halign { \strut\quad#\quad & \vrule# & \quad#\quad &\vrule#& % PAI4 \quad#\quad &\vrule#& % LSODA \quad#\quad \cr % RKF4RW Programm && PAI4 && LSODA && RKF4RW\cr \noalign{\hrule} Rechenzeit && $72.07$ && 77.24 && 127.38\cr }} $$

Getestet und einander gegenübergestellt wurden die drei Programmme PAI4, LSODA und RKF4RW an insgesamt 61 Differentialgleichungen, nämlich die 25 Gleichungen aus NSDTST (NA1--NE5), die 30 Gleichungen aus STDTST (A1--F5) und weitere 6 Gleichungen (G1--G6). Verglichen wurden natürlich, entsprechend der Konzeption der Formeln, nur die Rechenzeiten bei den Genauigkeiten $10^{-2}$, $10^{-4}$ und $10^{-6}$. % Die “exakten” Lösungen wurde mit mit dem Programm LSODA berechnet, bei einer Genauigkeitsanforderung von $10^{-8}$. Auffällig ist, daß LSODA bei der sehr geringen Genauigkeitsanforderung von $10^{-2}$ öfters länger rechnet, als bei einer Genauigkeitsanforderung von $10^{-4}$. Dieses Verhalten wurde ähnlich bei den zwei Programmen TENDLER und STINT beobachtet. Für solche niedrigen Genauigkeiten treten weitere Effekte auf. Bruder/Strehmel/Weiner (1988) stellen ihre Ergebnisse auch graphisch dar.

Es folgen nun die Gesamtrechenzeiten der drei Programme PAI4, LSODA und RKF4RW für die Genauigkeitsanforderung von $\varepsilon=10^{-4}$, siehe folgende Doppeltabelle. Die globalen Fehler sind annähernd vergleichbar. Beachtenswert ist hier die doppelte Rechenzeit von RKF4RW ggü. LSODA und die Verdreifachung ggü. dem schaltfähigen Programm PAI4. Obwohl das vom Programm RKF4RW verwendete Verfahren $A$-stabil ist, jedoch nicht $L$-stabil, ergeben sich für das Problem B5 ungewöhnlich lange Rechenzeiten. Die hohen Rechenzeiten des Programmes LSODA sind bekanntlich auf die nicht genügend großen Widlund-Winkel der BDF$i$, für $i>2$, zurückzuführen. Wüßte man von einer Differentialgleichung im voraus um die Lage und Größe der Eigenwerte der Jacobimatrix, so könnte man bei LSODA die Höchstordnung auf 2 begrenzen und damit auch das Problem B5 effizient lösen, wie Gaffney (1984)_{Gaffney, Patrick W.}, aber auch _{Tischer, Peter E.}Tischer (1989), deutlich machen.

Bei den Problemen G5 und G6 versagt das Programm PAI4: Rechenzeit und Genauigkeit sind nicht akzeptabel. Dies geschieht ohne sonstige Benutzerinformation oder Warnung. Benutzt wurde eine Genauigkeitsanforderung von $\varepsilon=10^{-4}$. Die Spalte $\varepsilon(t=1)$ gibt den globalen Fehler am Endpunkt $t=1$ an.

$$ \vbox{\offinterlineskip\halign{\strut# && \qquad# & \quad$#$\cr & \multispan2\qquad PAI4 & \multispan2\qquad LSODA\cr & CPU & \varepsilon(t=1) & CPU & \varepsilon(t=1)\cr \noalign{\hrule} G1 & 0.14 & 7\x6 & 0.22 & 2\x6\expstrut\cr G2 & 0.18 & 2\x5 & 0.22 & 2\x5\cr G3 & 0.18 & 2\x5 & 0.22 & 1\x5\cr G4 & 0.30 & 1\x5 & 0.36 & 3\x4\cr G5 & 5.20 & 3\x2 & 0.66 & 7\x5\cr G6 & 11.82 & 8\x2 & 1.04 & 1\x4\cr }} $$

In den nachfolgenden Tabellen bezeichnen $\sigma_{n-1}$ und $\sigma_n$ wie üblich die Standardabweichung und Varianz, gemäß der Formeln

$$ \sigma_{n-1} = \sqrt{{1\over n-1}\sum_{1\le i\le n}(x_i-\overline{x})^2} = \sqrt{{1\over n-1} \left\{\sum_{1\le i\le n} x_i^2 - {1\over n}\left(\sum_{1\le i\le n} x_i\right)^2 \right\}} $$

und

$$ \sigma_n = \sqrt{{1\over n}\sum_{1\le i\le n} (x_i-\overline{x})^2} = \sqrt{{1\over n}\left\{\sum_{1\le i\le n} x_i^2 - {1\over n}\left(\sum_{1\le i\le n} x_i\right)^2 \right\}}. $$

% % Für die nicht-steifen Testdifferentialgleichungen aus DETEST ergab sich nun das folgende, vergleichende Bild zwischen den drei schaltfähigen Programmen PAI4, LSODA und RKF4RW, wiederum bei einer Genauigkeitsanforderung von $\varepsilon=10^{-4}$. Erwähnenswert ist die geringe Gesamtrechenzeit des linearen Mehrschritt-Verfahrens LSODA\null. Zum Teil liegt dies an den “Ausreissern” bei den Runge-Kutta Verfahren. Hier wird deutlich, wie wichtig Robustheit ist. Auffällig ist hier, daß LSODA bei den relativ großdimensionalen Testgleichungen, wie NC4 und NC5, mit den jeweiligen Dimensionen von 51 und 30, deutliche Vorteile zeigt. Bei den kleindimensionalen Testgleichungen ist die Rechenzeit ohnehin nicht von tragender Bedeutung, obwohl dies natürlich von der verwendeten Rechenanlage abhängig ist. Auf kleinen Kompaktrechnern können diese kleindimensionalen Differentialgleichungen durchaus längere Rechenzeit beanspruchen.

\setbox0=\vbox{\offinterlineskip\halign { \strut\quad# && \quad$#$\cr Problem & \hbox{PAI4} & \hbox{LSODA} & \hbox{RKF4RW\quad}\cr \noalign{\hrule} NA1 & 0.06 & 0.24 & 0.14\cr NA2 & 0.06 & 0.22 & 0.10\cr NA3 & 0.36 & 0.76 & 0.62\cr NA4 & 0.06 & 0.14 & 0.12\cr NA5 & 0.04 & 0.18 & 0.16\cr \bisselPlatz NB1 & 0.44 & 1.24 & 0.92\cr NB2 & 0.20 & 0.38 & 0.36\cr NB3 & 0.16 & 0.36 & 0.28\cr NB4 & 0.52 & 0.86 & 1.20\cr NB5 & 0.34 & 0.72 & 0.76\cr \bisselPlatz NC1 & 0.54 & 0.62 & 0.88\cr NC2 & 1.90 & 1.14 & 2.28\cr NC3 & 0.80 & 0.76 & 1.06\cr NC4 & 9.56 & 3.72 & 4.46\cr NC5 & 5.16 & 1.38 & 4.66\cr \bisselPlatz ND1 & 0.98 & 0.82 & 1.50\cr ND2 & 1.24 & 1.14 & 1.54\cr ND3 & 1.44 & 1.50 & 1.86\cr ND4 & 2.24 & 2.08 & 2.62\cr ND5 & 4.90 & 3.20 & 4.16\cr \bisselPlatz NE1 & 0.20 & 0.40 & 0.44\cr NE2 & 0.52 & 1.74 & 1.00\cr NE3 & 0.70 & 1.04 & 1.62\cr NE4 & 0.06 & 0.08 & 0.14\cr NE5 & 0.12 & 0.14 & 0.46\cr \noalign{\hrule} Summe & 32.60 & 24.86 & 33.34\cr Mittel & 1.30 & 0.99 & 1.33\cr $\sigma_{n-1}$ & 2.19 & 0.92 & 1.36\cr $\sigma_n$ & 2.15 & 0.90 & 1.33\cr }}

\setbox1=\vbox{\offinterlineskip\halign { \strut\quad# && \quad$#$\cr Problem & \hbox{PAI4} & \hbox{LSODA} & \hbox{RKF4RW\quad}\cr \noalign{\hrule} A1 & 0.62 & 0.74 & 0.92\cr A2 & 2.14 & 1.78 & 3.46\cr A3 & 1.06 & 1.42 & 2.48\cr A4 & 4.10 & 3.08 & 12.60\cr \bisselPlatz B1 & 3.84 & 5.14 & 9.88\cr B2 & 0.74 & 0.96 & 1.02\cr B3 & 0.78 & 1.08 & 1.04\cr B4 & 1.02 & 1.36 & 1.62\cr B5 & 1.70 & 3.12 & 4.80\cr \bisselPlatz C1 & 0.76 & 1.00 & 0.96\cr C2 & 0.70 & 1.00 & 1.28\cr C3 & 0.72 & 1.02 & 1.92\cr C4 & 1.56 & 1.14 & 3.22\cr C5 & 2.04 & 1.66 & 4.00\cr \bisselPlatz D1 & 1.10 & 0.88 & 4.58\cr D2 & 0.66 & 0.78 & 1.60\cr D3 & 0.72 & 0.92 & 2.10\cr D4 & 0.22 & 0.30 & 0.58\cr D5 & 0.26 & 0.42 & 0.42\cr D6 & 0.26 & 0.52 & 0.36\cr \bisselPlatz E1 & 0.36 & 0.46 & 0.48\cr E2 & 0.48 & 1.10 & 0.98\cr E3 & 0.62 & 0.82 & 3.10\cr E4 & 2.10 & 2.50 & 8.06\cr E5 & 0.36 & 0.34 & 1.06\cr \bisselPlatz F1 & 0.36 & 0.56 & 0.78\cr F2 & 0.36 & 0.60 & 0.78\cr F3 & 3.96 & 4.78 & 12.46\cr F4 & 2.72 & 2.00 & 3.14\cr F5 & 3.10 & 10.90 & 4.36\cr \noalign{\hrule} Summe & 39.42 & 52.38 & 94.04\cr Mittel & 1.31 & 1.75 & 3.13\cr $\sigma_{n-1}$ & 1.17 & 2.11 & 3.38\cr $\sigma_n$ & 1.15 & 2.10 & 3.33\cr }} \twounequalboxes

Insgesamt gesehen scheint es hervorhebenswert, daß das Programm LSODA als Vertreter für ein Programm basierend auf linearen Mehrschrittverfahren, bei den nicht-steifen Testdifferentialgleichungen vergleichsweise geringere Rechenzeiten aufweist, als diejenigen Programme, die Runge-Kutta-Verfahren verwenden.