, 9 min read

Das Schalten im Programm PAI4

Das Programm PAI4, welches man beschrieben findet in dem Aufsatz von Bruder/Strehmel/Weiner (1988), 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.

1. Kurzbeschreibung des Programmes PAI4

1. Im Rahmen seiner Dissertation entwickelte Jürgen Bruder das schaltfähige Programm PAI4. Dieses Programm benutzt für den nicht-steifen Teil das Runge-Kutta Verfahren 3.ter Ordnung mit dem Parameterschema (Butcher-Schema)

$$ \begin{array}{c|ccc} \displaystyle\frac{1}{2} & \displaystyle\frac{1}{2} \\[6pt] 1 & -1 & 2 \\[6pt] \hline 1 & \displaystyle\frac{1}{6} & \displaystyle\frac{4}{6} & \displaystyle\frac{1}{6} \end{array} $$

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

$$ \begin{array}{c|ccc} \displaystyle\frac{1}{2} & \displaystyle\frac{1}{2}R_1^2 \\[6pt] 1 & -40R_1^3+78R_2^3 & \quad 41R_1^3-78R_2^3 \\[6pt] \hline & \displaystyle\frac{1}{2}R_1^4-R_2^4+R_3^4 & \displaystyle\frac{4}{3}R_1^4-2R_3^4 & \displaystyle-\frac{2}{3}R_1^4+R_2^4+R_3^4 \end{array} $$

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), 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:

  1. 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.
  2. 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.
  3. Nach jedem Schritt wird geprüft, ob ein Wechsel der Integrationsart möglich ist.
  4. Als Matrixnorm, für die Berechnung der für das Schalten so wichtigen Norm der Jacobimatrix, wird sowohl die 1-Norm $\Vert{}\cdot{}\Vert_1$ als auch die Maximumnorm $\Vert{}\cdot{}\Vert|_\infty$ benutzt, und es wird dann das Minimum dieser beiden Werte genommen.
  5. 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.

2. 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 (1983a) 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) beträgt die Rechenzeitersparnis jedoch schon 43%.

Programm PAI4 LSODA RKF4RW
Rechenzeit 72.07 77.24 127.38

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), aber auch 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.

  PAI4
CPU		 PAI4
$\varepsilon(t=1)$		 LSODA
CPU		 LSODA
$\varepsilon(t=1)$
G1 0.14 $7\cdot10^{-6}$ 0.22 $2\cdot10^{-6}$
G2 0.18 $2\cdot10^{-5}$ 0.22 $2\cdot10^{-5}$
G3 0.18 $2\cdot10^{-5}$ 0.22 $1\cdot10^{-5}$
G4 0.30 $1\cdot10^{-5}$ 0.36 $3\cdot10^{-4}$
G5 5.20 $3\cdot10^{-2}$ 0.66 $7\cdot10^{-5}$
G6 11.82 $8\cdot10^{-2}$ 1.04 $1\cdot10^{-4}$

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

Nicht-steife Probleme.

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

Steife Probleme.

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

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.

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. Gaffney, Parick W.: “A Performance Evaluation of Some FORTRAN Subroutines for the Solution of Stiff Oscillatory Ordinary Differential Equations”, ACM Transactions on Mathematical Software, Vo. 10, No. 1, March 1984, pp. 58-72
  3. Rentrop, Peter (*1948)
  4. 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.
  5. Rentrop, Peter: “Partitioned Runge-Kutta Methods with Stiffness Detection and Stepsize Control”, Numerische Mathematik, Vol 47, 1985, pp.545–564
  6. Tischer, Peter E.: “A New Order Selection for Ordinary Differential Equation Solvers”, SIAM Journal on Scientific and Statistical Computing, Vol 10, No 5, September 1989, pp.1024–1037