, 25 min read

TENDLER: 5. Die Fehlerkontrolle

Fortsetzung der TENDLER Programmbeschreibung.

  1. TENDLER: 1. Grobaufbau und prinzipielle Überlegungen
  2. TENDLER: 2. Benutzung des Programmes
  3. TENDLER: 3. Der Prädiktor
  4. TENDLER: 4. Die Korrektoriteration

Inhalt.

However, we can say that if you are using a 12-line solver for differential equations on anything bigger than a hand calculator, you should consider using one of the cited packages instead. Recently, we have noted that there is commercially available “software” for differential equations with no error control, a user-specified fixed step size, no warning messages, and so on. We strongly advise against using such programs, even on a personal computer. The reasons are straightforward. For all but trivial problems, such programs cannot be sufficient reliable for accurate computational results.

G.D. Byrne, A.C. Hindmarsh (1987)

In diesem Segment wird nun versucht zu überprüfen, ob die berechneten Näherungen die geforderte Genauigkeitsbedingung, die der Programmbenutzer angibt, erfüllen. Es wird an dieser Stelle deutlich werden, warum es günstig ist, eine Darstellung der berechneten Lösungen in Form von rückwärtsgenommenen Differenzen zu wählen. Sie dienen als Näherungen für entsprechend hohe Ableitungen, welche seinerseits zur Schätzung des lokalen Fehlers herangezogen werden.

Es folgt also die Beschreibung der Fehlerkontrolle. Die Fehlerkontrolle wird nach der Korrektoriteration ausgeführt. Dabei werden die verwendete Heuristik und weitere Einzelheiten angegeben. Die Fehlerkontrolle ist in zwei Phasen aufgeteilt. Jede Phase für sich fällt Entscheidungen über das Fortrechnen. Die Fehlerkontrolle kann entweder nicht bestanden oder passiert werden. In Abhängigkeit dieser Entscheidung müssen weitere Entscheidungen getroffen werden, wie weiterverfahren werden soll.

Im Falle, daß die Fehlerkontrolle nicht bestanden wurde, stehen einem zwei massgebliche Parameter zur Verfügung. Dies ist zum einen die Ordnung $p$ und zum anderen die Schrittweite $h$.

Zuerst wird versucht durch Verändern der Schrittweite dem Fehlertestversagen zu begegnen. Erst nach mehrmaligen Nichtbestehen wird die Ordnung ebenfalls verändert.

Bei einem Programm basierend auf zyklischen Formeln, sind nach einem Fehlertestversagen mehrere Fortsetzungsmöglichkeiten denkbar. Zum einen könnte nur eine einzige Stufe wiederholt werden, oder es könnte der gesamte Zyklus wiederholt werden. Beide Möglichkeiten haben ihre Vor- und Nachteile. Im Programm TENDLER wird die letztere Möglichkeit gewählt. Statistiken belegen, daß die Furcht vor dem Wiederholen ganzer Zyklen nicht gerechtfertigt ist. Erfahrungsmaterial spricht nicht eindeutig gegen vielstufige Zyklen, allerdings auch nicht eindeutig für sie.

Ein paar Worte zum Begriff “Fehlerkontrolle” und “Konvergenztest”. Die weiter unten zu beschreibende “Fehlerkontrolle” kann nicht wirklich globale Fehler kontrollieren, zumindestens nicht in jedem Falle. Sie kann also nicht Fehler, die entstehen durch Ersetzung des kontinuierlichen Problems (Differentialgleichungsanfangswertproblem) durch ein diskretes Problem (implizite Differenzengleichung) genau quantifizieren. Etwas korrekter wäre von einem Versuch der Schätzung zu sprechen. Den Ersetzungsfehlern überlagern sich weitere Fehler, wie z.B. Rundungsfehler, Iterationsfehler, Startfehler, etc. Mit dem Begriff Konvergenztest verhält es sich ähnlich. Der Konvergenztest ist kein Test im Sinne eines Majorantenkriteriums, Verdichtungskriteriums oder dergleichen, sondern eine Entscheidung über den Abbruch der Iteration, also die ausdrückliche Finitisierung eines unendlichen Vorganges.

Bibliographisch.

  1. George Dennis Byrne
  2. Alan C. Hindmarsh (*1942)
  3. Byrne, George D. und Hindmarsh, Alan C.: “Stiff ODE Solvers: A Review of Current and Coming Attractions”, Journal of Computational Physics, Vol 70, No 1, May 1987, pp.5—62

1. Abbruch- und Diskretisierungsfehler

Die Definition des lokalen Fehlers sind in Hairer/Wanner/Nørsett (1987), Albrecht (1985) und in Reimer (1982) jeweils verschieden, obwohl natürlich letztendlich sinnentsprechend. Gelegentlich wird auch zwischen den beiden Begriffen Abbruch- und Diskretisierungsfehler nicht unterschieden, oder nur einer der beiden eingeführt. Die hier beiden festgelegten Begriffe sind konzeptionell verschieden, führen jedoch unter milden Voraussetzungen zu gleichen Ergebnissen.

1. Sei

$$ \varphi\colon\mathbb{R}^{n\ell}\times(\mathbb{R}^{n\ell})^k \to \mathbb{R}^{n\ell} $$

eine vorerst nicht weiter eingeschränkte Vektorfunktion mit $k+1$ Vektorargumenten. $n$ erinnere hier an eine von außen kommende Dimension, und $\ell$ erinnere an eine ebenfalls von außen induzierte Stufenzahl. Betrachtet werde die Differenzengleichung

$$ \varphi(u_{\nu+1},u_\nu,\ldots,u_{\nu-k+1}) = 0, \qquad \nu=k,k+1,\ldots, $$

wobei die $u_0,\ldots,u_k$ vorgegeben gedacht seien. $u_k$ berechnet sich als Lösung von $\varphi(u_k,u_{k-1},\ldots,u_0)=0$, $u_{k+1}$ berechnet sich als Lösung von $\varphi(u_{k+1},u_k,\ldots,u_1)=0$, u.s.w. Die eindeutige Auflösbarkeit sei von $\varphi$ im folgenden stets, ohne ausdrückliche Betonung, vorausgesetzt. Dies ist beispielsweise erfüllt, wenn $\varphi$ bzgl. der ersten Vektorkomponente kontrahierend ist; ein häufig auftretender Spezialfall.

2. Beispiele: (1) Implizites Eulerverfahren. Hier ist

$$ \varphi(u_{m+1},u_m) = u_{m+1} - u_m - h f(u_{m+1}). $$

(2) Zyklische lineare Mehrschrittverfahren. Hier wäre

$$ \varphi(u_{m+1},u_m,\ldots,u_{m-k+1}) = \sum_{i=0}^k A_i u_{m+1+i-k} - h B_i \dot u_{m+1+i-k}. $$

3. Sei $y$ eine beliebige Kurve des $\mathbb{R}^{n\ell}$. Die Werte $u_m$ sollen $y(t_m)$ möglichst gut annähern. Die Differenz von $u_m$ und $y(t_m)$ ist der globale Fehler, oder, wenn man betonen will, daß man $y$ diskret annähert, der globale Diskretisierungsfehler, wobei man dann $\varphi$ als Diskretisierungsverfahren oder Differenzenverfahren bezeichnet. Der globale Fehler soll klein sein. Diese Forderungen führt zu Bedingungen an $\varphi$. Jedoch ist dieser globale Fehler schwer bestimmbar, mindestens so schwer bestimmbar wie die Kurve $y(t)$ selber. Einfacher zu bestimmen sind lokaler Abbruchfehler und lokaler Diskretisierungsfehler, zumindestens Anschmiegungen entsprechend hoher Ordnungen.

4. Definition: Sei $y(t)$ eine beliebige aber feste Kurve, z.B. Lösung des Anfansgswertproblems $\dot y=f(t,y)$, $y(t_0)=y_0$.

(1) Sei $u_{m+1}$ Lösung von $\varphi(u_{m+1},u_m,\ldots,u_{m-k+1})=0$ ($\forall m\ge k$). Dann heißt

$$ d(t_{m+1},y,h) := u_{m+1} - y(t_{m+1}) $$

globaler Diskretisierungsfehler an der Stelle $t_{m+1}$ für das Differenzenverfahren $\varphi$.

(2) Das Residuum

$$ r(t_{m+1},y,h) := \varphi(y(t_{m+1}),y(t_m),\ldots,y(t_{m-k+1})) \qquad (m\ge k) $$

heißt lokaler Abbruchfehler an der Stelle $t_{m+1}$ für das Differenzenverfahren $\varphi$.

(3) Sei $u_{m+1}$ Lösung von $\varphi(u_{m+1},y(t_m),\ldots,y(t_{m-k+1}))=0$. Dann heißt die Differenz

$$ s(t_{m+1},y,h) := u_{m+1} - y(t_{m+1}) $$

lokaler Diskretisierungsfehler an der Stelle $t_{m+1}$ für das Differenzenverfahren $\varphi$.

5. Bemerkung: Beim lokalen Abbruchfehler werden in alle $k+1$ Argumente von $\varphi$, die exakten Werte $y(t_i)$ eingesetzt und das entsprechende Residuum heißt lokaler Abbruchfehler. Beim lokalen Diskretisierungsfehler nimmt man für die letzten $k$ Argumente von $\varphi$ die exakten Kurvenwerte und der sich dann ergebende Lösungswert der Differenzengleichung wird mit dem exakten Wert der Kurve verglichen. Der lokale Diskretisierungsfehler ist ein spezieller globaler Diskretisierungsfehler, nämlich derjenige, der zusätzlich die Bedingung $\varphi(u_{m+1},y(t_m),\ldots,y(t_{m-k+1}))=0$ erfüllt, also ein globaler Fehler unter Nebenbedingungen.

Bei einfachen Zusatzvoraussetzungen an $\varphi$, sind die Bedingungen (2) und (3) zueinander gleichwertig.

6. Satz: Voraussetzungen: Für

$$ \varphi(\xi,v_1,\ldots,v_k) = \pmatrix{ \varphi^1(\xi,v_1,\ldots,v_k)\cr \vdots\cr \varphi^{n\ell}(\xi,v_1,\ldots,v_k)\cr } \in \mathbb{R}^{n\ell}, \qquad \xi,v_1,\ldots,v_k\in\mathbb{R}^{n\ell}, $$

sei $\varphi_1'$ die Jacobimatrix von $\varphi$ bzgl. des ersten Vektorargumentes, wobei jede Zeile der Jacobimatrix an verschiedenen Stellen ausgewertet wird, also

$$ \varphi_1'(\boldsymbol{\eta},v_1,\ldots,v_k) = \pmatrix{ D_1 \varphi^1(\eta_1,v_1,\ldots,v_k)\cr \vdots\cr D_1 \varphi^{n\ell}(\eta_{n\ell},v_1,\ldots,v_k)\cr } \in \mathbb{R}^{n\ell\times n\ell}, \qquad \eta_1,\ldots,\eta_{n\ell},v_1,\ldots,v_k\in\mathbb{R}^{n\ell}, $$

wobei $D_1$ Ableitung nach dem ersten Vektorargument bezeichnet, also

$$ D_1 \varphi^i(\xi,v_1,\ldots,v_k) = \left( {\partial\varphi^i\over\partial\xi_1}, \ldots, {\partial\varphi^i\over\partial\xi_{n\ell}} \right) = \nabla_\xi \varphi^i, \qquad \xi_1,\ldots,\xi_{n\ell}\in\mathbb{R}. $$

Für $u_{m+1}$ gelte $\varphi(u_{m+1},y(t_m),\ldots,y(t_{m-k+1}))=0$, und $\varphi_1'$ sei beschränkt und invertierbar für $h\to0$.

Behauptung:$\qquad r={\cal O}(h^{p+1}) \quad\iff\quad s={\cal O}(h^{p+1})$.

Beweis: Es gilt

$$ \eqalign{ 0 &= \varphi(u_{m+1},y(t_m),\ldots,y(t_{m-k+1}))\cr &= \varphi(y(t_{m+1}),y(t_m),\ldots,y(t_{m-k+1})) + \varphi_1'(\boldsymbol{\eta},y(t_m),\ldots,y(t_{m-k+1})) \left(u_{m+1} - y(t_{m+1})\right),\cr } $$

woraus wegen der Beschränktheit und Invertierbarkeit für $\varphi_1'$ die Behauptung folgt.     ☐

7. Corollar: Für das zyklische lineare Mehrschrittverfahren

$$ \varphi(u_{m+1},u_m,\ldots,u_{m-k+1}) = \sum_{i=0}^k A_i u_{m+1+i-k} - h B_i \dot u_{m+1+i-k} $$

gilt für beliebige feste $\mu$:

$$ \eqalign{ & r(t_{m+1},y,h) = (c_{p+1}\otimes I) y^{p+1}(t_\mu) h^{p+1} + {\cal O}(h^{p+1})\cr \iff& s(t_{m+1},y,h) = (A_k^{-1}c_{p+1}\otimes I) y^{p+1}(t_\mu) h^{p+1} + {\cal O}(h^{p+1}).\cr } $$

Deutet man $c_{p+1}$ als Vektor der Fehlerfaktoren der Stufen, so besagt das Resultat, wie man vom lokalen Abbruchfehler in den lokalen Diskretisierungsfehler umrechnet, nämlich einfach durch Multiplikation mit $A_k^{-1}$, zumindestens bis auf Anschmiegungen der Ordnung $p+2$ in $h$. Im Falle $\det A_k=0$, also $\varphi_1'$ nicht invertierbar, gilt der Satz nicht mehr so ohne weiteres.

8. Definition: Sei $y$ ein beliebige Kurve des $\mathbb{R}^{n\ell}$.

(1) Das Verfahren $\varphi$ heißt (polynomial) konsistent mit der (polynomialen) Konsistenzordnung $p$ genau dann, wenn

$$ r(t_{m+1},\mathbb{P}_\mu,h) = 0, \quad\hbox{ für } \mu=0,1,\ldots,p, \quad \forall m\ge k. $$

Äquivalent hierzu ist unter der Bedingung $\varphi_1'$ beschränkt und invertierbar, daß $s(t_{m+1},\mathbb{P}_\mu,h)=0$.

(2) Das Verfahren $\varphi$ heißt (polynomial) konvergent der (polynomialen) Konvergenzordnung $p$ genau dann, wenn $d(t_{m+1},\mathbb{P}_\mu,h)=0$, für $\mu=0,1,\ldots,p$, $h\to0$.

Die Betonung auf polynomialer Konsistenz und polynomialer Konvergenz beruht darauf, daß es auch andere Konsistenz- und Konvergenzarten gibt, z.B. basierend auf trigonometrischen Funktionen oder anderen Funktionensystemen.

Der nächste Abschnitt ist im großen und ganzen nichts weiter als eine längliche Fassung des obigen Corollars.

2. Einfache Beziehung zwischen exakten und berechneten Werten

1. Zu einem $\ell$-stufigen Zyklus mit $k$ Startwerten sei betrachtet

$$ L(t_{m\ell},y,h) := \sum_{i=0}^\kappa \left[(A_i\otimes I)Y_{m+i} - h(B_i\otimes I)\dot Y_{m+i}\right] $$

mit

$$ A_i = \pmatrix{ \alpha_{1,(i-\kappa)\ell+1} & \ldots & \alpha_{1,(i-\kappa)\ell+\ell}\cr \vdots & \ddots & \vdots\cr \alpha_{\ell,(i-\kappa)\ell+1} & \ldots & \alpha_{\ell,(i-\kappa)\ell+\ell}\cr }, \qquad i=0,\ldots,\kappa $$

und

$$ Y_{m+i} = \pmatrix{ y(t_{(m+i-\kappa)\ell+1})\cr \vdots\cr y(t_{(m+i-\kappa)\ell+\ell})\cr }, \qquad i=0,\ldots,\kappa. $$

$\dot Y_{m+i}$ enthält die entsprechenden Ableitungen. Die $B_i$ sind ganz analog wie die $A_i$. Hierbei sind

$$ Y_{m+i}\in\mathbb{R}^{n\ell}, \qquad A_i,B_i\in\mathbb{R}^{\ell\times\ell}, \qquad I = \mathop{\rm diag}(1,\ldots,1)\in\mathbb{R}^{n\times n}. $$

Es ist $\kappa = \lceil k/\ell \rceil$, mit der Anzahl der Startwerte $k$. Es ist

$$ \alpha_{ij} = \beta_{ij} = 0, \qquad j \lt -k+1, \quad i=1,\ldots,\ell . $$

Die Assoziation von $L$ gemäß obiger Festlegung mit $t_{m\ell}$ besitzt eine gewisse Freizügigkeit. Ein fester Wert $t$ innerhalb des Intervalles $[t_{m\ell+1}, t_{m\ell+\ell}]$ beispielsweise, wäre ebenfalls genauso brauchbar.

Sei $L=(L_1,\ldots,L_\ell)$, mit $L_i\in\mathbb{R}^n$ $(i=1,\ldots,\ell)$, so ist hiermit

$$ L_i(t_{m\ell},y,h) = \sum_{j=-k+1}^{\ell} \left[ \alpha_{ij}y(t_{m\ell+j}) - h\beta_{ij}\dot y(t_{m\ell+j}) \right]. $$

2. Ist die Funktion $y$ Lösung der Differentialgleichung

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

so drückt $L_i$ im Falle $L_i(t_{m\ell},y,h)\ne0$ aus, in wie weit die exakte Lösung der Differentialgleichung die $i$-te Stufe nicht erfüllt. Man findet jetzt durch Taylorentwicklung aller $y(t_{m\ell+j})$ mit der Entwicklungsstelle $t_{m\ell+\hat\imath}$, daß

$$ L_i(t_{m\ell},y,h) = c_{i,p_i+1} y^{(p_i+1)}(t_{m\ell+\hat\imath})h^{p_i+1} + {\cal O}(h^{p_i+2}), $$

wobei die $c_{ij}$ gegeben sind durch

$$ c_{ij} = \sum_{\nu=-k+1}^{\ell} \left[ {(\nu-\hat\imath)^j \over j!}\alpha_{i\nu} - {(\nu-\hat\imath)^{j-1} \over (j-1)!}\beta_{i\nu} \right], \qquad j\gt 0, $$

und für $j=0$ durch

$$ c_{i0} = \sum_{\nu=-k+1}^{\ell} \alpha_{i\nu}, \qquad j=0. $$

Die Koeffizienten $c_{ij}$ nehmen also Vorzeichen und Fakultäten der Taylorentwicklung in sich auf. Man beachte, daß alle Stufen der Tendlerschen Zyklen jeweils die gleiche Konsistenzordnung haben, daher erscheint letztendlich für jedes $i\in\{1,\ldots,\ell\}$ die gleiche $h$-Potenz, also schließlich $p\equiv p_i$.

Die $c_{ij}$ für $0\le j\le p_i$ hängen nicht von $\hat\imath$ ab, und $\hat\imath$ muß nicht notwendigerweise ganzzahlig sein. Die $c_{ij}$ für $j>p_i$ können von $\hat\imath$ abhängen.

Für $L=(L_1,\ldots,L_\ell)$ erhält man damit

$$ L(t_{m\ell},y,h) = (c_{p+1}\otimes I) y^{(p+1)}(t_{m\ell+\hat\imath})h^{p+1} + {\cal O}(h^{p+2}), $$

wobei

$$ c_{p+1} = \pmatrix{ c_{1,p+1}\cr \vdots\cr c_{\ell,p+1}\cr }. $$

Der Vektor $c_{p+1}$ hängt nicht von der Taylor-Entwicklungsstelle $\hat\imath$ ab.

Die Konsistenzordnung ist immer gleich der Konvergenzordnung bei sämtlichen Tendlerschen Zyklen, wie man durch längere, aber elementare Rechnung nachweist. Bei zyklischen Formeln, selbst wenn alle Stufen die gleiche Konsistenzordnung aufweisen, muß die Konvergenzordnung nicht grundsätzlich mit der Konsistenzordnung übereinstimmen. Diese Überlegungen findet man in allgemeiner Form dargestellt bei Albrecht (1985)._{Albrecht, Peter}

Seien $y_{m\ell+i}\in\mathbb{R}^n$ im folgenden die vom Programm TENDLER berechneten Näherungen für $y(t_{m\ell+i})\in\mathbb{R}^n$. Die Komponenten von $y_{m\ell+i}$ bzw. von $y(t_{m\ell+i})$ seien mit $y_{\nu,m\ell+i}$ bzw. mit $y_\nu(t_{m\ell+i})$ bezeichnet, $\nu=1,\ldots,n$.

3. Unter der vereinfachenden Annahme, daß

$$ y(t_{m\ell+j}) = y_{m\ell+j}, \qquad j=-k+1,\ldots,0, $$

erhält man durch komponentenweises Anwenden eines Mittelwertsatzes unter Beachtung von beschränkten, geeigneten Inversen ($h$ ausreichend klein), daß

$$ \pmatrix{ y(t_{m\ell+1})\cr \vdots\cr y(t_{m\ell+\ell})\cr } - \pmatrix{ y_{m\ell+1}\cr \vdots\cr y_{m\ell+\ell}\cr } = (d_{p+1}\otimes I) y^{(p+1)}(t_{m\ell+\hat\imath})h^{p+1} + {\cal O}(h^{p+2}) \tag{1} $$

mit

$$ d_{p+1} = A_\kappa^{-1} c_{p+1}. $$

Die vereinfachende Annahme liegt hier darin, daß Gleichheit für “vergangene” Werte $y(t_{m\ell+j})$ und $y_{m\ell+j}$ für geeignete $j$, angenommen wird. Diese Größen können jedoch grundverschieden sein und dürfen nicht miteinander verwechselt werden.

4. Anstatt zyklusweise vorzugehen, kann man auch stufenweise vorgehen. Nimmt man vereinfachend an, daß

$$ y(t_{m\ell+j}) = y_{m\ell+j}, \qquad j=-k+1,\ldots,i-1, $$

für ein festes $i\in\{1,\ldots,\ell\}$, so erhält man

$$ y(t_{m\ell+j}) - y_{m\ell+j} = {c_{i,p+1}\over\alpha_{ii}} y^{(p+1)}(t_{m\ell+i})h^{p+1} + {\cal O}(h^{p+2}), \qquad j=i+1. \tag{2} $$

Die Vereinfachung liegt hier wieder in der vorausgesetzten Exaktheit vergangener Werte. Gleichungen $(1)$ und $(2)$ sind die jeweilige Grundlage für Phase 1 und Phase 2 der Fehlerkontrolle im Programm TENDLER. In der Phase 2 wird lediglich der Fall $i=3$ und ggf. der Fall $i=4$ betrachtet.

5. Dem Term $y^{(p+1)}(t_{m\ell+\hat\imath}){\mskip 3mu} h^{p+1}$ mißt man jetzt in der Fehlerkontrolle große Bedeutung bei. Alle anderen einfliessenden Terme, also insbesondere die Terme höherer Ordnung, oder andere Störgrößen, die z.B. aus der Korrektoriteration stammen könnten, werden nicht weiter gesondert beachtet. Aber auch Effekte, die ihren Ursprung in einem Schrittweiten- oder Ordnungswechsel haben, oder Startfehler, alles dies wird völlig in den nachfolgenden Überlegungen ignoriert. Dies heißt nicht, daß diese Größen keinen Einfluß hätten. Im Vergleich zu den einfachen Tests, wie sie nachfolgend beschrieben werden, entziehen sich diese Einflüsse aber leider einer einfach zu handhabenden, kostengünstigen Überprüfung ihres quantitativen Einflusses. Allerdings wird dann an anderer Stelle versucht durch Sicherheitsfaktoren und Sicherheitsreserven hier einen gewissen Ausgleich zu erzielen.

Es ist nun notwendig $y^{(p+1)}(t_{m\ell+\hat\imath}) {\mskip 3mu} h^{p+1}$ durch leichter berechenbare Größen zu ersetzen. Hierzu wird $\nabla^{p+1}y_{m\ell+i}$ genommen, wobei wie üblich $\nabla$ die rückwärtsgenommene Differenz bezeichnet.

Es gilt nun, z.B. siehe Hildebrand, Francis B., Hildebrand (1956) [S.47, (2.6.6) und S.92, (4.2.7)], oder siehe Fröberg, Carl-Erik, Fröberg (1970) [S.170], daß

$$ \nabla^{p+1}y(t_{m\ell+i}) = y^{(p+1)}(\xi) {\mskip 3mu} h^{p+1}, % \quad\hbox{für eine Zwischenstelle}\quad \qquad \xi \in \left[t_{m\ell+i-p-1}{\mskip 3mu},{\mskip 5mu}t_{m\ell+i}\right]. $$

Für diese Zwischenstelle $\xi$ gilt nicht notwendig $\xi=t_{m\ell+i}$. Da die rückwärtsgenommenen Differenzen ehedem zu Interpolationszwecken berechnet werden müssen, erfordert die Bereitstellung von $\nabla^{p+1}y_{m\ell+i}$ nur eine Vektoroperation. Dieses Argument kann man natürlich auch anders herum betrachten: Da schon sowieso eine Fehlerkontrolle über $\nabla^{p+1}y_{m\ell+i}$ durchgeführt wird, kann man auch leicht mit rückwärtsgenommenen Differenzen interpolieren.

Nach diesen Vorbereitungen sollen nun die beiden wichtigen Glieder der Fehlerkontrolle vorgestellt werden.

Bibliographisch:

  1. Hildebrand, Francis Begnaud: “Introduction to Numerical Analysis”, McGraw Hill Book Company, New York Toronto London, 1956, x+511 S.
  2. Fröberg, Carl-Erik: “Introduction to Numerical Analysis&rdquo, Addison-Wesley, Reading (Massachusetts) Menlo Park London Sydney Manila, 2nd Edition, 1970, xii+433 S.

3. Die Phase 1 der Fehlerkontrolle

1. Der Fehlerkontroll-Teil ist in zwei Phasen eingeteilt. In der ersten Phase wird nach den ersten beiden Schritten $\nabla^{p+1}y_{m\ell+2}$ ausgerechnet. Die berechneten beiden Näherungen $y_{m\ell+1}$ und $y_{m\ell+2}$ werden nun akzeptiert, falls

$$ \left\|\nabla^{p+1}y_{m\ell+2}\right\|_w \: \le\: {\varepsilon\over\mathopen|d_{2,p+1}\mathclose|}. \tag{Phase 1} $$

Hierbei sind $d_{i,p+1}$, mit $i=1,\ldots,\ell$, die Komponenten des Vektors $d_{p+1}$. Der Nenner wurde gewählt als

$$ \max\left( \mathopen|d_{1,p+1}\mathclose|, \mathopen|d_{2,p+1}\mathclose|\right) = \mathopen|d_{2,p+1}\mathclose|. $$

Diese Beziehung gilt für alle Ordnungen $p=1,\ldots,7$.

Aus Sicherheitsgründen und um möglichst die vorgegebene Toleranzschwelle nicht zu überschreiten, wird natürlich dem Maximum der beiden Werte der Vorzug eingeräumt. Im Falle, daß Phase 1 oder Phase 2 einen Schritt zurückweisen, kommen dann anschliessend die gleichen Strategien zur Anwendung.

D.h. in der zweiten Stufe entscheidet man über die zweite berechnete Näherung und gleichzeitig aber auch über die erste Näherung. Dies liegt daran, daß nicht immer genügend viele vergangene Lösungswerte bekannt sind, um die rückwärtsgenommenen Differenzen entsprechend hoher Ordnung zu bilden. Dies ist beim Start der Integration der Fall.

Die Vektoren $d_{p+1}=A_\kappa^{-1}c_{p+1} \in \mathbb{R}^{\ell\times\ell}$, für die Ordnungen $p=1,\ldots,7$, sind wie nachfolgend, teilweise gerundet.

$i$ $d_2$ $d_3$ $d_4$ $d_5$ $d_6$ $d_7$ $d_8$
1 -0.5 $-0.\overline{2}$ $-0.1\overline{36}$ -0.096 -0.072992 -0.058309 -0.048209
2 -1 $-0.\overline{518}$ -0.359504 -0.28032 -0.232830 -0.203812 -0.178361
3 -1.5 -0.839506 1.680440 -0.704419 -0.427500 -0.396400 -0.468812
4 * * * * -1.174794 -1.295310 -1.185795

Hierbei sieht man also tatsächlich, daß von den ersten beiden Komponenten, die zweite betragsmässig immer die größere ist. Nur die zweite Zeile findet in dem Programm TENDLER Anwendung, die anderen Zeilen werden nicht benötigt. Abgespeichert werden jedoch nicht diese Werte sondern die Kehrwerte.

2. Es ist hierbei $|\cdot |_w$ dieselbe gewichtete Norm wie sie auch bei der Korrektoriteration auftrat. Eine Möglichkeit den Gewichtsvektor zu aktualisieren ist wie folgt. Man wählt für jede Komponente des Gewichtsvektor das Maximum von

  1. eins und
  2. dem Betrag der größten Näherung, die bisher berechnet wurde.

Diese Wahl der Gewichte würde die Fehlerkontrolle in der folgenden Weise beeinflussen: Die Fehlerkontrolle wäre

  1. relativ für diejenigen Komponenten die betragsmässig größer sind als eins,
  2. absolut für die Komponenten, die betragsmässig kleiner als 1 sind und so bleiben und
  3. absolut gefolgt von relativ für diejenigen Komponenten, die erst betragsmässig kleiner und anschließend größer sind als 1.

Die genaue Wahl des Gewichtsvektors kann der Benutzer beliebig vorgeben. Das Programm TENDLER erzwingt oder erfordert hier keinerlei Einschränkungen. Nichts desto trotz hat die Wahl der Norm Einfluß auf die Effizienz des Programmes.

Nicht unüblich ist hierbei auch eine Kopplung des Gewichtsvektors mit der geforderten Genauigkeit. Aber auch mit berechneten Werten für die Ableitungen kann man die Normroutine versehen. Diese Berechnungsvorschriften hängen z.T. von der zugrunde liegenden Zielsetzung des Programmbenutzers ab.

Eine Beschreibung zu dem Komplex der Normberechnung findet man z.B. bei Shampine/Gordon (1984). Wenn Geschwindigkeitsüberlegungen im Vordergrund stehen, ist eine auf der Maximumnorm aufbauende Norm kostengünstiger zu berechnen, als eine Routine basierend auf der euklidischen Norm. Mit beiderlei Arten von Normroutinen wurde das Programm TENDLER getestet.

Biographisch:

  1. Lawrence F. Shampine
  2. Marilyn K. Gordon

4. Die Phase 2 der Fehlerkontrolle

1. Wie vorher wird $y^{(p+1)}(t_{m\ell+\hat\imath}) {\mskip 3mu} h^{p+1}$ durch $\nabla^{p+1}y_{m\ell+i}$ angenähert, insbesondere mit $\hat\imath=i$. Die Werte $1/d_{2,p+1}$ liegen in etwa zwischen 1 und 6. Die Werte $\alpha_{ii}/c_{i,p+1}$ liegen größenordnungsmäßig zwischen 1 und 20.

Die maximale Stufenzahl $\ell$ ist hierbei $\ell=3$, für $p\in\{1,2,3,4\}$ und $\ell=4$, für $p\in\{5,6,7\}$. In der zweiten Phase der Fehlerkontrolle wird für die letzte $\ell$-te Stufe, bzw. letzten beiden (d.h. für $\ell-1$ und $\ell$) Stufen im Zyklus, der folgende Test angewendet, um auf Akzeptanz zu prüfen:

$$ \left\| \nabla^{p+1}y_{m\ell+i} \right\|_w \: \le\: {\varepsilon \over c_{i,p+1}/ \alpha_{ii}}, \qquad\hbox{für}\quad i=3,\ldots,\ell. \tag{Phase 2} $$

Wird dieser Test nicht bestanden, so wird der gesamte Zyklus zurückgewiesen, und es wird für den zu wiederholenden Zyklus eine geeignete neue Schrittweite bestimmt.

Je nach gewählter Option können die Tests während der Phase 2 zu einem einzigen Test zusammenschrumpfen. Testläufe zeigen, daß die maßgeblichste Überprüfung, diejenige der Phase 1 ist, also gemeint sind damit die Überprüfungen, die zu einem Wiederholen des Zyklus führen. Aufgrund der Größenverhältnisse der verwendeten Konstanten ist dies nicht verwunderlich. Dies darf natürlich nicht dahingehend interpretiert werden, daß dies bei allen zyklischen Formeln so der Fall ist, auch nicht bei den zyklischen Formeln von J.M. Tendler. Zu den Gründen warum gleich der gesamte Zyklus wiederholt wird, anstatt nur die betreffende Stufe, gelten die Überlegungen, die beim Konvergenztest angegeben wurden, entsprechend (Stummelzyklusproblem). Gelegentlich tritt der Fall ein, daß ein zurückgewiesener Schritt seine Ursachen in einer zu alten Iterationsmatrix $W$ hatte, jedoch ist dies selten. Es ist gerade die Aufgabe des Hindmarsh-Testes hier für eine entsprechend günstige Aktualisierung zu sorgen.

2. Zur Interpolation der Zwischengitterpunkte wird das Unterprogramm tdlintsglbase() aufgerufen, welches wiederum $(p-1)$-mal das Unterprogramm tdlnewtform() aufruft.

Da $y_{m\ell}$ schon bekannt ist, müssen nur $(p-1)$-mal, anstatt $p$-mal Werte berechnet werden. Die insgesamt $p$ Ordinatenwerte müssen dann umgeformt werden in rückwärtsgenommene Differenzen. Dies geschieht mit dem Unterprogramm tdlformbackdiff().

Schließlich muß noch für den Prädiktor der Wert $hf_{m\ell-1}$ auf die neue Schrittweite angepasst werden. Dies besorgt das Unterprogramm tdladapthz().

Der Arbeitsaufwand für einen Schrittweitenwechsel beträgt insgesamt ${\cal O}({3\over2}p^2+{1\over2}p)$ Vektoroperationen. Die Bereitstellung und die ständige Pflege und Aktualisierung der Differenzentabelle pro Stufe kostet $p$ Vektoroperationen. Das Programm DE/STEP benötigt hier $(p+2)$ Vektoroperationen. Bzgl. der Überlegungen, ob es günstiger ist den gesamten Zyklus zu wiederholen oder nur die zurückgewiesene Stufe, beachte man, daß innerhalb eines Zykluses dann eine Differenzentabelle für Ableitungswerte mitgeführt werden müsste. Dies würde den Speicherplatzbedarf auf jeden Fall anschwellen lassen und möglicherweise auch den Rechenzeitbedarf. Wie erwähnt treten auch theoretische Schwierigkeiten bzgl. des Stummelzyklus, also des nicht zu Ende geführten Zyklus, hinzu.

Nocheinmal soll darauf hingewiesen werden, daß in dem Programm TENDLER nicht alle Stufen einer Fehlerkontrolle unterzogen werden. Beachtet man gewisse Eigenheiten beim Start der Integration, so liesse sich diese Eigenschaft leicht verändern. Der Integrationsanfang ist insofern einer Sonderbehandlung zu unterziehen, als daß dort nicht gleich genügend viele rückwärtsgenommene Differenzen entsprechender Ordnung zur Verfügung stehen.

5. Die verwendeten Fehlerfaktoren im Programm TENDLER

Die Reziprokwerte der Beträge der skalierten Fehlerfaktoren, also

$$ \left|{c_{i,p+1} / \alpha_{ii}}\right|, \qquad i=1,\ldots,4, \quad p=1,\ldots,7 $$

sind in der fünfzeiligen Tabelle zusammengefaßt. Die Werte sind gerundet.

$p$ 1 2 3 4 5 6 7
$\left\vert{\alpha_{11}/c_{1,p+1}}\right\vert$ 2.0 4.5 $7.\overline{3}$ 10.4167 13.7 17.15 20.7429
$\left\vert{\alpha_{22}/c_{2,p+1}}\right\vert$ 2.0 4.5 $7.\overline{3}$ 10.4167 13.7 17.9504 20.7429
$\left\vert{\alpha_{33}/c_{3,p+1}}\right\vert$ 2.0 4.5 6.0 9.3 13.8687 17.4349 15.921
$\left\vert{\alpha_{44}/c_{4,p+1}}\right\vert$ * * * * 9.6904 9.472 14.7809

Man behalte im Auge, daß für die Fehlerkontrolle, nur die letzten beiden Zeilen der bzw. die letzten drei Zeilen der angegebenen fünfzeiligen Tabelle verwendet werden und zwar für die zweite Phase des Fehlerkontrollteils.

Die erste Phase benutzte die Werte $d_{p+1}=A_\kappa^{-1}c_{p+1}$. Die obigen Werte $\alpha_{ii}/c_{i,p+1}$ werden in transponierter Form in der Konstantenmatrix alfa_div_errfac[8][4] abgespeichert. Die erste Zeile dieser Tabelle wird also in dem Programm TENDLER nirgendwo benutzt. Dennoch sind diese Werte gespeichert, da sie u.U. bei einer Änderung des Konvergenztestes benutzt werden könnten.

Für die erste Phase werden ja die Reziprokwerte der Beträge $\mathopen|d_{2,p+1}\mathclose|$ für die Ordnungen $p=1$ bis $p=7$ benutzt. Diese lauten gerundet:

$p$ 1 2 3 4 5 6 7
$\mathopen\vert1/d_{2,p+1}\mathclose\vert$ 1.0 1.9286 2.7816 3.5674 4.295 4.9065 5.6066

Diese Konstanten sind in dem Konstantenvektor errconst[8] abgespeichert.

6. Programmnahe Beschreibung der Fehlerkontrolle

Since for a stiff problem it is stability and not accuracy that controls step-length, an error control mechanism based on accuracy requirements is unnecessary. The numerical results certainly bear this out. We are left with the interesting thought that, in general, algorithms for stiff systems, based on explicit methods used in some acceptably stable manner, do not need to have an error control mechanism.

Jack D. Lambert (1985)

1. Die erste Stufe des Zykluses wird immer durchgelassen. Im Programm TENDLER ist hier also $i=0$. In diesem Falle finden somit keinerlei weitere Überprüfungen statt, insbesondere keine Normberechnungen. Gesetzt wird nun die maximal zulässige Schranke $\varepsilon_{\rm max}$ für die weiteren Stufen, also für $i=1,2$ und gegebenenfalls noch für $i=3$. Man beachte hier die Indizierung von 0 an, anstatt von 1, im Gegensatz zu oben. Man setzt nun

$$ \varepsilon_{\rm max} = \varepsilon\cdot\cases{ \texttt{errconst[}p\texttt{]}, & falls $i=1$,\cr \texttt{alfa_div_errfac[}p\texttt{][}i\texttt{]}, & falls $i\gt 1$.\cr} $$

und testet jetzt hiermit auf Akzeptanz des Schrittes durch den Test

$$ \mathopen\|\nabla^{p+1}y\mathclose\| \le \varepsilon_{\rm max}. $$

Wird dieser Test nicht bestanden, so wird die Schrittweite $h$ auf jeden Falle verkleinert und zwar um den Faktor ratio, mit

$$ \texttt{ratio} = \cases{ \texttt{ecsaffac[}\xi\texttt{]}\cdot {\displaystyle\root p+1\of{\varepsilon_{\max}\over \mathopen\|\nabla^{p+1}y\mathclose\|}}, & beim ersten Fehlschlag,\cr 1/2, & beim zweiten Fehlschlag,\cr 1/10, & beim dritten Fehlschlag.\cr } $$

Hierbei ist

$$ \xi = \begin{cases} 0, & \text{falls Iterationsart} \in\left\{ \texttt{MNI}, \texttt{MNIFLOAT} \right\}, \\ 1, & \text{falls Iterationsart} \in\left\{ \texttt{FIX}, \texttt{FIXFLOAT}, \texttt{FIXPECE} \right\}. \\ \end{cases} $$

Schließlich ist ecsaffac[] = $\{0.7, 0.6\}$.

Beim dritten konsekutiven Fehlschlag, also drei Fehlschläge direkt hintereinanderfolgend bei demselben Zyklus, wird die Ordnung $p$ auf 1 gesetzt und dies unabhängig davon, welche Ordnung vorher auch immer benutzt wurde.

Erfahrung und Testläufe zeigten, daß ein Fehlschlag selten ist, circa 1 Fehlschlag auf 400 Schritte, als grober Hinweis. Dies ist natürlich von der Differentialgleichung und der angeforderten Genauigkeit abhängig, aber auch —und das ist wichtig— von der Schrittweiten- und Ordnungssteuerung.

Ein doppelter Fehlschlag ist noch seltener. Dreifache Fehlschläge in dem Programm TENDLER sind ebenfalls sehr selten. Für Aussagen solcher Art dienen u.a. die generierten Statistiken. Gleichzeitig ergab sich, daß bei vierstufigen Zyklen, also bei den Ordnungen $p\in\{5,6,7\}$, die Fehlerkontrolle für die dritte Stufe ($i=2$) nie versagte. Anteil hat maßgeblich die zweite Stufe und danach die letzte Stufe. Diese Beobachtung gilt ersteinmal nur für das Programm TENDLER, mit seinen Faktoren $\left|\alpha_{33}/c_{3,p+1}\right|$. Für andere Programme, basierend auf möglicherweise anderen zyklischen Formeln, kann das Bild vielleicht schon ganz anders aussehen.

Bei der Zurückweisung eines Schrittes werden die neu benötigten Werte auf dem neuen, engeren Gitter durch Interpolation bestimmt. Bei einem Einschrittverfahren entfällt diese Aufgabe. Wenn in dem Programm TENDLER mit dem impliziten Eulerverfahren gerechnet werden soll, ist also bei einem Schrittweitenwechsel keine Interpolation erforderlich. Die Schrittweite wird um den Faktor ratio verkleinert. Die Interpolation mit Hilfe rückwärtsgenommener Differenzen geschieht durch das Unterprogramm tdloldbdf().

2. Wird ein Schritt zurückgewiesen, so werden die folgenden Maßnahmen ergriffen:

  1. Die Ordnung wird niemals erhöht, höchstens drastisch erniedrigt.
  2. Falls die Fehlerkontrolle weniger als drei Mal nicht erfolgreich war, so wird die Schrittweite verkleinert, sofern sie die minimale Schrittweite noch überschreiten sollte.
    • Nach dem ersten Fehlschlag, wird die neue Schrittweite ähnlich bestimmt, wie im Abschnitt über die Schrittweiten- und Ordnungssteuerung.
    • Nach dem zweiten Fehlschlag wird die Schrittweite halbiert, soweit sie dieses zuläßt.
  3. Sollte die Schrittweite die kleinstzulässige Schrittweite nicht überschreiten, so wird die Ordnung auf 2 gesetzt, falls die Ordnung größer ist als 2, und sie wird auf 1 gesetzt, falls sie vorher nur noch 2 war.
  4. In den beiden Fällen (2) und (3) wird dann der gesamte Zyklus wiederholt.
  5. Sollte das nicht möglich sein, d.h. die Schrittweite ist minimal und die Ordnung ist gleich eins, dann wird der Zyklus zwangsweise akzeptiert. Je nach entsprechend gesetztem Modus erfolgt eine Fehlermeldung und eine entsprechende Markierung in den Rückgabeparametern.
  6. Nach drei Zurückweisungen wird der Zyklus gänzlich neu gestartet, und zwar fast so, wie bei dem allerersten Aufruf des Programmes TENDLER. Die Schrittweite wird, soweit sie dieses zuläßt, noch gezehntelt.
  7. Das Betreten der Schrittweiten- und Ordnungssteuerung wird durch Setzen der Sperrvariablen problem und wait[] unterbunden.

Die in dem Programm TENDLER verwendeten Taktiken sind denjenigen des Programmes STINT und DIFJMT ähnlich. Die verwendeten Strategien sind überdies ähnlich denjenigen Taktiken, wie sie auch verwendet werden in dem Programm DIFSUB, siehe Gear (1971b), siehe Gear (1971c) [9.3, S.156--166]. Ähnlichkeiten bestehen auch zu dem Programm DSTIFF, siehe Gupta, Gopal K., Gupta (1985), S.942–943, welches seinerseits sich anlehnt an das Programm DIFSUB. Es gibt jedoch auch eine Reihe von Unterschieden. Z.B. wird an keiner Stelle pauschal 10-mal gewartet bis die Schrittweiten- und Ordnungssteuerung erneut betreten wird.

Bibliographisch:

  1. Lambert, Jack D.: A stable sequence of steplengths for Euler's rule applied to stiff systems of differential equations, Computer & Mathematics with Applications, Vol. 12, Issues 5-6, Part 2, September-December 1986, pages 1141–1151
  2. Gupta, Gopal K.: “Description and Evaluation of a Stiff ODE Code DSTIFF”, SIAM Journal on Scientific and Statistical Computing, Vol 6, No 4, October 1985, pp.939–950
  3. Gear, Charles William (1935—2022)
  4. Gear (1968), Gear, Charles William: “The Automatic Integration of Stiff Ordinary Differential Equations”, IFIP Congress 68, Edinburgh, Scotland, $5^{\rm th}$–$10^{\rm th}$ August 1968, North Holland Publishing Company, Amsterdam, pp.A81–A85
  5. Gear (1971a), Gear, Charles William: “Experience and Problems with the Software for the Automatic Solution of Ordinary Differential Equations”, in “Mathematical Software”, ACM Monograph Series, Editor John Rice, Academic Press, New York London, 1971, pp.211–227
  6. Gear (1971b), Gear, Charles William: “The Automatic Integration of Ordinary Differential Equations”, CACM, Vol 14, No 3, March 1971, pp.176–179
  7. Gear (1971b), Gear, Charles William: “Algorithm 407. DIFSUB for Solution of Ordinary Differential Equations [D2]”, CACM, Vol 14, No 3, March 1971, pp.185–190
  8. Gear (1971c), Gear, Charles William: “Numerical Initial Value Problems in Ordinary Differential Equations”, Prentice-Hall, Englewood Cliffs (New Jersey), 1971, _xvii+253 S.
  9. Tendler, Joel Marvin: A Stiffly Stable Integration Process Using Cyclic Composite Methods
  10. Tendler, J.M., Bickart, T.A., Picel, Z.: A Stiffly Stable Integration Process Using Cyclic Composite Methods. ACM Trans. Math. Softw. 4(4), 339–368 (1978)