Teil 1 - Technische Universität Chemnitz

Skript zur Vorlesung
Numerik gewöhnlicher
Differentialgleichungen
Dr. M. Weise
gehalten im WS2015/16
Technische Universität Chemnitz
Website zur Vorlesung:
http://www-user.tu-chemnitz.de/~wmicha/w15num_ode/
Dieses Vorlesungsskript wurde im WS 2006/07 erstellt von Roland Herzog unter
Benutzung früherer Vorlesungen und Übungen von Prof. Peter Benner, Prof. Rüdiger
Weiner, Prof. Joachim Stöckler, Prof. Hans Josef Pesch und Prof. Rolf Rannacher,
denen ich herzlich für das Zurverfügungstellen der Unterlagen danke.
Im WS 2007/08, WS 2012/2013 und WS 2013/2014 wurde es durch Rene Schneider
teilweise überarbeitet und erweitert. Im WS 2010/11 hat Herbert Egger die Vorlesung gehalten und viele Hinweise zur Überarbeitung gegeben, die im WS 2011/12
zu einer wesentlichen Verbesserung geführt haben.
Material für: 20–21 Vorlesungen à 90 Minuten
Fehler und Kommentare bitte an: [email protected]
Stand: 9. Oktober 2015
Inhaltsverzeichnis
Kapitel 0. Einführung und Motivation
5
1
Aufgabenstellung, Beispiele und Ziele der Vorlesung
5
2
Wiederholung und Notation
8
Kapitel 1. Numerische Lösung von Anfangswertproblemen
15
3
Einschrittverfahren
3.1
Grundbegriffe und ein erstes Verfahren
3.2
Konvergenzanalysis allgemeiner Einschrittverfahren
3.3
Runge-Kutta-Verfahren
3.3.1
Allgemeine Ordnungsbedingungen
3.3.2
Explizite Runge-Kutta-Verfahren
3.3.3
Motivation impliziter Runge-Kutta-Verfahren
3.3.4
Implizite Runge-Kutta-Verfahren
3.3.5
Kollokationsverfahren und IRKV
3.3.6
Einige Klassen von IRKV
3.4
Stabilitätsbegriffe bei Einschrittverfahren
3.4.1
Stabilitätsgebiet und A-Stabilität
3.4.2
Steife Differentialgleichungen
3.4.3
L-Stabilität
3.5
Symplektische Verfahren
3.6
Praktische Aspekte von Einschrittverfahren
3.6.1
Schrittweitensteuerung
3.6.2
Linear implizite Runge-Kutta-Verfahren
3.6.3
Verfahren mit dichter Ausgabe
3.7
Einschrittverfahren für DGL zweiter Ordnung
15
15
22
33
38
42
46
48
54
57
61
61
67
70
74
79
79
91
96
98
4
Mehrschrittverfahren
4.1
Einleitung und Grundbegriffe
4.2
Konvergenzuntersuchung bei Mehrschrittverfahren
4.3
Stabilitätsbegriffe bei Mehrschrittverfahren
4.4
Praktische Aspekte bei Mehrschrittverfahren
103
103
117
134
141
5
Exponentielle Integratoren
145
6
Discontinuous-Galerkin-Verfahren
148
Kapitel 2. Weiterführende Themen
7
161
RWP mit gewöhnlichen Differentialgleichungen
8
Differentiell-algebraische Systeme
8.1
Einführung
8.2
Eigenschaften linearer DAE-Systeme
3
161
167
167
169
4
Kapitel 0. Inhaltsverzeichnis
8.3
9
Numerische Behandlung linearer DAE-Systeme
Anfangsrandwertprobleme
Literaturverzeichnis
182
188
195
KAPITEL 0
Einführung und Motivation
Inhalt
§1
1
Aufgabenstellung, Beispiele und Ziele der Vorlesung
5
2
Wiederholung und Notation
8
Aufgabenstellung, Beispiele und Ziele der Vorlesung
Wir betrachten in dieser Vorlesung hauptsächlich Anfangswertprobleme für Systeme gewöhnlicher Differentialgleichungen (Dgl) erster Ordnung:

Gesucht ist eine differenzierbare Funktion y : [0, T ] → Cn ,

sodass y 0 (t) = f (t, y(t)) für alle t ∈ [0, T ]
(AWP)


und y(0) = ya gilt.
Dabei sind die rechte Seite f : R × Cn → Cn , der Anfangswert (AW) ya ∈ Cn ,
die Anfangszeit t = 0 und die Endzeit t = T gegeben. Die unabhängige Variable
t spielt in vielen, aber nicht in allen Problemen die Rolle der Zeit. Die Ableitung
y 0 (t) wird dann oft mit ẏ(t) bezeichnet. Die gesuchte Funktion y(t) heißt auch der
Zustand der Dgl. Ist die rechte Seite f von t unabhängig, so heißt das Problem
autonom. Meist beschränken wir uns auf reelle Zustände, also
ya ∈ Rn
und
f : R × Rn → Rn .
Zahlreiche Fragestellungen führen auf Probleme von diesem Typ.
Beispiel 1.1 (Räuber-Beute-Modell)
[Heuser, 1991, § 59] Wir nehmen an, eine Population y1 (Beutetiere) lebe von einer
ausreichend vorhandenen Nahrungsquelle. Der einzige natürlich Feind ist eine Population y2 von Raubtieren, die auf y1 angewiesen sind. Ein einfaches Modell für die
Populationsentwicklung ist gegeben durch
Dabei sind α1 > 0 und α2 > 0 die natürliche Vermehrungs- bzw. Verminderungsraten
der Beute bzw. der Räuber in Abwesenheit der jeweils anderen Population. Der
5
6
Kapitel 0. Einführung und Motivation
Wechselwirkungsterm geht von einer Begegnungshäufigkeit von Räuber und Beute
aus, die über Konstanten β1 , β2 > 0 proportional zu beiden Beständen ist. Ein
typischer Lösungsverlauf ist in Abbildung 1.1 zu sehen.
Predator−Prey model with α1 = 1, α2 = 2, β1 = β2 = 0.05
120
prey
predator
100
population
80
60
40
20
0
0
2
4
6
8
10
time
12
14
16
18
20
Abbildung 1.1. Verlauf der Lösung aus Beispiel 1.1 (Räuber-BeuteModell) für die Parameterwahl α1 = 1, α2 = 2, β1 = β2 = 0.05 zum
Anfangswert y(0) = (100, 10)> .
Differentialgleichungen höherer Ordnung lassen sich durch Einführung zusätzlicher
unbekannter Funktionen auf Systeme erster Ordnung reduzieren:
Beispiel 1.2 (Newtonsches Bewegungsgesetz)
Nach dem Newtonschen Bewegungsgesetz gilt für einen Körper der Masse m die
Beziehung F = ma zwischen der auf ihn wirkenden Kraft F und der resultierenden
Beschleunigung a. Bezeichnen wir mit s(t) ∈ R3 den Ort des Körpers, so gilt s00 (t) =
a(t) und damit das Dgl-System zweiter Ordnung s00 (t) = F (t)/m. Durch Einführung
der Geschwindigkeit v(t) = s0 (t) ∈ R3 als Zwischengröße können wir ein äquivalentes
Dgl-System erster Ordnung in y(t) = (s(t), v(t))> ∈ R6 formulieren:
Beim Anfangswertproblem sind dann AW für Ort und Geschwindigkeit gegeben. Das (AWP) wird man nur in seltenen Fällen analytisch lösen können. Wir behandeln
deshalb in dieser Vorlesung numerische Verfahren für die näherungsweise Lösung
von (AWP). Um die Brauchbarkeit der numerischen Lösung beurteilen zu können,
benötigen wir Aussagen über ihre Genauigkeit. Ziele der Vorlesung sind
• das Kennenlernen von Verfahren zur numerischen Lösung von (AWP),
• die Analysis dieser Verfahren bzgl. Konvergenz und Stabilität,
• deren praktische Umsetzung in Computerprogramme sowie
§ 1. Aufgabenstellung, Beispiele und Ziele der Vorlesung
7
• Techniken der Schrittweitensteuerung zur effizienten numerischen Lösung.
Zwei weitere Beispiele zur Anwendung:
Beispiel 1.3 (Schwingung einer Membran)
wird modelliert durch die partielle DGL
0 = ρ(x)utt (x, t) + b(x) ut (x, t) − T0 ∆u(x, t) − f (x, t)
∀t ≥ 0, ∀x ∈ Ω,
(1.1)
wobei u(x, t) Auslenkung der Membran über dem Gebiet Ω ⊂ R2 und Zeit t, ρ(x)
Dichte der Membran, b(x) Dämpfungskoeffizient, T0 die Vorspannung der Membran
und f (x, t) eine von außen wirkende Kraftdichte bezeichnen. Hinzu kommen Randbedingungen
u(x, t) = d(x, t),
∀t ≥ 0, ∀x ∈ ∂Ω,
(1.2)
und Anfangsbedingungen
u(x, 0) = g0 (x)
ut (x, 0) = g1 (x)
∀x ∈ Ω
∀x ∈ Ω.
Eine Semidiskretisierung im Ort (z. B. mit Finite-Elemente-Methode) führt auf ein
System gewöhnlicher DGLn zweiter Ordnung
0 = M ü + C u̇ + Ku − f (t),
(1.3)
wobei u(t) Vektor der Funktionswerte u(xi , t) in den Stützstellen der Ortsdiskretisierung ist, M die Massematrix, C Dämpfungsmatrix, K Steifigkeitsmatrix der
Membran und f (t) aus den äußeren Kräften resultiert. Die Dimension n der Vektoren und Matrizen ist dabei i. A. sehr groß, 103 ≤ n ≤ 107 sind heute typische
Werte.
Ein Beispiel ist in Abbildung 1.2 gegeben.
Abbildung 1.2. Beispiel einer Lösung Beispiel 1.3 (Membran) für
eine FEM-Diskretisierung mit n = 24.576 Freiheitsgraden auf Gebiet
Ω = (−1, 1)2 \ (0, 1) × (−1, 0) (L-förmig), ρ ≡ 1, b ≡ 1, T0 = 1, f ≡ 0,
g0 ≡ g1 ≡ 0, d(x, t) „periodisches Kippen eines Teilrandes“.
Auch in der Umsetzung von Steuerungs- und Regelungsaufgaben für gewöhnliche
Differentialgleichungen werden die hier zu besprechenden Verfahren benötigt.
8
Kapitel 0. Einführung und Motivation
Abbildung 1.3. Legway: Modell eines selbst-balancierenden, einachsigen Fahrzeugs.
Beispiel 1.4 (Legway)
(Ilka Riedel) Siehe Abbildung 1.3.
§2
Wiederholung und Notation
Wir wiederholen kurz einige wichtige Ergebnisse aus der Theorie der gewöhnlichen
Differentialgleichungen sowie der Analysis. Hier und in Zukunft bedeutet kyk für
einen Vektor y ∈ Rn immer die Euklidische Norm, also kyk = (y > y)1/2 . (Wir könnten
auch jede andere Norm verwenden.)
Satz 2.1 (Picard-Lindelöf )
[Heuser, 1991, Satz 12.2, Aufgabe 12.6] Es sei
Q = [a1 , b1 ] × [a2 , b2 ] × · · · × [an , bn ] ⊂ Rn
ein kompakter Quader und ya ∈ int(Q). Weiterhin sei f : [0, T ] × Q → Rn stetig
und bzgl. y (gleichmäßig) Lipschitz-stetig, d. h., es gibt ein L > 0, sodass
kf (t, y1 ) − f (t, y2 )k ≤ L ky1 − y2 k
(2.1)
gilt für alle t ∈ [0, T ] und alle y1 , y2 ∈ Q.
Satz 2.2 (Stetige Abhängigkeit von den Anfangswerten)
Unter den Voraussetzungen von Satz 2.1 gilt für die Lösung y(t; s) des (AWP)
y 0 (t) = f (t, y(t)),
y(0) = s
zu verschiedenen AW s1 , s2 ∈ int(Q) die Abschätzung
ky(t; s1 ) − y(t; s2 )k ≤ eL t ks1 − s2 k
für alle t ≥ 0 aus dem gemeinsamen Intervall der Lösungen zu s1 und s2 .
(2.2)
Der Beweis von Satz 2.2 läuft über das Lemma von Gronwall. Wir geben hier eine
recht allgemeine Version dieses Resultats an, das später (in diskreter Form) bei der
Stabilitätsuntersuchung wieder eine Rolle spielt.
§ 2. Wiederholung und Notation
9
Exponential deviation of solutions for different initial values. y’(t) = y + 0.02
1.5
1
y(t)
0.5
0
−0.5
−1
0
0.5
1
1.5
time t
2
2.5
3
Abbildung 2.1. Illustration der exponentiellen Abweichung (2.2)
über die Zeit bei der Lösung eines (AWP) zu verschiedenen Anfangswerten.
Lemma 2.3 (Lemma von Gronwall (Thomas Hakon Grönwall, 1877–1932))
Es seien ϕ, L, C : [0, T ] → R nicht-negative stetige Funktionen. Die Funktion C sei
monoton wachsend. Wenn die Abschätzung
erfüllt ist, dann gilt:
0 ≤ ϕ(t) ≤ C(t) exp
Z
t
L(s) ds
0
für alle t ∈ [0, T ].
(2.3)
Insbesondere im Fall C(t) ≡ 0 ist also ϕ(t) ≡ 0. Falls L konstant ist, so vereinfacht
sich (2.3) zu
0 ≤ ϕ(t) ≤ C(t) eL t
für alle t ∈ [0, T ].
(2.4)
Beweis: in Übung 1, Aufgabe 2
Voraussetzung 2.4 (Generalvoraussetzung)
Wir betrachten das (AWP). Wir wollen im Folgenden stets voraussetzen, dass . . .
10
Kapitel 0. Einführung und Motivation
Bemerkung 2.5 (Hinreichende Bedingung für Lipschitz-Stetigkeit)
Es sei f : [0, T ]×Q → Rn stetig. Falls f außerdem stetige erste partielle Ableitungen
∂f /∂yi , i = 1, . . . , n besitzt, so ist die Lipschitz-Bedingung (2.1) mit
erfüllt.1
Die Fehlerabschätzungen „leben“ später von der Glattheit (Differentiationsordnung)
der Lösung. Diese erhält man aus folgendem Satz:
Satz 2.6 (Maximale Regularität)
Es sei y : [0, T ] → Q ⊂ Rn eine Lösung des (AWP) auf [0, T ]. Falls f ∈ C k ([0, T ]×Q)
ist für ein k ∈ N0 , dann gilt y ∈ C k+1 ([0, T ]; Rn ).
Beweis:
Wir verwenden später häufig Abschätzungen, die auf dem Satz von Taylor bzw. dem
daraus folgenden „Schrankensatz“ basieren. Da wir es mit vektorwertigen Funktionen
zu tun haben, formulieren wir den Satz von Taylor in Integralform.
Satz 2.7 (Satz von Taylor in Integralform)
Es sei G ⊂ Rn offen und F : G → Rm (k + 1)-mal stetig partiell differenzierbar
(diffbar). Falls die Punkte x und x + d und die gesamte Verbindungsstrecke in G
1Dabei
ist kfy (t, y)k die durch die Euklidische Vektornorm induzierte Matrixnorm (Spektralnorm) der Jacobimatrix fy (t, y).
§ 2. Wiederholung und Notation
liegen, dann gilt:
k=0:
k=1:
Z
11
1
F 0 (x + s d) d ds
0
Z 1
0
(1 − s) F 00 (x + s d) [d, d] ds
F (x + d) = F (x) + F (x) d +
F (x + d) = F (x) +
0
k ∈ N0 : F (x + d) =
k
X
F (j) (x) [d(j) ]
j=0
j!
+
Z
1
0
(1 − s)k (k+1)
F
(x + s d) [d(k+1) ] ds.
k!
(Siehe [Heuser, 2002, Satz 168.4] für den Fall k = 1.)
Dabei ist die Jacobimatrix F 0 (x) an der Stelle x als eine 1-Form (lineare Form)
zu verstehen, die Rn nach Rm abbildet. Dementsprechend ist die j-te Ableitung
F (j) (x) an der Stelle x eine multilineare Abbildung, die [Rn ]j nach Rm abbildet. Die
Koeffizienten von F (j) (x) werden gebildet von den partiellen Ableitungen von F der
Ordnung j an der Stelle x. Der Ausdruck F (j) (x)[d(j) ] bedeutet gerade das j-fache
Einsetzen der Richtung d in diese multilineare Abbildung, die übrigens (Satz von
Schwarz) symmetrisch ist.
Folgerung 2.8 (Schrankensatz)
Es seien die Voraussetzungen des Satzes 2.7 erfüllt. Weiter seien alle partiellen Ableitungen von F der Ordnung k + 1 in G beschränkt. Dann existiert eine Konstante
Ck+1 , sodass die folgende Abschätzung für alle x, x + d ∈ G erfüllt ist:
In die Konstante geht die Ableitung F (k+1) ein.
Wir werden diese Sätze oft mit x = (t, y) anwenden. Mit Hilfe des Schrankensatzes
für k = 0 kann man z. B. Bemerkung 2.5 beweisen. Wir schreiben auch
k = 0 : F (x + d) = F (x) + O(kdk)
k = 1 : F (x + d) = F (x) + ∇F (x)> d + O(kdk2 )
etc. mit Hilfe der folgenden Notation:
Definition 2.9 (O- und o-Notation)
(a) Man sagt, eine Funktion f : R+ → Rn verhält sich wie O(g(h)) mit einer
anderen Funktion g : R+ → R+ (sprich: „groß O von g von h“) für h & 0,
falls h0 > 0 und eine von h unabhängige Konstante c existieren, sodass gilt:
kf (h)k ≤ c g(h) gilt für alle 0 < h ≤ h0 .
12
Kapitel 0. Einführung und Motivation
(b) Eine Funktion f : R+ → Rn verhält sich wie o(g(h)) mit einer anderen
Funktion g : R+ → R+ (sprich: „klein o von g von h“) für h & 0, falls gilt:
kf (h)k
= 0.
h&0 g(h)
lim
Oft schreibt man auch kurz: f (h) = O(g(h)) oder f (h) ∈ O(g(h)) bzw. f (h) =
o(g(h)) oder f (h) ∈ o(g(h)).
Beispiel:
Satz 2.10 (Banachscher Fixpunktsatz)
[Heuser, 2002, Satz 111.11] Es sei X eine nichtleere, abgeschlossene Teilmenge des
Rn und A : X → X eine kontrahierende Selbstabbildung von X; für alle x, y ∈ X
gelte also
kA x − A yk ≤ q kx − yk mit einem festen 0 < q < 1.
Dann besitzt A genau einen Fixpunkt x∗ in X. Dieser kann iterativ gewonnen werden: Wählt man einen beliebigen Startwert x0 ∈ X und setzt
xn+1 := A xn
n = 0, 1, 2, . . . ,
so konvergiert xn → x∗ . Überdies gilt die Fehlerabschätzung
Der Satz bleibt richtig, wenn man statt Rn einen allgemeinen Banachraum mit der
Norm k·k zugrunde legt.
Bemerkung 2.11 (Autonomie-Transformation)
Eine Dgl erster Ordnung der Form y 0 (t) = f (y(t)) heißt autonom. Eine nichtautonome Dgl y 0 (t) = f (t, y(t)) kann durch Einführung einer zusätzlichen Zustandsvariablen immer auf ein autonomes Problem zurückgeführt werden: Man setzt
welches dann die Dgl
§ 2. Wiederholung und Notation
13
erfüllt.
KAPITEL 1
Numerische Lösung von Anfangswertproblemen
Inhalt
3
Einschrittverfahren
3.1
Grundbegriffe und ein erstes Verfahren
3.2
Konvergenzanalysis allgemeiner Einschrittverfahren
3.3
Runge-Kutta-Verfahren
3.3.1
Allgemeine Ordnungsbedingungen
3.3.2
Explizite Runge-Kutta-Verfahren
3.3.3
Motivation impliziter Runge-Kutta-Verfahren
3.3.4
Implizite Runge-Kutta-Verfahren
3.3.5
Kollokationsverfahren und IRKV
3.3.6
Einige Klassen von IRKV
3.4
Stabilitätsbegriffe bei Einschrittverfahren
3.4.1
Stabilitätsgebiet und A-Stabilität
3.4.2
Steife Differentialgleichungen
3.4.3
L-Stabilität
3.5
Symplektische Verfahren
3.6
Praktische Aspekte von Einschrittverfahren
3.6.1
Schrittweitensteuerung
3.6.2
Linear implizite Runge-Kutta-Verfahren
3.6.3
Verfahren mit dichter Ausgabe
3.7
Einschrittverfahren für DGL zweiter Ordnung
§3
15
15
22
33
38
42
46
48
54
57
61
61
67
70
74
79
79
91
96
98
4
4.1
4.2
4.3
4.4
Mehrschrittverfahren
Einleitung und Grundbegriffe
Konvergenzuntersuchung bei Mehrschrittverfahren
Stabilitätsbegriffe bei Mehrschrittverfahren
Praktische Aspekte bei Mehrschrittverfahren
103
103
117
134
141
5
Exponentielle Integratoren
145
6
Discontinuous-Galerkin-Verfahren
148
Einschrittverfahren
§ 3.1
Grundbegriffe und ein erstes Verfahren
Wir erinnern noch einmal an unsere Generalvoraussetzung 2.4. Die exakte Lösung
des (AWP) bezeichnen wir mit y(t). Alle in diesem Kapitel behandelten numerischen Verfahren bestimmen Näherungswerte der exakten Lösung auf einem Punktgitter.
15
16
Kapitel 1. Numerische Lösung von Anfangswertproblemen
Definition 3.1 (Gitter, Gitterfunktion)
(a) Ein Gitter auf [0, T ] ist eine endliche Punktmenge
T = {t0 , t1 , . . . , tN } mit 0 = t0 < t1 < · · · tN = T.
Die {ti } heißen auch die Stützstellen. Die Abstände hi = ti+1 − ti heißen
die Schrittweiten des Gitters T .
(b) Für ein gegebenes Gitter bezeichnet
die (maximale) Gitterweite oder den Diskretisierungsparameter.
Man sagt dann: Das Gitter T gehört1 zur Klasse Th .
(c) Ist hi = h für alle i, so heißt das Gitter äquidistant.
(d)
Die Differentialgleichung
y 0 (t) = f (t, y(t))
besagt, dass die gesuchte Lösung y(t) an der Stelle t die Tangentensteigung f (t, y(t))
hat. Mit anderen Worten, die Steigung der Lösung y(t) passt sich in das durch f (t, y)
definierte Richtungsfeld ein. Für skalare Differentialgleichungen, also im Fall n = 1,
lässt sich das einfach grafisch veranschaulichen (Abbildung 3.1 links).
Das einfachste numerische Verfahren besteht nun darin, die Lösungskurve y(t) im
Intervall [t, t + h] durch ihre Tangente im Punkt (t, y(t)) zu approximieren, deren
Steigung aus dem Richtungsfeld, d. h. aus der rechten Seite f (t, y(t)) abgelesen werden kann.
Auf einem Gitter T mit Stützstellen ti erhalten wir sukzessive die Näherungen yi =
yh (ti ) für die Werte der exakten Lösung y(ti ) gemäß
1T
gehört dann auch zur Klasse Th0 für alle h0 ≥ h.
17
§ 3. Einschrittverfahren
Slopefield for f(t,y) = 1 − t * y and solution for y(0) = −0.8
1
0.8
0.6
y(t+h)
0.4
Value y
0.2
y(t)
0
−0.2
−0.4
yh(t+h)
=yh(t)
−0.6
−0.8
−1
0
1
2
3
Time t
4
5
t
t+h
Abbildung 3.1. Richtungsfeld einer Dgl (links) und Konstruktionsidee des expliziten Euler-Verfahrens (Polygonzugmethode) (rechts).
Die Näherungsvorschrift (3.1) heißt auch Methode von Euler, explizites EulerVerfahren oder Polygonzugmethode. In der englischsprachigen Literatur wird
das Verfahren auch als Vorwärts-Euler-Verfahren (Forward Euler) bezeichnet.
Eine wichtige Fragestellung ist die Untersuchung der Konvergenz und Konvergenzgeschwindigkeit der numerischen Lösungen, wenn die Gitterweite h & 0 strebt. Wir
wollen dies zunächst exemplarisch für das explizite Euler-Verfahren untersuchen.
Dazu bezeichnen wir für ein gegebenes Gitter T mit
den (globalen) Fehler des Näherungsverfahrens am Gitterpunkt ti ∈ T . Auch eh
ist eine Gitterfunktion.
Die Abbildung 3.2 zeigt, dass sich der Fehler ei+1 an einem Gitterpunkt ti+1 aus
zwei Anteilen zusammensetzt:
18
Kapitel 1. Numerische Lösung von Anfangswertproblemen
=yh(t 2 )
=yh(t 1 )
z2
=yh(t 0 )
y(t 1 )
y(t 2 )
t1
t2
y(t 0 )
t0
Abbildung 3.2. Fehlerfortpflanzung bei einem Einschrittverfahren.
Wir schätzen beide Anteile ab.
(a) Zur Abschätzung des lokal hinzukommenden Fehlers definieren wir den sogenannten lokalen Diskretisierungsfehler di+1 des Euler-Verfahrens an
der Stelle ti+1 über
Man nennt di+1 auch den Konsistenzfehler (an der Stelle ti+1 ), da er sich
aus dem Einsetzen der exakten Lösung in das diskrete Schema (hier (3.1))
ergibt.
Um den lokalen Diskretisierungsfehler abzuschätzen, verwenden wir die Taylorentwicklung (Satz 2.7)
§ 3. Einschrittverfahren
19
Hierfür setzen wir voraus, dass y in C 2 ([0, T ]; Rn ) liegt.2 Wir setzen in (3.2)
ein und erhalten die Fehlerdarstellung
(3.3)
Wegen der Beschränktheit von y 00 erhalten wir die Abschätzung
(3.4)
(b) Zur Untersuchung der Fortpflanzung des Fehlers ei müssen wir am Zeitpunkt ti zwei explizite Eulerschritte betrachten zu den Startwerten yh (ti )
und y(ti ), die nach Definition gerade um ei = yh (ti ) − y(ti ) auseinander
liegen:
Deren Differenz lässt sich abschätzen als:
(3.5)
2Hinreichend
dafür ist nach Satz 2.6 f ∈ C 1 ([0, T ] × Q).
20
Kapitel 1. Numerische Lösung von Anfangswertproblemen
Hier setzen wir voraus, dass die rechte Seite f in einer geeigneten Menge
Lipschitz-stetig bzgl. y mit L-Konstante L ist.3
Zusammen erhalten wir mit (3.2) und (3.5) also die Abschätzung
für den globalen Fehler. Mit (3.6) haben wir eine Rekursion für die Fehlerfortpflanzung von Schritt zu Schritt gewonnen. Wir bringen diese durch ein Teleskopsummenargument 4 auf die Form
Zur expliziten Abschätzung benutzen wir das folgende diskrete Lemma von Gronwall
mit den Setzungen
Lemma 3.2 (Diskretes Lemma von Gronwall, vgl. Lemma 2.3)
Es seien {ϕn }, {Ln }, {Cn } nicht-negative Zahlenfolgen, n ∈ N0 . Die Folge {Cn } sei
monoton wachsend. Wenn die Abschätzung
0 ≤ ϕ n ≤ Cn +
3Genauer:
n−1
X
j=0
Lj ϕj
für alle n ∈ N0
Die diskreten Lösungen {yh (ti )} müssen auch in Q liegen.
Einsetzen kei k in kei+1 k, beginnend bei ke0 k
4wiederholtes
§ 3. Einschrittverfahren
21
erfüllt ist, dann gilt:
(3.8)
Insbesondere im Fall Cn ≡ 0 ist also ϕn ≡ 0. Falls Lj ≡ L ist, so vereinfacht sich
(3.8) zu
(3.9)
Beweis: eigene Übung
Durch die Anwendung des Lemmas ergibt sich aus (3.7) die Abschätzung:
Durch Einsetzen der Abschätzung (3.4) für den lokalen Fehler, also
ergibt sich unter Beachtung von
Pn
j=1
hj−1 = tn folgender Konvergenzsatz:
Satz 3.3 (Konvergenz des expliziten Euler-Verfahrens)
Die Lösung des (AWP) sei C 2 ([0, T ]; Rn ). Dann gilt für den Fehler der Näherungslösungen des expliziten Euler-Verfahrens auf einem beliebigen Gitter {0 =
t0 , t1 , . . . , tN = T } der Klasse Th die Abschätzung
(3.10)
22
Kapitel 1. Numerische Lösung von Anfangswertproblemen
für alle n = 0, 1, . . . , N und damit
Bemerkung 3.4
Wesentliche Bestandteile des Konvergenzbeweises am Beispiel des expliziten EulerVerfahrens sind:
Im nächsten Abschnitt werden wir dieses Vorgehen auf allgemeine Einschrittverfahren übertragen.
§ 3.2
Konvergenzanalysis allgemeiner Einschrittverfahren
Als Verallgemeinerung von (3.1) betrachten wir im Rest von § 3 Verfahren der Bauart
§ 3. Einschrittverfahren
23
(3.12)
die aus der Kenntnis eines Näherungswertes yi für y(ti ) und der Schrittweite hi einen
neuen Näherungswert yi+1 für y an der Stelle ti+1 = ti + hi bestimmen.
Definition 3.5 (Ein(zel)schrittverfahren)
Ein numerisches Verfahren zur Lösung des (AWP), das der Vorschrift (3.12) genügt, heißt Einschrittverfahren (ESV). Φ heißt die Verfahrensfunktion oder
Inkrementfunktion des Verfahrens.5
Wie wir bereits in § 3.1 gesehen haben, spielt der Konsistenzfehler eine wesentliche
Rolle. Er bezieht sich immer auf das Einsetzen der exakten Lösung in das Verfahren.
Definition 3.6 (Konsistenzfehler, Konsistenzordnung bei ESV)
Es sei y(·) die Lösung der Dgl y 0 (t) = f (t, y(t)) auf dem Intervall [0, T ].
(a) Die Größe
(3.13)
heißt der Konsistenzfehler oder Abschneidefehler des ESV (3.12) an
der Stelle t zur Schrittweite h.
Später werden wir auch andere (beliebige) Funktionen z in d einsetzen,
d(z, ., .), um Fehlerfortpflanzung zu untersuchen.
(b) Das ESV heißt konsistent mit dem (AWP), falls gilt:
(3.14)
5Natürlich
hängt die Verfahrensfunktion Φ auch von der rechten Seite f ab (sonst würde ja das
Verfahren die Dgl gar nicht berücksichtigen), aber diese Abhängigkeit drücken wir nicht explizit
aus.
24
Kapitel 1. Numerische Lösung von Anfangswertproblemen
(c) Das ESV besitzt die Konsistenzordnung p für das (AWP), falls
Die Konsistenz und Konsistenzordnung beziehen sich dabei jeweils nur auf die Differentialgleichung des (AWP) auf [t, t + h], es werden hier exakte Anfangswerte bei
t verwendet. Aus (3.13) folgt, dass ein ESV offenbar genau dann konsistent mit dem
(AWP) ist, wenn gilt:
lim Φ(t, y(t), h) = f (t, y(t)) gleichmäßig in [0, T − h].
h&0
Wie bereits beim Euler-Verfahren gesehen, bestimmt man die Konsistenzordnung
mittels Taylorentwicklung, vgl. (3.3). Dafür mussten wir hinreichende Glattheit der
exakten Lösung y annehmen, also hinreichende Glattheit der rechten Seite f . Allgemeiner geht maxt∈[0,T ] ky (p+1) (t)k in die Konstante c ein.
Fazit: Ein Verfahren erreicht bei einer konkreten Aufgabenstelle „seine“ maximale
Konsistenzordnung nur dann, wenn die rechte Seite hinreichend glatt ist.
Exemplarisch listen wir nun einige weitere explizite und implizite ESV auf.
Beispiel 3.7 (Einige explizite und implizite Einschrittverfahren)
(a) Das explizite Euler-Verfahren (Forward Euler) hat die von h
unabhängige Verfahrensfunktion
(b) Das implizite Euler-Verfahren (Backward Euler) hat die implizit definierte Verfahrensfunktion
§ 3. Einschrittverfahren
25
(Lösbarkeit im Kontext allgemeinerer impliziter Verfahren, § 3.3.4)
(c) Das verbesserte Euler-Verfahren (Collatz, 1960) hat die Verfahrensfunktion
(d) Das Verfahren von Heun (Heun, 1900) hat die Verfahrensfunktion
Alle obigen Verfahren sind Beispiele für sogenannte Runge-Kutta-Verfahren, die
systematisch in § 3.3 behandelt werden. Die Struktur der mehrfachen (rekursiven)
f -Auswertungen pro Zeitschritt ist typisch für diese Verfahren.
Beispiel 3.8 (Taylorreihen-Methode)
Weitere explizite ESV, die nicht zu den Runge-Kutta-Verfahren gehören, lassen sich
durch die Taylorreihen-Methode konstruieren.6 Aus der Darstellung
6Dies
ist das einzige Verfahren der ganzen Vorlesung, das Ableitungen von f verwendet, um
eine höhere Ordnung zu erzielen.
26
Kapitel 1. Numerische Lösung von Anfangswertproblemen
Ein Nachteil der Taylorreihen-Methode liegt darin, dass man entsprechend hohe
Ableitungen von f berechnen muss. Der Aufwand zur Verwendung dieser steigt
exponentiell mit p und polynomiell mit n.
Beachte: Bereits y 000 enthält den Term fyy , eine bilineare Abbildung von Rn nach
Rn .
Bemerkung 3.9 (Lokaler Diskretisierungsfehler)
Der in Definition 3.6 eingeführte Konsistenzfehler wird auch als lokaler Diskretisierungsfehler bezeichnet, da wir ihn — wie bereits in § 3.1 gesehen — wie folgt
interpretieren können, siehe auch Abbildung 3.3
27
§ 3. Einschrittverfahren
y(t+h)
h d(y(.),t+h,h)
y(t)
yh(t+h)
=yh(t)
t
t+h
Abbildung 3.3. Veranschaulichung lokaler Diskretisierungsfehler.
Das heißt: Wenn man bei t mit dem exakten Startwert yh (t) = y(t) startet, so ist der
Fehler nach einem einzelnen Integrationsschritt der Länge h gerade h d(y(·), t+h, h).
Achtung: In manchen Büchern wird der Konsistenzfehler als h d(y(·), t + h, h) definiert.
Definition 3.10 (Fehler, Konvergenz, Konvergenzordnung bei ESV)
Es sei y(·) die Lösung des (AWP) auf dem Intervall [0, T ]. Weiter sei yh (·) eine mit
einem ESV (3.12) erzeugte Näherungslösung auf einem Gitter T = {0 = t0 < t1 <
. . . < tN = T }.
(a) Die Größe
(b) Das ESV heißt konvergent, falls auf jeder Folge von Gittern {Th } mit
h & 0 für die erzeugten Näherungslösungen gilt:
(c) Das ESV hat die Konvergenzordnung p ∈ N, falls eine von h unabhängige
Konstante c > 0 sowie h0 > 0 existieren, sodass für die Näherungslösung
auf einem beliebigen Gitter der maximalen Gitterweite h ≤ h0 gilt:
28
Kapitel 1. Numerische Lösung von Anfangswertproblemen
Beim expliziten Euler-Verfahren hatten wir mit Hilfe einer Stabilitätsabschätzung
(diskretes Lemma von Gronwall) nachgewiesen, dass aus der Konsistenz des Schemas
die Konvergenz mit derselben Ordnung folgt. Dies wiederholen wir jetzt für ein
allgemeines ESV.
Lemma 3.11 (Diskrete Stabilität von ESV)
Die Verfahrensfunktion des ESV (3.12) erfülle die Lipschitz-Bedingung (Stabilitätsbedingung)
kΦ(t, w, h) − Φ(t, z, h)k ≤ K kw − zk
(3.16)
für alle t ∈ [0, T ], alle w, z ∈ Rn sowie alle h ≤ h0 . Es sei T = {0 = t0 , t1 , . . . , tN = T }
ein Gitter auf [0, T ] mit maximaler Schrittweite h ≤ h0 und wh und zh zwei beliebige
Gitterfunktionen auf T . Dann gilt die diskrete Stabilitätsabschätzung
Beweis:
§ 3. Einschrittverfahren
29
Bemerkung 3.12 (zu Lemma 3.11)
(a) Beim expliziten Euler-Verfahren galt Φ(t, y, h) = f (t, y) und daher K = L,
die Lipschitz-Konstante der rechten Seite f , und zwar für alle h > 0. Auch
bei allen anderen in dieser Vorlesung behandelten ESV (siehe § 3.3) ist die
Bedingung (3.16) erfüllt. Bei impliziten ESV verlangt sie tatsächlich eine
Beschränkung der Schrittweite h ≤ h0 , was natürlich für die Konvergenzbetrachtung keine Einschränkung darstellt. Im Allgemeinen hängt K von
der Lipschitz-Konstanten der rechten Seite f ab.
(b) Der Beweis von Lemma 3.11 zeigt, dass die (globale) Lipschitz-Bedingung
(3.16) nicht notwendig für alle w, z ∈ Rn benötigt wird, sondern nur für
solche, die in einer „hinreichend großen Menge“ liegen, die für h ≤ h0 alle
Näherungslösungen enthält. Dasselbe trifft auch auf Satz 3.13 zu.
Satz 3.13 (Konvergenz von ESV)
Die Verfahrensfunktion des ESV (3.12) erfülle die Lipschitz-Bedingung (3.16) für alle
h ≤ h0 . Besitzt das Verfahren die Konsistenzordnung p ∈ N und gehört die Lösung
des (AWP) zu C p+1 ([0, T ]; Rn ), dann gilt für die Näherungslösungen auf einem
beliebigen Gitter {0 = t0 , t1 , . . . , tN = T } der Klasse Th mit maximaler Schrittweite
h ≤ h0 die Abschätzung
(3.18)
(3.19)
30
Kapitel 1. Numerische Lösung von Anfangswertproblemen
Das Verfahren ist also (bei Verwendung exakter Anfangswerte) konvergent, und die
Konvergenzordnung entspricht der Konsistenzordnung p.
Beweis: Wir gehen vor wie im Beweis von Satz 3.3. Es sei y(·) die exakte Lösung
unseres (AWP). Wir verwenden Lemma 3.11 (diskrete Stabilität), um die Fehlerfortpflanzung abzuschätzen. Dazu setzen wir für zh die Näherungslösung des ESV
ein und für wh die exakte Lösung y (eingeschränkt auf das Gitter). Dann gilt nach
Voraussetzung für die jeweiligen lokalen Fehler
von der (3.19) eine direkte Konsequenz ist.
Folgerung 3.14
Bemerkung 3.15 (zur Fehlerabschätzung (3.18))
(a) Der globale Fehler kann exponentiell mit der Zeit t anwachsen. Dieses Verhalten entspricht dem der Differenz zweier Lösungen zu verschiedenen Anfangsdaten, siehe Satz 2.1.
(b)
§ 3. Einschrittverfahren
31
(c) (3.18) ist eine a-priori Abschätzung, da die Berechneten yi nicht in die
Schranke eingehen, man also diese Schranke (aufgrund erwarteter Eigenschaften der Lösung, hier Glattheit) bereits erhält bevor man die yi überhaupt berechnet.
Im Gegensatz dazu werden wir später a-posteriori Abschätzungen kennen lernen, welche die berechnete Näherungslösung verwenden um schärfere
Aussagen zu treffen.
(d) Da man die Konstante c des lokalen Fehlers (die maxt∈[0,T ] ky (p+1) (t)k enthält) und die Lipschitz-Konstante K des Verfahrens (beide abhängig von
der rechten Seite f ) im Allgemeinen nicht kennt, kann man (3.18) kaum
dazu verwenden, den globalen Fehler zu schätzen bzw. um eine Schrittweite
h zu bestimmen, um diesen unter einer vorgegebenen Schranke zu halten.
Hierfür sind a-posteriori Schranken besser geeignet.
Die Fehlerabschätzung (3.18) lässt sich leicht dahingehend erweitern, dass Rundungsfehler bei der numerischen Durchführung berücksichtigt werden. In dem Fall
ist also im Beweis von Satz 3.13 nicht mehr dzj = 0, sondern etwa hj−1 dzj = O(ε)
mit der Maschinengenauigkeit ε. Für äquidistante Gitter ergibt sich so die Fehlerabschätzung
In der Praxis zeigt sich allerdings, dass diese Abschätzung i. d. R. zu pessimistisch
ist, da sich Rundungsfehler und lokale Fehler gegenseitig aufheben können. Dennoch
gibt diese Überlegung den richtigen Hinweis, dass es für eine hohe Genauigkeit der
Lösung i. A. wesentlich zielführender ist, ein Verfahren höherer Ordnung mit
moderaten Schrittweiten h zu verwenden als ein Verfahren niedriger Ordnung
mit geringer Schrittweite h.
Kapitel 1. Numerische Lösung von Anfangswertproblemen
10
rounding error
truncation error
total error
9
8
7
6
errors
32
5
4
3
2
1
0
0
0.2
0.4
0.6
0.8
1
h
Abbildung 3.4. Einfluss von Konsistenz- und Rundungsfehlern in
der Fehlerabschätzung.

Download Report