• Keine Ergebnisse gefunden

Numerik großer nichtlinearer Systeme

N/A
N/A
Protected

Academic year: 2021

Aktie "Numerik großer nichtlinearer Systeme"

Copied!
290
0
0

Wird geladen.... (Jetzt Volltext ansehen)

Volltext

(1)

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

(2)

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

(3)

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

(4)

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

(5)

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

(6)

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

(7)

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

(8)

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

(9)

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.

(10)

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.

(11)

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) := nk=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

(12)

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 = nk=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Φ11+ (AΦ22+· · · + (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.

(13)

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

(14)

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

(15)

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.

(16)

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

(17)

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

(18)

−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

(19)

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.)

(20)

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

(21)

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.

(22)

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.

(23)

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

(24)

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

nk=1 ∫ 1 0 αkΦ′k(x)Φ′m(x)dx = ∫ 1 0 f (x)· Φm(x)dx, m = 1, . . . , n. (42)

Da das Integral ∑nk=10k(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 nk=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:

nk=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.

(25)

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. nk=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

i

u

i−1

u

i+1

x

1

x

2

x

i−1

x

i

x

i+1

b

i

b

i+1

V

i

Finite−Volumen−Diskretisierung

(26)

C. 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.

(27)

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]

(28)

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

(29)

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.

(30)

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

(31)

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: nk=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

(32)

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− x12

,

31Tatsächlich kann der Wert u

i auch noch andere Bedeutungen als „Wert an einer Stelle“ haben, etwas

Mittelwert über Si.

(33)

wobei die Vektoren xi die Positionen der Funktionsauswertungen ui bezeichnen. Damit

ergibt sich dann die Ersatzgleichung 6 ∑ i=2 ui− u1 ∥xi− x12 · 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≈ nj=1 Wjg(xj)

(34)

mit den Quadraturknoten a≤ x1 < x2 <. .. < xn≤ b und den Gewichten Wj ∈ R über in λu(x) + nj=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,

Referenzen

ÄHNLICHE DOKUMENTE

Die Bedeutung des Satzes ist nicht, dass das kontinuierliche und das diskrete Dynamische System sich f¨ur jeden Startwert gleich verhalten, sondern dass es zu jedem Startwert

Eine mögliche Notation für die Darstellung der (statischen)

Geben Sie außerdem die exakte L¨ osung der Gleichung an und vergleichen Sie die beiden L¨ osungen, indem Sie den globalen Fehler auf dem Gitter angeben.. Laden Sie sich dazu die

Selbst wenn Dennis achtmal mehr Flaschen hätte, würde Falk immer noch 12 mehr haben als Dennis und Bernd zusammen.. Stellen Sie ein lineares Gleichungssystem zur Berechnung

Beim zweiten Einkauf habe ich eine Packung Fischfutter und 2 Packungen Hundefutter gekauft. Dabei musste ich 2 Packungen Katzenfutter vom ersten Einkauf

Zeigen Sie, dass unter dieser Voraussetzung jeder Endomorphismus von V Summe zweier diagonalisierbarer Endomorphismen ist.. Laza: Lineare Algebra individuell Online-Version

Tutorium zur Linearen Algebra I Blatt 12..

Ubungen zur Linearen Algebra II ¨ Bergische Universit¨ at Wuppertal.. Blatt 0