Numerische Methoden für grosse nichtlineare
Gleichungssysteme
SoSe 2013
Wolfgang Mackens
Institut für Mathematik, TUHH
1. April 2013
Zusammenfassung
Grundziel des Kurses ist, die Teilnehmer in die Lage zu versetzen, große nichtlineare Gleichungssysteme unter MATLAB effizient behandeln zu können. Dies beinhaltet -im Gegensatz zur reinen Anwendung von Programmen in der Vorlesung „Numerische Software“ - auch ein Verstehen der mathematischen Grundlagen der Algorithmen. Zugleich wird eine Mindestkenntnis von MATLAB angestrebt, wie sie etwa auf den ersten 70 Seiten der „Einführung in MATLAB“ von Peter Arbenz [Arb] vermittelt wird. Schließlich sollen auch Grundkenntnisse über das Auftreten von großen Syste-men vermittelt werden, weil sich wesentliche Eigenschaften der Systeme (und damit die Art ihrer numerischen Behandlung) aus ihrer Herkunft erklären.
Inhaltsverzeichnis
1 MATLAB 9
2 Einführung: Entstehung „großer Systeme“ 10
2.1 Entstehung durch Diskretisierung . . . 11
2.1.1 Galerkin-Approximationen . . . 12
2.1.2 Gewöhnliche Randwertaufgaben . . . 17
2.1.3 Partielle Randwertaufgaben . . . 26
2.1.4 Integralgleichungen . . . 33
2.2 Entstehen von Größe aus Nichtlinearität . . . 34
2.3 Größe durch Parameterabhängigkeit . . . 36
2.4 Größe durch Zeitentwicklung . . . 37
2.4.1 Raum-Diskretisierung . . . 37
3 Lösung linearer und nichtlinearer Gleichungssysteme kleiner bis
mode-rater Dimension 40
3.1 Vorbereitung: Banachs Fixpunktsatz . . . 40
3.1.1 Banach-Räume . . . 40
3.1.2 Der Fixpunktsatz . . . 41
3.1.3 Splitting-Verfahren . . . 43
3.1.4 Iteration und Neumannsche Reihen . . . 56
3.1.5 Störungslemmata . . . 57
3.1.6 Bestimmung einer Lipschitz-Konstante im nichtlinearen Fall . . . . 60
3.1.7 Der Satz von Ostrowski . . . 61
3.1.8 Drei Anwendungen des Satzes von Ostrowski: . . . 63
3.1.9 Wichtige Anwendungen des Banachschen Fixpunktsatzes: Picard-Lindelöf, Implizite Funktionen, Lokale Umkehrbarkeit . . . . 63
3.1.10 Andere Fixpunktsätze . . . 74
3.2 Verfahren vom Newton-Typ . . . 79
3.2.1 Lokale Konvergenz des Newton-Verfahrens . . . 81
3.2.2 Affin-Invarianz . . . 83
3.2.3 Globales Verhalten des Newtonverfahrens . . . 86
3.2.4 Globalisierungsstrategien: Dämpfung, Homotopie, Trust-Region . . 88
3.2.5 „Halbglobale Konvergenz“ der Newton-Iteration . . . 104
3.2.6 Noch einmal Dämpfung . . . 111
3.3 Behandlung parameterabhängiger Systeme moderater Größe . . . 114
3.3.1 Prädiktor-Korrektor-Astverfolgungsprobleme. . . 116
3.3.2 Differentialgleichungsberechnung von Lösungsästen . . . 121
3.3.3 Berechnung von Umkehrpunkten: . . . 125
3.3.4 Berandete lineare Systeme . . . 133
3.3.5 Block-System . . . 135
3.4 Newton-Verfahren für nicht quadratische Systeme . . . 136
3.4.1 Lösungsmannigfaltigkeiten . . . 136
3.4.2 Gauß-Newton-Iterationen . . . 137
3.5 Differenzen Newton-Verfahren. . . 144
3.5.1 Differenzen-Approximation von Ableitungen . . . 144
3.5.2 Wahl der Differenzenschrittweite . . . 145
3.5.3 Ergänzungen zu Ableitungsapproximationen der ersten Ableitung . 148 3.5.4 Strukturen in der Jacobimatrix . . . 152
4 Direkte Lösung „großer linearer Systeme“ 155
4.1 Bandreduktion bei SPD-Systemen . . . 160
4.2 Bandstruktur-Erzwingung . . . 170
4.3 Sherman-Morrison-Woodbury Ansätze . . . 173
5 Iterative Verfahren zur Lösung großer Linearer Systeme 175 5.1 Stationäre Verfahren: Konvergenzaussagen für Matrixklassen . . . 179
5.1.1 Diagonaldominanz und schwache Diagonaldominanz . . . 181
5.1.2 SPD-Systeme . . . 188
5.1.3 Konsistente Ordnung . . . 195
5.2 Nichtstationäre Verfahren 1: Konvergenzverbesserung . . . 196
5.2.1 Polynomiale Beschleunigung . . . 197
5.2.2 Unterraum-Minimierung . . . 202
5.3 Nichtstationäre Verfahren 2: CG-Typ-Verfahren . . . 203
5.3.1 Herleitung des Verfahrens als Ritz-Galerkin-Verfahren . . . 203
5.3.2 Skizze einer Fehleranalyse für den CG-Algorithmus . . . 205
5.4 Nichtstationäre Verfahren 3: CG-artige Verfahren für nicht-symmetrische Probleme . . . 209
5.5 Nichtstationäre Verfahren 4: Kazcmarz-Typ-Verfahren . . . 210
5.5.1 Zugang 1: Zeilenorientierte Lösung riesiger Systeme . . . 211
5.5.2 Zugang 2: Fehler-Minimierung durch Liniensuche . . . 212
5.5.3 Zugang 3: Allgemeine projektive Fehlerreduktion . . . 218
5.5.4 Konvergenz allgemeiner Projektionsverfahren . . . 218
5.5.5 Verallgemeinerungen . . . 223
5.6 Weitere Beispiele für Projektionsverfahren . . . 227
5.7 Substrukturen und Iterationsverfahren . . . 234
5.7.1 Die Schwarzsche Alternierende Iteration . . . 234
5.7.2 CG im Schur-Komplement . . . 241
5.8 Neuere Konvergenzanalysen für Kazcmarz-Typ Iterationen . . . 246
6 Direkte Iterationsvarianten für große nichtlineare Systeme 246 6.1 Linearisierungs„freie“ Verfahren . . . 247 6.1.1 Nichtlineares Jacobischrittverfahren . . . 247 6.2 Newton-Mix Verfahren . . . 249 6.2.1 XN-Verfahren . . . 249 6.2.2 NX-Verfahren . . . 251 6.2.3 Inexakte Newton-Methoden . . . 253 6.3 Update-Methoden . . . 256
7 Reduktionsansätze 257
7.1 Reduktion durch Teillösung . . . 257
7.1.1 Shooting-Verfahren . . . 257
7.1.2 Master-Slave-Einteilungen . . . 259
7.2 Verbindung von Diskretisierungen verschiedener Feinheiten . . . 259
7.2.1 Aufsteigende Gitter . . . 259
7.2.2 Extrapolationsverbesserung . . . 261
7.2.3 Mehrgitterverfahren . . . 264
7.3 Informations-Wiederverwendung . . . 264
7.3.1 Verwendung bestehender Löser . . . 264
7.3.2 Verwendung bestehender Lösungen . . . 265
7.3.3 Verwendung ähnlicher Probleme . . . 265
7.4 Zerlegungen . . . 265
7.4.1 Multiple Shooting . . . 265
7.4.2 Gebietszerlegungsverfahren . . . 265
7.4.3 Newton-Kopplung von Unterproblemen . . . 266
7.4.4 Weitere Block-Newton-Verfahren . . . 270 7.4.5 Master-Slave-Zerlegung . . . 270 7.4.6 ABS-Methoden . . . 271 7.4.7 Spektral-Zerlegungsmethoden . . . 272 7.4.8 Mehrgitter-Methoden . . . 272 8 Rekursive Projektion 272 8.1 CNSP: Condensed Newton/Supported Picard . . . 272
8.2 CNSP für parameterabhängige Gleichungen . . . 277
8.3 Recursive Projektion nach Shroff und Keller . . . 278
9 Methoden der Reduzierten Basen. 279 9.1 Generelle Idee . . . 279
9.2 Tangenten-orientierte Reduzierte Basen . . . 279
9.2.1 Grundlagen . . . 279
9.2.2 Zulässige Testmatrizen . . . 281
9.2.3 Berechnung geeigneter Basen . . . 283
9.3 POD-Ansätze . . . 285
12 Literatur 286
Abbildungsverzeichnis
1 Brücke bei Hochdonn mit Last . . . 10
2 Verformung durch eine inkrementell aufgebrachte Last . . . 11
3 Diskrete „Lösung“ einer Randwertaufgabe . . . 18
4 Stückweise lineare Funktion . . . 23
5 Finite Volumen . . . 25
6 Gitterfunktion . . . 28
7 Numerierung . . . 29
8 Lösungsapproximation zu Aufgabe (48) . . . 30
9 Finite-Element Gitter und Basisfunktion . . . 31
10 Finite Volumen Diskretisierung . . . 32
11 Zwei Lösungen des Bratuproblems (11). . . 35
12 Kontinuum von Bratu-Lösungen . . . 36
13 Schrittweitenvergleich impliziter vers. expliziter Löser . . . 39
14 ρ abhängig von ω für SOR und SSOR . . . 56
15 Zwei Lösungsäste . . . 68
16 Auflösbarkeit in verschiedenen Punkten . . . 69
17 Abbildungsgrade −1, 0 und 1 . . . . 74
18 Veranschaulichung des Satzes von Sard . . . 75
19 Grad für nichtreguläre y-Werte durch nahe gelegene y-Werte . . . . 76
20 Grad für nichtdifferenzierbare Funktion durch glatte Approximation . . . . 76
21 Homotope Verbindung mit einer Gerade . . . 77
22 Newton-Iteration für f :R 7−→ R . . . . 80
23 Gleichungssystem „Zwei Kreise“ . . . 83
24 Äquivalent zu „Zwei Kreise“ . . . 84
25 Kreisbüschel . . . 84
26 Konvergenzprobleme bei Newton . . . 87
27 Julia-Menge der Newton-Iteration für f (x) = x3− 1. . . . 87
28 Zu großer Newton-Schritt . . . 89
29 Armijo klappt . . . 90
30 Armijo steuert Nullstelle in Zoom-Viereck an . . . 90
31 Zoom: Lokale Konvergenz . . . 91
32 L(f, x0), Z(f, x0) und obere Schranke für Testfunktion f (x)2 . . . 91
34 Fortsetzungsmethode . . . 94
35 Umkehrpunkt . . . 96
36 Vergleich von Newtonfluss und Davidenko . . . 97
37 Disjunkte Branin-Pfade . . . 99
38 Trust-Region-Schritt . . . 100
39 Dogleg-Schritt, Definition . . . 103
40 Dogleg-Version für Gleichungssystem (165) . . . 103
41 Newton-Fluss liefert keine Nullstelle auf unbeschränkter Level-Menge . . . 105
42 Newton-Fluss liefert keine Nullstelle bei offenem Rand der Levelmenge . . 106
43 Zusammenhangskomponenten der Levelmenge . . . 106
44 Zusammenhangskomponenten der Levelmenge . . . 107
45 Newton-Fluss liefert keine Nullstelle . . . 108
46 Kontinuierliches Newton-Verfahren für z3− 1 = 0 . . . 110 47 Rosenbrock-Level . . . 111 48 Rosennbrock-Level . . . 112 49 Natürliche Testfunktion . . . 113 50 Treppe: Escher . . . 113 51 Verzweigungsdiagramm . . . 115
52 Reine und gestörte „Pitchfork-Verzweigung“ . . . 116
53 Prädiktor/Korrektor-Grundstruktur . . . 116
54 Festlegung der Tangentialrichtung . . . 118
55 Tangentialfeld . . . 118
56 Lösung der Linearisierung . . . 119
57 Festlegung der Tangentialrichtung . . . 120
58 Schritt senkrecht zur Tangente am aktuellen Ort . . . 120
59 Tangentialfeld . . . 121
60 Ganzer - sauberer - Kreis . . . 123
61 Low quality integration . . . 123
62 Ast kehrt in (0,0) nicht um . . . 126
63 Einfacher Umkehrpunkt . . . 127
64 Umkehr- und Verzweigungspunkt . . . 128
65 Mannigfaltigkeiten singulärer Jacobischer . . . 128
66 L-Umkehrpunkte . . . 131
67 Newton-Fluss in (1,1)-Umkehrpunkte . . . 132
68 Fehlereinfluss bei Kreisberechnung . . . 138
71 Fluss zum Ausgleichskreis . . . 143
72 Beobachteter Fehler bei einseitiger Approximation der ersten Ableitung . . 145
73 Theoretischer Abbruchfehler und beobachteter Fehler . . . 146
74 Abbruch-, Rundungs- und Gesamtfehler der einsitigen Differenzenapproxi-mation . . . 147
75 Feiner . . . 148
76 Feinstruktur Blow up . . . 148
77 Zusätzlicher Effekt . . . 149
78 Verschiedene Approximationsordnungen . . . 151
79 Approximation der vierten Ableitung . . . 152
80 Spaltengestörte Diagonalmatrix . . . 154
81 Tridiagonalmatrix . . . 154
82 Adjazenzgraph zur Matrix B . . . 161
83 Diskretisierungsstruktur . . . 163
84 Erste Numerierung . . . 163
85 Matrix zur ersten Numerierung . . . 164
86 Zweite Numerierung . . . 164
87 Matrix zur zweiten Numerierung . . . 164
88 Zwei CM-Numerierungen . . . 165
89 Matrizen zu den beiden Numerierungen . . . 165
90 Dritte Numerierung . . . 167
91 Matrix zur dritten Numerierung . . . 168
92 Matrix zur dritten Numerierung,rückwärts . . . 168
93 CM-Matrix mit Fill-In . . . 169
94 RCM-Matrix mit Fill-In . . . 170
95 Zyklischer Graph zu periodischen Randwerten . . . 170
96 Gerichtete Graphen . . . 183
97 Konvergenzrate und Kondition . . . 189
98 MPE bei glattem und rauem Fehlervektor . . . 198
99 MPE mit Fehleraufschaukelung . . . 199
100 L-förmiges Gebiet . . . 234
101 Rechteckige Teilgebiete . . . 235
102 Schwarzsche alternierende Iteration für 1D-Beispiel . . . 237
103 Schwarzsche alternierende Iteration: Konvergenz . . . 237
104 Ungeeignete Aufteilung . . . 238
105 Drei Gebiete . . . 239
107 Schwarzsches Problem . . . 240
108 Ω3 als „Interface“ . . . 241
109 Substruktur-Numerierung . . . 242
110 Graph von Φ(s) . . . 258
111 Differenzen-Lösungsaproximationen für das Bratu-Problem . . . 260
112 Die Approximationen bei x = 0.5 . . . 261
113 Werte bei x = 0.5 . . . 261
114 Gauss-Seidel-Newton . . . 267
115 Tangential-Block-Newton . . . 268
116 Hierarchische Blockstruktur . . . 271
117 CNSP-Grundstruktur . . . 276
1
MATLAB
Zur Gewährleistung eines qualifizierten Umganges mit MATLAB arbeiten die Studierenden parallel zur Vorlesung und den Übungen (mindestens) die ersten 70 Seiten von [Arb] durch, und zwar nach dem folgenden Plan:
1. Vorlesungswoche: Kapitel 1 und 2 2. und 3. Vorlesungswoche: Kapitel 3 4. Vorlesungswoche: Kapitel 4 und 5 5. Vorlesungswoche: Kapitel 6 und 7
Zusätzlich zu den Übungsaufgaben mit mathematischen Inhalten werden in den ersten Wochen verstärkt Aufgaben zum Umgang mit MATLAB gestellt.
2
Einführung: Entstehung „großer Systeme“
Zunächst einmal wollen wir kurz erklären, dass wir uns hier nicht so sehr für Gleichungssy-steme interessieren, die einfach nur furchtbar viele Unbekannte und Gleichungen haben1. Große (nichtlineare) Gleichungssysteme entstehen in den Anwendungen meistens durch ganz bestimmte Prozesse2. Die Kenntnis dieser Prozesse ist hilfreich, denn durch den je-weiligen Prozess erhalten die Systeme bestimmte Strukturen, die sich nicht nur bei der Aufstellung der Gleichungen, sondern auch im Lösungsprozess oft arbeitsvermindernd aus-nutzen lassen.
Eine große Klasse großer Systeme entsteht z.B. durch sogenannte Diskretisierung von Pro-blemen mit unendlich vielen3 Freiheitsgraden. Dabei bedeutet „Diskretisierung“ die Appro-ximation der unendlichdimensionalen Probleme durch endlichdimensionale Näherungen, vgl. Abschnitt 2.1. Man kann sich gut vorstellen, dass solche Approximationen in der Re-gel umso besser werden, je mehr Freiheitsgrade in den Diskretisierungen verwendet werden. Damit ist dann auch klar, dass man beliebig große Probleme erzeugen kann.
Andererseits kann ein Problem auch dadurch „größer“ werden, dass es zwar relativ wenige Variablen enthält, dass dieses Problem aber für mehrere, meist sehr viele „Fälle“ gerechnet werden muss.
Solche Fälle können durch verschiedene „Lastfälle“ gegeben sein, bei der z.B. ein und das-selbe mechanische System4 durch mehrere verschiedene Sätze einwirkender Kräfte belastet wird5.
Abbildung 1: Brücke bei Hochdonn mit Last
Solche „Fälle“ können andererseits aber auch dadurch erzeugt sein, dass das zu behandelnde System von einem oder gleich mehreren Parametern abhängt6, und dass man das System in einer Parameterstudie untersuchen will. Dabei vergleicht man entweder die Lösung für verschiedene Parameter, oder man will die Parameter so anpassen, dass bestimmte weitere Bedingungen erfüllt werden. Vgl. hierzu Abschnitt 2.3.
1Zwar werden wir auch zu diesem allgemeinen Fall Aussagen machen, aber die Menge solch allgemeiner
Aussagen ist ziemlich schnell erschöpft, wenn die Systeme nicht noch zusätzliche Struktur aufweisen.
2Das ist eigentlich einleuchtend, weil kein Mensch Gleichungssysteme mit um eine Million Gleichungen
und Unbekannten eingeben könnte. Wenn ihnen aber ein bestimmter Bildungsprozess zugrunde liegt, so wird man diesen auch nutzen, um das System in den Rechner einzubringen.
3i.a. sogar überabzählbar vielen
4Etwa eine Brücke, siehe Abbildung 1.
Abbildung 2: Verformung durch eine inkrementell aufgebrachte Last
Weiterhin kann ein sehr oft wiederholtes Lösen von Gleichungssystemen dadurch nötig werden, dass man beim Lösen von Evolutionsgleichungen implizite Löser verwendet, die pro Zeitschritt die Lösung eines Systems erfordern, vgl. Abschnitt 2.4.
2.1
Entstehung durch Diskretisierung
Den meisten großen Systemen der Numerik liegen Gleichungen in Funktionenräumen zu-grunde, die zunächst einmal unendlich-dimensional sind. Dabei wird analog zum linearen Gleichungssystem
Ay = f
eine gesuchte Lösung y durch eine lineare oder nichtlineare allgemeinere Abbildung A auf einen Bildpunkt f abgebildet.
Während im endlichdimensionalen Fall A eine Matrix ist und y und f Vektoren sind (meist gleicher Dimension), sind y und f hier Funktionen. Die Abbildung A wird dann meist „Operator“ genannt, der aus Funktionen andere Funktionen macht.
A kann z.B. ein gewöhnlicher Differentialoperator sein, wie zum Beispiel
(Ay)(x) := y(n)(x) +
n−1
∑
k=0
ak(x)· y(k−1)(x), x∈ [0, 1]
wobei A Funktionen aus Cn[0, 1] mit vorgegebenen Randbedingungen auf Funktionen
aus C0[0, 1] abbildet7.
A kann bei Funktionen y = y(x1, . . . , xn) ein partieller Differentialoperator sein,
wie zum Beispiel der Laplace-Operator (Ay)(x) := n ∑ k=1 ∂2y ∂xk2 (x) = f (x), x ∈ [0, 1]n= Ω,
der die auf Ω zweimal stetig differenzierbaren Funktionen mit Null-Randwerten C2 0(Ω) in die auf Ω stetigen Funktionen abbildet.
7Die Koeffizientenfunktionen a
A kann ebenso ein Integraloperator sein: (Ay)(x) :=
∫ 1 0
A(x, s)y(s)ds
der (bei z.B. stetigem A(x, s)) stetige Funktionen auf stetige Funktionen abbildet, und es kann A - weitaus komplizierter - aus mehreren solchen Abbildungen
zusammen-gesetzt sein.
Es kann A aber auch einfach selbst eine Matrix sein, die aber so groß ist, dass man das zugehörige Gleichungssystem nicht durch die üblichen Methoden lösen kann. Wir wollen all diese Fälle hier kurz beleuchten. Für eine genauere Analyse der Gleichungen und der Approximationsgüte müssen wir auf entsprechende Vorlesungen verweisen. Da die Galerkin-Approximation in dieser Vorlesung von großer Bedeutung ist, wollen wir sie für den technisch einfachen Fall sehr großer Gleichungssysteme erklären.
2.1.1 Galerkin-Approximationen
Bei sogenannten Galerkin8-Diskretisierungen von Operatorgleichungen
Ay = f, A : U −→ V (1)
mit einem Operator, der den Vektorraum U in den Vektorraum V abbildet, wird eine Lösungsapproximation ˜y in einem endlichdimensionalen Unterraum ˜U von U gesucht. Mit einer Basis Φ1, . . . , Φn von ˜U macht man dazu den Ansatz
˜ y = n ∑ k=1 αkΦk.
Setzt man ˜y für y in (1) ein, so gelangt man im Falle der Linearität von A zur Aufgabe A˜y =
n
∑
k=1
αk(AΦk) = f,
Hier muss offenbar das Element f ∈ V aus den Bildern AΦ1, . . . , AΦk kombiniert werden.
Das wird selten möglich sein9.
Ist dies nicht möglich, so wird die Wunsch-Gleichung
(AΦ1)α1+ (AΦ2)α2+· · · + (AΦn)αn= f (2)
nur näherungsweise gelöst werden können.
8Sprich „Galjorkin“; im Russischen wird der Familienname des Erfinders des hier zu besprechenden
Verfahrens, „Boris Grigorjewitsch Galjorkin“, mit einem „e“ mit Doppelpunkten geschrieben, „ë“, das „ jo“
gespochen wird: ΓaΛëpκnH.
9Ein Beispiel, in dem so etwas klappen kann, ist die Gewinnung spezieller Lösungen von
Differential-gleichungen mit konstanten Koeffizienten (vgl. Mathematik III). Hier kann man häufiger die rechte Seite f
als zugehörig zu einem invarianten Unterraum von A identifizieren. In diesem Fall wird man die Φ1, . . . , Φn
als Basis dieses Unterraumes wählen.
Ist auf V ein inneres Produkt ⟨·, ·⟩ :
{
V × V −→ R
(f1, f2) 7−→ ⟨f1, f2⟩
gegeben, so bietet sich an, ˜y als Projektion von f auf span{AΦ1, . . . , AΦn} zu realisieren.
Aus der linearen Algebra ist bekannt, dass diese Projektion ˜y eindeutig bestimmt ist und dass jeder Koeffizienten-Vektor α = (α1, . . . , αn)T von Koeffizienten α1, . . . , αn für eine
Darstellung ˜y =∑nk=1αk(AΦk) das Gramsche System
Gα = F (3)
erfüllt, wobei
Fk =⟨AΦk, f⟩, k = 1, . . . , n
und
Gi,j =⟨AΦi, AΦj⟩, i, j = 1, . . . , n
gelten10.
Im Fall, dass A einfach nur eine große Matrix aus R(N,N ) ist, die von U =RN nach V = U abbildet, die Φk-Ansatzvektoren in eine Matrix
Φ := (Φ1, . . . , Φn)∈ R(N,n)
eingeordnet sind und das innere Produkt das normale Euklidische Innere Produkt ist, bekommt das Gramsche System (3) als Normalgleichung für das Ausgleichsproblem
⟨AΦα − f, AΦα − f⟩ = min
α (4)
die Form (
ΦTATAΦ)α = ΦTATf. (5)
Wie in den Normalgleichungen
ATAx = ATb zum Ausgleichsproblem
∥Ax − b∥2 = min führt ATA zu einer Quadratur der Kondition von A.
Während im Ausgleichsproblem diese Konditionsverschlechterung durch die Verwendung von Orthogonalisierungsverfahren vermieden wird, ist hier eine einfache Idee zur Vermei-dung eines zweiten As, einfach ein AT auf beiden Seiten fortzulassen.
Man gelangt dann zum System (
ΦTAΦ)α = ΦTf, (6)
welches entsteht, wenn man die Gleichung (2) mit Φ1, . . . , Φn innerlich multipliziert.
Im allgemeinen Fall A : U −→ V ist dabei natürlich Vorbedingung, dass die Elemente Φ1, . . . , Φn∈ U auch zu V gehören.
Die so entstehende Diskretisierung (6) heißt Ritz-Galërkin-Diskretisierung.
10Dabei ist G bekanntermaßen regulär und die Lösung α eindeutig, wenn die Elemente AΦ
1, . . . , AΦn
li-near unabhängig sind. Andernfalls ist die Approximation ˜y wohl existent und eindeutig, die zur Darstellung
Schon im Rn-Fall sieht man leicht, dass dies Verfahren nicht unbedingt zu einer Lösung führen muss:
Will man nämlich z.B. die Lösung des Gleichungssystems ( 1 0 0 −1 ) x = ( 1 1 )
im reduzierten Raum span {(
1 1
)}
lösen, so wird das Gleichungssystem (2.1.1) zu 0· α1 = 2,
was natürlich nicht hilft.
Problemlos wird das Verfahren (im Endlichdimensionalen) dann, wenn A symmetrisch und positiv definit ist. Dann ist die Matrix ΦTAΦ ebenfalls symmetrisch und positiv definit. Tatsächlich hat das reduzierte System dann auch noch gleich eine andere schöne und hilfreiche Interpretation. Wir schreiben diese im nächsten Unterparagraphen sowohl für den endlichdimensionalen Fall auf als auch für den Fall der approximativen Lösung gewisser gewöhnlicher Randwertaufgaben. Dazu starten wir ganz neu.
a. Ritz-Galërkin-Approximation In der Mathematik III haben Sie gelernt, dass man eine Minimierungsaufgabe
f (x, y) = min, f ∈ C2(R2,R) dadurch löst, dass man den Gradienten von f gleich Null setzt
∇f(x, y) = 0, (7)
und an den Lösungen dieses 2× 2-Gleichungssystems prüft, ob die Hessesche von f positiv definit ist, um zu sichern, dass es sich bei der Lösung um ein Minimum handelt.
In Vorlesungen zur numerischen Mathematik haben Sie vielleicht schon erfahren, dass das Lösen des Gleichungssystems (7) im Ernstfall durchaus keine Trivialität sein muss, und das man über die verschiedenen Verfahren hierfür ohne weiteres ein bis zwei Semester sprechen kann.
Und Vorlesungen über Optimierung haben Ihnen möglicherweise schon vermittelt, dass es keineswegs klar ist, dass man bei der Minimierung von f überhaupt den Weg über das Nullstellenproblem (7) gehen muss. Man kann auch andere Verfahren nutzen, die - ausge-hend von einer Start-Näherung für die Lösung - konsekutiv bessere Näherungen erzeugen, wobei in jedem Schritt der Funktionswert geeignet verkleinert wird.
In Vorlesungen über die numerische Lösung großer linearer Gleichungssysteme haben Sie unter Umständen schließlich gehört, dass die modernen Methoden zur Gewinnung einer Lösung u∈ Rn eines Gleichungssystemes
Au = b, (8)
im Falle einer symmetrischen und positiv definiten Systemmatrix A ∈ R(n,n) den Spieß sogar umkehren, indem Sie statt (8) das Minimierungsproblem
F (u) := 1 2u
Dabei gehen verschiedene Verfahren fast immer so vor, dass sukzessive Lösungsapproxima-tionen F (u) über wachsende Unterräume11 minimieren. Wie sich herausstellt, erfüllen die Lösungen dabei Gleichungen der Form (6).
Wenn u eine Lösung von (9) ist, so muss die Richtungableitung von F bei u für alle Richtungen v ∈ Rn verschwinden d dtF (u + tv)|t=0 = 0. Wegen F (u + tv) = 1 2(u + tv) T A(u + tv)− bT(u + tv) = 1 2u T Au + t uTAv + t21 2v T Av− bTu− tbTv heißt dies, dass
uTAv = bTv für alle v ∈ Rn. (10)
Wenn v in dieser Gleichung den ganzenRndurchläuft, folgt hieraus natürlich sofort wieder, dass u das System (8) lösen muss.
Interessant wird’s, wenn wir F einmal nicht über ganz Rn sondern nur über einen Unter-raum minimieren. Für diese Betrachtung wird die Bedingung (10) noch nützlich werden.
Fehlerabschätzung für eine Ritz-Lösung Die Idee, Minima quadratischer Funktio-nale12 durch Minima dieser Funktionale über niedrigdimensionaler Teilräume zu approxi-mieren, wird allgemein auf Walter Ritz13 zurückgeführt.
Ist Rh ein Teilraum14 des Gesamtraumes Rn und
uh = arg minu∈Rh ( 1 2u TAu− bTu ) , so ist völlig analog zu (10) die schwache Gleichung
uThAvh = bTvh für alle vh ∈h (11)
zu gewinnen.
Wegen Rh ⊂ Rn können in (10) natürlich auch die v-Elemente aus Rh eingesetzt werden,
so dass das Minimum u im Gesamtraum ebenfalls
uTAvh = bTvh für alle vh ∈ Rh (12)
erfüllt. Subtraktion der Gleichung (11) von der Gleichung (12) führt auf die wichtige Glei-chung
(u− uh)TAvh = 0 für alle vh ∈ Rh. (13)
Da A eine SPD-Matrix ist, ist durch
⟨v, w⟩A:= vTAw
11Wenn man nicht bei Null startet, über wachsende affine Räume
12Im vorliegenden endlichdimensionalen Fall das Funktional (8), in späteren Anwendungen aber auch
allgemeinere Funktionale, vgl. (33) unten.
13http://de.wikipedia.org/wiki/Walter_Ritz
14Das tiefgestellte h soll hier schon darauf hinweisen, dass Teilräume von Funktionenräumen, wie wir sie
später zur Minimierung von (33) verwenden werden, stets von einem Längenmaß h abhängen werden, welches - wie die Differenzenschrittweite bei den Differenzenverfahren - durch Verkleinerung zu einer Verbesserung der Approximation führt.
ein inneres Produkt definiert und durch ∥v∥A=
√ ⟨v, v⟩A
eine zugehörige Norm.
Die Gleichung (13) besagt mit dieser Interpretation, dass uh die beste Approximation von
u aus Rh bezüglich der durch
da(v, w) :=∥v − w∥A
definierten Metrik ist:
∥u − uh∥A≤ ∥u − vh∥A für alle vh ∈ Rh. (14)
Hiermit wird nun vielleicht klar, weshalb in der Linearen Algebra I gleich ganz allgemeine innere Produkte und nicht nur das euklidische Produkt behandelt wurden. Mit der dem A-inneren Produkt zugeordneten Norm berechnet die Methode die beste Approximation an die Lösung aus dem gewählten Unterraum.
Während die Verallgemeinerung der Sichtweise eine schöne und einfache Interpretation der Ritz-Approximation ergibt, hätte man aber natürlich trotzdem gern gewusst, wie gut denn diese Approximation in normalen Maßstäben ist, etwa in der euklidischen Norm, ∥ · ∥2. Hierfür nutzen wir die folgenden Abschätzungen
α∥v∥22 ≤ ∥v∥2A = vTAv, für alle v ∈ Rn (15)
und
|⟨v, w⟩A| = |vTAw| ≤ ∥A∥2∥v∥2∥w∥2 = C∥v∥2∥w∥2, (16) wobei α = λmin(A) der kleinste Eigenwert von A und C = ∥A∥2 = λmax(A) der größte Eigenwert von A ist.
Mit (15) und (16) sowie (14) können wir jetzt für alle vh ∈ Rh abschätzen.
α∥u − uh∥22 ≤ ∥u − uh∥2A ≤ ∥u − vh∥2A=|⟨u − vh, u− vh⟩A| ≤ C∥u − vh∥22 so dass ∥u − uh∥2 ≤ √ C α∥u − vh∥2 = √
cond(A)∥u − vh∥ für alle vh ∈ Rh (17)
ist.
Satz 2.1 (Fehlerabschätzung für die Ritz-Galerkin-Approximation)
Sei A∈ R(n,n) eine SPD-Matrix mit Spektralkondition cond(A) und seien für einen Unter-raum Rh ⊂ Rn definiert u = arg minu∈Rn ( 1 2u TAu− bTu )
und uh = arg minu∈Rh
( 1 2u TAu− bTu ) so gilt ∥u − uh∥2 ≤ √ C α∥u − vh∥2 für alle vh ∈ Rh (18)
Die Abschätzbarkeit (15) von ∥v∥A nach unten wird man im Rahmen der
Variationsme-thoden die „Koerzivität“ oder „Elliptizität“ der Bilinearform⟨u, v⟩A nennen, die
Abschätz-barkeit (16) nach oben seine „Stetigkeit“. Die Folgerung (17) aus beiden Eigenschaften wird später im allgemeinen Rahmen der partiellen Differentialgleichungen als Céa-Lemma
Anmerkungen 2.2 („Projektion eines Gleichungssystems“)
Wählt man in der Ritz-Galerkin-Approximation einer Lösung x∈ Rn von Ax = b aus dem linearen Unterraum Rh ⊂ R(n,n) eine Basis v1, . . . , vm von Rh, bildet mit diesen Vektoren
die Matrix V = (v1, . . . , vm) und stellt die Ritz-Lösung in der Form v
h = V y mit y ∈ Rm
dar, so werden die schwachen Gleichungen (11) alle zusammen erfasst durch
VT(AV y− b) = 0 (19)
oder (
VTAV)y = VTb. (20)
Die Matrix(VTAV)wird manchmal als „Projektion der Matrix A auf den Unterraum Rh“
bezeichnet15.
2.1.2 Gewöhnliche Randwertaufgaben A. Diskretisierung mit finiten Differenzen
Lineare Randwertaufgabe Wir erklären das Verfahren der Finiten Differenzen zu-nächst nur als besonders einfache Version für die lineare Differentialgleichung
−y′′(x) + p(x)y′(x) + q(x)y(x) = f (x) (21)
auf dem endlichen reellen Intervall (a, b) mit auf [a, b] stetigen Funktionen p, q und f sowie den einfachen Dirichletschen Randbedingungen
y(a) = A, y(b) = B. (22)
Wenn wir (der Einfachheit halber) zusätzlich annehmen, dass q(x)≥ 0, x ∈ [a, b], können wir davon ausgehen (vgl. [CL]), dass die Randwertaufgabe (21,22) eine eindeutige Lösung y∈ C2[a, b] besitzt16.
Bei der Methode der finiten Differenzen verzichtet man auf den Erwerb der Kenntnis von y(x) an allen Punkten des Intervalles [a, b] und bescheidet sich auf die Berechnung von Näherungen yi, i = 1, . . . , n von Funktionswerten y(xi) in Punkten xi eines Gitters
a = x0 < x1 < . . . < xn< xn+1 = b
in [a, b]. Üblicherweise wählt man die Gitterpunkte äquidistant xi+1− xi = h := (b− a)/(n + 1)
und nennt h die Schrittweite der Diskretisierung.
Um zu Gleichungen für die Näherungen yi ≈ y(xi) zu gelangen, schreibt man zunächst
einmal die Differentialgleichung (21) in den „inneren Punkten“ x1, . . . , xn des Gitters hin
−y′′(x
i) + p(x)y′(xi) + q(x)y(xi) = f (xi), i = 1, . . . , n. (23)
Mit Hilfe der per Taylorentwicklung für C4-Funktionen leicht zu bestätigenden Differen-zenapproximationsaussagen
y(x + h)− y(x − h)
2h = y
′(x) + O(h2), y(x− h) − 2y(x) + y(x + h)
h2 = y
′′(x) + O(h2)
15Obwohl man (n, n)-Matrizen natürlich überhaupt nicht auf Teilräume desRn projizieren kann.
16Es ist nützlich, vor dem Versuch einer numerischen Approximation zu wissen, dass eine Lösung
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
Exakte Lösung von −y’’+y=1, y(−1)=y(1)=0 und FD−Approximation
y 5 y(x
5)
Fehler wurde überzeichnet
Abbildung 3: Exakte Lösung (durchgezogen) + FD-Approximation (∗) ersetzt man nun die Ableitungen in (23) durch Differenzen mit dem Ergebnis17
−y(xi−1)+2y(xi)−y(xi+1)
h2 + p(xi)
−y(xi−1)+y(xi+1)
2h + q(xi)y(xi) = f (xi) + O(h 2), i = 1, . . . , n.
Die Approximationen yi ≈ y(xi) werden nun dadurch definiert, dass man von ihnen
ver-langt, dasjenige Gleichungssystem zu lösen, was sich ergibt, wenn man die (unbekannten) Differenzen-Fehler der Größe O(h2) einfach fortlässt.
−yi−1+ 2yi− yi+1
h2 + p(xi)
−yi−1+ yi+1
2h + q(xi)yi = f (xi); i = 1, . . . , n. (24) Setzt man hierin die Randwerte y0 = y(x0) = A und yn+1 = y(xn+1) = B ein, so sieht
man, dass es sich um ein lineares Gleichungssystem
Ahyh = fh+ Rh (25)
für den Vektor yh := (y1, . . . , yn)T handelt mit der Tridiagonalmatrix
Ah = 1 h2 2 + h2q(x 1) 12hp(x1)− 1 0 . . . 0 −1 2hp(x2)− 1 2 + h 2q(x 2) 12hp(x2)− 1 . .. ... 0 . .. . .. . .. 0 .. . . .. . .. . .. 12hp(xn−1)− 1 0 . . . 0 −12hp(xn)− 1 2 + h2q(xn) (26) dem Vektor der „diskretisierten“ rechten Seite
fh = (f (x1), . . . , f (xn))T
und dem Randwertevektor
Rh = (1 + hp(x1))Ah−2, 0, . . . , 0, (1− hp(xn))Bh−2)T.
Wie man sieht, können hier durch hinreichend kleine Wahl des Diskretisierungsparameters h Systeme beliebiger Größe entstehen.
Man sieht aber auch gleich eine durch die Diskretisierung erzeugte spezielle Form des Pro-blems:
Die Matrix ist tridiagonal, unter geeigneten Voraussetzungen kann das System durch Gaus-selimination ohne Zeilenvertauschungen gelöst werden, dann spielt sich der ganze Lösungs-prozess im Speicher für die drei Diagonalen ab (vgl. Lineare Algebra I, z.B. [MV]), und die Lösung benötigt nicht mehr als 5n Operationen und für die Matrixelemente weniger als 3n Speichervariablen.
Aufgabe 2.3
Zeigen Sie, dass unter der Voraussetzung q(x) ≥ m > 0 für alle x ∈ [a, b] ein h0 > 0 existiert, so dass die Matrix Ah aus (26) für 0 < h < h0 strikt diagonal dominant und damit regulär ist(vgl. Satz von Gerschgorin, Lineare Algebra II).
Nichtlineare Randwertaufgaben Randwertaufgaben sind oft auch nichtlinear. Als mathematisch wohl am einfachsten zu behandelnde Form kann man diejenigen Aufga-ben ansehen, bei denen nur die rechten Seiten nichtlinear in den Unbekannten sind. Ein Beispiel-Problem dieser Art, das wir später immer wieder aufgreifen werden ist das soge-nannte Bratu-Problem.
Beispiel 2.4 (Bratu-Problem)
Das Bratuproblem ist in der eindimensionalen Form die folgende von einem Parameter λ abhängende Randwertaufgabe
−y′′(x) = λ exp(y(x)), x∈ [0, 1], y(0) = y(1) = 0. (27)
Anmerkungen 2.5
1. Für λ∈ [0, λ∗) mit einem sog. kritischen Parameter λ∗ ≈ 33.513890719 hat es jeweils zwei Lösungen, die bei λ∗ zusammenfallen und wie folgt dargestellt werden können
u(x)m =−2 ln ( cosh((x− 12)θ2) cosh(θ4) ) .
Dabei sind die beiden Lösungen durch die zwei reellen Lösungen θ1 und θ2 der 1D-Fixpunktgleichung θ =√2λ cosh ( θ 4 ) bestimmt.
2. Rechts von λ∗ gibt es keine Lösung. Vgl. dazu die Abbildung 11. (Man beachte aber, dass der Abbildung 11 Randwertaufgabe mit Grundintervall [−1, 1] zugrunde liegt. Das erklärt die Verschiebung des λ-Umkehrwertes.)
3. Die Bratu-Aufgabe findet in vielen Bereichen der Natur- und Ingenieurwissenschaft Anwendung: Bei der Berechnung der Brennstoffentzündung in der Theorie der Ver-brennung, als Modell für exothermische Reaktionen, im Chandrasekhar Modell der Expansion des Weltalls, in der Theory allgemeiner chemischer Reaktionen, in der Theory der Strahlungswärmeübertragung und in der Nanotechnologie. In allen An-wendungen ist von großer Bedeutung, dass Lösungen recht des kritischen Punktes λ∗ nicht existieren.
4. Die Aufgabe (27) ist die Gleichung, die den stationären Zustand der Reaktionsdiffu-sionsgleichung (71) von Seite 37 charakterisiert.
Der Parameter λ wird uns erst später interessieren. Wir nehmen für ihn hier einmal den festen Wert 0.5 an.
Wenden wir das obige Diskretisierungsverfahren von Seite 17 für p = q = 0, A = B = 0 und f (x) := 0.5 exp(y(x)) an, so erhalten wir das nichtlineare endlichdimensionale Problem
F (y) := 2 −1 0 · · · 0 −1 2 . .. ... ... 0 . .. ... ... 0 .. . . .. ... 2 −1 0 · · · 0 −1 2 y1 .. . .. . .. . yn −h2 2 exp(y1) .. . .. . .. . exp(yn) = 0. (28)
Im Gegensatz zur Aufgabe (27) ist dieses System - wie die meisten nichtlinearen Systeme - nicht geschlossen lösbar. Wie aus den Vorlesungen Anaylysis I-III bekannt sein sollte, muss man zur (approximativen) Lösung von (28) iterative Verfahren anwenden. Aus den genannten Vorlesungen sind sogar schon zwei Verfahren bekannt: Fixpunktiterationen nach Banach und das Verfahren sukzessiver Linearisierung nach Newton.
Wir werden auf beide in diesem Skript intensiv eingehen. Eine mögliche Fixpunktiterati-on zu Erzeugung einer hoffentlich gegen eine Lösung kFixpunktiterati-onvergierenden Folge y0, y1, y2, . . . könnte aus (28) schnell wie folgt entwickelt werden
2 −1 0 · · · 0 −1 2 . .. ... ... 0 . .. ... ... 0 .. . . .. ... 2 −1 0 · · · 0 −1 2 yk+11 .. . .. . .. . yk+1n = h 2 2 exp(yk 1) .. . .. . .. . exp(yk n) , k = 0, 1, 2, . . . .
Im Newtonverfahren berechnet man ausgehend von einem Näherungsvektor y0 für eine Lösung sukzessive eine Folge y1, y2, . . . von hoffentlich besseren Näherungen, indem man das Systeme
F (y) = 0 (29)
bei yk linearisiert
F (yk) + F′(yk)(y− yk)≈ F (y)
und anstelle der nichtlinearen Gleichung (29) ersatzweise die linearisierte Gleichung F (yk) + F′(yk)(y− yk) = 0
löst und yk+1 durch deren Lösung definiert.
Formal schreibt man das18 als
Bei der konkreten Rechnung vermeidet man natürlich die Bildung der Inversen von F′(yk). Man versteht−F′(yk)−1F (yk) nur als eine Kurzschreibweise für die Lösung δkdes linearen
Gleichungssystems
F′(yk)δk =−F (yk). (31)
Algorithmus 2.6 (Newton, n-ter Schritt) 1. Bilde F′(yk) und F (yk).
2. Löse (31) nach δk
3. Setze yk+1= yk+ δk.
Speziell für die Aufgabe (28) bemerkt man, dass die Jacobimatrix des Systems auch wieder eine Tridiagonalmatrix ist, so dass bei der Lösung von (31) wieder die obigen Anmerkung zur Ausnutzung der Tridiagonalität zur Anwendung kommen können. Wenn die Jacobi-Matrix dichter besetzt ist, wird man gegebenenfalls andere Techniken ausnutzen müssen, um den Newton-Schritt auszuführen oder um einen ausführbaren Ersatz für ihn zu finden.
B. Finite Elemente Sehr häufig ergeben sich Differentialgleichungen als Euler-Lagrange-Gleichungen gewisser Funktionale, deren Minima die gesuchten Funktionen sind.
Wenn für x∈ [0, 1] die Funktion u ∈ C1[0, 1] mit U (0) = a und U (1) = b die Verformung einer durch eine Streckenlast f ∈ C[0, 1] belasteten Saite ist, so minimiert u das Funktional der potentiellen Energie
P (u) = 1 2 ∫ 1 0 (u′)2dx− ∫ 1 0 f · u dx, (32) also
u = arg min{u∈C1[0,1] mit u(0)=a und u(1)=b}P (u) (33)
(Dabei haben wir - wie man das in der Mathematik gern macht - alle auftretenden Kon-stanten gleich 1 angenommen und gleich noch das Grundintervall zu [0, 1] skaliert.) Wenn u Minimum von P ist, so muss für alle kleinen v ∈ C1[0, 1] mit v(0) = v(1) = 0 (dies sichert, dass u+tv für alle t∈ R auch die Randbedingungen erfüllt) die Richtungsableitung von P bei u in Richtung von v gleich Null sein. Wegen
0 = dtdP (u + tv)|t=0 = dtd ( 1 2 ∫1 0(u′ + tv′) 2dx−∫1 0 f· (u + tv)dx ) t=0 = (∫01(u′+ tv′)· v′dx−∫01f · vdx ) t=0 = ∫01u′· v′dx−∫01f· vdx (34) gilt also ∫ 1 0 u′· v′dx− ∫ 1 0 f · vdx = 0 für alle v ∈ C01[0, 1], (35) wobei die tiefgestellte 0 bei C1
0[0, 1] signalisiert, dass dies die C1-Funktionen v mit Null-randwerten v(0) = v(1) = 0 sind.
Wenn wir annehmen, dass u ein weiteres Mal stetig differenzierbar ist19, können wir im ersten Term von (35) partiell integrieren mit dem Resultat20
∫ 1 0 −u′′· vdx −∫ 1 0 f · vdx = 0 für alle v ∈ C01[0, 1], (36) also ∫ 1 0 (u′′(x) + f (x))· vdx = 0 für alle v ∈ C01[0, 1]. (37) Für stetige u′′ und f folgt hieraus die sogenannte Euler-Lagrange-Gleichung zum Variati-onsproblem (32):
−u′′(x) = f (x), x∈ (0, 1). (38)
Wäre nämlich R(x) := u′′(ˆx) + f (ˆx)̸= 0 für ein ˆx ∈ (0, 1), etwa R(ˆx) = w > 0,
so gäbe es eine ganze Umgebung [ˆx− ε, ˆx + ε] von ˆx in (0, 1) mit R(x)≥ w/2 für x ∈ [ˆx − ε, ˆx + ε]. Wählt man nun eine Funktion ˆv in C01[−1, 1] mit
ˆ v ≥ 0 in [ˆx − ε, ˆx + ε], ˆ v = 0 in [0, 1]\ [ˆx − ε, ˆx + ε], ∫x+εˆ ˆ x−ε v(x)dx = I > 0,ˆ , so ist offenbar ∫ 1 0 R(x)ˆv(x)dx = ∫ x+εˆ ˆ x−ε R(x)ˆv(x)dx≥ w/2 · I > 0, was (37) widerspricht.
Zusammen mit den Randbedingungen des Variationsproblems ergibt (38) die Randwert-aufgabe
−u′′(x) = f (x), x∈ (0, 1), u(0) = a, u(1) = b. (39)
Diese Randwertaufgabe steht im gleichen Verhältnis zum Variationsproblem (33) wie das Gleichungssystem (8) zum Variationsproblem (9), und genau so wie man die Lösung von (8) durch Minimierung von F (u) über einen niederdimensionalen Unterraum approximieren kann, versucht man, das Randwertproblem (39) durch Minimierung von P (u) aus (32) über einen endlichdimensionalen (niederdimensionalen) Unterraum zu approximieren. Indem man von u zunächst die Funktion R(x) = a∗ (1 − x) + bx abzieht, findet man, dass V = u− R die randhomogene Randwertaufgabe
−V′′(x) = f (x), x∈ (0, 1), V (0) = 0, V (1) = 0 (40)
erfüllt.
Eine numerische Approximation einer Lösung dieser Randwertaufgabe kann man nun zu finden versuchen, indem man das Funktional P (V ) aus (32) über einen endlichdimensio-nalen Unterraum von C01[0, 1] minimiert.
Wenn uh einmal das als existent angenommene Minimum von P (u) über einen
(endlich-dimensionalen) Lösungsansatzraum Lh bezeichnet, so leitet man ganz analog und formal
die schwache Gleichung ∫ 1 0 u′h· vh′dx− ∫ 1 0 f · vhdx = 0 für alle vh ∈ Lh, (41) her.
In (14) haben wir gesehen, dass der Approximationsfehler eine Güteabschätzung in einer mit dem Problem verwobenen Norm erfüllt.
Dasselbe leitet man hier auf ähnlichem Weg her, wobei die Norm wie oben mit der schwa-chen Gleichung verbunden ist und so aussieht:
∥uh∥a :=
{∫ 1 0
u′h· u′hdx }1/2
Wenn man im Funktionenraum mit dieser Norm minimieren will, sollte der Raum in dieser Norm abgeschlossen sein. Schließt man den C01[0, 1] dazu in dieser Norma ab, gelangt man zum sogenannten „Sobolev-Raum“ W1
0[0, 1], der auf [0, 1] „schwach differenzierbaren Funk-tionen mit einer schwachen Ableitung“, deren Quadrat über [0, 1] integrierbar ist21. Wir werden das Konzept der schwachen Ableitung in dieser Vorlesung nicht erklären können22 Wichtig ist für uns nur, dass die stückweise stetig differenzierbaren Funktionen zu dieser Menge gehören, so dass wir als Unterraum beispielsweise den Raum der bezüglich eines vorgegebenen Gitters
0 = x0 < x1 <· · · < xn+1 = 1
stückweise linearen randhomogenen Funktionen uh(x) =
n
∑
k=1
Φk(t)αk
wählen können, worin die kte Hutfunktion
Φk(x) = x−xk−1 xk−xk−1 für x∈ [xk−1, xk], x−xk+1 xk−xk+1 für x∈ [xk, xk+1], 0 sonst
einfach mit dem Wert αk der Funktion vh an der Stelle xk gewichtet wird.
Die folgende Skizze zeigt eine solche Funktion für
x = [0 0.2 0.4 0.5 0.7 1] und y = [0 0.6 2 1.5 1.8 0]. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
2Zusammensetzung einer stückweise linearen Funktion aus Hutfunktionen
Abbildung 4: Stückweise lineare Funktion
21genau genommen. Lebesgue-integrierbar, aber dies werden wir in dieser Vorlesung nicht exakt erklären.
22Informationen darüber kann man z.B. in der paralleleln Vorlesung „Numerik partieller
Geht man nun mit dem Ansatz uh(x) =
∑n
k=1Φk(t)αk
in die Gleichung (41) und lässt dort die vh alle Hutfunktionen Φm, m = 1, . . . , n
durch-laufen, so gelangt man zu dem linearen Gleichungssystem
n ∑ k=1 ∫ 1 0 αkΦ′k(x)Φ′m(x)dx = ∫ 1 0 f (x)· Φm(x)dx, m = 1, . . . , n. (42)
Da das Integral ∑nk=1∫01Φ′k(x)Φ′m(x)dx nur dann von Null verschieden ist, wenn sich die Indizes k und m höchstens um den Wert 1 unterscheiden23, entsteht wieder ein tridiagonales lineares Gleichungssystem.
Funktionen wie die Hutfunktionen, die einen kleinen Träger haben24, heißen Finite Ele-mente. Die eben demonstrierte Ritz-Galerkin Approximation mit Hut-Funktionen ist ein Beispiel einer Finite Element Methode (FEM).
Genauso wie in (14) ist die FEM-Approximation in einem gewissen25 Sinne eine beste Approximation aus dem gewählten Finite-Element-Unterraum.
Die „Kunst der Finiten Elemente“ liegt deshalb darin, den Unterraum und seine Basis so zu wählen, dass
1. Das entstehende Gleichungssystem (42)gut lösbar ist. 2. Der Unterraum die Lösung gut approximieren kann.
c. Petrov-Galërkin Wenn man in der schwachen Gleichung (41) mit dem Ansatz uh(x) =
∑n
k=1Φk(t)αk für die vh anstelle der den Ansatzraum aufspannenden Funktionen Φk(t)
an-dere Funktionen Ψk(x) als sogenannte „Testfunktionen“ wählt, mit dem Ergebnissystem n ∑ k=1 ∫ 1 0 αkΦ′k(x)Ψ′m(x)dx = ∫ 1 0 f (x)· Ψm(x)dx, m = 1, . . . , n. (43)
so spricht man von einem Petrov-Galerkin-Vorgehen anstelle eines Ritz-Galerkin-Verfahrens. Manchmal setzt man die Summe uh(x) =
∑n
k=1Φk(t)αkauch einfach in die
Differentialglei-chung der Randwertaufgabe (40) ein26, multipliziert die ganze Gleichung für m = 1, . . . , n mit den Testfunktionen Ψm, m = 1, . . . , n und integriert von 0 bis 1:
− n ∑ k=1 ∫ 1 0 αkΦ′′k(x)Ψm(x)dx = ∫ 1 0 f (x)· Ψm(x)dx, m = 1, . . . , n. (44)
Anmerkungen 2.7 (Petrov-Galerkin für den endlichdimensionalen Fall) Bei der Ritz-Approximation der Lösung eines (n× n)-Gleichungssystems
Ax = b
durch Approximationen aus einem Teilraum Rh mit Basisvektoren v1, . . . , vm wurde man
für die Näherung vh = V y auf das „projizierte Problem“
VT(AV y− b) = 0
23Andernfalls sind sie auf keinem gemeinsamen Teilintervall von Null verschieden.
24Der Träger einer Funktion ist die Menge, auf der sie nicht verschwindet.
geführt (vgl. Anmerkung 2.2).
Bei Matrizen A, die nicht SPD sind, kann man nicht erwarten, dass die „projizierte Ma-trix“ VTAV weiter SPD ist, ja noch nicht einmal die Regularität kann regelhaft ausgesagt
werden.
Beim Petrov-Galerkin-Verfahren wird das überbestimmte Gleichungssystem AV y− b = 0
durch „Testen der Gleichung“ mit von v1, .., vm verschiedenen Vektoren w1, . . . , wm in der
Bilddimension auf m reduziert:
(wk)T(AV y− b) = 0, k = 1, . . . , m oder - mit W = (w1, . . . , wm)
-WTAV y = WTb. (45)
Durch geeignete Wahl von W kann die Systemmatrix WTAV manchmal angenehme
Ei-genschaften erhalten.
d. Kollokationsverfahren Wählt man in (44) für die Ψm(x) Delta-Distributionen
Ψm(x) = δzm(x)
so dass „das Integral einer mit Ψm(x) multiplizierten Funktion f deren Wert bei zm ergibt“,
so erhält man. − n ∑ k=1 αkΦ′′k(zm) = f (zm), m = 1, . . . , n. (46)
Dies bezeichnet man als „Kollokationsverfahren“.
0 5 10 15 20 −2 −1 0 1 2 3 4 5 6
u
iu
i−1u
i+1x
1x
2x
i−1x
ix
i+1b
ib
i+1V
i
Finite−Volumen−DiskretisierungC. Finite Volumen Bei der Finite-Volumen-Diskretisierung der Randwertaufgabe (39), die wir hier27 noch einmal (mit anderem Grundbereich) wiederholen
−u′′(x) = f (x), x∈ (0, 1), u(0) = a, u(20) = b, (47)
werden zu einem nicht notwendig äquidistanten Gitter 0 = x1 < x2 < · · · < xi−1 <
xi < xi+1 <. .. < xN = 2 sogenannte Kontroll-Volumina Vi = [bi, bi+1], i = 2, . . . , N − 1
eingeführt, mit xi−1 < bi < xi. Zur Bestimmung von Approximationen ui der Lösungswerte
u(xi), i = 2, . . . , N − 1 werden N − 2 Gleichungen wie folgt erzeugt:
1. Es wird die Differentialgleichung über das Kontrollvolumen Vi integriert:
− ∫ bi+1 bi u′′(x)dx =− ∫ bi+1 bi f (x)dx
2. Die rechte Seite dieser Gleichung wird entweder exakt ausgerechnet oder durch Qua-dratur approximiert mit dem Ergebnis Fi.
3. Auf der linken Seite wird einmal integriert, so dass sich
u′(bi)− u′(bi+1) = Fi, i = 2, . . . , N− 1
ergibt.
4. Die Ableitungen in den letzten Gleichungen werden nun durch Approximationen er-setzt, die man mit den ui-Werten ausdrückt. Da die Punkte bi in der Abbildung
gerade den Mitten des einschließenden Intervalle [xi−1, xi] liegen, bieten sich die
Nä-herungen u′(bi)≈ ui− ui−1 xi− xi−1 , i = 2, . . . , N− 2 an.
Dies führt wieder auf ein lineares Gleichungssystem für die ui-Werte, wobei für die
Rand-punkte natürlich keine Gleichungen aufgestellt zu werden brauchen.
Es können auch ganz andere Approximationen für die u′(xi)-Werte verwendet werden. Für
das Verfahren wichtig ist, dass an der Grenze zweier aneinanderstoßender Kontrollvolu-mina dieselbe Approximation verwendet wird. Die Ableitung am Rand entspricht einem Fluss von einem Volumen in das andere. Was aus dem einen herausfließt, sollte aus Erhal-tungsgründen in das anschließende hineinfließen.
Im Falle partieller Differentialgleichungen wird dies noch klarer werden.
2.1.3 Partielle Randwertaufgaben
Partielle Differentialgleichungen werden auf ganz analoge Weise wie die gewöhnlichen dis-kretisiert.
A. Finite Differenzen Als einfaches Beispiel der FD-Diskretisierung einer elliptischen Randwertaufgabe28 betrachten wir die Diskretisierung des Poisson-Problems
−(∆u)(x, y) = 1, (x, y) ∈ Ω = [0, 1]2, u
|∂Ω = 0. (48)
In Absatz A. von Abschnitt 2.1.2 wurde erklärt, wie aus einer gewöhnlichen Randwertauf-gabe
−y′′ = f (x, y); y(a) = A, y(b) = B
mit Hilfe der Differenzenapproximation
y(x− h) − 2y(x) + y(x + h)
h2 = y
′′(x) + O(h2) (49)
ein Gleichungssystem für Näherungen
yi ≈ y(a + ih), i = 1, . . . , n; h = (b − a)/(n + 1)
der Lösung y(x) eingeschränkt auf ein äquidistantes Gitter a = x0 < x1 < . . . < xn< xn+1 = b
erhalten werden kann.
Zur numerischen Behandlung der Aufgabe (48) ersetzen wir die partiellen Ableitungen des Laplace-Operators ∆u(x, y) = ∂ 2u(x, y) ∂x2 + ∂2u(x, y) ∂y2
auf einem zweidimensionalen Gitter (vergleiche Abbildung 6) ganz analog zu (49) durch Differenzen
u(x− h, y) − 2u(x, y) + u(x + h, y)
h2 = ∂2u(x, y) ∂x2 + O(h 2 ) und
u(x, y− h) − 2u(x, y) + u(x, y + h)
h2 =
∂2u(x, y)
∂y2 + O(h 2
)
und schreiben für jeden der inneren Punkte (xi, yj) = (ih, jh), i, j = 1, . . . , 4; h = 0.2
mit Approximationen ui,j ≈ u(xi, yj), i, j = 1, . . . , 4 eine Ersatzgleichung für die
Differen-tialgleichung aus (48) hin29:
−ui−1,j− ui+1,j + 4ui,j − ui,j−1− ui,j+1
h2 = 1, i, j = 1, 2, 3, 4, (50)
Mit dem „diskreten Laplaceoperator“ wird üblicherweise ein sogenannter Differenzenstern assoziert. 1 h2 −1 −14 −1 −1 (51)
Seine Anwendung an der Stelle (xi, yj) führt auf die linke Seite von (50).
Dabei werden die Randwerte eingesetzt, wenn Punkte auf dem Rand ∂Ω von Ω liegen.
28Der Unterschied zwischen elliptischen, parabolischen und hyperbolischen Differentialgleichungen sollte
schon in den einführenden Mathematik-Vorlesungen des Bachelor-Studiums erklärt worden sein. Sollte dies nicht der Fall sein, konsultuiere man z.B. das Skript [Urb]
Für den Punkt (x3, y2) liegen alle zu berücksichtigenden Nachbarn in Ω, so dass die Glei-chung
−u2,2− u4,2+ 4u3,2− u3,1− u3,3 = h2
lautet. Für (x2, y1) liegt der Nachbarpunkt (x2, y0) auf dem Rand, so dass die Gleichung hier die folgende Form hat:
−u1,1− u3,1+ 4u2,1− u2,2 = h2
Abbildung 6: Gitterfunktion
Nun handelt es sich bei den 4× 4 Gleichungen für die 16 Unbekannten ui,j, i, j = 1, .., 4
sicher um ein lineares Gleichungssystem. Um dies in der Form Matrix mal Vektor = Vektor ⇐⇒ Au = f
aufschreiben zu können, muss die Matrix aus Unbekannten (ui,j)i,j=1,...,4 in einen Vektor
(uk)k=1,...,16 verpackt werden. Da die rechten Seiten alle gleich h2sind, macht das Aufstellen
des Vektors f keine Probleme.
Das Einordnen der ui,j in den Vektor u ∈ R16 ist nichts anderes als ein Numerieren der
0 2 4 6 8 10 0 2 4 6 8 10 u 1 u 2 u 3 u4 u 5 u6 u7 u8 u 12 u 11 u 10 u 9 u 13 u14 u 15 u 16 Abbildung 7: Numerierung so bekommt die Matrix A die folgende Gestalt:
A = 4 −1 0 0 −1 0 0 0 0 0 0 0 0 0 0 0 −1 4 −1 0 0 −1 0 0 0 0 0 0 0 0 0 0 0 −1 4 −1 0 0 −1 0 0 0 0 0 0 0 0 0 0 0 −1 4 0 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 0 4 −1 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −1 4 −1 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −1 4 −1 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −1 4 0 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 0 4 −1 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −1 4 −1 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −1 4 −1 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −1 4 0 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 0 4 −1 0 0 0 0 0 0 0 0 0 0 0 −1 0 0 −1 4 −1 0 0 0 0 0 0 0 0 0 0 0 −1 0 0 −1 4 −1 0 0 0 0 0 0 0 0 0 0 0 −1 0 0 −1 4 .
Die Linien wurden dabei eingezogen, um die Block-Struktur der Matrix deutlich zu machen. Die Diagonalblöcke B = 4 −1 0 0 −1 4 −1 0 0 −1 4 −1 0 0 −1 4
verbinden die Unbekannten in jeweils derselben Reihe des ui,j-Gitters. Der in A ganz
oben links stehende B-Block regelt die Differenzenbeziehung zwischen der untersten Reihe der ui,j. Die negativen Einheitsmatrizen links von den Diagonalblöcken B greifen auf die
vorhergehende Schicht zu. In der ersten Blockzeile gibt es die noch nicht. Hier würde diese negative Einheitsmatrix auf die Randwerte zugreifen, und wenn diese nicht Null wären, würden sie als bekannte Werte mit auf die rechte Seite zum f -Vektor kommen.
Genauso greift die negative Einheitsmatrix rechts der Diagonalblöcke auf die Folgeschicht zu.
Wenn wir die Matrix SO aufstellten, könnten wir die Approximationen in MATLAB einfach durch
u = A\f erhalten30
Wenn man die Lösungskomponenten wieder in eine Gitterform gebracht hat (in die man auch noch die Randwerte geeignet einbringt) kann man die Lösung im MATLAB leicht plotten mit dem Ergebnis:
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Lösung von −∆ u = 1, u ∂Ω=0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
Abbildung 8: Lösungsapproximation zu Aufgabe (48)
B. Finite Elemente Wie zur Differentialgleichung (38) die schwache Gleichung (35) und die Variationaufgabe (33) zum Funktional (32) gehörten, gehört zur partiellen Rand-wertaufgabe
−(∆u)(x, y) = f(x, y), (x, y) ∈ Ω = [0, 1]2, u
|∂Ω= 0. (52)
die schwache Gleichung ∫ Ω ∇u · ∇v = ∫ Ω f vd(x, y), für alle v∈ H01(Ω) (53)
und das Variationsproblem
u = arg minu∈H1 0(Ω)P (u) (54) mit P (U ) := ∫ Ω ( 1 2∥∇u(x, y)∥ 2 2− f(x, y) · u(x, y) ) d(x, y). (55)
Der Übergang zwischen der Differentialgleichung und der schwachen Form, die sich im Wesentlichen durch Anwendung der partiellen Integration auf den ersten Summanden im
Integral ∫
1 0
(−u′′(x)− f(x)) v(x)dx = 0
ergab, erhält man hier aus dem Gaussschen Integralsatz: ∫ Ω f vd(x, y) = − ∫ Ω div(∇u)vd(x, y) (56) = − ∫ Ω div(v∇u)d(x, y) + ∫ Ω ∇u · ∇vd(x, y) (57) = ∫ Ω ∇u · ∇vd(x, y). (58)
Analog zum oben diskutierten eindimensionale Fall sucht man wieder ein Minimum von P (u) in einem endlichdimensionalen Unterraum.
Die sich dabei aus der schwachen Gleichung (53) für eine Ansatzfunktion uh(x, y) =
∑n
k=1Φk(x, y)αk ergebenden Gleichungen sehen strukturell genau wie (42) aus: n ∑ k=1 ∫ Ω αk∇Φk(x, y)· ∇Φm(x, y)d(x, y) = ∫ Ω f (x, y)· Φm(x, y)dx, m = 1, . . . , n. (59)
Um zu den Hutfunktionen bei eindimensionalem Grundgebiet analoge Funktionen definie-ren zu können, teilt man das Gebiet Ω möglichst gut in Dreiecke ein, von denen je zwei höchstens eine Seite oder eine Ecke gemeinsam haben. Zu gegebenen Funktionswerten in den Ecken wird als Ansatzfunktion die auf den Dreiecken bilinearen Funktionen gewählt, die an den Dreiecksecken die gegebenen Werte haben. Eine Hutfunktion zu einem inneren Dreieckspunkt erhält man dann, indem man diese Funktion in diesem Punkt gleich eins und in allen anderen Punkte gleich Null wählt.
Noch stärker als im eindimensionale Fall haben die Finiten Elemente gegenüber den Finiten Differenzen den Vorteil, dass sie durch nichtgleichmäßige Gitter von Funktionswerten -besser möglichen lokalen schnellen Änderungen der Lösungen anpassen können.
In Figur 9 ist ein Finite-Element-Gitter dargestellt sowie die Hutfunktion, welche in allen Gitterpunktem den Wert 0, im Gitterpunkt Nr 1 jedoch den Wert 1 hat.
0 1 2 3 4 5 6 7 8 9 10 −1 0 1 2 3 4 5
Trinangulierung mit vier inneren Punkten
1 2 3 4 0 2 4 6 8 10 0 2 4 −0.2 0 0.2 0.4 0.6 0.8 1 Basisfunktion P 1
C. Finite Volumen Die Gleichung ∆u = f diskretisiert man mit der Finite Volumen Methode nach der folgenden an Figur 10 orientierten Vorgehensweise:
Man teilt das Gebiet in Kontrollvolumina Si ein, wovon jedes einen Punkt mit
Funktions-wert uienthält31. Nützlich ist es, wenn - wie in der Skizze - die Ränder der Kontrollvolumina
auf den Verbindungslinien der Punkte senkrecht stehen32.
−1 0 1 2 3 4 5 6 −1 0 1 2 3 4 5 6 Finite Volumen im R2 S 3 S 4 S 6 u 1 u 2 u 3 u 4 u S5 5 u 6 S 2 L 12 L 14 L 16 L 13 L 15 S 1
Abbildung 10: Finite Volumen Diskretisierung, Erkl. siehe Text
Nun gewinnt man eine Gleichung, die die über Ihre zusammanstoßenden Volumina benach-barten Werte u1 bis u6 verbindet, wie folgt: Man integriert die Differentialgleichung über
S1: ∫ S1 div(grad u)dx = ∫ S1 f dx. Dann wendet man den Gaussschen Integralsatz an und erhält
∫ ∂S1 (grad u)· ν dx = ∫ S1 f dx,
wobei ν die äußere Normale des Rands ∂S1 (an den glatten Stellen des Randes) bezeichnet. Das Randintegral wird aufgeteilt in die Teilrandbereiche ∂S1,i, die die Bereiche S1 und Si
trennen für i = 2, . . . , 6. 6 ∑ i=2 ∫ ∂S1,i grad u· νi dx = ∫ S1 f dx.
Nun wird auf den Rändern die Werte grad u· νi mit Hilfe der gegebenen Funktionswerte
genähert, etwa durch
(gradu· νi)|S1,i ≈
ui− u1 ∥xi− x1∥2
,
31Tatsächlich kann der Wert u
i auch noch andere Bedeutungen als „Wert an einer Stelle“ haben, etwas
Mittelwert über Si.
wobei die Vektoren xi die Positionen der Funktionsauswertungen ui bezeichnen. Damit
ergibt sich dann die Ersatzgleichung 6 ∑ i=2 ui− u1 ∥xi− x1∥2 · L1i= QuadS1f,
wobei QuadS1 eine Kubatur-Formel bezeichnet, die das Integral von f über S1 nähert. Die Approximationen für die Randwerte können natürlich auch noch genauer gemacht werden. Wichtig (und charakteristisch für FV-Ansätze) ist, dass in den Gleichungen für verschiedene aneinanderstoßende Kontrollvolumina für die Berührungskanten stets diesel-ben Näherungen genommen werden.
2.1.4 Integralgleichungen
Eine „lineare Integralgleichung“ für eine Funktion u :R ⊃ [a, b] −→ R hat die Gestalt λ(x)u(x) +
∫ b a
k(x, y)u(y) dy = f (x), x∈ [a, b]. (60)
Dabei heißt die Funktion k : [a, b]2 −→ R „Kern“ der Integralgleichung. Für festes x ist k(x, y) bezüglich y hinreichend „nett“, so dass das Integral für - sagen wir - stetige Funktio-nen u existiert und möglichst wieder auf eine stetige Funktion K(u)(x) :=∫abk(x, y)u(y) dy führt33.
Definition 2.8 (Integralgleichungs-Typen)
Je nach dem Aussehen der Funktion λ(x) erhält die Gleichung (60) Zusatznamen. Sie heißt Integralgleichung 1. Art wenn λ(x) = 0 für alle x
Integralgleichung 2. Art wenn λ(x) = λ∈ R (oder in C) Integralgleichung 3. Art sonst
In Abhängigkeit von den Integrationsgrenzen wird weiter aufgeteilt in Fredholmsche Gleichungen für x-unabhängige Grenzen und Volterrasche Gleichungen bei x-Abhängigkeit
mindestens einer Integrationsgrenze.
Wir werden als Beispiele nur Fredholmsche Gleichungen 2. Art betrachten. Solche Aufgaben λu(x) +
∫ b a
k(x, y)u(y) dy = f (x), x∈ [a, b]. führt man durch Anwendung einer Quadraturformel
∫ b a g(x)dx≈ n ∑ j=1 Wjg(xj)
mit den Quadraturknoten a≤ x1 < x2 <. .. < xn≤ b und den Gewichten Wj ∈ R über in λu(x) + n ∑ j=1 Wjk(x, xj)u(xj) dy = f (x), x∈ [a, b]. (61)
Setzt man für x hierin die xi-Werte ein, gelangt man zu einem linearen Gleichungssystem
zur Bestimmung von Näherungen ui von u(xi) für i = 1, . . . , n über:
λui+ n
∑
j=1
Wjk(xi, xj)uj = f (xi), i = 1, . . . , n. (62)
Nystrøm (vgl. [Atk], Seite 88ff) hat bemerkt, dass man nach Berechnung der Werte ui, i =
1, . . . , n die Gleichung (61) als systemangepasste Interpolationsformel benutzen kann, um eine kontinuierliche Näherung uh(x) für die Lösung u(x) zu erhalten:
uh(x) =−λ−1 n
∑
j=1
Wjk(x, xj)uj+ f (x), x∈ [a, b].
Da die Integralkerne k(x, y) normalerweise an fast allen Punkten von Null verschieden sind, ist für die Diskretisierung (62) im Gegensatz zu Diskretisierungen von Differentialgleichun-gen typisch, dass die Systemmatrizen voll besetzt sind.
2.2
Entstehen von Größe aus Nichtlinearität
Abgesehen vom Problem (27) waren alle bisherigen Beispielprobleme lineare Aufgaben und damit nach der Diskretisierung dem Aufgabengebiet der numerischen linearen Algebra zugehörig. Das Problem (28) zeigt eine strukturell relativ harmlose Nichtlinearität. Die einzelnen Gleichungen hängen nichtlinear jeweils nur von einer einzigen Komponente ab. Wenn man die Jakobimatrix von
F (y) := 2 −1 0 · · · 0 −1 2 . .. ... ... 0 . .. ... ... 0 .. . . .. ... 2 −1 0 · · · 0 −1 2 y1 .. . .. . .. . yn − h2 2 exp(y1) .. . .. . .. . exp(yn) (= 0) (63)
bildet, trägt der nichtlineare Term mit den Exponentialfunktionen nur zur Diagonale der Matrix bei34, so dass die Struktur der Matrix nicht zerstört wird.
Aber auch dieses einfache nichtlineare Diagonalfeld weist35 ein typisches Phänomen bei nichtlinearen Problemen auf: Die Aufgabe hat mehr als eine Lösung, von der aber oft nur eine von Interesse ist. In Figure 11 ist nur die untere Lösung von praktischem Interesse, weil die obere Lösung in Bezug auf die der Aufgabe zugrunde liegende Zeitentwicklungsauf-gabe36 instabil und daher in der Wirklichkeit nicht beobachtbar ist37. Hat man die letzte
34Man nennt dies dann auch ein Diagonalfeld.
35Wie das approximierte kontinuierliche Bratu-Problem
36Vgl. Seite 37).
37Bitte wiederholen Sie: Wann werden stationäre Punkte von Differentialgleichungssystemen als stabil,
