79
Vorlesung Modellierung und Simulation I Prof. Dr. Gillian Queisser 1. Oktober 2013
Vorlesung Modellierung und Simulation I
Prof. Dr. Gillian Queisser
1. Oktober 2013
Inhaltsverzeichnis
1 Gewöhnliche Differentialgleichungen 31.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Anfangswertproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Einschrittverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.1 Explizites Eulerverfahren . . . . . . . . . . . . . . . . . . . . . . . 31.3.2 Implizites Eulerverfahren . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Taylorreihenverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.5 Diskretisierungsfehler und Fehlerordnung . . . . . . . . . . . . . . . . . . . 5
1.5.1 Lokaler und globaler Diskretisierungsfehler . . . . . . . . . . . . . . 51.5.2 Abschätzen des lokalen Diskretisierungsfehlers . . . . . . . . . . . . 61.5.3 Fehlerordnung und Konsistenz . . . . . . . . . . . . . . . . . . . . 8
1.6 Verbesserte Methoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.6.1 Verbesserte Polygonzugmethode . . . . . . . . . . . . . . . . . . . . 81.6.2 Trapezmethode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.6.3 Verfahren von Heun . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.7 Stabilität . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2 Partielle Differentialgleichungen 132.1 Gängige Operatoren in mehrdimensionaler Analysis . . . . . . . . . . . . . 132.2 Beispiele partieller Differentialgleichungen . . . . . . . . . . . . . . . . . . 14
2.2.1 Die Diffusionsgleichung . . . . . . . . . . . . . . . . . . . . . . . . 152.2.2 Die Wellengleichung . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2.3 Poisson- und Potential-Gleichung . . . . . . . . . . . . . . . . . . . 162.2.4 Poisson-Nernst-Planck Gleichungen . . . . . . . . . . . . . . . . . . 17
3 Diskretisierung I: Differenzenverfahren für partielle Differentialgleichungen 183.1 Gebietsdiskretisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193.2 Approximationseigenschaften im R1 . . . . . . . . . . . . . . . . . . . . . . 193.3 Erstellen eines linearen Gleichungssystems . . . . . . . . . . . . . . . . . . 213.4 Finite Differenzen in R2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4.1 Matrix-Vektor-Schreibweise in R2 . . . . . . . . . . . . . . . . . . . 233.4.2 Sternoperatoren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.4.3 Eigenschaften von Differenzensternen . . . . . . . . . . . . . . . . . 27
3.5 M-Matrizen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.5.1 Wiederholung von speziellen Matrixeigenschaften . . . . . . . . . . 283.5.2 Eigenschaften von M-Matrizen . . . . . . . . . . . . . . . . . . . . 28
2
3.5.3 Abschätzen der Eigenwertbereiche einer Matrix . . . . . . . . . . . 293.5.4 Zusammenhang zwischen M-Matrix und Spektralradius . . . . . . 32
3.6 Eigenschaften der Systemmatrix der Poisson-Gleichung . . . . . . . . . . . 353.6.1 Gebräuchliche Matrixnormen und positiv definite Matrizen . . . . 363.6.2 Matrixeigenschaften von Kh . . . . . . . . . . . . . . . . . . . . . . 39
3.7 Konvergenzuntersuchung für das Finite Differenzen Verfahren . . . . . . . 403.7.1 Stetige Abhängigkeit von den Randdaten . . . . . . . . . . . . . . 413.7.2 Konvergenz, Konsistenz und Stabilität . . . . . . . . . . . . . . . . 42
3.8 Das Neumann-Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 443.8.1 Diskretisierung des Neumann-Problems . . . . . . . . . . . . . . . 453.8.2 Lösen des Neumann-Problems . . . . . . . . . . . . . . . . . . . . . 46
3.9 Differenzenverfahren für allgemeine Probleme zweiter Ordnung . . . . . . 48
4 Diskretisierung II: Finite Elemente Verfahren 494.1 Funktionalanalytische Grundlagen . . . . . . . . . . . . . . . . . . . . . . 49
4.1.1 Normierte Räume . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.1.2 Banach-Räume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.1.3 Der Sobolev-Raum L2(Ω) . . . . . . . . . . . . . . . . . . . . . . . 524.1.4 Schwache Differenzierbarkeit . . . . . . . . . . . . . . . . . . . . . . 524.1.5 Die Hilbert-Räume Hk(Ω) und Hk
0 (Ω) . . . . . . . . . . . . . . . . 534.1.6 Dualräume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.2 Variationsformulierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564.2.1 Untersuchung des elliptischen Differentialoperators zweiter Ordnung 574.2.2 Existenz und Eindeutigkeit für das Variationsproblem . . . . . . . 594.2.3 Schwache Lösung des Randwertproblems . . . . . . . . . . . . . . . 624.2.4 Variationsproblem der Neumann-Randwertaufgabe . . . . . . . . . 63
4.3 Galerkin-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634.4 Finite Elemente Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.4.1 Beispiel von Courant . . . . . . . . . . . . . . . . . . . . . . . . . . 664.4.2 Triangulierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.4.3 Finite Elemente im R1 . . . . . . . . . . . . . . . . . . . . . . . . . 704.4.4 Finite Elemente im R2 . . . . . . . . . . . . . . . . . . . . . . . . . 734.4.5 Konvergenzaussagen zu Finite Elemente Verfahren . . . . . . . . . 76
5 Ausblick 775.1 Finite Volumen Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . 775.2 Lösen von linearen Gleichungssystemen . . . . . . . . . . . . . . . . . . . . 78
3
1 Gewöhnliche Differentialgleichungen
Dieses Kapitel beschäftigt sich mit den Grundlagen zur Lösung von gewöhnlichen Diffe-rentialgleichungen. Die behandelten numerischen Verfahren lassen sich später im Rahmender zeitlichen Diskretisierung bei partiellen Differentialgleichungen einsetzen.
1.1 Motivation
Bespiel überlegen
1.2 Anfangswertproblem
Definition 1. Anfangswertproblem (AWP) Sei A : D ⊂ R× RN → RN dann heißt
∂u
∂t= A(t, u(t)) t ∈ [T0, T ] (1.1)
u(T0) = u0 (1.2)
1.3 Einschrittverfahren
Definition 2. Einschrittverfahren Ein Einschrittverfahren berechnet iterativ eine dis-krete Approximation u∆t(t+ ∆t) für u(t+ ∆t) bei bekanntem Startzeitpunkt T0 und be-kannten Anfangswerten u0 nach folgender Vorschrift:
u∆t(T0) = u0 (1.3)u∆t(t+ ∆t) = u∆t(t) + ∆t · φ(t, u∆t(t+ ∆t), u∆t(t)) (1.4)
hierbei heißt φ(t, u∆t(t+ ∆t), u∆t(t)) Verfahrensfunktion.
1.3.1 Explizites Eulerverfahren
Wir betrachten eine gewöhnliche Differentialgleichung (ODE: ordinary differential equa-tion) der Art
y′(x) = f(x, y(x)) (1.5)
Gesucht ist die Lösung y(x) mit vorgegebenem Startwert y(x0) = y0, d.h. die Steigungam Startpunkt ist gegeben.
4
Idee 1. Lokale lineare Approximation der gesuchten Funktion. Gegeben sei xk, yk undsuche yk+1 für xk+1:
yk+1 = yk + hf(xk, yk) (1.6)
Dieses Verfahren heißt Eulerverfahren und ist ein explizites Einschrittverfahren.
Definition 1. Ein Einschrittverfahren heißt explizit, wenn f nur von (xk, yk) abhängtund nicht von yk+1. Ein Verfahren heißt implizit, wenn f von (xk, yk, yk+1) abhängt.
1.3.2 Implizites Eulerverfahren
1.4 Taylorreihenverfahren
Idee 2. Verwende mehr Summanden aus der Taylorreihe.
y(x) = y(x0)+(x− x0)
1!y′(x0)+
(x− x0)2
2!y′′(x0)+ . . .+
(x− x0)p
p!y(p)(x0)+Rp+1. (1.7)
Mit xk+1 − xk =: h folgt
yk+1 = yk +h
1!y′k +
h2
2!y′′k +
h3
3!y
(3)k + . . .+
hp
p!y
(p)k +Rp+1. (1.8)
Der spätere Abbruch der Taylorreihe liefert eine bessere lokale Approximation der gesuch-ten Funktion, die Schwierigkeit besteht jedoch in der Berechnung zusätzlicher höhererAbleitung bis zur p-ten Ordnung.
Beispiel 1. Betrachte die Funktion
y′ = −2xy2 (1.9)
mit y(0) = 1. Die Taylorentwicklung von y um x liefert
y(x) = yk + c1(x− xk) + c2(x− xk)2 + c3(x− xk)3 + c4(x− xk)4 + . . . (1.10)
mit unbekannten Koeffizienten ci. Setze in y′ = −2xy2 ein, mit x = xk + h:
c1 + 2c2h+ 3c3h2 + 4c4h
3 + . . . = −2(xk + h)(yk + c1h+ c2h
2 + c3h3 + c4h
4 + . . .)2
=
= −2(xk + h)(y2k + 2c1ykh+ (c2
1 + 2c2yk)h2+
+ (2c1c2 + 2c3yk)h3 + . . .
)= −2xky
2k + (−2y2
k − 4c1xkyk)h+(−4c1yk − 2xk(c
21 + 2c2yk)
)h2 +
+(−2(c2
1 + 2c2yk)− 4xk(c1c2 + c3yk))h3 + . . .
5
Ein Koeffizientenvergleich liefert
c1 = −2xky2k (1.11)
c2 = −(yk + 2c1xk)yk (1.12)
c3 =−(4c1yk + 2xk(c
21 + 2c2yk))
3(1.13)
c4 = −(
1
2c2
1 + c2yk + xk(c1c2 + c3yk)
)(1.14)
Bei Betrachtung des Approximationsfehlers ek := y(xk)− yk weist dieses Verfahren einebessere Näherung als das explizite Eulerverfahren auf.
1.5 Diskretisierungsfehler und Fehlerordnung
Frage 1. Wie gut sind unsere Approximationen und wie schnell wird die Approximationbesser für h→ 0?
Betrachteyk+1 = yk + hΦ(xk, yk, yk+1, h) (1.15)
Die Funktion Φ wird Verfahrensfunktion genannt. Für unsere bisherigen Verfahren ist Φgegeben durch
1. Explizites Eulerverfahren: Φ(xk, yk, h) = f(xk, yk)
2. Implizites Eulerverfahren: Φ(xk, yk, yk+1, h) = f(xk, yk)
3. Tayloralgorithmus: Φ(xk, yk, yk+1, h) = c1 + c2h+ c3h2 + . . .+ cph
(p−1)
1.5.1 Lokaler und globaler Diskretisierungsfehler
Seinen im Folgenden xk = x0 + k · h äquidistante Stützpunkte, y(xk) die exakte Lösungin xk und yk die approximierte Lösung in xk.
Definition 2. Unter dem lokalen Diskretisierungsfehler in xk+1 versteht man den Wert
dk+1 := y(xk+1)− y(xk)− h · Φ(xk, y(xk), y(xk+1), h) (1.16)
Die Bezeichnung lokal wird deshalb gewählt, weil dk+1 der Fehler aus einem enzelnenSchritt des Einschrittverfahrens ist und gibt Auskunft über die Qualität der Verfahrens-funktion. In der Praxis ist die folgende Frage wichtig.
Frage 2. Wie gut ist die berechnete Approximation nach mehreren Schritten?
Definition 3. Als globalen Fehler gk in xk bezeichnet man die Differenz
gk := y(xk)− yk (1.17)
6
Idee 3. Bringe dk und gk in Beziehung und schätze gk mit Hilfe von dk ab.
Voraussetzung: Φ(x, y, z, h) erfülle die lokale Lipschitz-Bedingung
|Φ(x, y, z, h)− Φ(x, y∗, z, h)| ≤ L |y − y∗| (1.18)|Φ(x, y, z, h)− Φ(x, y, z∗, h)| ≤ L |z − z∗| (1.19)
(1.20)
für x, y, y∗, z, z∗, h ∈ B und 0 > L < ∞. Es folgt aus dk+1 := y(xk+1) − y(xk) − h ·Φ(xk, y(xk), y(xk+1), h):
y(xk+1) = y(xk) + hΦ(xk, y(xk), y(xk+1), h) + dk+1 (1.21)
Subtrahiereyk+1 = yk + hΦ(xk, yk, yk+1, h) (1.22)
⇒ gk+1 = gk + h (Φ(xk, y(xk), y(xk+1), h)− Φ(xk, yk, y(xk+1), h) (1.23)+Φ(xk, yk, y(xk+1), h)− Φ(xk, yk, yk+1, h) + dk+1) (1.24)
Wegen der Lipschitzbedingung und hL < 1 folgt
|gk+1| ≤ |gk|+ h (L · |y(xk)− yk|+ L · |y(xk+1 − yk|) + dk+1 (1.25)
⇒ |gk+1| ≤1 + hL
1− hL|gk|+
|dk+1|1− hL
(1.26)
Für den expliziten Fall gilt:
|gk+1| ≤ (1 + hL)|gk|+ |dk+1| (1.27)
Es existiert ein K > 0, sodass 1+hL1−hL = 1 + hK. Dann können wir schreiben
|gk+1| ≤ (1 + a)|gk|+ b (1.28)
mit geeigneten Konstanten a, b.
1.5.2 Abschätzen des lokalen Diskretisierungsfehlers
Sei D ≥ 0 so definiert, dassmaxk|dk| ≤ D (1.29)
Hilfssatz 1. Erfülle gk die Bedingung
|gk+1| ≤ (1 + a)|gk|+ b. (1.30)
Dann gilt:
|gn| ≤(1 + a)n − 1
b+ (1 + a)n|g0| ≤
b
a(ena − 1) + ena|g0| (1.31)
7
Beweis.
|gn| ≤ (1 + a)|gn−1|+ b ≤ (1 + a)2|gn−2|+ ((1 + a) + 1)b
...≤ (1 + a)n|g0|+
((1 + a)n−1 + . . .+ (1 + a) + 1
)b
=(1 + a)n − 1
a· b+ (1 + a)n|g0|
Ferner gilt: et konvex ⇒ (1 + t) ≤ et, ∀t.
⇒ (1 + a)n ≤ ena
Mit g0 = y(x0)− y0 = 0 ergibt sich aus dem Hilfssatz folgender
Satz 1. Für den globalen Fehler gn an der Stelle xn = x0 + nh gilt für eine expliziteMethode
|gn| ≤D
hL
(enhL − 1
)≤ D
hL· enhL. (1.32)
Für den impliziten Fall gilt
|gn| ≤D
hK(1− hL)
(enhK − 1
)≤ D
hK(1− hL)enhK (1.33)
Der gobale Fehler hängt also stark von D ab, sowie von K und L.
Beispiel 2. Für das explizite Eulerverfahren gilt
dk+1 = y(xk+1)− y(xk)− hf(xk, y(xk))
Entwickle y(xk+1) mit der Taylorexpansion
y(xk+1) = y(xk) + hy′(xk) +1
2h2y′′(xk + Θh)
mit 0 < Θ < 1. Wegen f(xk, y(xk)) = y′(xk) folgt
dk+1 = y(xk) + hy′(xk) +1
2h2y′′(xk + Θh)− y(xk)− hy′(xk) =
1
2h2y′′(xk + Θh)
SeiM2 := max
x0≤ξ≤xn
∣∣y′′(ξ)∣∣⇒ max |dk+1| ≤
1
2h2M2
Eingesetzt:
|gn| ≤D
hLenhL =
hM2
2LenhL
D.h. gn nimmt proportional zu h ab. Das Eulerverfahren besitzt die Fehlerordnung 1.
8
1.5.3 Fehlerordnung und Konsistenz
Definition 4. Ein Einschrittverfahren besitzt die Fehlerordnung p, falls für seinen lokalenDiskretisierungsfehler dk gilt:
max |dk| ≤ D = C · hp+1 = O(hp+1) (1.34)
Dabei ist C eine Konstante.
Folgerung. Der globale Fehler gn einer expliziten und impliziten Methode ist beschränktdurch
|gn| ≤C
LenhLhp = O(hp) (1.35)
Beispiel 3. Für die Taylorreihenmethode bis Ordnung p gilt:
dk+1 = y(xk+1)− y(xk)− hy′(xk)−h2
2!y′′(xk)− . . .−
hp
p!y(p)(xk)
=hp+1
(p+ 1)!y(p+1)(xk + Θh)
(1.36)
für 0 < Θ < 1 Daraus ergibt sich für die Taylorreihenmethode die Fehlerordnung p.
Definition 5. Ein Einschrittverfahren heißt mit der gewöhnlichen Differentialgleichungkonsistent, falls ihre Fehlerordnung mindestens 1 ist.
1.6 Verbesserte Methoden
Im Folgenden werden verbesserte Diskretisierungsverfahren für ODEs hergeleitet.
1.6.1 Verbesserte Polygonzugmethode
Idee 4. Kombiniere zwei Schrittweiten.
y(1)k+1 = yk + hf(xk, yk) (1.37)
y(2)
k+ 12
= yk +h
2f(xk, yk) (1.38)
y(2)k+1 = y
(2)
k+ 12
+h
2f(xk +
h
2, y
(2)
k+ 12
) (1.39)
Durch Anwendung der Richardson-Extrapolation folgt:
yk+1 = 2y(2)k+1 − y
(1)k+1 = 2y
(2)
k+ 12
+ hf
(xk +
h
2, y
(2)
k+ 12
)− yk − hf(xk, yk)
= yk + hf
(xk +
h
2, yk +
h
2f(xk, yk)
)Daraus folgt die
9
Verbesserte Polygonzugmethode
k1 := f(xk, yk) (1.40)
k2 := f
(xk +
h
2, yk +
1
2k1
)(1.41)
yk+1 = yk + hk2 (1.42)
1.6.2 Trapezmethode
Durch eine äquivalente Integraldarstellung des ursprünglichen ODE-Problems gelangenwir zu der Trapezmethode.
Äquivalente Integraldarstellung
Im Folgenden sei die Gleichung
y′(x) = f(x, y(x))
zugrunde gelegt. Auf dem Intervall [xk, xk+1] wird auf beiden Seiten integriert:∫ xk+1
xk
y′(x)dx =
∫ xk+1
xk
f(x, y(x))dx (1.43)
⇔ y(xk+1)− y(xk) =
∫ xk+1
xk
f(x, y(x))dx (1.44)
Da y(x) unbekannt ist, wird das Integral auf der rechten Seite approximiert durch∫ xk+1
xk
f(x, y(x))dx ≈ f(xk, yk)(xk+1 − xk) +
+1
2(xk+1 − xk) (f(xk+1, yk+1)− f(xk, yk))
=h
2(f(xk, yk) + f(xk+1, yk+1))
Definition 6. Die so konstruierte Trapezmethode ist durch die Verfahrensfunktion
yk+1 = yk +h
2(f(xk, yk) + f(xk+1, yk+1)) (1.45)
definiert.
Da die Trapezmethode ein implizites Verfahren ist, muss eine Näherung für yk+1 gefun-den werden. Dies kann durch eine Fixpunktiteration geschehen.
10
Fixpunktiteration.
y(0)k+1 = yk + hf(xk, yk) (1.46)
y(n+1)k+1 = yk +
h
2
(f(xk, yk) + f(xk+1, y
(n)k+1
)(1.47)
Diese Iteration konvergiert gegen einen Fixpunkt yk+1, falls die Lipschitz-Bedingungerfüllt ist und hL
2 < 1, bzw. h < 2L .
Fehlerordnung der Trapezmethode
Die Trapezmethode besitzt die Verfahrensfunktion
Φ(xk, yk, yk+1, h) :=1
2(f(xk, yk) + f(xk+1, yk+1)) .
Der lokale Diskretisierungsfehler lässt demnach berechnen durch
dk+1 = y(xk+1)− y(xk)−h
2(f(xk, y(xk)) + f(xk+1, y(xk+1)))
= y(xk+1)− y(xk)−h
2
(y′(xk) + y′(xk+1)
)= hy′(xk) +
h2
2y′′(xk) +
h3
6y′′′(xk) +O(h4)
−h2
(y′(xk) + y′(xk) + hy′′(xk) +
h2
2y′′′(xk) +O(h3)
)= − 1
12h3y′′′(xk) +O(h4)
Daraus folgt, dass der Hauptteil des lokalen Diskretisierungsfehlers proportional zu h3 istund damit hat die Trapezmethode Fehlerordnung 2. Dies ist identisch mit der verbesser-ten Polygonzugmethode, wir werden aber später sehen, dass die Trapezmethode bessereStabilitätseigenschaften hat.
1.6.3 Verfahren von Heun
In obigem Ansatz bleibt zu entscheiden wie viele Iterationsschritte in der Fixpunktiterati-on ausgeführt werden. In der Praxis wird oft nur ein Schritt ausgeführt, aufgrund dessen,dass durch das numerische Verfahren ohnehin Näherungen yk+1 ≈ y(xk+1) berechnetwerden. Das Verfahren von Heun besteht aus einem Eulerschritt als Prädiktor und einemKorrekturschritt und wird deshalb auch als Prädiktor-Korrektor-Methode bezeichnet:
y(p)k+1 = yk + hf(xk, yk) (1.48)
yk+1 = yk +h
2(f(xk, yk) + f(xk+1, y
(p)k+1) (1.49)
11
Das algorithmische Vorgehen ist dabei
k1 = f(xk, yk) (1.50)k2 = f(xk + h, yk + hk1) (1.51)
yk+1 = yk +1
2h(k1 + k2) (1.52)
Dies entspricht einer Mittelung der Steigungen k1 und k2 in den Punkten (xk, yk) und(xk+1, y
(p)k+1. Die Fehlerordnung den Heun-Verfahrens ist 2, der Beweis ist analog zur
Polygonzugmethode zu führen.
1.7 Stabilität
Neben der Feststellung der Konsistenz eines Einschrittverfahrens, stellen wir uns dieFrage nach der Konvergenz eines Verfahrens, speziell die
Frage 3. Welche Anforderung muss an die Zeitschrittweite gestellt werden, um Konver-genz des Verfahrens zu garantieren?
Wir betrachten folgendes
Beispiel 4.u′(t) = λ · u(t), λ < 0
Anwendung des expliziten Eulerverfahrens. Die Verfahrensfunktion lautet
Φ(t, uh(t), uh(t+ h), h) = f(t, uh(t), h) = λ · uh(t)
=⇒ u(t+ h) = u(t) + λ · h · u(t) = (1 + λh)u(t)
=⇒ u((k + 1) · h) = (1 + λh)k · u0
Dies bedeutet, dass das Verfahren nur dann konvergiert, wenn
|1 + λh| ≤ 1
⇒ 1 ≥ 1 + λh > −1
⇒ h < − 2
λ
Diese Anforderung an die Zeitschrittweite heißt Courant-Friedrichs-Levy-Bedingung undh · λ die CFL-Zahl.
Anwendung des impliziten Eulerverfahrens. Die Verfahrensfunktion lautet
Φ(t, uh(t), uh(t+ h), h) = f(t, uh(t+ h)) = λuh(t+ h)
⇒ uh(t+ h) = uh(t) + hf(t, uh(t+ h))
= uh(t) + hλuh(t+ h)
⇒ uh(t+ h) =uh(t)
1− hλ=
1
1− hλ)ku0
12
Dies zeigt, dass implizite Eulerverfahren ist stabil für alle λ < 0.
Fazit. Senkenterme in der Differentialgleichung sollten implizit behandelt werden, Quel-lenterme explizit.
Definition 7. Folgende Aussagen definieren die Steifigkeit eines Problems:
1. Ein Problem heißt steif, wenn explizite Einschrittverfahren zu kleine Schrittweitenbenötigen.
2. Ein Problem heißt steif, wenn für die Eigenwerte λ der Jacobi-Matrix Dλf dergewöhnlichen Differentialgleichung gilt:
minλi<0
λi maxλi<0
λi
Fazit. Explizite Verfahren sind nicht geeignet für steife Probleme.
13
2 Partielle Differentialgleichungen
In diesem Kapitel werden Herleitungen für partielle Differentialgleichungen (PDEs) un-terschiedlicher physikalischer Phänomene gezeigt. Die so motivierten PDEs werden inden späteren Kapiteln numerisch behandelt.
2.1 Gängige Operatoren in mehrdimensionaler Analysis
Im Folgenden werden einige grundlegenden Begriffe und Operatoren definiert, die späterin den klassischen Anwendungsfällen, d.h. für klassische partielle Differentialgleichun-gen, benötigt werden. Dieser Abschnitt liefert also nur einen marginalen Einblick in diemehrdimensionale Analysis, die ohnehin eine Grundvoraussetzung für die Analyse höher-dimensionaler Probleme ist.
Definition 8. Sei U ⊂ Rn eine offene Menge und f : U → R eine reelle Funktion. fheißt im Punkt x ∈ U partiell differenzierbar bzgl. der i-ten Koordinatenrichtung, falls
Dif(x) := limh→0
f(x+ hei)− f(x)
h(2.1)
existiert. Dabei ist ei ∈ Rn der i-te Einheitsvektor
ei = (0, . . . , 0, 1, 0, . . . , 0)
Definition 9. Sei U ⊂ Rn offen und f : U → R partiell differenzierbar. Dann heißt
grad f(x) :=
(∂f
∂x1(x), . . . ,
∂f
∂xn(x)
)(2.2)
der Gradient von f im Punkt x ∈ U .
Satz 2. Seien f, g : U → R zwei partiell differenzierbare Funktionen. Dann gilt
grad (f · g) = g · grad f + f · grad g (2.3)
Bemerkung 1. Statt grad(f) wird häufig auch ∇f geschrieben mit
∇ :=
(∂
∂x1, . . . ,
∂
∂xn
)als vektorwertiger Differentialoperator.
14
Definition 10. Sei U ⊂ Rn eine offene Menge und
V = (v1, . . . , vn) : U → Rn
eine partiell differenzierbares Vektorfeld (d.h. alle vi sind partiell differenzierbar). Dannheißt die Funktion
div v :=n∑i=1
∂vi∂xi
(2.4)
die Divergenz des Vektorfeldes v.
Definition 11. Sei U eine offene Menge im R3. Für ein partiell differenzierbares Vek-torfeld v : U → R3 bezeichnet man
rot v :=
(∂v3
∂x2− ∂v2
∂x3,∂v1
∂x3− ∂v3
∂x1,∂v2
∂x1− ∂v1
∂x2
)(2.5)
als Rotation von v.
Bemerkung 2. Die Rotation von v lässt sich auch als Vektorprodukt
rot v = ∇v
schreiben.
Definition 12. Sei U → Rn offen und f : U → R zweimal stetig partiell differenzierbar.Dann ist der Laplace-Operator definiert als
∆f := div · grad f = div · ∇f =∂2f
∂x21
+ . . .+∂2f
∂x2n
(2.6)
∆ :=∂2
∂x21
+ . . .+∂2
∂x2n
=n∑i=1
∂2
∂x2i
(2.7)
2.2 Beispiele partieller Differentialgleichungen
Im vorigen Kapitel wurden gewöhnliche Differentialgleichungen vom Typ
dudt
= f(t, u)
betrachtet. Die numerische Behandlung dieser Gleichung führte zu unterschiedlichen Zeit-schrittverfahren, d.h. der Operator d
dt wurde diskret behandelt, wobei f(t, u) eine kon-tinuierliche Funktion war. Im jetzigen Kontext sei f jedoch keine explizite Funktion,sondern nur über ihre Ortsableitungen definiert. Dies verlangt nach einer diskreten Be-schreibung im Ort und führt somit zu Gleichungen mit Differentialoperatoren in Zeit undOrt, also partiellen Differentialgleichungen.
15
2.2.1 Die Diffusionsgleichung
Sei c(x, t) eine Funktion in Raum und Zeit, welche die Konzentrationsverteilung einerSpezies, z.B. Kalziumionen in einer Nervenzelle, oder die Wärmeverteilung in einemWärmeleiter beschreibt. Dann ist der Fluss F von c(x, t)
F = −D · grad c, (2.8)
mit Diffusionskoeffizient D. Durch Forderung der Massenerhaltung und Anwendung desGauss’schen Integralsatzes erhält man die Diffusionsgleichung.
Massenerhaltung. Sei V ein beliebiges Volumenelement und c die Konzentration in V .Die Änderung von c in V ist beschrieben durch∫
V
∂c
∂td~x.
Unter der Voraussetzung, dass keine Quellen oder Senken in V enthalten sind gilt:
−∫∂VF · ~ndS =
∫V
∂u
∂td~x. (2.9)
Satz 3. Der Gauss’sche Integralsatz besagt∫Vdiv F (~x)d~x =
∫∂VF (~x) · ~ndS (2.10)
Daraus folgt
−∫∂VF · ~ndS = −
∫Vdiv F (~x)d~x =
∫V
∂c
∂td~x (2.11)
=⇒ −div F =∂c
∂t(2.12)
=⇒ ∂c
∂t= −div (D∇c) (Diffusionsgleichung) (2.13)
2.2.2 Die Wellengleichung
Wir betrachten für die Wellengleichung die Größen Geschwindigkeit v, Dichte ρ undDruck p.
1. Es gilt∂ρ
∂t= −ρ0div v (2.14)
wobei ρ0 eine feste Dichte definiert. Die Herleitung obiger Gleichung geschieht ana-log zur Diffusionsgleichung also über die Annahme der Massenerhaltung, gefolgtvon der Anwendung des Gauss’schen Integralsatzes.
16
2. Newton’sches Gesetz: Es gilt
ρ0∂v
∂t= −grad p (2.15)
und bedeutet, dass eine räumliche Änderung des Druckfeldes eine Beschleunigungbewirkt.
3. Zustandsgleichung: Der Druck p ist bei konstanter Temperatur proportional zurDichte
⇒ p = c2 · ρ (2.16)
⇒ ∂2
∂t2ρ = −ρ0div
(∂v∂t
)= −div
(ρ0∂v
∂t
)(2.17)
⇒ 1
c2
∂2
∂t2p = div (grad p) (2.18)
⇔ ∂2
∂t2p = c2 · div (grad p) = c2∆p (2.19)
Beispiel 5. Wellengleichung in R1 und R2.
1. Im R1 beschreibt die Wellengleichung eine schwingende Saite:
utt = uxx
2. In R2 beschreibt die Wellengleichung eine schwingende Membran:
utt = c2∆u
2.2.3 Poisson- und Potential-Gleichung
Sei Ω ⊂ Rd, d = 2, 3 und f : Ω → R eine bekannte Ladungsdichteverteilung in Ω. Fürdie Spannung u gilt
−∆u = f(x) in Ω (2.20)
Ein Spezialfall der Poissongleichung ist die Potentialgleichung
∆u = 0 in Ω ⊂ R2 (2.21)
und beschreibt einen Ladungsfreien Raum.
17
Beispiel 6. Lösung der Potentialgleichung auf einer Kreisscheibe mit Radius 1. BetrachteΩ := (x, y) ∈ R2;x2 + y2 < 1 und transformiere x und y in Polarkoordinaten um:
x := r · cosφ
y := r · sinφ
Transformiere die Potentialgleichung in Polarkoordinaten um:
⇒ ∆u =∂2u
∂r2+
1
r
∂u
∂r+
1
r2
∂2u
∂φ2(2.22)
Dann erfüllen rk cos kφ und rk sin kφ die Potentialgleichung. Zu wählen ist noch dieRandbedingung für Radius r = 1:
u|Rand = u(cosφ, sinφ) = a0 +
∞∑k=1
(ak cos kφ+ bk sin kφ) (2.23)
Dann lautet die Lösung im Innern (r < 1):
u(x, y) = a0
∑rk · (ak cos kφ+ bk sin kφ) (2.24)
2.2.4 Poisson-Nernst-Planck Gleichungen
Als ein Beispiel für gekoppelte und nichtlineare PDEs können die Poisson-Nernst-PlanckGleichungen erwähnt werden. Diese Gleichungen beschreiben den Prozess der Elektrodif-fusion, also die Diffusion von Teilchen gekoppelt mit einem Konvektionsterm, der durchdas zu berechnende elektrische Feld bestimmt ist:
∂ci∂t
= ∇(Di∇ci +Di
ziF
RTci∇Φ
)(2.25)
∇(ε∇Φ) = −
(ρf +
∑i
ciziF
)(2.26)
Dabei sind die ci verschiedene Ionenspezies, Di die zugehörigen Diffusionstensoren, zi dieTeilchenladung, F die Faraday-Konstante, R die allg. Gaskonstante, T die Temperaturund Φ das Potential mit spezifischer Konstante ε und statischer Ladungsverteilung ρf .
18
3 Diskretisierung I: Differenzenverfahrenfür partielle Differentialgleichungen
In diesem Kapitel werden wir ein Approximationsverfahren für das kontinuierliche PDE-Problem herleiten, welches auf der Approximation der Ortsableitungen fundiert. Dasals Differenzenverfahren bezeichnete Diskretisierungsverfahren wird für ein- und zwei-dimensionale Fälle hergeleitet. Weiter werden wir Eigenschaften der Systemmatrix deshergeleiteten Gleichungssystems analysieren, die uns am Ende des Kapitels zu einemKonvergenzbeweis für das Differenzenverfahren führen. Die sichtbar werdenden Vor- undNachteile dieses Verfahrens leiten über in das folgende Kapitel alternativer und allgemei-nerer Diskretisierungsverfahren.
Betrachte die Poissongleichung
−∆u = f auf Ω (3.1)
!
Abbildung 3.1: Kontinuierliches Rechengebiet Ω
Die Poissongleichung ist zunächst auf einem kontinuierlichen Gebiet Ω definiert, d.h. anunendlich vielen Punkten. Dies ist für ein numerisches Verfahren nicht zugänglich.
Idee 5. Wähle endlich viele Punkte in Ω aus, in denen −∆u = f erfüllt ist. Dies führt zueiner Approximation des kontinuierlichen Gebiets, sowie der kontinuierlichen Gleichung.
19
3.1 Gebietsdiskretisierung
Betrachte beispielsweise das Einheitsquadrat als kontinuierliches Gebiet
Ω = (x, y) : 0 < x < 1, 0 < y < 1. (3.2)
Das Vorgehen zur Approximation des kontinuierlichen Gebiets, d.h. Diskretisierung vonΩ ist das Folgende:
1. Überziehe Ω mit einem gleichmäßigen Gitter. Die Menge der Gitterknoten wird alsΩh bezeichnet. Eine feste Schrittweite h liefert
Ωh =
(x, y) ∈ Ω :x
h,y
h∈ Z
.
2. Erfülle die Differentialgleichung in jedem Punkt aus Ωh.
• ersetze u(x) durch uh(x). Dabei ist u(x) die kontinuierliche und uh(x) diediskrete Lösung.
• Approximiere die Ableitungen:
limh→0
u(x+ h)− u(x)
h≈ u(x+ h)− u(x)
h
!h
Abbildung 3.2: Diskretisiertes Gebiet Ωh. Dieses entsteht durch Überziehen des konti-nuierlichen Gebiets Ω mit einem Gitter. Die Gitterknoten definieren dasendliche diskretisierte Gebiet
3.2 Approximationseigenschaften im R1
Betrachte das eindimensionale Problem
20
u′′(x) = f(x) inΩ = (0, 1)
u(0) = ϕ0
u(1) = ϕ1
Die Approximation der Ableitung kann auf verschiedene Arten geschehen:
1. rechtsseitig: δ+u(x) = u(x+h)−u(x)h
2. linksseitig: δ−u(x) = u(x)−u(x−h)h
3. symmetrisch: δ0u(x) = u(x+h)−u(x−h)2h
Für die zweite Ableitung können links- und rechtsseitige Differenzen δ+ und δ− kombi-niert werden:
u′′(x) = δ+δ−u(x) =u(x+h)−u(x)
h − u(x)−u(x−h)h
h=u(x+ h)− 2u(x) + u(x− h)
h2(3.3)
Lemma 1. Sei [x− h, x+ h] ⊂ Ω. Es gilt
1. δ±u(x) = u′(x) + hR mit |R| ≤ 12‖u
′′‖∞
2. δ0u(x) = u′(x) + h2R mit |R| ≤ 16‖u
′′′‖∞
3. δ+δ−u(x) = u′′(x) + h2R mit |R| ≤ 112‖u
(4)‖∞
Beweis.zu 1.:
u(x± h) = u(x)± hu′(x) +h2
2u′′(x) + . . .
= u(x)± hu′(x) +h2
2u′′(ξ), mit x ≤ ξ ≤ x+ h
⇔ u(x+ h)− u(x)
h= u′(x) +
h
2u′′(ξ)
≤ u′(x) +h
2‖u′′‖∞
⇒ 1.
zu 2.: Es gelten die Taylorentwicklungen um x± h:
u(x+ h) = u(x) + hu′(x) +h2
2u′′(x) +
h3
6u′′′(ξ)
u(x− h) = u(x)− hu′(x) +h2
2u′′(x)− h3
6u′′′(ξ)
21
Subtraktion ergibt:
δ0u(x) =2hu′(x) + h3
6 (u′′′(ξ) + u′′′(ξ))
2h
= u′(x) +h2
12(u′′′(ξ) + u′′′(ξ)
≤ u′(x) +h2
12· 2‖u′′′‖∞ = u′(x) +
h2
6‖u′′′‖∞
⇒ 2.
zu 3.: Betrachte die Entwicklung um x± h bis zur Ordnung 4:
u(x+ h) = u(x) + hu′(x) +h2
2u′′(x) +
h3
6u′′′(x) +
h4
24u(4)(ξ)
u(x− h) = u(x)− hu′(x) +h2
2u′′(x)− h3
6u′′′(x) +
h4
24u(4)(ξ)
Addition der obigen Gleichungen, Subtraktion von 2u(x) und Division durch h2 liefert
δ+δ−u(x) =h2u′′(x) + h4
24
(u(4)(ξ) + u(4)(ξ)
)h2
⇒ δ+δ−u(x) = u′′(x) +h2
24
(u(4)(ξ) + u(4)(ξ)
)≤ u′′(x) +
h2
12‖u(4)‖∞
⇒ 3.
3.3 Erstellen eines linearen Gleichungssystems
Durch die Approximation der Ableitungen auf einem diskreten Gebiet entsteht ein endli-ches lineares Gleichungssystem, dass es schlußendlich zu lösen gilt. Wir betrachten weiterdie Poisson-Gleichung
u′′(x) = ∆u(x) = f(x),
und approximieren ∆ ≈ ∆h = δ+δ−. Wir erhalten also eine diskretisierte Gleichung
δ+δ−u(x) = f(x) +O(h2)
In Matrix-Vektor-Schreibweise lässt sich obige Gleichung zu dem Beispiel darstellen als
22
10 x3x2x1
Abbildung 3.3: Beispieldiskretierung eines eindimensionalen Gebiets mit drei innerenKnoten
δ+δ−u(x) =1
h2
−2 1 01 −2 10 1 −2
· uh(x1)
uh(x2)uh(x3)
=
f(x1)− 1h2u(0)
f(x2)f(x3)− 1
h2u(1)
Im allgemeinen Fall erhalten wir
1
h2·
−2 1 0 . . . 0 0 01 −2 1 0 . . . 0 00 1 −2 1 0 . . . 0
. . . . . .. . . . . .
0 0 . . . 0 1 −2 10 0 . . . 0 0 1 −2
u1
...
un
=
f1 − u0h2
f2
...
fn−1
fn − unh2
Wir erhalten also ein Gleichungssystem von der Form
Kh · uh = fh
Bemerkung 3. Kh ist dünnbesetzt, d.h.
#Ki,j 6= 0; i, j = 1 . . . n = O(n)
3.4 Finite Differenzen in R2
Wir betrachten nun ein zweidimensionales Gebiet Ω und zwar das kontinuierliche, wieauch diskretisierte Einheitsquadrat
Ω := (x, y) : 0 < x < 1, 0 < y < 1Ωh := (x, y) : (x, y) ∈ Ω;x/h, y/h ∈ Z
mit den Gebietsrändern
Γ := (x, y) : x ∈ 0, 1, y ∈ 0, 1Γh := (x, y) ∈ Γ : x/h, y/h ∈ Z
Wir betrachten nun das Randwertproblem
23
−∆u = f in Ω
u = ϕ auf Γ
mit der diskreten Form
−∆huh := (−δ−x δ+x − δ−y δ+
y )uh(x) (3.4)
Definition 13. Mit uh wird die Gitterfunktion von u bezeichnet, also die Reduktion vonu auf das Gitter Ωh.
Wendet man den diskreten Laplace-Operator ∆h auf die Gitterfunktion uh an, so erhältman
−∆huh = (−δ−x δ+x − δ−y δ+
y )uh(x)
= − 1
h2(u(x− h, y) + u(x+ h, y) + u(x, y − h) + u(x, y + h)− 4u(x, y))
Dies zeigt, dass u an 5 Gitterpunkten ausgewertet wird. Deshalb wird obige Darstellungauch als Fünfpunktformel bezeichnet.
3.4.1 Matrix-Vektor-Schreibweise in R2
Im eindimensionalen Fall existierte eine natürliche Nummerierung der Knoten, welchedie Matrixstruktur festlegte (die Ordnung der Knoten wurde stillschweigend verwendet).In R2 sind unterschiedliche Ordnungen der Knoten denkbar.
1 2
4
3
5 6
7 8 9
Abbildung 3.4: Lexikographische Nummerierung der inneren Knoten
24
Lexikographische Knotennummerierung
Die lexikographische Knotennummerierung ist eine zeilenweise Nummerierung der Kno-ten und definiert eine Matrix der Form
1
h2
4 −1 0 −1 0 0 0 0 0−1 4 −1 0 −1 0 0 0 0
0 −1 4 −1 0 −1 0 0 0−1 0 0 4 −1 0 −1 0 0
0 −1 0 −1 4 −1 0 −1 00 0 −1 0 −1 4 0 0 −10 0 0 −1 0 0 4 −1 00 0 0 0 −1 0 −1 4 −10 0 0 0 0 −1 0 −1 4
· uh = fh (3.5)
wobei mit fh die rechte Seite f versehen mit den Einträgen die aus der Randbedingungentstehen bezeichnet wird.Obige Matrix hat eine Blocktridiagonalstruktur der Form
Kh =1
h2
D −I 0 0−I D −I 0
. . .0 0 −I D
(3.6)
mit
D =
4 −1 0 · · · 0−1 4 −1 0 · · ·
. . .. . .
0 · · · 0 4 −1
−I =
−1 0 0 · · · 0
0 −1 0 · · · 0. . .
. . .0 · · · 0 0 −1
Schachbrettnummerierung
Die Schachbrettnummerierung nummeriert die Knoten in der Schwarz/Weiß-Abfolge ei-nes Schachbretts.Daraus entsteht folgende Matrixstruktur
25
1 2
4
3
5
6
7 8
9
Abbildung 3.5: Schachbrettnummerierung der inneren Knoten
1
h2
4 0 0 0 0 −1 −1 0 00 4 0 0 0 −1 0 −1 00 0 4 0 0 −1 −1 −1 −10 0 0 4 0 0 −1 0 −10 0 0 0 4 0 0 −1 −1−1 −1 −1 0 0 4 0 0 0−1 0 −1 −1 0 0 4 0 0
0 −1 −1 0 −1 0 0 4 00 0 −1 −1 −1 0 0 0 4
· uh = fh (3.7)
Dies definiert eine Matrix mit der Struktur
Kh =
(D1 LLT D2
)(3.8)
mit den aus Kh ersichtlichen Submatrizen Di, L und LT .
Bemerkung 4. Die Kontennummerierung ändert nichts an den algebraischen Eigen-schaften des Systems, kann aber technische Aspekte bei der Implementierung beeinflussen.
3.4.2 Sternoperatoren
Die Fünfpunktformel definiert die Rechenvorschrift zur Approximation des Laplace-Operators in R2 über finite Differenzen. Eine vereinfachte Schreibweise hierfür ist dieDefinition eines Sternoperators.
Definition 14. Die Fünfpunktformel wird über einen Fünfpunktstern folgendermaßendefiniert
26
−∆huh =1
h2(−u(x− h, y)− u(x+ h, y)− u(x, y − h)− u(x, y + h) + 4u(x, y))
=:1
h2
−1−1 4 −1
−1
= −∆h
Bemerkung 5. Der Fünfpunktstern ist keine Matrix, sondern lediglich eine Schablonedie auf Ωh aufgelegt die Rechenvorschrift der Fünfpunktformel vorgibt und damit amEnde auch die Matrixstruktur von Kh. Die Nummerierung der Knoten geht nicht in denFünfpunktstern ein, da er direkt auf das Gitter aufgelegt wird.
Weitere Sternoperatoren und Rechenvorschriften
1. in R1
a) δ+ = 1h ·[
0 −1 1]
b) δ− = 1h ·[−1 1 0
]c) δ0 = 1
2h ·[−1 0 −1
]2. Allgemeiner Differenzenstern:
1
hk
...c−1,1(x, y) c0,1(x, y) c1,1(x, y)
· · · c−1,0(x, y) c0,0(x, y) c1,0(x, y) · · ·c−1,−1(x, y) c0,−1(x, y) c1,−1(x, y)
...
=
1
hk·∑i,j
cijuh(x+ ih, y + jh)
3. Multiplikation von Sternen: Am Beispiel von zwei eindimensionalen Sternenwird die Multiplikation von Sternen demonstriert.
[a b c
] [d e f
]uh =
[a b c
]· (d · uh(x− h) + e · uh(x) + f · uh(x+ h))
= a(duh(x− 2h) + euh(x− h) + fuh(x))
+b(duh(x− h) + euh(x) + fuh(x+ h))
+c(duh(x) + euh(x+ h) + fuh(x+ 2h))
=[ad ae+ bd af + be+ cd bf + ce cf
]
27
3.4.3 Eigenschaften von Differenzensternen
Wir wollen einige typische (aber nicht allgemein gültige) Eigenschaften von Differen-zensternen zusammenstellen. Basierend darauf definieren wir danach eine Klasse vonMatrizen, die abschließend wichtige Eigenschaften der Systemmatrix der diskretisiertenModellgleichung liefern wird.Folgende gängige Eigenschaften von Differenzensternen, bzw. der Systemmatrix werdenbetrachtet:
1. Die Zeilensumme∑n
j=1 cij ist Null.
n∑j=1
cij = 0, ∀i = 1 . . . n
2. Vorzeichenmuster: Für ∆h gilt:
cii > 0, cij ≤ 0 (i 6= j)
Man beachte das Vorzeichenmuster gilt nicht immer, z.B. für ∆2h ist das Muster
nicht erfüllt.
3. Diagonaldominanz:
a) Schwache Diagonaldominanz
aii ≥n∑
j = 1j 6= i
|aij|
.
b) Starke Diagonaldominanz
aii >
n∑j = 1j 6= i
|aij|
4. Kh ist symmetrisch
Aus diesen Eigenschaften lässt sich eine Klasse von speziellen Matrizen, den sogenanntenM-Matrizen, aufbauen.
3.5 M-Matrizen
In diesem Abschnitt wird die Definition einer M-Matrix eingeführt. Anschließend zeigenwir, dass die Systemmatrix Kh M-Matrix Eigenschaften besitzt. Diese werden bei derAnalyse von Kh nützlich sein, um abschließend wichtige Matrixeigenschaften zu spezifi-zieren.
28
3.5.1 Wiederholung von speziellen Matrixeigenschaften
Definition 15. Eine Matrix A ∈ M(n × n,K) heißt invertierbar, wenn es eine MatrixA ∈M(n× n, k) gibt, mit
A · A = A ·A = En
wobei En die Einheitsmatrix ist.
Definition 16. Eine Matrix A ∈M(n× n,K) heißt regulär, wennn∑j=1
λiAj = 0⇔ λi = 0, ∀i = 1 . . . n
wobei Aj die Spaltenvektoren von A sind.
Lemma 2. Eine reguläre, bzw. nicht-singuläre Matrix A ∈M(n×n,K) ist invertierbarund lässt sich als endliches Produkt von Elementarmatrizen darstellen.
Lemma 3. Für A ∈M(n× n,K) sind folgende Bedingungen äquivalent:
1. A ist invertierbar
2. AT ist invertierbar
3. Spaltenrang von A = n, Zeilenrang von A = n
Insbesondere gilt: (AT )−1 = (A−1)T .
3.5.2 Eigenschaften von M-Matrizen
Im Folgenden werden Matrizen A,B ∈ M(n × n,K) mit Einträgen (aij)i,j=1...n bzw.(bij)i,j=1...n betrachtet.
Definition 17. A ≥ B wird komponentenweise definiert, d.h. es gilt aij ≥ bij , ∀i, j =1 . . . n. Analog lässt sich A ≤ B, A < B und A > B definieren.
Definition 18. Eine n×n-Matrix heißt M-Matrix, wenn folgende Eigenschaften erfülltsind:
1. Vorzeichenbedingung: aii > 0, aij ≤ 0 ∀i, j = 1 . . . n, i 6= j
2. A ist regulär und A−1 ≥ 0.
Frage 4. Gelten für die Matrix Kh aus dem Modellproblem
−∆u = f
Khuh = fh
(3.9)
die M-Matrix Eigenschaften?
29
Die Vorzeichenbedingung lässt sich direkt an Kh ablesen. Zu zeigen bleibt, dass A regulärist und A−1 ≥ 0.
Definition 19. Eine Matrix A heißt irreduzibel, falls jeder Index i ∈ 1, . . . , n mitjedem Index j ∈ 1, . . . , n verbunden ist.
1. i ist mit j direkt verbunden, wenn aij 6= 0.
2. i ist mit j verbunden, wenn eine Kette aus direkten Verbindungen ik, k = 1 . . . pexistiert, sodass
i = i1, i2, . . . , ip = j
mit aik−1ik 6= 0.
3. Eine Alternative Definition ist
A irreduzibel ⇔ Es existiert keine Permutation Π
mit
ΠTAΠ =
(A1 00 A2
)
3.5.3 Abschätzen der Eigenwertbereiche einer Matrix
Wir werden das Kriterium von Gerschgorin in zwei Fassungen beweisen. Dieses Kriteriumwird uns Auskunft über die Eigenwertbereiche geben.
Kriterium von Gerschgorin
Gegeben sei eine Matrix A = (aij) ∈ Cn×n. Betrachte dazu
Ki := z ∈ C : |z − aii| ≤n∑
j = 1j 6= i
|aij |, i = 1 . . . n
Satz 4. Die n Eigenwerte von A liegen in der Vereinigung
n⋃i=1
Ki.
Beweis.Sei λ Eigenwert von A und v Eigenvektor von A.oBdA. kann ‖v‖∞ := maxni=1 |vi| = 1 = |vr| angenommen werden. Es gilt
(A− λEn)v = 0
30
Für die r-te Zeile folgt
(arr − λ)vr = −n∑
i=1,i 6=rarivi
Anwendung der Dreiecksungleichung ergibt
|arr − λ| = |(arr − λ)vr| =
∣∣∣∣∣∣n∑
i=1,i 6=jarivi
∣∣∣∣∣∣≤
n∑i=1,i 6=j
|ari||vi| ≤n∑
i=1,i 6=j|ari|
⇒ |arr − λ| ≤n∑
i=1,i 6=j|ari|
⇔ λ ∈ Kr = z ∈ C : |z − arr| ≤n∑
i=1,i 6=r|ari|
Damit ist gezeigt, dass jeder Eigenwert in der Vereinigung aller Kreisscheiben Ki, i =1 . . . n liegt. Wendet man das Kriterium von Gerschgorin auf die Systemmatrix Kh an,dann folgt:
1.∑n
i=1,i 6=j aij = 4, ∀j ⇒ Radius für alle Kreise ist 4.
2. ajj = 4
=⇒ λ ∈[0, 8h−2
]Im nächsten Schritt wollen wir zeigen, dass die Eigenwerte von Kh echt größer Null sind,d.h. wir wollen zeigen
λ ∈(0, 8h−2
)Dazu beweisen wir die
Verschärfung des Kriteriums von Gerschgorin
Unter der Voraussetzung, dass A irreduzibel ist, gilt:
λ ∈
(n⋃r=1
Kr
)∪
(n⋂r=1
∂Kr
)mit
31
Kr := z ∈ C : |z − arr| <n∑
i=1,i 6=r|ari|
∂Kr := z ∈ C : |z − arr| =n∑
i=1,i 6=r|ari|
(3.10)
Beweis. λ sei ein Eigenwert von A mit dem zugehörigen Eigenvektor v, mit ‖v‖∞ = 1.Falls λ ∈
⋃nr=1Kr ⇒ Verschärfte Fassung von Gerschgorin bereits erfüllt. Sei also
λ /∈⋃nr=1Kr sondern in der Vereinigung der Ränder.
Hilfsbehauptung. Sei arj 6= 0 (also r direkt verbunden mit j). Dann folgt
|vr| = 1 und |λ− arr| =n∑
i=1,i 6=r|ari|
⇒ |vj | = 1 und |λ− ajj | =n∑
i=1,i 6=j|aji|
Der Satz von Gerschgorin beweist die Existenz eines r ∈ 1, . . . , n mit |vr| = 1 und|λ − arr| ≤
∑ni=1,i 6=r |ari|. Da λ auf dem Rand liegt (siehe Annahme) lässt sich der
Hilfssatz auf λ anwenden.Da A irreduzibel =⇒ Für jedes j ∈ 1, . . . , n existiert eine Verbindung
r = i0, i1, . . . , ik = j mit aip−1ip 6= 0
Die Hilfsbehauptung liefert
|vip | = 1 und |λ− aipip | =n∑
i=1,i 6=ip
|aipipi| ∀p = 0, . . . , k.
Insbesondere gilt |vj | = 1 und |λ− ajj | ∈ ∂Kj . j wurde dabei beliebig gewählt.=⇒ λ muss auf jedem Rand ∂Ki ∀i = 1, . . . , n liegen.=⇒ λ ∈
⋂ni=1 ∂Ki.
=⇒ Behauptung
Beweis der Hilfsbehauptung. Aus |λ−arr| ≤∑n
i=1,i 6=r |ari| folgt wegen der Annahme|λ− arr| =
∑ni=1,i 6=r |ari| und |vr| = 1. Aus der Dreiecksungleichung folgt
n∑i=1,i 6=r
|ari||vi| =n∑
i=1,i 6=r|ari|
32
Dies gilt insbesondere für i = j und da |vj | ≤ 1 muss summandenweise gelten:
|arj ||vj | = |arj | und |arj | 6= 0
⇒ |vj | = 1
Frage 5. Was folgt mit den Kriterien von Gerschgorin für die Matrix Kh?
Es treten die folgenden 3 Fälle auf:
aii =4
hund
ri =n∑
j=1,j 6=i|aij | =
2/h2 für die Eckpunkte3/h2 für die Seitenpunkte4/h2 für die inneren Punkte
Da Kh symmetrisch und reell ist, sind alle Eigenwerte von Kh reell und liegen auf dem re-ellen Zahlenstrahl. Die Kreisränder ∂Kj haben den gleichen Mittelpunkt mit unterschied-lichen Radien und sind deshalb disjunkt aus dem verschärften Kriterium von Gerschgorinfolgt deshalb
λ ∈(
0,8
h2
)3.5.4 Zusammenhang zwischen M-Matrix und Spektralradius
Zunächst definieren wir die Begriffe Diagonaldominanz, Irreduzibilität und den Spektral-radius. Anschließend wird der Zusammenhang zwischen M-Matrix und Spektralradiusanalysiert.
Definition 20. (Diagonaldominanz)
1. A heißt diagonaldominant, falls∑j,j 6=i|aij | < |aii| ∀i = 1, . . . , n (3.11)
2. A heißt irreduzibel diagonaldominant, falls A irreduzibel ist und∑j,j 6=i|aij | < |aii| für ein i (3.12)
∑j,j 6=i|aij | ≤ |aii| ∀i = 1, . . . , n (3.13)
Bemerkung 6. Aus irreduzibel und diagonaldominant folgt irreduzibel diagonaldomi-nant, die Rückrichtung gilt jedoch nicht.
33
Definition 21. (Spektralradius)Der Spektralradius %(A) einer Matrix A ist definiert als
%(A) := max|λ| : λ ist Eigenwert von A. (3.14)
Satz 5. Folgende Aussagen gelten für die Matrix D−1B mit A = D−B, D := diagaii :i = 1, . . . , n und B := D −A:
1. A sei diagonaldominant oder irreduzibel diagonaldominant
⇒ %(D−1B) < 1.
2. A erfülle die Vorzeichenbedingung. Dann gilt
A ist M-Matrix ⇔ %(D−1B) < 1.
Beweis von 1. Betrachte C := D−1B mit
cij = −aijaii, (i 6= j)
cii = 0
1. Sei Diagonaldominanz vorausgesetzt. Dann gilt:
ri :=n∑
j=1,j 6=i|cij | < 1 ∀i = 1, . . . , n.
Aus dem Kriterium von Gerschgorin folgt:
λ ∈n⋃i=1
Kri(cii) =n⋃i=1
Kri(0)
⇒ |λ| ≤ maxi=1...n
ri < 1
⇒ %(C) = %(D−1B) < 1
2. Sei nun die irreduzible Diagonaldominanz vorausgesetzt:
A irreduzibel diagonaldominant ⇒ rj ≤ 1 ∀j = 1, . . . , nrj < 1 für mindestens ein i
Aus der scharfen Version von Gerschgorin folgt
λ ∈n⋃j=1
Krj (0) ∪
n⋂j=1
∂Krj (0)
34
Zu zeigen ist nun⋃nj=1Krj (0) ∪
(⋂nj=1 ∂Krj (0)
)⊂ K1(0).
Fall 1: Alle rj = r sind gleich. Da es mindestens ein i gibt mit ri < 1, folgt
n⋂j=1
∂Krj (0) = ∂Kr(0) ⊂ K1(0)√
Fall 2: Falls nicht alle rj gleich sind, dann ist⋂nj=1 ∂Krj (0) = ∅ (da alle Kreise
den Ursprung 0 haben).⇒ Behauptung.
Beweis von 2. A ist genau dann eine M-Matrix, wenn %(D−1B) < 1 und die Vorzei-chenbedingung gilt. Zu zeigen ist also: A ist nicht-singulär und A−1 ≥ 0⇔ %(D−1B) < 1.Wir zeigen zunächst die
Rückrichtung. Sei %(D−1B) = %(C) < 1. Dann konvergiert die Neumann-Reihe
S :=∞∑ν=0
Cν = (I − C)−1
⇒ S(I − C) = I
⇔ SD−1(D −B) = SD−1A = I
⇒ A−1 = SD−1
Da D−1 ≥ 0, B ≥ 0⇒ C ≥ 0, Cν ≥ 0 und S ≥ 0.⇒ A−1 ≥ 0⇒ 2. Bedingung für M-Matrix.⇒ A ist eine M-Matrix.
Hinrichtung. Sei A eine M-Matrix. λ sei Eigenwert zu dem Eigenvektor u 6= 0 vonD−1B. Es gilt
|λ| · |u| = |λu| = |D−1Bu| ≤ D−1B|u|
Ferner gilt: A−1 ≥ 0 und D ≥ 0⇒ A−1D ≥ 0.
⇒ −A−1DD−1B|u| ≤ −A−1D|λ||u|⇒ |u| = A−1(D −B)|u| = A−1D(I −D−1B)|u|
≤ A−1D|u| −A−1D|λ||u| = (1− |λ|)A−1D|u|
Wäre nun |λ| ≥ 1⇒ |u| − (1− |λ|)A−1D|u| ≤ 0. Da
1− (1− |λ|)A−1D ≥ 0
folgt |u| ≤ 0⇒ u = 0. Dies ist ein Widerspruch zur Annahme.
35
Satz 6. Eine irreduzible M-Matrix A hat eine positive Inverse
A−1 > 0.
Beweis. Seien α, β ∈ 1, . . . , n. Da A irreduzibel ist, existiert eine Verbindung
α = α0, α1, . . . , αk = β
C := D−1B und cαi−1αi > 0
⇒ (Ck)αβ =∑
γ1,...,γk
cαγ1cγ1γ2 . . . cγk−1β ≥ cαα1cα1α2 · . . . · cαk−1β > 0
Wie oben bewiesen gilt %(C) < 1.⇒ S :=
∑∞ν=0C
ν konvergiert. Da Sαβ ≥ (Ck) > 0 ist S durch Ck > 0 beschränkt.⇒ S > 0. Mit A−1 = SD−1 > 0 folgt A−1 > 0
3.6 Eigenschaften der Systemmatrix der Poisson-Gleichung
Nachdem verschiedene Matrixnormen eingeführt werden, wenden wir uns in diesem Ab-schnitt den Eigenschaften der Matrix Kh aus der Poisson-Gleichung zu.
Definition 22. (Vektornorm)V sei ein Vektorraum über K (= R oder C). ‖.‖ heißt Norm in V , falls gilt
‖u‖ = 0⇐⇒ u = 0 (3.15)‖u+ v‖ ≤ ‖u‖+ ‖v‖ ∀u, v ∈ V (3.16)‖λu‖ = |λ|‖u‖ λ ∈ K,u ∈ V (3.17)
Definition 23. (Matrixnorm)V sei ein Vektorraum versehen mit einer Vektornorm ‖.‖. Dann ist
‖A‖M := supu∈V
‖Au‖‖u‖
: 0 6= u ∈ V
(3.18)
die der Vektornorm ‖.‖ zugeordnete Matrixnorm.
Lemma 4. Es gilt‖A‖M ≥ %(A) (3.19)
Beweis. Mit ‖A‖M = supu∈V
‖Au‖‖u‖ : 0 6= u ∈ V
und einem Eigenvektor v von A gilt:
‖Av‖‖v‖
=‖λv‖‖v‖
= |λ|
⇒ sup
‖Au‖‖u‖
≥ max
v∈V
‖Av‖‖v‖
= |λ∗| = %(A)
36
3.6.1 Gebräuchliche Matrixnormen und positiv definite Matrizen
Zeilensummennorm
Die zur Maximumsnorm ‖.‖∞ zugehörige Matrixnorm
‖A‖∞ := maxα∈1,...,n
∑β∈1,...,n
|aαβ|
(3.20)
heißt Zeilensummennorm.
=⇒ ‖Au‖∞ = ‖∑j
aijuji‖∞ (3.21)
für ein i so, dass∑
j |aij | maximal.
Bemerkung 7. Aus A ≥ B folgt ‖A‖∞ ≥ ‖B‖∞, da A ≥ B komponentenweise definiertist.
Für einen späteren Beweis werden wir folgenden Satz benötigen:
Satz 7. Sei A eine M-Matrix. Existiert ein Vektor w mit Aw ≥ 1, dann gilt
‖A−1‖∞ ≤ ‖w‖∞ (3.22)
Beweis. Es gilt|u| ≤ ‖u‖∞ · 1 ≤ ‖u‖∞ ·Aw
für
|u| :=
|u1||u2|...|ui|
Da A eine M-Matrix ist, gilt A−1 ≥ 0.
|A−1u| ≤ A−1|u| ≤ A−1‖u‖∞Aw= ‖u‖∞A−1Aw = ‖u‖∞ · w
⇒ |A−1u|‖u‖∞
≤ w und speziell‖A−1u‖∞‖u‖∞
≤ ‖w‖∞
⇒ ‖A−1‖∞ ≤ ‖w‖∞
Satz 8. Seien A und B M-Matrizen mit B ≥ A. Dann gilt:
0 ≤ B−1 ≤ A−1 und ‖B−1‖∞ ≤ ‖A−1‖∞ (3.23)
37
Beweis. Es giltA−1 −B−1 = A−1(B −A)B−1.
Da A,B M-Matrizen sind, folgt A−1 ≥ 0 und B−1 ≥ 0 und mit B ≥ A folgt B −A ≥ 0.
⇒ A−1 −B−1 = A−1(B −A)B−1 ≥ 0
⇔ A−1 ≥ B−1 und ‖B−1‖∞ ≤ ‖A‖∞
Spektralnorm
Zur euklidischen Vektornorm ‖u‖2 =√c∑n
i=1 |ui|2, c = const > 0 lässt sich eine zuge-hörige Matrixnorm, die Spektralnorm definieren.
Lemma 5. Die zur Vektornorm gehörige Matrixnorm ‖.‖2 heißt Spektralnorm und lässtsich schreiben als
‖A‖2 =√%(AtA) (3.24)
Beweis.‖A‖2 = sup
u∈V
‖Au‖2‖u‖2
: u 6= 0
= max
u∈V,‖u‖2=1‖Au‖2
⇒ ‖A‖22 = maxu∈V,‖u‖2=1
‖Au‖22 = max‖u‖2=1
〈Au,Au〉 = max‖u‖2=1
〈AtAu, u〉
Da AtA symmetrisch ist, liefert eine Hauptachsentransformation
P tAtAP = diag(λi) = D
mit P = (v1, . . . , vn), mit den Eigenvektoren vi von AtA. D enthält die zugehörigenEigenwerte auf der Diagonalen. Für P gilt PP t = 1.
⇒ ‖A‖22 = max‖u‖2=1
〈AtAu, u〉 = max‖u‖2=1
〈AtAPP tu, PP tu〉
=︸︷︷︸u=P tu
max‖u‖2=1
〈AtAPu, P u〉 = max‖u‖2=1
〈P tAtAPu, u〉 = max‖u‖2=1
〈Du, u〉
= max‖u‖2=1
〈n∑i=1
λiui, ui〉 = max‖u‖2=1
(n∑i=1
λi|ui|2)
= λmax = %(AtA)
⇒ ‖A‖2 =√%(AtA)
Bemerkung 8. Falls A symmetrisch ist, kann man obiges Vorgehen direkt auf A an-wenden.
⇒ ‖A‖2 = %(A), für A symmetrisch
38
Positiv definite Matrizen
Bevor wir im folgenden Abschnitt die wichtigsten Eigenschaften der Systemmatrix zurdiskreten Poisson-Gleichung zusammenstellen, benötigen wir noch den Begriff der posi-tiven Definitheit.
Definition 24. Eine Matrix A heißt positiv definit, wenn A symmetrisch ist und
〈Au, u〉 > 0 ∀u 6= 0. (3.25)
Lemma 6. A ist positiv definit, genau dann wenn alle Eigenwerte von A positiv sind.
Beweis. A symmetrisch ⇒ ∃P : P tAP = diag(λi). Dabei sind λi die Eigenwerte von Aund P die Matrix zusammengesetzt aus den Eigenvektoren.
⇒ 〈Au, u〉 = 〈APu, Pu〉 = 〈P tAPu, u〉
= 〈n∑i=1
λiu, u〉 =
n∑i=1
〈λiu, u〉 > 0, gdw λi > 0
Daraus folgt folgendes
Lemma 7. Eine positiv definite Matrix ist regulär und besitzt eine positiv definite Inverse.
Lemma 8. Für A symmetrisch und diagonaldominant (oder irreduzibel diagonaldomi-nant) mit aii > 0, so ist A positiv definit.
Beweis.n∑
j=1,j 6=i|aij | < aii ⇒ Gerschgorinkreise liegen in (0,∞)
⇒ λi positiv.⇒ A ist positiv definit.
Lemma 9. λmin und λmax seien minimaler bzw. maximaler Eigenwert von A (positivdefinit). Dann gilt:
‖A‖2 = λmax
‖A−1‖2 =1
λmin
Beweis. A ist positiv definit. ⇒ A symmetrisch ⇒ ‖A‖2 = %(A)
⇒ ‖A‖2 = λmax
‖A−1‖2 = %(A−1) =1
λmin
39
3.6.2 Matrixeigenschaften von Kh
Die Summe der in den vorigen Abschnitten zusammengefassten Sätze und Matrixeigen-schaften ermöglichen den Beweis folgender Aussagen:
Satz 9. Die Matrix Kh, definiert durch den Stern 11 −4 1
1
besitzt folgende Eigenschaften:
1. Kh ist eine M-Matrix.
2. Kh ist positiv definit.
3. ‖Kh‖∞ ≤ 8h2
und ‖K−1h ‖∞ ≤
18 .
4. ‖Kh‖2 ≤ 8h2
cos2(πh2
)< 8
h2und
5. ‖K−1h ‖2 ≤
18h
2 sin−2(πh2
)= 1
2π2 +O(h2) < 116
Beweis.
1. Kh erfüllt die Vorzeichenbedingung und ist irreduzibel diagonaldominant. In Satz5 haben wir gezeigt, dass Kh dann auch eine M-Matrix ist.
2. Kh ist symmetrisch und irreduzibel diagonaldominant ⇒ Kh ist positiv definit(siehe Lemma 8).
3. a) ‖Kh‖∞ = maxi=1,...,n
(∑nj=1 |Kij |
)= 1
h2max 6, 7, 8 = 8
h2
b) Zu zeigen ist: ‖K−1h ‖∞ ≤
18 . Nutze dazu Satz 7 mit
w(x, y) =x(1− x)
2
und betrachte
Khwh = −(x− h)(1− (x− h))
2h2− (x+ h)(1− (x+ h))
2h2+
2 · x(1− x)
2h2= 1
Khw ≥ 1. Für ‖w‖∞ gilt:
d
dxw(x, y) = −x+
1
2⇒ max
x,yw(x, y) =
1
2
Wegen Satz 7 gilt: ‖K−1h ‖∞ ≤ ‖w‖∞ = 1
8 .
40
4. Folgt aus folgendem
Lemma 10. Die Eigenvektoren und Eigenwerte von Kh sind
uνµ(x, y) = sin(νπx) sin(µπy), (x, y) ∈ Ω (3.26)
mitλνµ =
4
h2
(sin2
(νπh
2
)+ sin2
(µπh
2
)), 1 ≤ ν, µ ≤ n− 1 (3.27)
Beweis.
Khuνµ =
1
h2(− sin(νπ(x− h)) sin(µπy) + 4 sin(νπx) sin(µπy)
− sin(νπ(x+ h)) sin(µπy)− sin(νπx) sin(µπ(y − h))− sin(νπx) sin(µπ(y + h)))
Wende die Regelsin(x± y) = sinx cos y ± sin y cosx
mit x := νπx und y := νπh auf obigen Ausdruck an.
⇒ Khuνµ = λνµu
νµ
⇒ ‖Kh‖2 = %(A) = λmax ≤ 8h2
sin2(π(n−1)h
2
)= 8
h2cos2
(πh2
)< 8
h2
5. Es gilt
‖K−1h ‖2 = %(K−1
h ) =1
λminmit
λmin = λ11 =8
h2sin2
(πh
2
)=
8
h2
(πh
2+O(h2)
)2
=8
h2
((πh
2
)2
+O(h4)
)= 2π2 +O(h2)
⇒ ‖K−1h ‖2 ≤ 1
2π2+O(h2)
3.7 Konvergenzuntersuchung für das Finite DifferenzenVerfahren
Dass ein numerisches Approximationsverfahren mit zunehmender Gitterfeinheit gegen diekontinuierliche Lösung konvergiert, ist Grundvoraussetzung für den effizienten Einsatzdes Verfahrens. Wir werden im Folgenden, dass das behandelte Differenzenverfahren fürdie Poissongleichung stabil und konvergent ist.
41
3.7.1 Stetige Abhängigkeit von den Randdaten
Maximumsprizip
Lemma 11. Sei uh 6= 0 eine Lösung von −∆huh = fh mit fh = 0 und uh|∂Ωh= ϕh. Dann
nimmt uh sein Maximum bzw. Minimum auf dem Rand ∂Ωh an, oder uh ist konstant.
Beweis. Für −∆huh = 0 gilt
uh(x, y) =1
4(uh(x− h, y) + uh(x+ h, y) + uh(x, y − h) + uh(x, y + h))
Angenommen uh ist nicht konstant und uh(x, y) sei Maximum in Ωh.⇒ uh(x± h, y) und uh(x, y ± h) müssen maximal sein. Da Kh aber irreduzibel ist, setztsich diese Bedingung auf ganz Ωh fort.⇒ uh ist konstant. Dies ist ein Widerspruch zu obiger Annahme.
Vergleichsprinzip
Seien uh, vh Lösungen von
−∆huh = fh mit uh|∂Ωh= ϕuh
−∆hvh = fh mit vh|∂Ωh= ϕvh
Dann gilt:
1. ‖uh − vh‖∞ ≤ maxx∈∂Ωh|ϕu(x)− ϕv(x)|
2. uh ≤ vh, falls ϕu ≤ ϕv auf ∂Ωh.
Beweis. Betrachte wh := vh−uh. wh ist Lösung der Poisson-Gleichung −∆hwh = 0 undwh ≥ 0 auf ∂Ωh, falls
ϕu ≤ ϕv.
Aus dem Maximumsprinzip folgt wh ist konstant oder wh > 0. ⇒ Behauptung 2.Mit uh = ϕu und vh = ϕv auf ∂Ωh gilt
−max |ϕu − ϕv| ≤ wh ≤ + max |ϕu − ϕv| auf ∂Ωh
Aufgrund des Maximumsprinzips gilt dies auf ganz Ωh ⇒ Behauptung 1.
42
3.7.2 Konvergenz, Konsistenz und Stabilität
Wir wollen nun das kontinuierliche mit dem numerisch approximierten Problem verglei-chen. Wir betrachten dazu
−∆u = f in Ω
u|Γ = ϕ
und
−∆huh = fh in Ωh
uh|Γh= ϕh
Bemerkung 9. Das kontinuierliche und diskrete Problem existieren in unterschiedlichenGebieten. Um beide Probleme vergleichen zu können, müssen sie in einem Raum definiertsein.
Definition 25. (Restriktion)
Rh : C(Ω) −→ Uh
u 7→ Rhu
mit (Rhu)(~x) = u(~x) für alle ~x ∈ Ωh ist die Restriktion der Lösung in Ω auf Ωh.Dabei sei h ∈ H ⊂ R+ und H eine Menge ohne Häufungspunkt (in unserem Fall istH =
1n : n ∈ N
).
Bemerkung 10. Durch die Restriktion der kontinuierlichen Lösung auf das diskreteGebiet können die Lösungen u und uh nun verglichen werden.
Konvergenz
Die diskrete Lösung uh ∈ Uh konvergiert bzgl. ‖.‖h, h ∈ H gegen u, falls
‖uh −Rhuh‖h −→ 0 (3.28)
Definition 26. uh −Rhu heißt Diskretisierungsfehler des Diskretisierungsverfahrens.
Definition 27. Die Diskretisierung Kh heißt stabil bzgl. ‖.‖∞ für h ∈ H ⊂ R+, falls
suph∈H‖K−1
h ‖∞ ≤ C <∞ (3.29)
Bemerkung 11. Betrachte die zwei Systeme
Kh(uh) = fh
Kh(uh) = fh + ε
43
Dann folgt
uh = K−1h (fh)
uh = K−1h (fh + ε)
⇒ ‖uh − uh‖ ≤ C · ‖ε‖
Stabilität bedeutet also, dass kleine Änderung auf der rechten Seite nur kleine Änderungenin der Lösung bewirken.
Bemerkung 12. Der Fünfpunktstern und somit die Diskretisierung Kh zur Poissonglei-chung hat die Eigenschaft
‖K−1h ‖∞ ≤
1
8,
ist also stabil.
Konsistenz
Sei Khuh = fh die Diskretisierung von Ku = f . K sei ein Differentialoperator derOrdnung m. Weiter seinen Rh und Rh Restriktionsoperatoren für u und f . Die Diskreti-sierung (Kh, Rh, Rh) des Differetialoperators K hat die Konsistenzordnung k bzgl. ‖.‖∞,falls gilt
‖KhRhu− RhKu‖∞ ≤ C · hk · ‖u‖Ck+m(Ω) ∀u ∈ Ck+m(Ω) (3.30)
Beispiel 7. Sei Rh = Rh gegeben durch
(Rhu)(~x) = u(~x) ∀~x ∈ Ωh.
Dann ist (Kh, Rh, Rh) = (∆h, Rh, Rh) konsistent mit Ordnung 2 bzgl. ‖.‖∞.
Beweis. Wir erinnern uns an
(∂−∂+u)(x) = u′′(x) + h2R, |R| ≤ 1
12‖u‖C4(Ω)
Im R2 wenden wir diesen Ansatz in x und y-Richtung an. Taylorentwicklung liefert
−∆hRu(x, y) = −∆u(x, y) + h2(Rx +Ry)
mit |Rx| ≤1
12‖u(4)‖C0(Ω) ≤
1
12‖u‖C4(Ω)
Analog dazu folgert man: |Ry| ≤ 112‖u‖C4(Ω).
⇒ C = 16 mit ‖KhRhu−RhKu‖ ≤ C · h2‖u‖C4(Ω)
Satz 10. (Satz über die Konvergenz)Sei die Diskretisierung (Kh, Rh, Rh) konsistent von der Ordnung k. Sei die zu Kh gehö-rende Matrix Kh stabil bzgl. ‖.‖∞. Dann ist das Verfahren konvergent von der Ordnungk, falls u ∈ Ck+m(Ω). Dabei bezeichne m die Ordnung des Differentialoperators Kh.
44
Beweis. Betrachte wh = uh −Rhu. Unser Ziel ist es zu zeigen, dass
wh → 0 für h→ 0⇒ uh → u für h→ 0.
Es gilt:
Khwh = Khuh −KhRhu = fh −KhRhu = Rhf −KhRhu = RhRu−KhRhu
Da wh|Γh= 0 gilt Khwh = Khwh.
⇒ wh = K−1h (RhKu−KhRhu)
⇒ ‖wh‖∞ = ‖uh −Rhu‖∞ ≤ ‖K−1h ‖∞ · ‖RhKu−KhRhu ‖∞
⇒ ‖uh −Rhu‖∞ ≤ Chk‖u‖Ck+m(Ω)
Frage 6. Was ergibt sich dann für den 5-Punkt-Stern?
Sei u ∈ C4(Ω). Dann gilt:
‖K−1h ‖∞ <
1
8, c =
1
6⇒ ‖uh −Rhu‖∞ ≤
h2
48· ‖u‖C4(Ω)
3.8 Das Neumann-Problem
Bisher wurden die Funktionswerte der gesuchten Funktion auf dem Rand des Gebietesvorgegeben mit u|Γ = ϕ. Stattdessen können auch die Ableitungen von u in Normalen-richtung auf dem Rand vorgegeben werden, d.h. wir betrachten das Problem
−∆u = f in Ω = (0, 1)× (0, 1)
∂u
∂n= ϕ auf Γ
Bemerkung 13. Dadurch, dass ∂u∂n vorgegeben wird, gibt es potentiell unendlich viele
Lösungen. Denn falls u Lösung des Problems ist, so erfüllt auch u+ c das Problem.
Eine eindeutige Lösung erfordert also eine zusätzliche Bedingung um die Variable c ein-deutig festzulegen. Diese Zusatzbedingung nennt man Kompatibilitätsbedingung.
Lemma 12. (Kompatibilitätsbedingung)Es muss gelten ∫
∂Ωϕds = −
∫Ωfdx (3.31)
Beweis.
−∫
Ωfdx =
∫Ω
∆udx =
∫Ωdiv(∇u)dx
=
∫∂Ω∇u · ~nds =
∫∂Ω
∂u
∂~nds =
∫∂Ωϕds
45
3.8.1 Diskretisierung des Neumann-Problems
Wir betrachten das Neumann-Problem auf dem Gebiet Ω = [0, 1] × [0, 1]. Dabei tretendrei Typen von Knoten auf:
1. innere Knoten
2. Knoten auf Kanten
3. Knoten in den Eckpunkten von Ω
Die Diskretisierung des Laplace-Operators führt in diesen drei Fällen zu:
Innere Punkte Für Punkte im Innern, die keine Randknoten zum Nachbarn haben gilt:−∆uh(x, y) = 1
h2(4uh(x, y)− uh(x− h, y)− uh(x+ h, y)− uh(x, y − h)− uh(x, y + h))
Randnahe Punkte Für innere Knoten am rechten Rand von Ω ergibt sich:−∆uh(x, y) = 1
h2(3uh(x, y)− uh(x− h, y)− uh(x, y − h)− uh(x, y + h))
Ecknahe Punkte Für den Knoten rechts unten ergibt sich:−∆uh(x, y) = 1
h2(2uh(x, y)− uh(x− h, y)− uh(x, y + h))
Behandlung der Neumann Randbedingung
Die Ableitung von u in Normalenrichtung ~n kann über einseitige Differenzenquotientenapproximiert werden. Im eindimensionalen Fall lässt sich der rechte Rand folgendermaßendarstellen:
∂uh∂~n
(x) ≈ (∂−n uh)(x) =1
h(uh(x)− u(x− h~n)) = ϕ(x) (3.32)
Für unser Modellbeispiel folgt damit:
1
h(uh(x, 0)− uh(x, h)) = ϕ(x, 0) (unten)
1
h(uh(x, 1)− uh(x, 1− h)) = ϕ(x, 1) (oben)
1
h(uh(0, y)− uh(h, y)) = ϕ(0, y) (links)
1
h(uh(1, y)− uh(1− h, y)) = ϕ(1, y) (rechts)
Das resultierende diskrete Problem
Khuh = fh = fh −1
hϕh
ϕh =
ϕ(x, y) (x, y) ∈ Γ0 sonst
ist nicht notwendigerweise regulär. Deshalb muss zusätzlich die diskrete Kompatibilitäts-bedingung erfüllt sein.
46
Diskrete Kompatibilitätsbedingung
Die diskrete Kompatibilitätsbedingung lautet
−h2∑x∈Ωh
f(x) = h∑x∈Γ
ϕ(x). (3.33)
Satz 11. Ist neben Khuh = fh auch die diskrete Kompatibiltätsbedingung erfüllt, dannist das Problem lösbar und zwei Lösungen unterscheiden sich nur in einer Konstante c.
Beweis. Anwendung des Laplace-Operators auf eine konstante Funktion ergibt Null, d.h.
Kh · c · 1 = 0⇒ c · 1 ∈ ker(Kh)
und
ker(Kh) = span(1)
⇒ dim(ker(Kh)) = 1
Ferner gilt: Khuh = fh lösbar ⇒ fh ∈ Bild(Kh) und
Bild(Kh) = ker(Kh)⊥ = span(1)⊥
⇒ fh ∈⊥ 1 ⇔ fh · 1 = 0⇐⇒∑x∈Ωh
fh(x) = 0
⇒∑x∈Ωh
fh(x) =∑Ωh
fh −1
hϕh ⇐⇒
∑Ωh
fh =1
h
∑∂Ωh
ϕh
3.8.2 Lösen des Neumann-Problems
Die Lösbarkeit des Neumann-Problems setzt, wie oben gezeigt, die Erfüllung der Kom-patibilitätsbedingung voraus. Wir werden zur Vereinfachung der Notation im Folgendenmit fh die rechte Seite mit den Beiträgen auf dem Rand bezeichnen, d.h. fh − 1
h ; fh.
Lemma 13. Man wähle x0 ∈ Ωh beliebig und normiere uh durch
uh(x0) = 0.
Dann ist das um eine Zeile und eine Spalte reduzierte System
Khuh = fh
lösbar.
Beweis.
1. Kh ist symmetrisch (da Kh symmetrisch)
47
2. Kh ist eine M-Matrix, da Kh irreduzibel diagonaldominant ist
Daraus folgt, dass Kh regulär und somit das reduzierte System lösbar ist.
Zur Lösung des Problems wird also eine Korrektur vorgenommen. In obigem Fall durchElimination eines Eintrags aus uh und fh.⇒ Die Korrektur ist konzentriert in fh(x0). Eine Alternative dazu ist, die Korrekturdurch Erweiterung des Systems zu verteilen.
Verteilung der Korrektur
BetrachteKhuh = fh
mit
Kh =
(Kh 1
1T 0
), uh =
(uhλ
), fh =
(fhσ
)mit beliebigem σ.
Lemma 14. Khuh = fh ist eindeutig lösbar.
Beweis. Rang(Kh,1) = Rang(Kh) + 1⇒ Kh ist regulär ⇒ Lösbarkeit.
Lemma 15. Ist uh Lösung von Khuh = fh, so impliziert dies uh ist Lösung von Khuh =fh.
Beweis.
1. λ = 0 impliziert eine Normierungsbedingung für uh mit
1T · uh =
∑x∈Ωh
uh(x) = σ.
2. Sei nun λ 6= 0. uh ist Lösung von Khuh = fh mit fh := fh − λ · 1.
Dadurch wird eine Verteilung der Korrektur auf alle Vektoreinträge erreicht.
48
3.9 Differenzenverfahren für allgemeine Probleme zweiterOrdnung
In diesem Abschnitt betrachten wir allgemeine Probleme zweiter Ordnung und die darausresultierenden diskreten Gleichungen. Wir betrachten nun folgendes Problem
Ku = f in Ω
K =n∑
i,j=1
aij(x)∂2
∂xixj+
n∑i=1
ai(x)∂
∂xi+ a(x) (3.34)
Bemerkung 14. Es gilt aij(x) = aji(x), da ∂2
∂xixj= ∂2
∂xjxifür zweifach differenzierbare
Funktionen.⇒ A(x) = (aij(x))i,j=1...n ist symmetrisch. Beachte, dass aij(x) ortsabhängig sein kön-nen. Der oben definierte Differentialoperator zweiter Ordnung ist in der allgemeinstenOperatorform notiert.
Definition 28. (Elliptizität)Ein Differentialoperator heißt elliptisch, falls alle Eigenwerte des Hauptteils des Diffe-rentialoperators das gleiche Vorzeichen besitzen.
Für den Operator K zweiter Ordnung (siehe oben) gilt aij = aji, d.h. A = AT ⇒ A istpositiv definit. Dies hat zur Folge, dass alle Eigenwerte positiv sind, d.h. K ist elliptisch.
Bemerkung 15. −∆u = F ist eine elliptische partielle Differentialgleichung.
Bemerkung 16. Ein allgemeiner Differentialoperator ist elliptisch, wennn∑
i,j=1
aij(x)ξiξj > 0 ∀x ∈ Ω, 0 6= ξ ∈ Rn
Definition 29. Ein Differentialoperator K heißt gleichmäßig elliptisch in Ω, falls giltn∑
i,j=1
aij(x)ξiξj ≥ c(x)|ξ|2, c(x) > 0 ∀x ∈ Ω, 0 6= ξ ∈ Rn
Im zweidimensionalen Fall auf dem Gebiet Ω = (0, 1)× (0, 1) erhalten wir für randfernePunkte den diskreten Operator
a11(x, y)∂+x ∂−x + 2a12∂
0x∂
0y + a22(x, y)∂+
y ∂−y + a1(x, y)∂0
x + a2(x, y)∂0y + a(x, y)
= h−2
−12a12(x, y) a22(x, y) 1
2a12(x, y)a11(x, y) −2(a11(x, y) + a22(x, y)) a11(x, y)
12a12(x, y) a22(x, y) −1
2a12(x, y)
+
+1
2h
0 a2(x, y) 0−a1(x, y) 0 a1(x, y)
0 a2(x, y) 0
+
0 0 00 a(x, y) 00 0 0
Upwind-Methods missing (10.-11.12.2012)
49
4 Diskretisierung II: Finite ElementeVerfahren
Während das Diskretisierungsverfahren der Finiten Differenzen einen einfachen mathe-matischen Zugang über die direkte Approximation der Differentialoperatoren bietet, inder Implementierung verhältnismäßig leicht zu handhaben ist, so gibt es jedoch gewis-se Einschränkungen aufgrund des fundamentalen Vorgehens in diesem Verfahren. ZumBeispiel erfordert die Behandlung komplexer Gebiete mit krummen Rändern besonderenumerische Behandlung, was wiederum die Gesamtkomplexität beeinflusst. Ferner sahenwir in den Konvergenzbeweisen für das Finite Differenzen Verfahren sehr hohe Regulari-tätsannahmen für die Lösungsfunktion, die eine starke Einschränkung darstellen kann.Um eine gebietsunabhängige Diskretisierungsherangehensweise aufzubauen, sowie einVerfahren das substantiell schwächere Regularitätsvoraussetzungen stellt, werden wireinen funktionalanalytischen Weg finden. Anstelle der direkten Approximation von Diffe-rentialoperatoren tritt die Approximation der Lösungsräume in denen die Lösungsfunk-tion liegt durch endlich dimensionale, approximierende Funktionenräume. Ein darausresultierendes numerisches Verfahren ist das Verfahren der Finite Elemente, welchem wiruns in diesem Kapitel widmen werden.
4.1 Funktionalanalytische Grundlagen
Wie oben schon erwähnt, werden wir über die Funktionalanalysis eine neue Kategorie derDiskretisierung herleiten. Dazu sind einige funktionalanalytische Grundlagen notwendig,die in diesem Abschnitt zusammengestellt werden.
4.1.1 Normierte Räume
Definition 30. Sei X ein Vektorraum über R oder C und ‖.‖ : X → [0,∞) eine Norm.Dann wird
(X, ‖.‖X)
normierter Raum genannt.
Beispiel 8. Die stetigen Funktionen auf Ω bilden den normierten Raum C0(Ω) mit derSupremumsnorm ‖.‖∞.
Definition 31. Zwei Normen ‖.‖(1) und ‖.‖(2) auf X heißten äquivalent, wenn 0 < C <∞ existiert, mit
1
C‖x‖(1) ≤ ‖x‖(2) ≤ C‖x‖(1) ∀x ∈ X (4.1)
50
Operatoren
Definition 32. (Operatornorm)Seien X und Y normierte Räume mit Normen ‖.‖X und ‖.‖Y . Die Operatornorm einesOperators T : X → Y ist definiert als
‖T‖Y←X := supx∈X
‖Tx‖Y‖x‖X
: 0 6= x ∈ X
(4.2)
Bemerkung 17. Ist ‖T‖Y←X endlich, dann ist T beschränkt.
Bemerkung 18. Die beschränkten Operatoren bilden einen linearen Raum L(X,Y ) mit(T1 + T2)x = T1x+ T2x.
Offene Räume
Definition 33. (X, ‖.‖) sein ein normierter Raum. A ⊂ X heißt offen, falls für allex ∈ A ein ε > 0 existiert, sodass
Kε(x) := y ∈ X : ‖x− y‖ < ε
in A enthalten ist.
4.1.2 Banach-Räume
Sei (X, ‖.‖) ein normierter Raum. X heißt vollständig, wenn jede Cauchy-Folge konver-giert. Ein normierter und vollständiger Raum heißt Banach-Raum
Definition 34. (Cauchy-Folge)Eine Folge xh ∈ X : n ≥ 1 heißt Cauchy-Folge, wenn gilt
sup‖xn − xm‖X : n,m ≥ k → 0, für k →∞ (4.3)
oder alternativ∀ε > 0 ∃n0 ∈ N : ∀n,m ≥ n0 : ‖xn − xm‖ < ε (4.4)
Der Raum(L∞(D), ‖.‖L∞(D)
)L∞(D) bezeichnet den Raum der auf D beschränkten lokal integrierbaren Funktionen.L∞(D) besteht aus Äquivalenzklassen, wobei
f = g, falls f = g fast überall
‖u‖L∞(D) := infA
sup
x∈D\A|u(x)| : x ∈ D \A : µ(A) = 0
Wir werden im Folgenden erklären, was mit obiger Terminologie gemeint ist.
51
Definition 35. (fast überall)Zwei Funktionen f und g heißen fast überall gleich, falls
f = g bis auf Nullmenge A giltf = g auf Ω \A
gilt.
Definition 36. (Nullmenge)Eine Nullmenge besitzt das Maß µ(A) = 0.
Definition 37. (Maß)Sei J n die Menge der beschränkten Intervalle des Rn. Eine Intervallfunktion ist eineAbbildung
ϕ : J n → R
1. ϕ ist monoton, wenn gilt
I1, I2 ∈ J n und I1 ⊂ I2 ⇒ ϕ(I1) ≤ ϕ(I2).
2. ϕ heißt additiv, wenn gilt
I1, I2, I ∈ J n mit I1 ∩ I2 = ∅, I1 ∪ I2 = I
⇒ ϕ(I) = ϕ(I1) + ϕ(I2)(⇒ ϕ monoton und additiv ⇒ ϕ(∅) = 0 und ϕ(I) ≥ 0, ∀I ∈ J n).
3. ϕ heißt regulär, falls gilt: Für jedes ε > 0 gibt es zu jedem I ∈ J n ein offenesIntervall I∗ mit I ⊂ I∗, sodass gilt:
ϕ(I) ≤ ϕ(I∗) < ϕ(I) + ε
Eine monotone, additive, reguläre Intervall-Funktion heißt Maß.
Beispiel 9. Folgende gängige Maße können definiert werden:
1. Sei I = [a1, b1]× . . .× [an, bn] ∈ J n. Dann ist
v : J n → R mit v(I) = Πnj=1(bj − aj)
ein Maß (Volumenmaß).
2. Für treppenfunktion-approximierbare Funktionen f ist∫fdϕ =
∫Rn
fdϕ =m∑i=1
ϕ(Ij)f(Ij)
auf den Intervallen I1, . . . , Im das Lebegue-Maß.
Definition 38. (meßbar)f : Rn → R heißt meßbar, wenn eine Folge von Treppenfunktionen snn=1,2,... existiert,mit
f = lim sn
.
52
4.1.3 Der Sobolev-Raum L2(Ω)
Definition 39. (Skalarprodukt)Die Abbildung (., .) : X ×X → K heißt Skalarprodukt auf X, falls gilt
(x, x) > 0 ∀x ∈ X,x 6= 0 (4.5)(λx+ y, z) = λ(x, z) + (y, z) ∀λ ∈ K, x, y, z ∈ X (4.6)
(x, y) = (y, x) (4.7)
Das Skalarprodukt definiert eine Norm durch
‖x‖ :=√
(x, x).
Definition 40. (Hilbertraum)Ein Banach-Raum X heißt Hilbert-Raum, wenn ein Skalarprodukt (., .)X auf X existiert.
Der Sobolev-Raum L2(Ω)
Sei Ω eine offene Teilmenge von Rn. L2(Ω) definiert den Raum der Äquivalenzklassenmeßbarer und quadrat-integrabler Funktionen:
L2(Ω) := f : Ω→ R : f meßbar, |f |2 integrierbar (4.8)
Zwei Funktionen f und g sind gleich, wenn sie bis auf in A mit µ(A) = 0 gleich sind. Mitdem Skalarprodukt
(u, v)0 = (u, v)L2(Ω) :=
∫Ωu(x)v(x)dx ∀u, v ∈ L2(Ω) (4.9)
und Norm
‖u‖0 = ‖u‖L2(Ω) :=
√∫Ω|u(x)|2 (4.10)
ist L2(Ω) ein Hilbertraum.
4.1.4 Schwache Differenzierbarkeit
In L2(Ω) können keine Aussagen über die Differenzierbarkeit von f ∈ L2(Ω) im klassi-schen Sinne gemacht werden. Aus der partiellen Integration folgernd gilt jedoch:∫
Ωf ′(x)ϕ(x)dx = −
∫Ωf(x)ϕ′(x)dx
für Funktionen ϕ mit ϕ|∂Ω = 0.
Definition 41. Falls für f ∈ L2(Ω) eine Funktion g ∈ L2(Ω) existiert, sodass gilt
−∫
Ωf ′(x)ϕ(x)dx = −
∫Ωg(x)ϕ(x)dx =
∫Ωf(x)ϕ′(x)dx (4.11)
dann heißt g die schwache Ableitung von f .
53
Bemerkung 19. Folgende Zusammenhänge bestehen zwischen klassischer und schwacherDifferenzierbarkeit:
1. Klassisch differenzierbare Funktionen sind auch schwach differenzierbar.
2. Die schwache Ableitung ist durch die integrale Definition nicht punktweise definiert.
3. Hinreichend oft schwach differenzierbar impliziert klassische Differenzierbarkeit (sie-he Sobolevscher Einbettungssatz).
Höhere schwache Differenzierbarkeit
Sei α ein Multiindex α = (α1, . . . , αn) mit
|α| =n∑i=1
αi
Dα :=∂|α|
∂α1x1 . . . ∂αnxn
und f ∈ L2(Ω). Dann heißt g α-fache schwache Ableitung von f , wenn gilt∫Ωg(x)ϕ(x)dx = (−1)|α|
∫Ωf(x)Dαϕ(x)dx (4.12)
(g, ϕ)0 = (−1)|α|(f,Dαϕ) ∀ϕ ∈ C∞(Ω), ϕ|∂Ω = 0 (4.13)
4.1.5 Die Hilbert-Räume Hk(Ω) und Hk0 (Ω)
Die Menge aller Funktionen u aus L2(Ω), die schwache Ableitungen Dαu ∈ L2(Ω) besit-zen, bilden den Sobolev-Raum
Hk(Ω) := u ∈ L2(Ω) : Dαu ∈ L2(Ω) für |α| ≤ k (4.14)
mit k ∈ N0. Hk(Ω) wird in der Literatur auch als W k2 (Ω) bzw. W k,2(Ω) bezeichnet.
Hk(Ω) ist ein Hilbert-Raum mit dem Skalarprodukt
(u, v)k := (u, v)Hk(Ω) :=
√∑|α|≤k
‖Dαu‖2L2(Ω)
(4.15)
Satz 12. (Sobolevscher Einbettungssatz)Es gilt
Hs(Rn) ⊂ Ck(Rn) falls s > k +n
2, k ∈ N0 (4.16)
54
4.1.6 Dualräume
Sei X ein normierter, linearer Raum über R. Mit X ′ wird der Dualraum bezeichnet,bestehend aus allen beschränkten, linearen Abbildungen von X nach R:
X ′ = L(X,R)
Bemerkung 20. Mit der Norm
‖x′‖ := ‖x′‖R←X := sup
|x′(x)|‖x‖X
: 0 6= x ∈ X
ist X ′ ein Banachraum. x′ ∈ X ′ heißen lineare Funktionale auf X:
x′(x) = 〈x, x′〉X×X′
Lemma 16. Seien X und Y normiert und T ∈ L(X,Y ). Für jedes y′ ∈ Y ′ definiert
〈Tx, y′〉Y×Y ′ = 〈x, x′〉X×X′ ∀x ∈ X (4.17)
ein eindeutiges x′ ∈ X ′. Der zugehörige Dualoperator ist durch
T ′ : Y ′ −→ X ′
y′ 7−→ x′
definiert. D.h. 〈Tx, y′〉Y×Y ′ = 〈x, T ′y′〉X×X′ = 〈x, x′〉X×X′.
Lemma 17. Es gilt‖T ′‖X′←Y ′ = ‖T‖Y←X (4.18)
Beweis.
‖T ′‖X′←Y ′ = supy′ 6=0
‖T ′y′‖X′‖y′‖Y ′
= sup
x,y 6=0
〈x, T ′y′〉X×X′‖x‖X‖y′‖Y ′
= sup
x,y 6=0
〈Tx, y′〉Y×Y ′‖x‖X‖y′‖Y ′
= sup
x′ 6=0
‖Tx‖Y‖x‖X
= ‖T‖Y←X
Adjungierte Operatoren
Sei X ein Hilbert-Raum über R. Jedes y ∈ X definiert durch
fy(x) := (x, y)X (4.19)
ein lineares Funktional fy ∈ X ′ mit ‖fy‖X′ = ‖y‖X . Umgekehrt definiert jedes Funktionalfy ein y ∈ X, wie folgender Satz zeigt:
Satz 13. (Darstellungssatz von Riesz)X sei ein Hilbert-Raum und f ∈ X ′ ein Funktional. Dann existiert genau ein yf ∈ X,sodass
f(x) = (x, y)X ∀x ∈ X, ‖f‖X′ = ‖yf‖X′ (4.20)
55
Beweis. X ist abgeschlossen. Sei N = x ∈ X : f(x) = 0 der Kern von f .
1. N = X: Wähle y = 0 ⇒ Behauptung.
2. Sei N ein abgeschlossener Teilraum von X. Urbilder abgeschlossener Mengen unterstetigen Funktionen sind abgeschlossen und auf abgeschlossenen Mengen existiertein z ∈ X mit
z /∈ N und 〈z, n〉 = 0 ∀n ∈ N.
Ferner gilt:
f
(x− f(x)z
f(z)
)= f(x)− f
(f(x)z
f(z)
)= f(x)− f(x)
f(z)f(z) = 0
⇒ x− f(x)z
f(z)∈ N
Eingesetzt liefert
〈x− f(x)
f(z)z, z〉 = 0⇔ 〈x, z〉 − f(x)
f(z)〈z, z〉 = 0
⇒ f(x) =〈x, z〉‖z‖2
f(z) = 〈x, f(z)z
‖z‖2〉
Mit y := f(z)‖z‖2 z folgt die Behauptung, bis auf die Eindeutigkeit.
Eindeutigkeit: Sei y ein weiteres Element mit
f(x) = 〈x, y〉 = 〈x, y〉⇒ 〈x, y − y〉 = 0 ∀x ∈ X⇒ y − y = 0⇒ Eindeutigkeit
Bilinearformen
Definition 42. Sei V ein Hilbertraum. Die Abbildung a(., .) : V × V −→ R heißt Bili-nearform, falls
a(x+ λy, z) = a(x, z) + λa(y, z) (4.21)a(x, y + λz) = a(x, y) + λa(x, z) ∀λ ∈ R, x, y, z ∈ V (4.22)
Definition 43. a(., .) heißt stetig, falls ein Cs existiert, sodass
|a(x, y)| ≤ Cs‖x‖V ‖y‖V ∀x, y ∈ V (4.23)
56
Lemma 18. Zu einer stetigen Bilinearform existiert ein eindeutiger Operator A ∈L(V, V ′) mit
a(x, y) = 〈Ax, y〉V ′×V ∀x, y ∈ V (4.24)‖A‖V ′←V ≤ Cs (4.25)
Beweis. Man halte x ∈ V fest.ϕx(y) := a(x, y)
ist ein Funktional ϕx ∈ V ′ mit ‖ϕx‖V ′ ≤ Cs‖x‖V . Da ϕx über die Bilinearform definiertist, ist ϕx linear.
Ax := ϕx
für x ∈ V und es gilt mit
‖Ax‖V ′ ≤ Cs‖x‖V ⇔‖Ax‖V ′‖x‖V
≤ Cs
auch‖A‖V ′←V := sup
‖Ax‖V ′‖x‖V
≤ Cs
Analog dazu: Halte y ∈ V fest.
Definition 44. (V-Elliptizität)Eine Bilinearform heißt V-elliptisch, falls sie auf V × V stetig ist und CE > 0 existiert,sodass
a(x, x) ≥ CE‖x‖2V ∀x ∈ V (4.26)
4.2 Variationsformulierung
Das Ziel in diesem Abschnitt ist es, über die Funktionalanalysis einen neuen Zugang zumModellproblem zu finden. Ein äquivalentes Problem zum klassischen Problem definie-ren wir über eine Variationsformulierung. Diese wird Lösungen des Modellproblems inHilberträumen angeben, die bis auf Nullmengen auch klassische Lösungen sein werden.
Satz 14. (Charakterisierungssatz)Sei V ein linearer Raum und
a : V × V −→ R
symmetrisch und positiv. f : V −→ R sei ein lineares Funktional. Die Abbildung
J(v) :=1
2a(v, v)− f(v)
nimmt in V ihr Minimum genau dann bei u an, wenn gilt:
a(u, v) = f(v) ∀v ∈ V
Die Lösung u ist eindeutig.
57
Beweis. Seien u, v ∈ V, t ∈ R. Betrachte
J(u+ tv) =1
2a(u+ tv, u+ tv)− f(u+ tv) =
=1
2
(a(u, u) + 2ta(u, v) + t2a(v, v)
)− f(u)− tf(v)
= J(u) + t(a(u, v)− f(v)) +1
2t2a(v, v)
1. Gelte für u ∈ V : a(u, v) = f(v)
t=1⇒ J(u+ v) = J(u) + (f(v)− f(v)) +1
2a(v, v)
= J(u) +1
2a(v, v) > J(u)
⇒ u ist Minimalpunkt.
2. Sei u ∈ V Minimalpunkt.
J(u+ tv) = J(u) + t(a(u, v)− f(v)) +1
2t2a(v, v)
⇒ 0 =dJ(u+ tv)
dt= a(u, v)− f(v) + ta(v, v)
⇒ J ′(u+ tv)|t=0 = a(u, v)− f(v)
⇔ a(u, v) = f(v)
⇒ a(u, v) = f(v)⇐⇒ u ist Minimalpunkt.
4.2.1 Untersuchung des elliptischen Differentialoperators zweiter Ordnung
Wir betrachten den elliptischen Differentialoperator
Lu := −n∑
i,k=1
∂i(aik∂ku) + a0u
mit a0(x) ≥ 0 (x ∈ Ω). Das zugehörige Problem
Lu = f in Ω
u = g auf ∂Ω
kann ohne Beschränkung der Allgemeinheit als homogenes Problem aufgefasst werden:
w := u− g⇒ Lw = f1 in Ω
w = 0 auf ∂Ω
Folgender Satz stellt den Zusammenhang zwischen Randwertaufgabe und Variationspro-blem her:
58
Satz 15. (Minimaleigenschaft)Jede klassische Lösung der Randwertaufgabe
−∑i,k
∂i(aik∂ku) + a0u = f in Ω (4.27)
u = 0 auf ∂Ω (4.28)
ist Lösung des Minimierungsproblems
J(v) :=
∫Ω
1
2
∑i,k
aik∂iv∂kv +1
2a0v
2 − fv
dx→ min (4.29)
unter allen Funktionen aus C2(Ω) ∩ C0(Ω) mit Nullrandwerten.
Beweis. Die Greensche Formel besagt∫Ωv∆u+∇v∇udx =
∫∂Ωv∂u
∂~nds
Angewandt auf∫
Ω v∂iwdx mit w := aik∂ku liefert∫Ωv∂iw + ∂ivwdx =
∫∂Ωv∂w
∂~nds
Mit v|∂Ω = 0 folgt ∫Ωv∂iw = −
∫Ω∂ivw
und mit w = aik∂ku ∫Ωv∂i(aik∂ku)dx = −
∫Ωaik∂iv∂kudx.
Setze
a(u, v) :=
∫Ω
∑i,k
aik∂iu∂kv + a0uvdx
f(v) :=
∫Ωfvdx
und summiere über i und j:
⇒∫
Ω
∑i,j
v∂i(aik∂ku)dx = −∫
Ω
∑i,j
aik∂iv∂kudx︸ ︷︷ ︸=a(u,v)−a0uv
Eine Erweiterung um∫
Ω fvdx liefert schließlich
a(u, v)− f(v) =
∫Ωv(−∑
∂i(aik∂ku) + a0u− f)dx
=
∫Ωv(Lu− f)dx =
falls Lu=f0
59
Der Charakterisierungssatz besagt also
u ist Lösung des Variationsproblems und
u ∈ C2(Ω) ∩ C0(Ω)
⇓u ist klassische Lösung
Wir klären nun die Frage nach der Existenz einer Lösung.
4.2.2 Existenz und Eindeutigkeit für das Variationsproblem
Betrachte das Dirichlet-Integral (Energiefunktional)
J(u) :=
∫Ω|∇u|2dx =
∫Ω
n∑i=1
u2xi(x)dx (4.30)
Aus der Minimaleigenschaft folgt
J(u)→ min⇐⇒ ∆u = 0 in Ω
u = ϕ auf Γ
Dirichlet-Prinzip
Dirichlet argumentierte, da J(u) ≥ 0 muss für ein u das Minimum angenommen werden.⇒ Es existiert eine Lösung zu
∆u = 0 in Ω
u = ϕ auf Γ
Einen Widerspruch zu dieser Aussage demonstrierte Weierstraß (1870):Betrachte dazu
J(u) =
∫ 1
0u2(x)dx −→ min für u ∈ C0([0, 1])
und u(0) = 1, u(1) = 0.Wir definieren nun die Funktionenfolge
un(x) =
1− nx 0 ≤ x ≤ 1
n0 x > 1
n
(4.31)
Es gilt limn→0 un = 0, d.h. das Infimum von J(u) ist
inf J(u) = 0,
wird aber für stetige Funktionen nicht angenommen. Dies ist ein Widerspruch zumDirichlet-Prinzip. Dieses Beispiel zeigt, um eine eindeutige Lösung des Variationspro-blems zu garantieren, muss der Funktionenraum richtig gewählt werden. Dies besagt derfolgende Satz:
60
1
1
u1(x)
1
n
Abbildung 4.1: Folgenfunktionen un(x)
Existenz und Eindeutigkeit
Satz 16. (Lax-Milgram)Sei V eine abgeschlossene, konvexe Menge in dem Hilbertraum H und a : H ×H −→ Reine elliptische Bilinearform. Für jedes l aus H ′ hat das Variationsproblem
J(v) :=1
2a(v, v)− 〈l, v〉 −→ min (4.32)
genau eine Lösung in V .
Beweis. J ist nach unten beschränkt:
J(v) ≥ 1
2CE‖v‖2 − ‖l‖‖v‖
=1
2CE‖v‖2 − ‖l‖‖v‖+ ‖l‖2 − ‖l‖2
=1
2CE(C2
E‖v‖2 − 2CE‖l‖‖v‖+ ‖l‖2)− 1
2CE‖l‖2
=1
2CE(CE‖v‖ − ‖l‖)2 − 1
2CE‖l‖2 ≥ −‖l‖
2
2CE
Wir setzen nun c1 := infJ(v) : v ∈ V . Sei (vn) eine Minimalfolge, d.h.
limn→∞
vn = c1.
Es gilt: CE‖vn − vm‖2 ≤ a(vn − vm, vn − vm), da a(., .) elliptisch ist. Mit Hilfe derParallelogramgleichung
‖vn + vm‖2 + ‖vn − vm‖2 = 2(‖vn‖2 + ‖vm‖2) (4.33)
61
folgt
a(vn + vm, vn + vm) + a(vn − vm, vn − vm) = 2(a(vn, vn) + a(vm, vm))
⇔ a(vn − vm, vn − vm) = 2a(vn, vn) + 2a(vm, vm)− a(vn + vm, vn + vm)
Angewandt auf
CE‖vn − vm‖2 ≤ a(vn − vm, vn − vm)
= 2a(vn, vn) + 2a(vm, vm)− a(vn + vm, vn + vm)
Mit J(v) = 12a(v, v)− 〈l, v〉 erhalten wir
4J(vn) = 2a(vn, vn)− 4〈l, vn〉
= 4J(vn)− 4〈l, vn〉+ 4J(vm)− 4〈l, vm〉 −(
8 · J(vn + vm
2
)− 4〈l, vn + vm〉
)= 4J(vn) + 4J(vm)− 8J
(vn + vm
2
)≤ 4J(vn) + 4J(vm)− 8c1
Die Forderung, dass V konvex ist, liefertvn + vm
2∈ V
(vn) ist eine Minimalfolge ⇒ limn→∞ 4J(vn,m) = 4c1.
⇒ CE‖vn − vm‖2 −→ 0 für n,m→∞,
also ‖vn − vm‖2 −→ 0.⇒ (vn) ist eine Cauchy-Folge in H. H ist ein Hilbertraum, deshalb existiert ein u ∈ Hund mit V abgeschlossen auch u ∈ V mit limn→∞ vn = u und J(u) = limn→∞ J(vn) =infv∈V J(v).⇒ J(v) = 1
2a(v, v)− 〈l, v〉 −→ min hat eine Lösung u ∈ V .
Eindeutigkeit: Seien u1 und u2 Lösungen, also Grenzwerte von Minimalfolgen. Dannist auch u1, u2, u1, u2, . . . eine Minimalfolge.⇒ u1, u2, u1, u2, . . . ist Cauchy-Folge. Dies ist aber nur der Fall wenn u1 = u2 ⇒ Eindeu-tigkeit.
Bemerkung 21. Folgende Schlussfolgerung können gemacht werden:
1. Mit V = H folgt: Zu jedem l ∈ H ′ gibt es ein u ∈ H mit
a(u, v) = 〈l, v〉 ∀v ∈ H
2. Sei a(u, v) := (u, v). Der Rieszsche Darstellungssatz liefert: Zu jedem l ∈ H ′ gibtes ein u ∈ H mit
(u, v) = 〈l, v〉 ∀v ∈ HDamit erhält man eine Einbettung
H ′ −→ H
l 7−→ v
62
4.2.3 Schwache Lösung des Randwertproblems
Bei Finite Differenzenverfahren zeigte sich in den Konvergenzbeweisen die Forderungnach starker Regularität der Lösungsfunktionen. Differenzierbarkeit der Lösung wur-de im klassischen Sinne aufgefasst. Um die funktionalanalytischen Eigenschaften vonHilbert-Räumen mit Integrabilitätsbedingungen (statt klassische Differenzierbarkeit) zunutzen führen wir hier den Begriff der schwachen Differenzierbarkeit ein. Wir werdensehen, dass dieser Begriff in großen Teilen von Ω mit der klassischen Differenzierbarkeitübereinstimmen kann, Ausnahmemengen jedoch erlaubt sind, in denen Lösungsfunktio-nen nicht differenzierbar sein müssen. Dies ermöglich es uns später numerisch einfachhandhabbare Funktionenräume zu definieren.
Definition 45. Eine Funktion u ∈ H10 (Ω) heißt schwache Lösung der Randwertaufgabe
2. Ordnung
Lu = f in Ω
u = 0 auf Γ
wenn gilt:
a(u, v) = (f, v)0 ∀v ∈ H10 (Ω) (4.34)
a(u, v) :=
∫Ω
∑i,k
aik∂iu∂kv + a0uvdx (4.35)
Satz 17. Sei L ein Differentialoperator 2. Ordnung. Dann hat das Randwertproblem
Lu = f in Ω
u = 0 auf Γ
stets eine schwache Lösung in H10 (Ω) und ist Minimum des Problems
1
2a(v, v)− (f, v)0 −→ min in H1
0 (Ω)
Beispiel 10. Für das Modellproblem
−∆u = f in Ω
u = 0 auf Γ
ist die zugehörige Bilinearform
a(u, v) =
∫Ω∇u∇vdx.
Das Variationsproblem lässt sich dann folgendermaßen stellen:Finde u ∈ H1
0 (Ω), sodass
(∇u,∇v)0 = (f, v)0 ∀v ∈ H10 (Ω)
Dieses u findet man als Lösung des Minimierungsproblems1
2
∫Ω∇u∇vdx− (f, v)0 −→ min
63
4.2.4 Variationsproblem der Neumann-Randwertaufgabe
Wir betrachten das Problem
Lu = f in Ω∑i,k
niaik∂ku = g auf Γ
wobei ni der i-te Anteil des lokalen Normalenvektors ist. Mit f ∈ L2(Ω) und g ∈ L2(Γ)ist das lineare Funktional
〈l, v〉 :=
∫Ωfvdx+
∫Γgvdx
definiert. Das Variationsproblem hat dann die Gestalt:Finde u ∈ H1(Ω), sodass
a(u, v) = (f, v)0,Ω + (g, v)0,Γ ∀v ∈ H1(Ω)
gilt.
Satz 18. Die Variationsaufgabe auf Ω mit stückweise glattem Rand und erfüllter Kegel-bedingung
J(v) :=1
2a(v, v)− (f, v)0,Ω − (g, v)0,Γ −→ min
ist äquivalent zu
Lu = f in Ω∑i,k
niaik∂ku = g auf Γ
falls u ∈ C2(Ω) ∩ C1(Ω).
Beachte: Die Variationsaufgabe hat in H1(Ω) eine eindeutige Lösung (nicht notwendi-gerweise in C2(Ω) ∩ C1(Ω)).
4.3 Galerkin-Verfahren
Um die schwache Lösung der Randwertaufgabe zu finden muss ein Minimierungsproblemgelöst werden, und zwar über einen unendlich dimensionalen Hilbertraum. Um einennumerischen Ansatz für dieses Problem zu finden, ist die Idee der Galerkin-Verfahren,den Funktionenraum endlich-dimensional zu approximieren. Wir werden also z.B. H1(Ω)durch einen approximierten Raum Vh ersetzen.Wir betrachten dazu die VariationsaufgabeFinde u ∈ H1(Ω) mit
a(u, v) = (f, v).
64
Beispiel 11. Sei die Randwertaufgabe
−∆u = f in Ω
u = u0 auf Γ
gegeben. Die schwache Formulierung lautet:Suche u ∈ H1
0 (Ω) mit ∫Ω∇u∇vdx =
∫Ωfvdx ∀v ∈ H1
0 (Ω)
mit
a(u, v) :=
∫Ω∇u∇vdx
(f, v) :=
∫Ωfvdx
Problem. Hier muss ein Minimierungsproblem in einem unendlich-dimensionalen RaumHm(Ω) bzw. Hm
0 (Ω) gelöst werden. Für die numerische Behandlung des Minimierungs-problems benötigen wir einen endlich-dimensionalen Raum.
Idee 6. Ersetze den Lösungsraum V durch Vh, wobei Vh endlich-dimensional ist, alsodurch eine endliche Basis darstellbar ist.
Damit wird das obige Minimierungsproblem zu
J(v) :=1
2a(v, v)− 〈l, v〉 −→ min
Vh(4.36)
uh ∈ Vh ist Lösung des Minimierungsproblems, wenn gilt
a(uh, v) = 〈l, v〉 ∀v ∈ Vh.
Das endliche Problem
Sei ψ1, ψ2, . . . , ψn eine Basis von Vh. Dann ist
a(uh, v) = 〈l, v〉 ∀v ∈ Vh
äquivalent zua(uh, ψi) = 〈l, ψi〉 i = 1, 2, . . . , N.
uh ∈ Vh lässt sich als Linearkombination der ψi darstellen:
uh =N∑k=1
zkψk (4.37)
65
mit zu berechnenden Koeffizienten zk. Dies führt zu einem GleichungssystemN∑k=1
a(ψk, ψi)zk = 〈l, ψi〉 i = 1, 2, . . . N
durch Einsetzen von uh =∑N
k=1 zkψk in a(uh, ψi) = 〈l, ψi〉 ∀i = 1, . . . , N . Mit Aik :=a(ψk, ψi) und bi := 〈, ψi〉 lässt sich das Gleichungssystem schreiben als
Az = b
Bemerkung 22. Falls a(., .) V -elliptisch ist, so ist die Matrix A positiv definit:
ztAz =∑i,k
ziAikzk = a
(∑k
zkψk,∑i
ziψi
)= a(uh, uh) ≥ CE‖uh‖2V
Frage 7. Wie gut approximiert uh ∈ Vh die Lösung u ∈ V ?
Lemma 19. (Céa-Lemma)Die Bilinearform a sei V -elliptisch und Hm
0 (Ω) ⊂ V ⊂ Hm(Ω). Dann gilt:
‖u− uh‖m ≤CSCE
infvh∈Vh
‖u− vh‖m (4.38)
Beweis. Es gilt:
a(u, v) = 〈l, v〉 ∀v ∈ Va(uh, v) = 〈l, v〉 ∀v ∈ Vh
Da Vh ⊂ V können wir schreiben
a(u, v)− a(uh, v) = a(u− uh, v) = 0 für v ∈ Vh.
Mit vh ∈ Vh können wir v ∈ Vh schreiben als
v = vh − uh ∈ Vh.
⇒ a(u− uh, vh − uh) = 0. Ferner ist a(., .) elliptisch, d.h.
CE‖u− uh‖2m ≤ a(u− uh, u− uh) = a(u− uh, u− uh) + a(u− uh, vh − uh)
≤ CS‖u− uh‖m‖u− vh‖mDaraus schließen wir
‖u− uh‖m ≤ CSCE‖u− vh‖m ∀vh ∈ Vh
⇒ ‖u− uh‖m ≤ CSCE
infvh∈Vh
‖u− vh‖m
Bemerkung 23. Aufgrund von a(u − uh, v) = 0 ist der Approximationsfehler u − uhorthogonal zu V .
66
Folgerung aus dem Céa-Lemma
Je besser Vh den Raum V approximiert, desto kleiner wird der Abstand zwischen u unduh
‖u− uh‖m.
Frage 8. Wie wählen wir Vh, bzw. die Basis von Vh?
4.4 Finite Elemente Verfahren
Die Frage nach einer geeigneten Basis für Vh beantwortet das Finite Elemente Verfahrendurch
1. Basisfunktionen niedriger Ordnung
2. lokale Funktionen, d.h. Funktionen mit kompaktem Träger
Die Approximationsgüte von Funktionen mit niedriger Ordnung ist zwar nur in kleinenBereichen vertretbar, durch die Einfachheit der Funktionen kann eine bessere Approxi-mation jedoch durch Gitterverfeinerung, auch in Kombination mit Erhöhung der Poly-nomordnung, erreicht werden. Durch Funktionen mit kompaktem Träger lassen sich dieEinträge der Steifigkeitsmatrix und der rechten Seite durch lokale Integralapproximationeinfach lösen und so das finale lineare Gleichungssystem aufstellen.Dieses kann mitunter sehr groß werden, da die Gitterfeinheit aufgrund der obigen Grun-dannahmen für die Basisfunktionen hoch werden kann. Courant demonstrierte 1943 fol-gendes grundlegende Beispiel für das Vorgehen bei der Finite Elemente Methode.
4.4.1 Beispiel von Courant
Betrachte das Poisson-Problem
−∆u = f in Ω = (0, 1)× (0, 1)
u = 0 auf Γ
Ω wird dabei in 8 Dreiecke I-VIII zerlegt, mit den Knoten L, R, O, U, LO, RU und Z.Wir wählen Vh folgendermaßen
Vh = v ∈ C(Ω) : v linear und v|Γ = 0.
v lässt sich also in jedem Gitterdreieck darstellen als
v(x, y) = a+ bx+ cy.
Da wir die Werte in den Gitterknoten kennen, lassen sich a, b und c eindeutig bestimmen.Bei N inneren Gitterknoten ist dimVh = N . Wir benötigen also N Basisfunktionen fürVh. Die Basisfunktionen ψiNi=1 seien definiert durch
ψi(Kj) = δij mit Kj = Knoten(j), j = 1, . . . , N.
Die Basisfunktion über dem inneren Knoten Z ist damit
67
I
II
III
IV
V
VI
VII
VIII
L R
O
U
LO
RU
Z
Abbildung 4.2: Unterteilung des Gebiets in Dreiecke
I II III IV V VI VII VIII
∂1ψZ−1h 0 1
h 0 0 1h 0 −1
h
∂2ψZ−1h 0 0 −1
h 0 1h
1h 0
Tabelle 4.1: Ableitungen von ψZ in die erste und zweite Raumrichtung
1. linear in jedem Dreieck
2. Null auf den umliegenden Knoten
erfüllt also die Finite Elemente Richtlinien niedrige Ordnung und lokal. Wir können nundie Ableitungen der Basisfunktion ψZ in den Dreiecken I-VIII berechnen (siehe Tabelle4.4.1)Zu lösen ist das Gleichungssystem
Au = b
mitAij = a(ψi, ψj).
In unserem Beispiel benötigen wir also a(ψZ , ψZ), a(ψZ , ψO), a(ψZ , ψU ), a(ψZ , ψL),a(ψZ , ψR), a(ψZ , ψLO) und a(ψZ , ψRU ).
68
Für a(ψZ , ψZ) gilt:
a(ψZ , ψZ) =
∫Ω
(∇ψZ)2dxdy =
∫I−V III
(∇ψZ)2dxdy
= 2 ·∫I+III+IV
((∂1ψZ)2 + (∂2ψZ)2
)= 2 ·
∫I+III
(∂1ψZ)2dxdy + 2 ·∫I+IV
(∂2ψZ)2dxdy
=2
h2
∫I+III
dxdy +2
h2
∫I+IV
dxdy
⇒ a(ψZ , ψZ) = 4
Für a(ψZ , ψO) gilt:
a(ψZ , ψO) =
∫I−V III
∇ψZ∇ψOdxdy
=
∫I+IV
∇ψZ∇ψOdxdy =
∫I+IV
∂1ψZ∂1ψO + ∂2ψZ∂2ψOdxdy
=
∫I+IV
∂2ψZ∂2ψOdxdy =
∫I+IV
−1
h· 1
hdxdy
= −1
h·∫I+IV
dxdy = −1
Aus Symmetriegründen folgt:
a(ψZ , ψO) = a(ψZ , ψU ) = a(ψZ , ψL) = a(ψZ , ψR) = −1
Durch Nachrechnen erhält man zudem:
a(ψZ , ψRU ) = a(ψZ , ψLO) = 0
Damit erhalten wir einen Matrix-Stern der Form 0 −1 0−1 4 −1
0 −1 0
Bemerkung 24. Die Identität der Sterne und damit der Systemmatrizen zwischen FiniteDifferenzen und Finite Elemente gilt im Allgemeinen nicht.
Das Beispiel von Courant zeigt uns:
1. Wir benötigen eine Triangulierung mit bestimmten Eigenschaften
2. Wir benötigen Ansatzfunktionen über dem diskreten Gittergebiet
69
Abbildung 4.3: Triangulierung eines Gebiets
4.4.2 Triangulierung
Ein Gebiet mit gekrümmten Rändern kann lokal linear approximiert werden.
Definition 46. (Zulässige Triangulierung)Eine Zerlegung T = T1, T2, . . . , TM von Ω in Dreiecks- bzw. Viereckselemente heißtzulässig, wenn folgende Eigenschaften erfüllt sind:
1. Ω =⋃Mi=1 Ti
2. Ti ∩Tj besteht aus genau einem Punkt, dann ist dieser ein Eckpunkt von Ti und Tj
3. Ti ∩ Tj mehr als ein Punkt, dann ist Ti ∩ Tj eine Kante von Ti und Tj
Die Punkte 2 und 3 definieren konforme Gitter.
A B
Abbildung 4.4: A: konformes Gitter, B: nicht konformes Gitter
In drei Raumdimensionen kommen verschiedene Gitterelemente zum Einsatz (siehe Abb.4.4.2).
Definition 47. (Finite Elemente Raum)Für ein Gitter Ωh über Ω ⊂ Rd ist
V ph (T ) = u ∈ H1 : für alle T ∈ T : u|T ∈ Pp
der konforme Finite Elemente Raum der Ordnung p.
Aus dem Beispiel von Courant können wir ein allgemeines Verfahren im Rd, d = 1, 2, 3,ableiten.
70
A DCB
Abbildung 4.5: A: Tetraeder, B: Quader, C: Prisma, D: Pyramide
4.4.3 Finite Elemente im R1
Wir betrachten als Modellproblem die Helmholtz-Gleichung
−∆u+ u = f mit a(u, v) =
∫Ω
(∇u∇v + uv)dx (4.39)
Sei N = a = x0, x1, x2, . . . , xN+1 = b die Diskretisierung des Gebiets [a, b] mit Gitter-weiten hi = xi+1 − xi und lokalen Ansatzfunktionen ϕi für die gilt
ϕi(xj) = δij
Damit hat die Lösung uh die Gestalt
uh =
N∑i=1
aiϕi.
Berechnung von uh auf Referenzintervallen
Die Ansatzfunktionen ϕi können durch Formfunktionen Φi dargestellt werden.Durch eine bijektive affine Transformation wird das Intervall [xi, xi+1] auf ein Referenz-intervall [0, 1] abgebildet. Die Lösung ui kann auf diese Weise immer auf dem Intervall[0, 1] über die Formfunktionen gelöst und dann zurücktransformiert werden. Dies hattechnische Vorteile gegenüber der Berechnung auf dem Originalelement.
Transformationsabbildung
Sei Ii = [xi, xi+1] und ξ ∈ [0, 1]. Dann definieren
xIi : [0, 1] −→ Ii
ξ 7−→ xi + hiξ
ξIi : Ii −→ [0, 1]
x 7−→ (x− xi)hi
eine Bijektion zwischen [0, 1] und [xi, xi+1]. Auf dem Referenzintervall können wir unsereLösung darstellen als
uh(ξ) = α1 + α2ξ.
71
1
1
1
x i x i+1x i ! 1
A B
Abbildung 4.6: A: Formfunktionen über Referenzintervall, B: Ansatzfunktionen überGitter
Dabei gilt ui = uh(0) = α1 und ui+1 = uh(1) = α1 + α2.
⇒ uh(ξ) = α1 + α2ξ = ui + (ui+1 − ui)ξ= (1− ξ)ui + ξui+1 = uiΦ1(ξ) + ui+1Φ2(ξ)
Bemerkung 25. Für die Formfunktionen gilt:
∀ξ ∈ [0, 1] : Φ1(ξ) + Φ2(ξ) = 1
Für die Ansatzfunktionen gilt
ϕi(x) =
Φ2(ξ(x)) x ∈ Ii−1
Φ1(ξ(x)) x ∈ Ii0 sonst
Wir können nun auf dem Referenzintervall und abhängig von den Formfunktionen dieMatrixeinträge der Systemmatrix berechnen.
Berechnung der Systemmatrix-Einträge
Zu berechnen sind die Einträge a(ϕi, ϕj) der Systemmatrix A (in unserem Fall für dasHelmholtz-Problem)
a(ϕi, ϕj) =N∑k=1
∫Ik
∇ϕi(x)∇ϕj(x) + ϕi(x)ϕj(x) dx.
72
Die ϕi, ϕj können durch Φn und Φm (n,m ∈ 1, 2) ersetzt werden.
⇒ (AIk)nm =
∫Ik
∇xΦn(ξ(x))∇xΦm(ξ(x)) + Φn(ξ(x))Φm(ξ(x))dx
= (xk+1 − xk)∫ 1
0∇ξΦn(ξ)ξ′(x(ξ))ξ′(x(ξ))∇ξΦm(ξ) + ΦnΦmdξ
= hk
∫ 1
0
1
h2k
∇Φn∇Φm + ΦnΦmdξ
Man sieht, die Systemmatrix ist zusammengesetzt aus 2× 2 Submatrizen AIk .Die Matrixeinträge von AIk lauten:
(AIk)11 =
∫ 1
0
1
hk∇Φ1∇Φ1 + Φ1Φ1 · hkdξ
=
∫ 1
0
1
hk+ hk(1− ξ)2dξ =
1
hk+
1
3hk
(AIk)12 =
∫ 1
0
1
hk∇Φ1∇Φ2 + Φ1Φ2 · hkdξ
=
∫ 1
0− 1
hk+ ξ(1− ξ) · hkdξ
= − 1
hk+
1
6hk
(AIk)21 = − 1
hk+
1
6hk
(AIk)22 =1
hk+
1
3hk
⇒ AIk = 1hk
[1 −1−1 1
]+ hk
[1/3 1/61/6 1/3
]Quadratische Elemente
Statt einer linearen Approximation können wir auch einen quadratischen Ansatz wählen.Dies erhöht die lokale Approximationsgüte. Wir stellen die Lösung uh(ξ) deshalb dar als
uh(ξ) = α1 + α2ξ + α3ξ2 auf I = [0, 1].
Um dieses Polynom eindeutig zu bestimmen, benötigen wir 3 Stützstellen, zusätzlich zuui und ui+1 können wir z.B.
ui+ 12
= uh
(xi + xi+1
2
).
73
An den Stützstellen gilt:
ui = uh(0) = α1
ui+1 = uh(1) = α1 + α2 + α3
ui+ 12
= uh
(1
2
)= α1 +
1
2α2 +
1
4α3
Berechnung der αi, i = 1, . . . , 3 liefert
uh(ξ) = uiΦ1(ξ) + ui+1Φ2(ξ) + ui+ 12Φ3(ξ)
mit
Φi(ξ) = 2
(ξ − 1
2
)(ξ − 1)
Φ2(ξ) = 2ξ
(ξ − 1
2
)Φ3(ξ) = 4ξ (1− ξ )
1
1
!3(!)
!2(!)!1(!)
Abbildung 4.7: Ansatzfunktionen für quadratische Elemente
Die Submatrizen AIk sind jetzt 3× 3-Matrizen und lassen sich analog zum linearen Fallausrechnen.
Bemerkung 26. Die Gradienten ∇Φi(ξ) sind nun keine Konstanten über dem Refe-renzintervall. Man benötigt also noch ein numerischen Verfahren zur Approximation derIntegrale.
4.4.4 Finite Elemente im R2
Wir übertragen nun die Vorgehensweise der Finiten Elemente im R1 auf den zweidi-mensionalen Fall. Gitterelemente, die im Eindimensionalen Intervalle waren, sind nunDreiecke oder Vierecke.
74
y
x
1
1
Abbildung 4.8: Affine Transformation von Dreiecken zwischen diskretisiertem Raum undReferenzelement
Die affine Transformation auf ein Referenzdreieck (bzw. Referenzviereck) ist deshalb einebilinear, bijektive Abbildung mit der Eigenschaft (im Fall des Dreieckselements):
x = x1 + (x2 − x1)ξ + (x3 − x1)η
y = y1 + (y2 − y1)ξ + (y3 − y1)η
⇔(xy
)=
(x1
y1
)+
(x2 − x1 x3 − x1
y2 − y1 y3 − y1
)(ξη
)Die Rücktransformation hat damit die Gestalt:(
ξη
)=
(x2 − x1 x3 − x1
y2 − y1 y3 − y1
)−1
︸ ︷︷ ︸A−1
(x− x1
y − y1
)
=1
detA
(y3 − y1 x1 − x3
y1 − y2 x2 − x1
)(x− x1
y − y1
)mit detA = (x2 − x1)(y3 − y1)− (x3 − x1)(y2 − y1). Für die partiellen Ableitungen gilt:
ux = uξξx + uηηx
uy = uξηy + uηηy
und
ξx =y3 − y1
detA
ηx =y2 − y1
detA
ξy =x3 − x1
detA
ηy =x2 − x1
detA
75
Bemerkung 27. detA beschreibt dabei die Volumenänderung des Dreiecks unter dervorgegebenen Transformation:
dxdy = detAdξdη
Mit ∇Φi, i = 1, 2, 3 und ξx, ξy, ηx, ηy sind die Komponenten der Matrix AIk vollständigbestimmbar. Dazu ist ein numerisches Integrationsverfahren notwendig.
Wahl der Φi: Lineare Elemente
Wir stellen die Lösung dar als
uh(ξ, η) = α1 + α2ξ + α3η
uj := uh(Pj), j = 1, 2, 3
Pj sind die Eckpunkte des Referenzdreiecks. Es gilt
u1 = uh(0, 0) = α1
u2 = uh(1, 0) = α1 + α2
u3 = uh(0, 1) = α1 + α3
Eingesetzt liefert dies
uh(ξ, η) = u1 + (u2 − u1)ξ + (u3 − u1)η = (1− ξ − η)u1 + ξu2 + ηu3.
Mit Φ1 = 1− ξ − η, Φ2 = ξ, Φ3 = η folgt
uh(ξ, η) = u1Φ1 + u2Φ2 + u3Φ3.
mit
Φi(Pj) = δij i, j = 1, 2, 33∑i=1
Φi(ξ, η) = 1 ξ, η ∈ T
Wahl der Φi: Quadratische Elemente
Wir benötigen für den Fall quadratischer Ansatzfunktionen 3 weitere Stützstellen, da wirdie Lösung uh nun schreiben als
uh(ξ, η) = α1 + α2ξ + α3η + α4ξ2 + α5ξη + α6η
2
Die Formfunktionen dazu lauten
Φ1 = (1− ξ − η)(1− 2ξ + 2η)
Φ2 = ξ(2ξ − 1)
Φ3 = η(2η − 1)
Φ4 = 4ξ(1− ξ + η)
Φ5 = 4ξη
Φ6 = 4η(1− ξ − η)
76
Bemerkung 28. Im R3 geht man analog vor. Der lineare Fall definiert
uh(ξ, η, ζ) = α1 + α2ξ + α3η + α4ζ.
Das Tetraederelement im dreidimensionalen Fall liefert die nötigen 4 Knoten für dieBestimmung von α1, . . . , α4.
4.4.5 Konvergenzaussagen zu Finite Elemente Verfahren
Abschätzung des Energiefehlers
Die Energienorm war definiert als
‖u‖a :=√a(u, u)
Wir wissen aus der Minimaleigenschaft
a(u− uh, u− uh) = minvh∈Vh
a(u− vh, u− vh)
⇒ ‖u− uh‖a = minvh∈Vh
‖u− vh‖a
Satz 19. Sei Ω ⊂ Rd, d ≤ 3 konvex und u ∈ H2(Ω) die schwache Lösung der Poisson-Gleichung. Sei uh ∈ Vh die Finite Elemente Lösung von
a(uh, vh) = 〈f, vh〉, vh ∈ Vh ⊂ H10 (Ω)
mit linearen, konformen finiten Elementen. Dann gilt für den Diskretisierungsfehler
‖u− uh‖a ≤ c · h · ‖f‖L2(Ω), f ∈ L2(Ω) (4.40)
Bemerkung 29. Falls für alle f ∈ L2(Ω) gilt
‖u‖H2(Ω) ≤ c‖f‖L2(Ω)
dann heißt das Problem H2-regulär.
Abschätzung des L2-Fehlers
Satz 20. Sei Ω ⊂ Rd, d ≤ 3 und u ∈ H2(Ω) schwache H2-reguläre Lösung der Poisson-Gleichung. Dann gilt
‖u− uh‖L2(Ω) ≤ c · h‖u− uh‖a (4.41)
‖u− uh‖L2(Ω) ≤ c · h2‖f‖L2(Ω) (4.42)
77
5 Ausblick
5.1 Finite Volumen Verfahren
Die Idee der Finite Volumen Diskretisierung ist sehr ähnlich zu den Finite ElementeVerfahren. Das Finite Elemente Konzept wird jedoch nicht auf das bisherige Gitter Ωh
angewandt, sondern auf ein verschobenenes, duales Gitter. Dies hat zur Folge, dass dasVerfahren lokal Flusserhaltend ist. Wir erinnern uns an Konzentrationsänderungen ineinem beliebigen Volumenelement die, ohne Quellen oder Senken im Volumen, gleichdem Fluss über den Elementrand sind:∫
B
∂u
∂td~x =
∫B−∆ud~x =
∫∂BF · ~nds
Wählen wir ein solches Volumenelement über jedem Knoten in Ωh, z.B. durch Verbindenvon Kantenmittelpunkten und Flächenschwerpunkte, so erhält man sogenannte Voronoie-Elemente die zusammengesetzt ein konformes duales Gitter über dem Gebiet Ω erzeugen.
⇒ −∫
Ωdiv(∇u)dx = −
M∑i=1
∫Bi
div(∇u)dx =M∑i=1
∫Bi
fdx
Für festes i gilt:
−∫Bi
div(∇u)dx =
∫Bi
fdx
⇔ −∫∂Bi
∇u · ~ndx =
∫Bi
fdx
⇔ −∫∂Bi
∂u
∂~ndx =
∫Bi
fdx
Dies führt zu einem Gleichungssystem
KFVh uFVh = fFVh
mit (KFVh
)Bi
= −∫∂Bi
∂ϕFVBi
∂~nds
78
5.2 Lösen von linearen Gleichungssystemen
Aus allen besprochenen Diskretisierungsverfahren gewinnen wir ein, mitunter sehr großes,lineares Gleichungssystem
K · uh = b,
wobei wir am Ende uh als Lösung des Problems berechnen wollen. D.h. wir benötigendie Darstellung
uh = K−1b,
müssen also K invertieren. Das Problem in realen Fällen ist jedoch, dass eine exakteInvertierung vom Rechen- und somit Zeitaufwand nicht vertretbar ist. Ein Ansatz zurLösung des Gleichungssystems ist, die Inverse von A zu approximieren und ein iterativesVerfahren anzusetzen. Welche Vor- und Nachteile solche Verfahren haben und wie maneffizient große lineare Systeme löst, ist Thema einer Folgevorlesung.
79
![футболвысокогорец-нт.рф/images/doc/sport podgotovka futbol.pdf · F b g b k l _ j k l s _ ] j h n _ k k b h g Z e v g h b y K \ _ j ^ h [ e Z k l b F b g b k](https://static.unterlagen.site/doc/80x56/5ea06bb051820068b4767b45/f-imagesdocsport-podgotovka-f-b.jpg)
![KD d / K >/D/d KE^K>/ d &/E E / > ^d d D Ed^ &KZ d, z Z E ... › CLIENT_CL_Corp_3513EFE5_AF37_1D84_E3… · KD d / K >/D/d ^dZ d '/ Z WKZd î d Z ] } v Z ] ^ P ] } ( } } u ] } >](https://static.unterlagen.site/doc/80x56/60ba3e73c531c90c56690c14/kd-d-k-dd-kek-d-e-e-d-d-d-ed-kz-d-z-z-e-a.jpg)
