UNIVERSITÄT KARLSRUHE





FAKULTÄT MATHEMATIK






Hauptseminar: Differentialgleichungen
bei Prof. Dr. Plum
WS 1997/1998




Einschließungen gewöhnlicher Differentialgleichungen
nach Dipl. Math. Rudolf Lohner

Michael Ralph Pape

Inhalt

• Inhalt
• 1 Einleitung
  • 1.1 Die Anfangswertaufgabe
  • 1.2 Lösungsmethoden
  • 1.3 Beschreibung des Verfahrens
• 2 Mathematische Betrachtungen
  • 2.1 Begriffe
  • 2.2 Grundproblem
  • 2.3 Umwandlung in ein Fixpunktproblem
  • 2.4 Taylorentwicklung $(p-1)$-ten Grades mit Restglied
• 3 Einschließung bei gewöhnlichen AWA
  • 3.1 Einschrittverfahren
  • 3.2 Entwicklung nach lokalen Fehlern
  • 3.3 Berechnung der $A_j$
  • 3.4 Gliederung des Verfahrens
  • 3.5 Behandlung des globalen Fehlers
  • 3.6 Ergebnisvergleich der verschiedenen Auswertungsmethoden

1 Einleitung

Wir haben folgendes Problem.

Wir haben eine Anfangswertaufgabe AWA, von der wir nur die Lösung an einer Stelle kennen. Wir wollen ein Verfahren entwickeln, daß uns eine Näherungslösung an einer nachfolgenden Stelle angibt, genauer, eine Einschließung der Näherungslösung an dieser Stelle.

Dieses Verfahren soll später auf einem Rechner impementiert werden.

1.1 Die Anfangswertaufgabe

Eine Anfangswertaufgabe ist gegeben durch

$$[1]\quad \left\{ \begin{array}{l} u'=f(u), f:D\to {\mathrm{I... ...m}R}}^n, u(t_0) = u_0,\ u_0 \in D . \end{array}\right.$$

Die Bedingung $u(t_0)=u_0$ heißt Anfangsbedingung. Eine Funktion $u:J\to {\mathrm{I\hspace*{-0.2em}R}}^n$ heißt Lösung der AWA auf einem Intervall $J$, wenn $u(t)$ dort differenzierbar ist, $t_0\in J$ ist, für alle $t\in J$ auch $(t,u(t)) \in D$ gilt, und wenn $u(t)$ auf $J$ die AWA [1] (1.2.2.2) erfüllt.

1.2 Lösungsmethoden

Man kann versuchen eine AWA, mit den Verfahren die wir in ANA III kennengelernt haben, exakt zu lösen. Dies geht leider nur, wenn man es schafft die DGL in eine dort aufgeführte Form umzuformen.

Man kann aber auch versuchen die Lösung durch Einschrittverfahren zu approximieren und sie dann einzuschließen.

1.3 Beschreibung des Verfahrens

Nun taucht die Frage auf, was das vorhin erwähnte Verfahren alles können soll?

Das Verfahren soll den Wert der Lösung $u$ an einer Stelle $t$ approximieren. Dabei soll ausgehend von $u$ mittels Betrachtung einer Taylorentwicklung mit Restglied auf die Lösung an der nachfolgenden Stelle $t_{j+1}$ geschlossen werden. Der dabei auftretenden lokale Fehler, das Restglied, wird dann approximiert und durch ein Intervall eingeschlossen.

Durch iteratives Anwenden des Verfahrens kann die Lösung z.B. an äquidistanten Punkten bestimmt werden. Dazu definieren wir ein Gitter $t_j=t_0+h\cdot j$ mit (konstanter) Schrittweite $h>0$. Wir betrachten obige AWA und bezeichnen dabei deren Lösungen $u(t_j)$ in den Gitterpunkten $t_j$ mit $u_j$.

Des weiteren gibt es zu sagen, daß das Verfahren durchweg Intervallrechnung verwenden soll.

2 Mathematische Betrachtungen

Nun wollen wir das Verfahren von der mathematischen Seite betrachten.

Dazu folgt erste eine kurze Einführung von Begriffen.

2.1 Begriffe

Intervallarithmetik:

Die Intervallaritmetik ist eine Arithmetik, die für abgeschlossene reelle Intervalle definiert ist. Derartige Intervalle werden folgendermaßen bezeichnet:

$$[a]:=[\underline a, \overline a]:= \{a\in {\mathrm{I\hspace*{-0.2em}R}} \vert \underline a \leq a\leq \overline a\} .$$

Sind Untergrenze und Obergrenze gleich, so liegt der Sonderfall eines entarteten Intervalls oder Punktintervalls vor. Puntintervalle werden mit den entsprechenden reellen Zahlen identifiziert, d.h. für eine reelle Zahl benutzten wir die Schreibweise ${a=[a]=[a,a]}$.

Beispiel:

$$\underbrace{[x]}_{\in {\mathrm{I\hspace*{-0.2em}R}}}+\unde... ...ine{x}+\overline{y}]} _{\in {\mathrm{I\hspace*{-0.2em}R}}}$$

Dabei ist ${\mathrm{I\hspace*{-0.2em}R}}$ die Menge der reellen Intervalle.

Des weiteren kann man Intervalle in Funktionen einsetzten. Die intervallmäßige Auswertung einer Funktion schließt dann den Wertebereich einer Funktion auf einem Intervall ein. Es gilt also:

$$W(f,[a,b])\subseteq f([a,b])$$

Das heißt nichts anderes, als das wir in der Intervallarithmetik mit Abschätzungen rechnen können. ohne das wir einen Fehler machen.

Maschinenintervallarithmetik:

Die Maschinenintervallarithmetik ist die Anwendung der Intervallarithmetik auf einem Computer. Wir gehen davon aus, daß im Rechner eine Fließkommaarithmetik verwendet wird. Die Menge der von einem Recher darstellbaren Zahlen, nennen wir die Menge der Maschinenzahlen und wollen sie mit R bezeichnen.

Beispiel:

Addieren wir zwei reelle Intervall im Rechner, so liefert uns die Intervallarithmetik wie wir oben gesehen haben $[\underline{x}+\underline{y},\overline{x}+\overline{y}]$. Diese Zahlen müssen aber für der Rechner nicht unbedingt darstellbar sein. Damit unser Ergebnis richtig bleibt, muss der Rechner dann runden. Und wie?

$$\underbrace{[x]}_{\in {\rm {\mathrm{I\hspace*{-0.2em}R}}}}... ... \in {\rm R}} , \bigtriangleup(\overline{x}+\overline{y})]$$

Dabei bedeutet $\bigtriangledown$ die Rundung zur nächst kleineren Gleitkommazahl und $\bigtriangleup$ die Rundung zur nächst größeren Gleitkommazahl.

Wir haben folgendes Grundproblem:

2.2 Grundproblem

Sei $J=[t_i,t_{i+1}]$ ein Intervall. Könnten wir die Lösung $u$ auf dem Intervall $J$ einschließen, so könnten wir auch deren Ableitung $u'=f(u)$ auf dem Intervall $J$ einschließen.

Wir suchen also eine Einschließung der Lösung $u$ auf dem Intervall $J$.

2.3 Umwandlung in ein Fixpunktproblem

Dazu wandeln wir zuerst das AWP in ein Fixpunktproblem um, bei dessen Lösung uns der Banachsche Fixpunktsatz behilflich ist.

Wegen der Stetigkeit von $f$ in der AWA [1] (1.2.2.2) kann die AWA durch komponentenweise Integration in eine äquivalente Integralform umgewandelt werden. D.h. es gilt

$$u(t)=u_0+\int\limits_{t_0}\limits^t f(u(s))ds.$$

Diese Integralgleichung schreiben wir mit $u=Tu$ als Operatorgleichung.

$T$ ist dabei der Integraloperator, der aus der ANA III Vorlesung bekannt ist.

Definition 2.1
Sei $[u^0]:=[\underline{u}^0, \overline{u}^0] \subseteq D\subseteq {\mathrm{I\hspace*{-0.2em}R}}^n$, $J$ ein Intervall. Dann sei

$$U:= C(J,[u^0])$$

die Menge der auf $J$ stetigen Funktionen, deren Bild in $[u^0]$ ist.

Gelten die Vorraussetzungen,

• wir haben einen Banachraum (vollständiger, metrischer Raum) auf einer Funktionen-Menge $U$
• der Operator $T$ sei auf $U$ erklärt und genüge auf $U$ einer Lipschitzbedingung, mit Lipschitzkonstanten $k<1$, d.h. der Operator sei kontrahierend. (Dies ist bei passender Normwahl kein Problem.)
• der Operator $T$ bilde $U$ in sich ab, d.h. $T(U)\subseteq U$,

dann sagt uns der BFPS, daß die Gleichung $T(u)=u$ auf $U$ genau eine Lösung hat.

Also, jede Lösung des AWP ist ein Fixpunkt des Operators $T$, d.h. $T(u)=u$.

Und umgekehrt ist jeder Fixpunkt des Operators $T$ Lösung der AWA.

Wir wenden diesen BFPS auf unser Problem an.

Satz 2.1 ( I )
Es seien die AWA [1] (1.2.2.2) und der Operator $T$ auf $U$ über dem Intervall $J=[t_0,t_0+h]$ gegeben. Gilt für ein Intervall $[u^0]\subseteq D$, daß

$$[2]\quad [u^1]:=u_0+[0,h]\cdot f([u^0])\subseteq [u^0]$$

ist, so besitzt die AWA auf ganz J eine eindeutige Lösung, die ganz in dem Intervall $[u^0]$ verläuft. Hierbei ist $f([u^0])$ eine intervallmäßige Auswertung von $f$ auf $[u^0]$.

Dieser Satz ist eine Folgerung aus dem BFPS. Zu beweisen ist noch die Gleichung $TU\subseteq U$. Die anderen Vorraussetzungen sind erfüllt.

Beweis 1
Wir zeigen nun, daß der Operator $T$ die zu $[u^0]$ gehörige Funktionenmenge $U$ in sich abbildet.

$U$ ist eine abgeschlossene Teilmenge des Banachraumes $C(J,{\mathrm{I\hspace*{-0.2em}R}}^n)$. Die Wertebereiche aller Funktionen $u\in U$ liegen in $[u^0]$. Sei ${u \in U}$ beliebig, dann ist also ${u(t) \in [u^0]}$ für alle $t\in J$ und es gilt:

$$\begin{eqnarray*} (Tu)(t)&=&u_0+\int\limits_{t_0}\limits^t f(u(s))ds \stackr... ...)\\ &=& [u^1]\stackrel{\rm nach Vor.}\subseteq [u^0]\quad . \end{eqnarray*}$$

Daraus folgt mit dem BFPS, daß die AWA auf ganz J genau eine Lösung besitzt die ganz in $[u^0]$ verläuft.

Das Gleichheitszeichen am Ende der ersten Zeile gilt wegen

$$\int\limits_a \limits^b [\underline{c}, \overline{c}] dx ... ...\overline{c} (b-a)]= [\underline{c}, \overline{c}] (b-a)$$

für $a\leq b$ und beliebige Intervalle $[c]=[\underline{c},\overline{c}]\subseteq {\mathrm{I\hspace*{-0.2em}R}}$.

Die Inklusion in der zweiten Zeile gilt wegen $t\in[t_0,t_0+h]$. Setzt man die Grenzen $t_0$ sowie $t_0+h$ für $t$ ein, so erkennt man, daß sich $(t-t_0)$ nur in dem Intervall $[0,h]$ bewegen kann.

Für eine beliebige Funktion $u\in U$ liegt also der Wertebereich von $Tu$ wieder in dem Intervall $[u^0]$, d.h. es ist $Tu\in U$. Da $T$ auf $U$ erklärt ist, dort einer Lipschitzbedingung mit Lipschitzkonstanten ${k=1/2<1}$ genügt und ferner $T$ die abgeschlossenen Teilmenge $U$ des Banachraumes $C(J,{\mathrm{I\hspace*{-0.2em}R}}^n)$ in sich abbildet, folgt aus dem Banachschen Fixpunktsatz, daß der Operator $T$ in $U$ genau einen Fixpunkt $u^*$ besitzt.

Wegen $(Tu)(t) \in [u^1]$ für alle Funktionen $u\in U$ und $t\in J$ verläuft die Fixpunktfunktion $u^*$ sogar ganz in dem Intervall $[u^1]$. Damit hat auch die AWA genau eine Lösung auf dem Intervall $J=[t_0,t_0+h]$, die ganz in dem Intervall $[u^0]$, sogar ganz in dem Intervall $[u^1]$ verläuft.

Die Bedingung [2] (1.2.2.9) $[u^1]\subseteq[u^0]$ garantiert zwar die Existenz und die Eindeutigkeit einer Lösung auf dem Intervall $J$, ist aber vom numerischen Standpunkt in vielen Fällen unbefriedigend, da meist sehr kleine Schrittweiten $h$ gewählt werden müssen!

Insgesamt können wir also mit Hilfe des Banachschen Fixpunktsatz ein grobes Einschließungsintervall $[u^0]$ bzw. $[u^1]$ der Lösung $u$ auf $J=[t_0,t_0+h]$ finden.

Im nächsten Abschnitt soll gezeigt werden wie man die Ableitungen und die Koeffizienten der Taylorentwicklung rekursiv berechnen kann.

2.4 Taylorentwicklung $(p-1)$-ten Grades mit Restglied

Mit Hilfe der rekursiven Berechnung der Ableitungen und der Taylorkoeffizienten soll durch entwickeln der Lösung $u_j$ auf die Lösung $u_{j+1}$ zum Zeitpunkt $t_{j+1}$ geschloßen werden.

2.4.1 Abkürzung für Taylorkoeffizienten

Im Folgenden wird die Abkürzung

$$\begin{array}{lc} [3]\quad& (u)_k:= \frac 1 {k!}\cdot u^... ...ad \textrm{f\uml {u}r } k=0,1,2,\ldots .\\ \end{array}$$

für die Taylorkoeffizienten verwendet.

2.4.2 Rechenregeln für Taylorkoeffizienten

Für Funktionen, die durch arithmetische Operationen aus anderen Funktionen zusammengesetzt sind, lassen sich Rechenregeln für die Taylorkoeffizienten angeben, wie z.B.

$$\begin{array}{ccll} (u \pm v)_k&=&(u)_k \pm (v)_k & \qua... ...v)_k \right) & \quad \textrm{Cauchyprodukt}\\ \end{array}$$

Diese Regeln sind lediglich die Regeln für das Rechnen mit Potenzreihen.

Berücksichtigt man noch die Beziehung

$$\begin{array}{cclll} (c)_0&=&c, &(c)_k=0 \quad \textrm{... ...geq 1, & \textrm{f\uml {u}r Konstanten } c, \end{array}$$

sowie die Beziehung

$$\begin{array}{cclll} (t)_0&=&t_0, &(t)_1=1, (t)_k=0 ... ...l {u}r die unabh\uml {a}ngige Variable } t,\\ \end{array}$$

so kann man für beliebige rationale Funktionen die Taylorkoeffizienten jeder Ordnung $k$ berechnen, indem man rekursiv für alle Teilausdrücke zunächst diese Koeffizienten für $k=0$, dann $k=1$, usw. berechnet.

Nun wollen wir die Taylorkoeffizienten und die Ableitungen des Taylorpolynoms rekursiv berechnen.

2.4.3 Taylorentwicklung

Der Einfachkeit halber nehmen wir an, die betrachteten Funktionen, sind reelle Funktionen von $t$ in einer Umgebung von $t_0$. Sie sollen noch $p$-mal stetig differenzierbar sein, mit $p$ genügend groß.

Mit der obigen Abkürzung haben wir z.B. eine Taylorreihe von $u$ um $t_0$ in der Form

$$[4]\quad u(t)=\sum\limits_{k=0}\limits^\infty (u)_k (t-t_0)^k .$$

Wir schreiben nun zur Funktion $u$ aus [4] (1.1.4.2) die Reihe für die Ableitung $u'$ auf. Diese ist

$$u'(t) \stackrel{\rm Def.}{=} \sum\limits_{k=1}\limits^\infty k\cdot (u)_k(t-t_0)^{k-1}$$

Mit einer Indextransformation erhalten wir dann

$$u'(t)= \sum \limits_{k=0} \limits^\infty (k+1)\cdot (u)_{k+1} (t-t_0)^k$$

Potenzreihe darf nur innerhalb ihres Konvergenzkreises gliedweise differenziert werden, 1.Ableitung einer PR ANA I-Stoff

Verwenden wir nun die Abkürzung (Beziehung [3]) entsprechend für $u'$, so erhalten wir die Koeffizienten

$$(u')_k= \frac 1 {k!}\cdot (u')^{(k)} (t_0)= \frac 1 {k!}\cdot u^{(k+1)}(t_0)$$

und damit die Reihe

$$u'(t)= \sum \limits_{k=0} \limits^\infty (u')_k (t-t_0)^k .$$

Ein Koeffizientenvergleich liefert

$$[5] (u')_k=(k+1)\cdot (u)_{k+1} \quad \textrm{bzw.} \quad (u)_{k+1}=\frac 1 {k+1} \cdot (u')_k .$$

Mit dieser Rekursionsformel lässen sich die Ableitung von $u$ und damit auch die Koeffizienten rekursiv berechnen.

2.4.4 Beispiel

An einem sehr kleinen Beispiel will ich nun die rekursive Berechnung der Taylorkoeffizienten demonstrieren. Gegeben sei die AWA

$$u'=-u^2, \quad u(0)=1$$

An der Stelle $t_0=0$ ist also

$$(u)_0=1 \quad \textrm{und} \quad (u)_{k+1} = \frac 1 {k+1}... ...um \limits_{j=0} \limits^{k} (u)_j (u)_{k-j}, k\geq 0 .$$

Schreiben wir uns die Taylorkoeffizienten explizit hin, so sehen wir

$$\begin{eqnarray*} \Rightarrow \quad (u)_1&=& - (u)_0(u)_0=-1, (u)_2&=& - \... ..._k&=& (-1)^k \quad \textrm{durch vollständige Induktion}.\\ \end{eqnarray*}$$

Die letzte Beziehung bestätigt sich auch sofort, wenn man die exakte Lösung der AWA betrachtet. Sie lautet

$$u(t)=\sum \limits_{k=0}\limits^\infty (-1)^kt^k=\frac 1{1+t} .$$

$\frac 1{1+t}$ ist der Reihenwert der Summe $\sum \limits_{k=0} \limits^\infty (-1)^kt^k$. Für $t<1$ ist ist die Summe die geometrische Reihe.

Nun kommen wir zum eigentlichen Verfahren.

3 Einschließung bei gewöhnlichen AWA

In diesem Abschnitt entwickeln wir unser Einschließungsverfahren für die Lösung einer gewöhnlichen AWA. Als Grundlage dient uns das explizite Einschrittverfahren.

3.1 Einschrittverfahren

Das explizite Einschrittverfahren ist ein Verfahren mit dem man aus der Lösung $u$ zum Zeitpunkt $t_j$ eine Lösung zum Zeitpunkt $t_{j+1}$ approximieren kann.

Wir betrachten das o.a. Gitter $t_j=t_0+h\cdot j$ mit konstanter Schrittweite $h>0$. Sei $\Phi=\Phi(u)$ die Verfahrensfunktion des Einschrittverfahrens, so ist

$$[6]\quad u_{j+1}=u_j+h\cdot \Phi(u_j)+z_{j+1},\quad j\geq 0, \textrm { mit } u_0 \textrm{ aus AWA.}$$

$z_{j+1}$ ist dabei er lokale Fehler der im $(j+1)$-ten Schritt gemacht wird.

Die Gleichung [6] ist eine exakte Darstellung der Lösung $u(t)$ an den Gitterpunkten. Klar, es wird ja der exakte lokale Fehler $z_{j+1}$ verwendet. Dieser ist natürlich nicht bekannt.

Wir erhalten eine Approximation für den lokalen Fehler, in dem wir $\Phi$ so wählen, daß $u_j+h\Phi(u_j)$ gerade ein Taylorpolynom um $u_j$ mit Entwicklungsstelle $t_j$, ausgewertet an der Stelle $h$ ist.

Haben wir also ein Taylorpolynom vom Grad $(p-1)$, so hat der lokale Fehler eine Darstellung in der Form

$$[7]\quad z_{j+1}= \frac {h^p} {p!}\cdot u^{(p)} (\tau_{j+1}), \hspace{5cm}$$

mit einer unbekannten Zwischenstelle $\tau_{j+1}\in (t_j,t_{j+1})$.

Da wir die p-te Ableitung über $[t_j,t_{j+1}]$ aus unserer Abschätzung $[u^0]$ bzw. $[u^1]$ für $u$ über dem Intervall $[t_j,t_{j+1}]$ rekursiv berechnen können, erhalten wir eine Einschließung für den lokalen Fehler.

Im nächsten Abschnitt wollen wir nach dem lokalen Fehler entwickeln.

3.2 Entwicklung nach lokalen Fehlern

Hierzu betrachten wir die Beziehung [6] völlig isoliert, d.h. lediglich als eine Vorschrift, die den Vektoren $z_0, z_1, \ldots, z_{j+1}$ (zunächst noch unabhängige Variablen) einen Vektor $u_{j+1}$ zuordnet:

$$u_{j+1}=u_{j+1}(z_0, z_1,\ldots,z_{j+1}).$$

Da die Verfahrensfunktion $\Phi$ stetig differenzierbar sein soll, ist $u_{j+1}$ eine stetig differenzierbare Funktion der $j+2$ Variablen $z_0, z_1, \ldots, z_{j+1}$. Deshalb können wir $u_{j+1}$ nach dem MWS komponentenweise um (zunächst) willkürlich gewählte $s_0, s_1, \ldots, s_{j+1}$ entwickeln. Wir erhalten:

$$[8]\quad u_{j+1}=\tilde u_{j+1}+\sum \limits_{k=0} \limits^{j+1} \frac {\partial u_{j+1}(\hat z)}{\partial z_k}(z_k-s_k).$$

Dabei ist $\tilde u_{j+1}=u_{j+1}(s_0, s_1,\ldots,s_{j+1})$ der Funktionswert an der Entwicklungsstelle, d.h. es gilt

$$\tilde u_{j+1}=\tilde u_j + h \Phi(\tilde u_j) + s_{j+1}$$

wobei $\tilde{u}_0=u_0=z_0=s_0$ ist.

Das Argument $\hat z$ in [8] (2.2.2) soll die bei der Entwicklung auftretenden Zwischenstellen andeuten, welche für die $j+2$ Komponenten von $u_{j+1}$ i.a. paarweise verschieden sind.

Die Ableitungen in [8] (2.2.2) können aus [6] (2.1.1) rekursiv berechnet werden. Mit der Kettenregel ist nämlich

$$\frac{\partial u_{j+1}} {\partial z_k}= \frac{\partial u_j... ...}{\partial z_k} + \frac{\partial z_{j+1}}{\partial z_k} .$$

Daraus ergibt sich dann

$$[9]\quad \frac {\partial u_{j+1}} {\partial z_k}= \left\{ \... ...I,& \textrm { f\uml {u}r } k=j+1\\ \end{array} \right .$$

$I$ ist dabei die Einheitsmatrix und $\Phi '$ ist die Ableitungsmatrix von $\Phi$ ist.

Mit den Abkürzungen

$$[10] (2.2.5) \quad A_k:= I+h\Phi '(u_k) \textrm{ und }\ ... ...:=\frac{\partial u_{j+1}(\hat z)}{\partial z_k}, k\leq j+1$$

sieht man sofort aus [9](2.2.4), daß gilt:

$$A_{j+1,k}= A_{j} A_{j,k}=\ldots= A_j A_{j-1}\cdot \ldots \cdot A_k, \quad k\leq j+1,$$

Das leere Produkt soll im Fall $k=j+1$ wie üblich den Wert $I$ haben.

Anmerkung:

Bei diesen Abkürzungen so wie auch in [9] (2.2.4) ist das Argument $\hat z$ weggelassen worden. Da sich die Zwischenstellen von einem Zeitschritt zum nächsten ändern, hängen die Matrizen $A_k$ wegen diesen Zwischenstellen auch noch von dem Index $j$ ab.

Auf den weiteren Index wird verzichtet, weil die Matrizen $A_k, k\leq j+1$ später in Intervallmatrizen eingeschloßen werden. Im Einschließungsalgorithmus werden dann ausschließlich diese Intervallmatrizen verwendet. Die zusätzliche Indizierung erübrigt sich, da zur Berechnung dieser Intervallmatrizen für die Fehler $z_k$ später Einschließungen $[z_k]$ eingesetzt werden. Diese Einschließungen (Intervalle) hängen nicht vom Zeitindex $j$ ab. Die Entwicklungspunkte $s_k \in [z_k]$ werden ebenfalls unabhängig von $j$ gewählt.

Damit sind also die Intervallmatrizen unabhängig von $j$.

Also können wir [8] (2.2.2) wie folgt schreiben:

$$\begin{eqnarray*}[11]\ (2.2.7)\quad u_{j+1} &=& \tilde u_{j+1}+\sum \limits_{... ... \sum \limits_{k=0} \limits^j A_{j,k}(z_k-s_k)+z_{j+1}.\\ \end{eqnarray*}$$

Daran sieht man, daß die Entwicklungsstelle $s_{j+1}$ im Zeitschritt $t_j\to t_{j+1}$ noch nicht benötigt wird. D.h sie kann später festgelegt werden.

Nun wollen wir uns der Berechnung der Matrizen $A_j$ zuwenden.

3.3 Berechnung der $A_j$

Die Berechnung der Matrizen $A_j$ bzw. deren intervallmäßige Berechnung ist ein relativ aufwendiger Teil des Verfahrens. Der Autor formuliert das entsprechende Resultat als Satz.

Satz 3.1 ( II )
Das verwendete Einschrittverfahren sei die Taylorentwicklung mit Restglied der Ordnung $p$. Werten wir dann das Taylorpolynom vom Grad $p-1$ der Lösung der (Matrix-)Variationsgleichung

$$[12]\quad U'(t)=\frac{\partial f}{\partial u} \cdot U(t),\quad t_j\leq t \leq t_{j+1}, \quad U(t_j)=I,$$

mit Entwicklungsstelle $t_j$ an der Stelle $t_{j+1}$ aus, so erhalten wir gerade die Matrix $A_j$.

Beweis bleibt weg.

Dieser Satz soll also zeigen, daß die Matrizen $A_j$ auf die gleiche Weise berechnet werden können wie die Lösung $u$ selbst, nämlich durch Anwendung der rekursiven Formeln auf die obige Variationsgleichung [12].

Nun haben wir alles zusammen um das Verfahren zu implementieren.

Das Gerüst dieses Verfahrens kann in den folgenden fünf Schritten festgehalten werden. Wir nehmen an, daß zum Zeitpunkt $t_j$ die Lösungseinschließung $[u_j]$ bereits berechnet worden ist, was ja für $j=0$ erfüllt ist.

3.4 Gliederung des Verfahrens

2.2.8 a)

berechne $[u_{j+1}^0]$

2.2.8 b)

berechne mit $[u_{j+1}^0]$ .. $[z_{j+1}]$

2.2.8 c)

berechne mit $[u_j]$ .. $[A_j]$

berechne damit und mit $[z_{j+1}]$ .. $[u_{j+1}]$

2.2.8 d)

iteriere a)-d)

2.2.8 e)

Vorbereitung auf den nächsten Zeitschritt

Nach Beendigung des letzten Schrittes sind alle notwendigen Daten bereitgestellt, die die Fortpflanzung des Verfahrens auf den nächsten Zeitschritt erlauben, so daß wir das Verfahren so lange fortsetzen können, bis wir das Ende des Integrationsintervalles erreicht haben oder bis einer der Teilschritte nicht mehr durchgeführt werden kann. (meist 2.2.8 a) , wenn die Lösung z.B eine Singularität besitzt.)

Nun will ich die Grobgliederung nochmals ausführlicher erklären.

zu 2.2.8 a)

Hier muß ein Intervallvektor berechnet werden welcher die Lösung auf dem ganzen Intervall $[t_j,t_{j+1}]$ enthält. Dazu verwendet man Satz I . Um die dort notwendige Inklusion zu erreichen, hat es sich in der Praxis herausgestellt, daß man das Intervall mittels $\varepsilon$-Inflation aufblähen sollte, bis die Inklusion in [2] gilt.

Kann dieser Schritt nicht erfolgreich beendet werden, so kann man dies erzwingen indem man die Schrittweite $h$ genügend klein wählt.

Dies macht man in der folgenden Weise.

Man wähle ein Startintervall$[u_{j+1}^0]$ welches $[u_j]$ im Inneren enthält. Dann iteriert man

$$\begin{eqnarray*}[u_{j+1}^0]&:=&(1+\varepsilon)[u_{j+1}^0]-\varepsilon[u_{j+1}^0]\\ {[u_{j+1}^1]}&:=&[u_j]+[0,h]\cdot f([u_{j+1}^0]) . \\ \end{eqnarray*}$$

Dies macht man solange bis $[u_{j+1}^1]\subseteq [u_{j+1}^0]$ ist. Das Ergebnisintervall dieser Iteration wird dann wieder $[u_{j+1}^0]$ genannt. Wird die Inklusion erreicht, so ist dieser Schritt beendet, da Satz I (1.2.2.8) die Existenz der Lösung auf ganz $[t_j,t_{j+1}]$ im Intervall $[u_{j+1}^0]$ garantiert.

zu 2.2.8 b)

Die Einschließung des lokalen Fehlers berechnet sich nach der vorhin gezeigten Methode, d.h. wir werten den lokalen Fehler über $[u_{j+1}^0]$ aus Punkt a) aus.

Hier muß eine Einschließung für den lokalen Fehler berechnet werden, der sich durch die Lösung $u$ bzw. deren Ableitungen ausdrücken lassen muß. Im Fall der Taylorentwicklung werten wir die Darstellung des lokalen Fehlers [7] (2.1.2) über dem Intervall $[u_{j+1}^1]$ aus Punkt a) aus, was uns eine Einschließung $[z_{j+1}]$ von $z_{j+1}$ liefert.

zu 2.2.8 c)

Die Matrizen $A_j$ können nach Satz II auf die gleiche Weise berechnet werden wie die Lösung $u$ selbst, nämlich durch Anwendung der rekursiven Formeln aus dem Abschnitt über die Taylorentwicklung auf die obige Variationsgleichung [12].

Die intervallmäßige Berechnung der Einschließungen $[A_j]$ wird dann durch intervallmäßige Auswertung dieser Rekursionsformeln erreicht.

Also die Sache mit den $[A_j]$ ist nichts anderes, als die Auswertung des Taylorpolynoms vom Grad $p-1$.

zu 2.2.8 d)

Wie man sofort aus der zur DGL äquivalenten Integralform erkennen kann, liegt, wenn man die Integration bei $t_j$ beginnt, $u(t)$ in dem Streifen

$$u(t) \in [u_j]+(t-t_j)\cdot f([u_{j+1}^1]),\quad t\in[t_j,t_{j+1}]$$

und beginnt man mit $t_{j+1}$ liegt $u(t)$ ebenso in

$$u(t) \in [u_{j+1}]+(t-t_{j+1})\cdot f([u_{j+1}^1]),\quad t\in[t_j,t_{j+1}] .$$

Wir können nun den Durchschnitt dieser beiden Streifen bilden, welcher die Lösungen auch enthält. Liegt dieser Durchschnitt im Intervall $[u_{j+1}^0]$ aus Schritt a), so können wir dieses Intervall als neues Startintervall in Schritt a) verwenden. Dies kann in Schritt b) zu einer engeren Einschließung des lokalen Fehlers führen und damit wiederum zu einer besseren Einschließung $[u_{j+1}]$ von $u_{j+1}$ in c).

Ein solcher Schnitt liefert leider nicht sehr viel bessere Einschließungen als das der Aufwand für einen solchen Schnitt vertretbar wäre.

zu 2.2.8 e)

Übergang zur nächsten Entwicklungsstelle. In diesem abschließenden Schritt des Verfahrens berechnen wir die Entwicklungsstelle $s_{j+1}$, mit der die Näherung $\tilde u_{j+1} = \tilde u_j + h \Phi(\tilde u_j)+ s_{j+1}$ berechnet werden muß, die dann im nächsten Zeischritt verwendet wird.

$s_{j+1}$ kann als beliebiger Punkt aus dem lokalen Fehler-Intervall $[z_{j+1}]$ das ihn Schritt b) berechnet wurde, gewählt werden (z.B. der Mittelpunkt).

WICHTIG: $s_{j+1}$ muß aus diesem Intervall gewählt wird, da mit diesem Intervall in c) die Einschließung der Matrix $A_j$ berechnet wird.

3.5 Behandlung des globalen Fehlers

In diesem letzten Abschnitt will ich nur kurz verschiedene Methoden erwähnen, mit denen man den globalen Fehler mittels der berechneten Näherungslösung $\tilde u_{j+1}$ gut einschließen kann. Der globale Fehler ist vom Autor definiert als die Differenz $r_{j+1}$ zwischen exakter Lösung $u_{j+1}$ und Näherungslösung $\tilde u_{j+1}$.

$$r_{j+1} = u_{j+1} + \tilde u_{j+1}$$

Dieser Teil des Einschließungsverfahrens ist wohl am schwierigsten allgemein beherrschbar, da bei verschiedenen DGL oft nur verschiedene Methoden gute Einschließungen liefern und auch bei ein und derselben DGL ist oftmals das Verhalten der speziell untersuchten Lösung für die Wahl der Einschließungsmethode entscheidend.

Wie man aus der Gleichung [8] sehen kann, kann man $u_{j+1}$ folgendermaßen schreiben:

$$\begin{eqnarray*} u_{j+1} &=& \tilde u_{j+1} + r_{j+1}, \quad {\rm wobei}\\ ... ...imits ^{j+1} A_{j+1,k} (z_k-s_k) = A_j r_j+(z_{j+1}-s_{j+1}) \end{eqnarray*}$$

Die Differenz $r_{j+1}$ zwischen exakter Lösung $u_{j+1}$ und Näherungslösung $\tilde u_{j+1}$ ist der globale Fehler, den wir genauer betrachten wollen. Wie man deutlich am letzten Gleichheitszeichen in der Darstellung des Fehlers sehen kann, haben wir es hier mit dem Problem der Fehlerfortpflanzung zu tun. D.h. mit der Frage, wie der im Intervall $[t_{j-1},t_j]$ vorliegende globale Fehler $r_j$ sich auf den globalen Fehler $r_{j+1}$ im Intervall $[t_j,t_{j+1}]$ überträgt, bzw. wie diese Übetragung am besten rechnerisch in den Griff zu bekommen ist.

Hierfür gibt der Autor verschiedene Auswertungsmethoden an, die ich nur mit Namen erwähnen will.

• 1) Intervallvektor-Auswertung
• 2) Parallelepiped-Auswertung
• 3) QR-Zerlegung-Auswertung

In der nun folgenden Tabelle sind die drei verschieden Auswertungen verglichen worden.

3.6 Ergebnisvergleich der verschiedenen Auswertungsmethoden

Die AWA lautet

$$\begin{eqnarray*} u'_1 &=& u_1-2u_2\quad,\quad u_1(0)=1,\\ u'_2 &=& 3u_1-4u_2\quad,\quad u_2(0)=-1\\ \end{eqnarray*}$$

und besitzt die exakte Lösung

$$\begin{eqnarray*} u_1(t) &=& 5 e^{-t}-4 e^{-2t},\\ u_1(t) &=& 5 e^{-t}-6 e^{-2t} \end{eqnarray*}$$

Die in der unten stehenden Tabelle angegebenen Ergebnisse wurden mit der Schrittweite $h=0.1$ und der Ordnung $p=14$ gerechnet. Es ist jeweils der größere der Durchmesser der Einschließungen beider Lösungskomponenten angegeben.

t	1	2	3
10	2.2E-05	5.9E-10	6.6E-14
20	7.8e+01	6.8E-10	6.7E-18
27	---	4.0E+28	8.3E-21
40	---	---	2.8E-26
10	---	---	6.2E-52
150	---	---	1.8E-73
200	---	---	4.7E-95
250	---	---	8.0E-98

1, 2, 3 bezeichnen die verschiedenen Auswertungen.