finite elemente i - tu-freiberg.deernst/lehre/fem_i/folien12/...finite elemente i 2 1 einleitung...
TRANSCRIPT
-
Institut für Numerische Mathematik und Optimierung
Finite Elemente IWintersemester 2012/13
Erste von zwei Vorlesungen im Modul
”Finite-Element Methoden für Mathematiker“
Hörerkreis: 5. Mm, 7. Mm, 9. Mm, 1. MWM
Oliver ErnstInstitut für Numerische Mathematik und Optimierung
-
Finite Elemente I 1
0 Vorbemerkungen
Vorgesehener Inhalt:
1. Einleitung:
2. Variationstheorie
3. Die FE-Methode anhand der Poisson-Gleichung
4. Kontruktion von FE-Räumen
5. Konvergenztheorie
6. Gleichungslöser
0 Vorbemerkungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 2
1 Einleitung
Bevor wir in die Details von Variationstheorie und Konstruktion von finite-Element-Räumen einsteigen wollen wir uns einen Eindruck verschaffen,wozu die FEM in der Praxis eingesetzt werden kann.
1 Einleitung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 3
1.1 Die Poisson-Gleichung
Viele physikalische Größen erfüllen eine elliptische Differentialgleichungzweiter Ordnung der Form
−∇·(k∇u) = f, u : Ω→ R (1.1)
wobei k : Ω → R eine positive Koeffizientenfunktion, f : Ω → R einensog. Quellterm darstellt und das Definitionsgebiet Ω ⊂ Rd, d = 2, 3, als be-schränkt, offen und zusammenhängend angenommen wird. Der Koeffizientk kann skalar oder auch ein positiv-definiter Tensor (d× d Matrix) sein undbeschreibt meist Materialeigenschaften.
Die typische Problemstellung besteht darin, u zu gegebenen f und k zubestimmen. (Hinzu kommen noch Randbedingungen auf ∂Ω.)
1.1 Die Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 4
Folgende Tabelle stellt einige physikalische Größen zusammen, welcheder Gleichung (1.1) genügen:
Anwendung u k f
Elektrostatik elektr. Potential Permittivität Ladungsdichte
Magnetostatik magn. Potential Permeablilität Ladungsdichte
Wärmetransport Temperatur Leitfähigkeit Wärmequelle
Grundwasserströmung Spiegelhöhe hydraul. Leitf. GW-Neubildung
elastische Membran Auslenkung M.-spannung Last
ideales Fluid Geschw.-Potential Dichte Zu/Abfluss
Gravitation Grav.-Potential 1/Grav.-Konst. Massendichte
Bemerkung 1.1 Oft ist u nur eine Hilfsgröße und der Flussvektor ∇u die wirklichinteressierende Größe. Letzterer wird oft durch (numerische) Differentiation der (nu-merisch berechneten) Lösung u gewonnen, ein instabiler Vorgang. Sog. gemischteFE-Formlierungen gestatten die direkte numerische Berechnung des Flusses.
1.1 Die Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 5
Ist k konstant in (1.1), so kann dies in f zusammengefasst werden und esergibt sich die Poisson-Gleichung
−∆u(x) = f(x), x ∈ Ω. (1.2a)
Den Gebietsrand Γ := ∂Ω zerlegen wir in ΓD und ΓN , Γ = ΓD ∪ ΓN ,ΓD ∩ ΓN = ∅ und stellen die Randbedingungen
u(x) = g(x) ∀x ∈ ΓD, kurz: u|ΓD = g, (1.2b)∂u
∂n(x) = h(x) ∀x ∈ ΓN , kurz: ∂nu|ΓN = h (1.2c)
mit zwei gegebenen, auf ΓD bzw. ΓN definierten Funktionen g und h. DasRandwertproblem (1.2) besitzt – unter geeigneten Voraussetzungen an dasGebiet Ω und die Daten f, g, h – eine eindeutig bestimmte klassische (d.h.in Ω zweimal stetig differenzierbare) Lösung.
1.1 Die Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 6
Wir betrachten als Beispiel das elektrische Feld im Außenraum zweier Elek-troden in einem Gehäuse. Die Anordnung sei uniform in einer Richtung: dieElektroden mögen recheckigen Querschnitt besitzen, das Gehäuse sei einKreiszylinder.
Auf den Elektroden Γ1 und Γ2 legen wir jeweils den konstanten Potential-wert φ = 1 bzw. φ = −1 fest, am Gehäuse Γ0 den Wert φ = 0. Ansonstensei die Anordnung ladungsfrei.
Es ergibt sich die elliptische Randwertaufgabe
∆φ = 0, in Ω,
φ = 0, auf Γ0,
φ = 1, auf Γ1,
φ = −1, auf Γ2.
1.1 Die Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 7
Surface: u Contour: u Arrow: grad(u)
-6
-4
-2
0
2
4
6
-6 -4 -2 0 2 4 6-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Max: 1.00
Min: -1.00
-0.95
-0.85
-0.75
-0.65
-0.55
-0.45
-0.35
-0.25
-0.15
-0.05
0.05
0.15
0.25
0.35
0.45
0.55
0.65
0.75
0.85
0.95
Max: 0.950
Min: -0.950
FEM-Approximation des Potentials im Außenraum der Kondensatoren.
1.1 Die Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 8
1.2 Die Helmholtz-Gleichung
Monochromatische Wellen mit Kreisfrequenz ω > 0 besitzen die Darstel-lung
w(x , t) = u(x )eiωt, x ∈ R, t > 0.
Einsetzen in die lineare Wellengleichung wtt − c2∆u = 0 mit Ausbreitungs-geschwindigkeit c > 0 liefert für den räumlichen Anteil u die Helmholtz-Gleichung (auch reduzierte Wellengleichung oder Schwingungsgleichung)
∆u+ k2u = 0, k =ω
c> 0
mit Wellenzahl k.
Man kann mit der Helmholtz-Gleichung z.B. die Akustik eines Konzertsaalsmodellieren.
1.2 Die Helmholtz-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 9
Hierzu sei der Querschnitt des Konzertsaals gegeben durch ein Polygon.Auf den Rändern stellen wir die sogenannte schallharte Randbedingung
∂u
∂n= 0.
Als Quellterm der Gleichung nehmen wir die exponentiall abklingendeFunktion
f(x, y) = e−10((x−1/2)2+(y−1/2)2).
Mit der Wellenzahl k =√
0.8 erhalten wir die Randwertaufgabe
∆u+ 0.8u = e−10((x−1/2)2+(y−1/2)2), in Ω,
∂u
∂n= 0, auf ∂Ω.
1.2 Die Helmholtz-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 10
-1
-0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
-0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5
Konzertsaal-Beispiel, Triangulierung
Surface: u Contour: u
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
-1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5
-0.2
-0.1
0
0.1
0.2
0.3
Max: 0.305
Min: -0.250
-0.236
-0.208
-0.181
-0.153
-0.125
-0.097
-0.07
-0.042
-0.014
0.014
0.042
0.069
0.097
0.125
0.153
0.18
0.208
0.236
0.264
0.292
Max: 0.292
Min: -0.236
Konzertsaal-Beispiel, Lösung
1.2 Die Helmholtz-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 11
1.3 Minimalflächen
Die durch eine Funktion u : Ω → R, Ω ⊂ R2 ein beschränktes Gebiet,definierte Fläche
{z = u(x, y) : (x, y) ∈ Ω}
ist gegeben durch
A(u) =
∫Ω
√1 + |∇u|2 dx .
Schreibt man noch die Randwerte der Funktion u vor, so kann man zeigen,dass die Fläche A(u) minimal wird, wenn u folgende Randwertaufgabelöst:
−∇·
(1√
1 + |∇u|2∇u
)= 0, in Ω,
u = g, längs ∂Ω,
wobei g die vorgeschriebenen Randwerte darstellt.
1.3 Minimalflächen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 12
Wir betrachten konkret das Gebiet Ω := {(x, y) ∈ R2 : x2 + y2 < 1} undschreiben auf dessen Rand für u die Funktionswerte u(x, y) = x2 vor.
Surface: u Height: u Contour: u Max: 1.00
Min: 00
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1Max: 0.975
Min: 0.02500.025
0.075
0.125
0.175
0.225
0.275
0.325
0.375
0.425
0.475
0.525
0.575
0.625
0.675
0.725
0.775
0.825
0.875
0.925
0.975
Minimalfläche
1.3 Minimalflächen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 13
1.4 Die Navier-Stokes-Gleichungen
Die Strömung eines inkompressiblen homogenen Newtonschen Fluids inzwei Raumdimensionen wird beschrieben durch ein Vektorfeld
u = u(x , t) =
[u(x , t)
v(x , t)
], u : Ω ⊂ R2 → R2
sowie eine skalare Funktion
p = p(x , t), p : Ω→ R,
die den Geschwindigkeitsvektor bzw. den Druck am Ortspunkt x ∈ Ω zurZeit t angeben.
1.4 Die Navier-Stokes-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 14
Die Größen u und p sind eindeutig gegeben als Lösung des Systemspartieller Differentialgleichungen
∂tu + (u · ∇)u − ν∆u +∇p = f ,∇·u = 0 ,
(ν bezeichnet die kinematische Viskosität des Fluids) zusammen mit ge-eigneten Anfangsbedingungen
u0 = u(x , 0), p0 = p(x , 0), x ∈ Ω,
sowie Randbedingungen längs ∂Ω.
1.4 Die Navier-Stokes-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 15
Als Besipiel betrachten wir das sog. DFG-Benchmarkproblem der Strömungin einem Kanal um einen Kreiszylinder:
Am linken Rand wird ein parabolisches Einströmprofil[u
v
]=
[1.2y(0.41− y)/0.412
0
]vorgeschrieben, am rechten Rand Neumann-Randbedingungen und anden übrigen Rändern u = 0 . Ferner ist ν = 0.001.
1.4 Die Navier-Stokes-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 16
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
Gitter aus 3608 Dreiecken (16 836 Freiheitsgrade)
1.4 Die Navier-Stokes-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 17
Time=6.98 Surface: Velocity field [m/s] Arrow: Velocity field
-0.6
-0.4
-0.2
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 2 2.2
Max: 2.174
Min: 00
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Strömungsgeschwindigkeitsfeld bei t = 7s.
1.4 Die Navier-Stokes-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 18
1.5 Die Maxwell-Gleichungen
Makroskopische elektromagnetische Phänomene werden beschriebendurch die Maxwell-Gleichungen
∇×E + Ḃ = 0 , (1.3a)
∇×H − Ḋ = J , (1.3b)∇·D = ρ, (1.3c)∇·B = 0, . (1.3d)
Hierbei bezeichnen E die elektrische Feldstärke, D die Verschiebungs-stromdichte, H die magnetische Feldstärke, B die magnetische Fluß-dichte, J die Leitungsstromdichte sowie ρ die elektrische Ladungsdichte.Gleichung (1.3a) ist das Faradaysche Induktionsgesetz, (1.3b) die Maxwell-AmpËre Gleichung, (1.3c) das Gaußsche Gesetz und (1.3d) besagt dassdas magnetische Feld stets rotationsfrei ist.
1.5 Die Maxwell-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 19
In linearen und isotropen Materialien gelten die konstitutiven Gleichungen
B = µH , D = �E , J = σE + J i, (1.4)
mit der elektrischen Permittivität �, magnetischen Permeabilität µ und elek-trischer Leitfähigkeit σ als skalare Proportionalitätsfaktoren.
Wir betrachten ein Beispiel aus der geoelektrischen Erkundung, bei demdurch die Messung elektrischer Felder an der Erdoberfläche Rückschlüsseauf die örtliche Verteilung der Leitfähigkeit σ im Untergrund und damit aufdessen Zusammensetzung gezogen werden.
1.5 Die Maxwell-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 20
• In geophysikalischen Anwendungen kann der Verschiebungsstrom Ḋvernachlässigt werden. Ferner kann µ = µ0 als konstant angenommenwerden.
• Um nur mit einem Feld rechnen zu müssen eliminiert man z.B. aus(1.3a) und (1.3b) das magnetische Feld H .
• Das verbleibende Feld wird auf einem Teilgebiet Ω ⊂ R3 berechnet.
• An dessen Rand ∂Ω sind geeignete Randbedingungen zu stellen; eineeinfache RB in diesem Fall ist n×E = 0 , n die äußere Einheitsnormalelängs ∂Ω.
• Zur Zeit t = t0 wird eine Anfangsbedingung E(x , 0) = E0(x ), x ∈ Ω,benötigt.
1.5 Die Maxwell-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 21
Man erhält so die Anfangs-Randwertaufgabe:
µ0σĖ +∇×∇×E = −µ0J̇ i in Ω, (1.5a)n ×E = 0 längs ∂Ω, (1.5b)E(x , 0) = E0(x ). (1.5c)
1.5 Die Maxwell-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 22
2 Variationstheorie
Die FEM wird angewandt auf die sog. Variationsformulierung von Rand-und Anfangswertaufgaben für Differentialgleichungen.
In diesem Abschnitt stellen wir die wichtigsten Existenz- und Eindeutig-keitssätze bereit und beleuchten die Beziehung zwischen klassischer undVariationsformulierung von Differentialgleichungen.
Die Charakterisierung von Funktionen, für die ein Funktional stationär wird,als Lösung von Differentialgleichungen ist das Thema der Variationsrech-nung, einer mathematischen Disziplin, deren Ausgangspunkt üblicherweiseauf das Jahr 1696 datiert wird, in dem Johann Bernoulli das Brachistochro-nenproblem der mathematischen Weltöffentlichkeit stellte.
2 Variationstheorie TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 23
2.1 Einstimmung: das Brachistochronenproblem
JOHANN BERNOULLI legte 1696 folgendes Problem allen ”tieferdenkendenMathematikern“ der Welt vor:a
Gegeben seien ein Punkt A und ein tiefer lie-gender Punkt B. Man bestimme diejenige Amit B verbindende Kurve, entlang der ein rei-bungsfrei gleitender Körper allein unter demEinfluss der Schwerkraft in kürzester Zeit vonA nach B gelangt.
{ LEIBNIZ UND NEWTON
Abb. 3 . . . und die Lösung
Johann Bernoullis ,,Kurve kürzester Fallzeit“,eine Zykloide, die Brachistochrone, wie sie genanntwurde, brachte den Stein ins Rollen, im wahrstenSinne des Wortes. Den Stein? Eine Lawine löste sieaus!
Der Streit beginntLeibniz hatte seiner Lösung der Bernoullischen Auf-gabe hinzugefügt, das Problem der Kurve kürzesterFallzeit sei so recht geeignet, die Vorzüge seiner Dif-ferentialrechnung ins rechte Licht zu setzen. Dennnur in dieser Bewanderte können diese Aufgabe glattlösen. Außer denen traue er nur Newton, Huygensund Hudde eine Lösung zu. Das war nicht falsch,aber Leibniz hätte besser nichts geschrieben, dennauf der anderen Seite des Ärmelkanals fühlte mansich tief beleidigt. Nicolas Fatio, ein in Englandlebender Schweizer, Mitglied der Royal Society, er-öffnete als erster den Angriff auf Leibniz. Das kamnicht von ungefähr. Fatio, ein enger Freund Newtons,war vor Jahren von Leibniz einmal herablassend be-handelt worden, und jetzt sprach ihm Leibniz auch
noch die Fähigkeit ab, das Problem der Brachisto-chrone lösen zu können. Fatio nahm Rache. In seiner,,Linea brevissimi descensus investigatio geometricaduplex“ von 1699 schrieb er: Newton sei der ersteund um mehrere Jahre älteste Erfinder dieses Kal-küls. Ob Leibniz, der zweite Erfinder, von Newtonetwas entlehnt hat, sollen andere beurteilen. ... Aberniemanden, der die Dokumente durchstudiert, wirddas Schweigen des allzu bescheidenen Newton oderLeibnizens vordringliche Geschäftigkeit täuschen“.
Der Prioritätsstreit war entbrannt. Den Frontal-angriff Fatios beantwortete Leibniz ein Jahr später,1700, in den Acta Eruditorium. Er schrieb ,,Newtonhabe in den Principia selber die Unabhängigkeit derLeibnizschen Entdeckungen anerkannt“. Möglicher-weise hätten sich die Wogen wieder geglättet, aberLeibniz beging wieder einen strategischen Fehler.1705 in einer Rezension von Newtons ,,Optik“, eineanonyme Rezension zwar, die aber nur von Leibnizstammen konnte, war zu lesen:
. . . Statt der Leibnizschen Differenzen benutztnun Herr Newton, und hat er immer benutzt, Flu-xionen, welche sich so nahe wie möglich wie diein gleichen kleinstmöglichen Zeitteilchen hervor-gebrachten Vermehrungen der Fluenten verhalten.Er hat davon in seinen mathematischen Prinzipiender Naturlehre und in anderen später veröffentlich-ten Schriften einen eleganten Gebrauch gemacht,wie auch später Honoratus Fabri in seiner SynopsisGeometrica .. .
Insgesamt war es eine wohlwollende Rezension,aber Newton und Fabri in einem Atemzug zu nen-nen, das war zu viel! Jetzt war endgültig Feuer imHaus!
Newtons Zorn und die Royal Society1708 ließ Newton durch den Oxforder Professor JohnKeill, der gerade Mitglied der Royal Society gewor-den war, Leibniz der Fälschung zeihen. Keill: Allediese Dinge folgen aus der .. . berühmten Methodeder Fluxionen, deren erster Erfinder Sir Isaac New-ton war, wie jeder leicht feststellen kann ... die selbeArithmetic wurde dann später von Leibniz in denActa Eruditorum veröffentlicht, der dabei nur denNamen und die Art Weise der Bezeichnung wechselte.
Das war eine bösartige Anschuldigung, undLeibniz beging, verlassen von seinem diploma-tischen Geschick und in Überschätzung seinerPosition, den dritten Fehler. Er rief seine heimli-chen Feinde in England in eigener Sache als Richter
aJohann Bernoulli: Problema ad cujus solutionem Mathematici invitantur. Acta Eruditorum,1696, pp. 269
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 24
Galileo Galilei
GALILEO zeigte 1638, dass ein Körperauf dem Weg von A nach B schnel-ler über den Umweg ADB rutscht alsauf direktem (geraden) Wege AB. Eben-so zeigt er, dass ACDB, ACDEB undACDEFB jeweils zu kürzeren Laufzei-ten führen, und schliesst, dass ein Kreis-bogen die Lösung sei.
A
B
C
D
EF
Falsche Brachistochrone desGALILEO
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 25
Gottfried Wilhelm Leibniz Jacob Bernoulli Leonhard Euler
• 1684, Acta Eruditorum: LEIBNIZens Abhandlung über Differential- undIntegralrechnung erscheint.
• Die Brüder JAKOB UND JOHANN BERNOULLI bauen diese Ideen syste-matisch zu praktischen und universell einsetzbaren Lösungtechnikenaus.
• Später im 18. Jahrhundert entwickelt dann LEONHARD EULER formaleLösungstechniken für Variationsprobleme.
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 26
Mathematisches Modell
−1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5 4−2.5
−2
−1.5
−1
−0.5
0
0.5
A
ByB
xB
∆x
∆y
∆s
∆s =√
∆x2 + ∆y2
=
√1 +
(∆y
∆x
)2∆x
v =∆s
∆t=√
2g√−y
∆x = ∆x(∆t)
∆y = ∆y(∆t)
∆t =1√2g
√1 +
(∆y∆x
)2√−y
∆x
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 27
Aufsummieren aller Teilzeiten ∆t zur Gesamtzeit T gefolgt vom Grenzüber-gang ∆t→ 0 führt auf den Ausdruck
T =1√2g
∫ xB0
√1 + y′(x)2√−y(x)
dx,
der die Laufzeit T in Abhängigkeit von der Bahnkurve y angibt.
Variationsproblem: Bestimme unter allen auf dem Intervall [0, xB ] stetigdifferenzierbaren Funktionen y = y(x) mit Randwerten
y(0) = 0, y(xB) = yB
diejenige(n), welche
T (y) :=
∫ xB0
√1 + y′(x)2√−y(x)
dx
minimieren.
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 28
Das Brachistochronenproblem besitzt folgende allgemeinere Struktur: ge-sucht ist eine zulässige Funktion, (Glattheit, Randbedingungen)
x : [0, tE ]→ Rn, t 7→ x (t),
welche ein Funktional (”Zielfunktional“)
ψ(x ) :=
∫ tE0
L(t,x (t), ẋ (t)) dt
minimiert, mit einer gegebenen Funktion L : R2n+1 → R.
Beim Brachistochronenproblem:
n = 1, x (t) = x(t), L(t, x, ẋ) =
√1 + ẋ2√−x
.
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 29
Lösungsansatz: Auf EULER und LAGRANGE geht folgender Ansatz zurück,der notwendige Bedingungen für das Vorliegen einer Minimallösung liefert:
Zu einer Minimallösung x ?(t) erhalten wir mit einer beliebigen, hinreichendglatten Funktion φ : [0, tE ]→ Rn mit φ(0) = φ(tE) = 0 durch
x�(t) := x?(t) + �φ(t)
Variationen von x ?, welche ebenfalls für das Variationsproblem zulässigsind. Insbesondere entspricht x ? = x0.
Einsetzen in das Zielfunktional ψ führt auf eine Funktion J = J(�) := ψ(x�)die an der Stelle � = 0 ein lokales Minimum besitzt, für die also
0 =dJ(�)
d�
∣∣∣∣�=0
=
∫ tE0
[∂L
∂x(t,x ?(t), ẋ ?(t))− ∂
2L
∂t∂ẋ(t,x ?(t), ẋ ?(t))
]φ(t) dt .
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 30
Dieses Integral kann nur dann für eine hinreichend große Klasse (z.B. C∞)Funktionen verschwinden, wenn auch
∂L
∂x(t,x ?(t), ẋ ?(t))− ∂
2L
∂t∂ẋ(t,x ?(t), ẋ ?(t)) = 0 . (2.1)
Gleichung (2.1) heißt Euler-Lagrange-Gleichung des Variationsproblemsund liefert eine notwendige Bedingung für das Vorliegen einer minimieren-den Funktion x ?(t).
Beim Brachistochronenproblem gilt ∂tL(t, x, ẋ) = 0 (autonomes Problem).
Folge: Für die Hamilton-Funktion
H(t, x, ẋ) := L(t, x, ẋ)− ẋ∂ẋL(t, x, ẋ)
gilt
d
dtH(t, x?(t), ẋ?(t)) ≡ 0, d.h. H(t, x?(t), ẋ?(t)) = const ∀t.
Man nennt H daher auch ein erstes Integral von L.
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 31
Speziell beim Brachistochronenproblem erhält man:
H(t, x, ẋ) =1
√−x√
1 + ẋ2 ≡ C(> 0),
oder äquivalent−x(1 + ẋ2) =: 2r > 0
mit einer neuen Konstanten r.Auflösen nach ẋ unter Beachtung von ẋ(t) < 0 führt auf die AWA
ẋ(t) = −
√2r + x(t)
−x(t), x(0) = 0,
wofür die Methode der Trennung der Veränderlichen auf
t = t(x) = −∫ x
0
√−ξ
2r + ξdξ
führt.
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 32
Die geschickte Substitution
ξ = −2r sin2 ϕ2, ϕ ∈
[0,π
2
], dξ = −2r sin ϕ
2cos
ϕ
2dϕ
ξ = 0⇔ ϕ = 0, ξ = x⇔ ϕ = ϑ
führt schließlich auf die parametrisierte Lösungsdarstellung
t(ϑ) = r(ϑ− sinϑ), x(ϑ) = −r(1− cosϑ)
bzw., nach Rückbenennung auf die alten Variablen x und y, auf
x(ϑ) = r(ϑ− sinϑ), y(ϑ) = −r(1− cosϑ), 0 ≤ ϑ ≤ ϑB .
Der Parameterwert ϑB am Endpunkt lässt sich durch Lösen der Gleichung
yB(ϑB − sinϑB) + xB(1− cosϑB) = 0 (2.2)
bestimmen; diese ergibt sich durch Elimination von r aus
y(ϑB) = yB und x(ϑB) = xB .
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 33
Die Lösungskurve (2.2) ist eine Zykloide, die Bahn eines Punktes auf demRand eines Kreises, der auf der x-Achse in positiver Richtung abrollt.
−1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5 4−3
−2.5
−2
−1.5
−1
−0.5
0
0.5
1
x
y
ϑ
2.1 Einstimmung: das Brachistochronenproblem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 34
2.2 Bilinearformen
Definition 2.1 Sei V ein reeller normierter Vektorraum. Eine Bilinearformist eine Abbildung
a : V × V → R,
welche linear in beiden Argumenten ist. Die Bilinearform a heißt stetig, fallses eine Konstante C gibt mit
|a(u, v)| ≤ C‖u‖‖v‖ ∀u, v ∈ V . (2.3a)
Die Bilinearform a heißt symmetrisch, falls
a(u, v) = a(v, u) ∀u, v ∈ V . (2.3b)
Die Bilinearform a heißt koerziv, falls es eine Konstante α > 0 gibt mit
a(u, u) ≥ α‖u‖2 ∀u ∈ V . (2.3c)
2.2 Bilinearformen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 35
Bemerkungen 2.2(a) Im Fall eines komplexen Vektorraums V fordert man Antilinearität
im zweiten Argument und spricht dann von einer Sesquilinearform.Anstelle der Symmetrie fordert man hier a(u, v) = a(v, u) und sprichtvon einer Hermiteschen Sesquilinearform.
(b) Eine koerzive symmetrische [Hermitesche] Bilinearform [Sesquilinear-form] definiert ein Innenprodukt auf dem reellen [komplexen] Vektor-raum V . Oft wird es das zur Bilinearform a(·, ·) gehörende Energie-Innenprodukt genannt.
(c) Die Eigenschaft der Stetigkeit (2.3a) zieht die Stetigkeit der Bilinearforma(·, ·) als Funktion ihrer beiden Argumente nach sich.
2.2 Bilinearformen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 36
Satz 2.3 Sei a : H ×H → R eine stetige, symmetrische und koerziveBilinearform auf dem reellen Hilbert-Raum H sowie ` : H → R einestetige Linearform. Dann gilt
(a) Die Minimierungsaufgabe
Minimiere J(u) := 12a(u, u)− `(u) unter allen u ∈H (2.4)
besitzt eine eindeutige Lösung.
(b) Die Minimierungsaufgabe (2.4) ist äquivalent zur Variationsaufgabe
Bestimme u ∈H sodass a(u, v) = `(v) ∀v ∈H . (2.5)
2.2 Bilinearformen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 37
2.3 Die Dirichletsche RWA für die Poisson-Gleichung
2.3.1 Die Poissongleichung
Ist k konstant in (1.1), so kann dies in f zusammengefasst werden und es ergibtsich die Poissongleichung
−∆u(x) = f(x), x ∈ Ω. (2.6a)
Den Gebietsrand Γ := ∂Ω zerlegen wir in ΓD und ΓN , Γ = ΓD ∪ ΓN , ΓD ∩ ΓN = ∅und stellen die Randbedingungen
u(x) = g(x) ∀x ∈ ΓD, kurz: u|ΓD = g, (2.6b)∂u
∂n(x) = h(x) ∀x ∈ ΓN , kurz: ∂nu|ΓN = h (2.6c)
mit zwei gegebenen, auf ΓD bzw. ΓN definierten Funktionen g und h. Das Rand-wertproblem (2.6) besitzt – unter geeigneten Voraussetzungen an das Gebiet Ω unddie Daten f, g, h – eine eindeutig bestimmte klassische (d.h. in Ω zweimal stetigdifferenzierbare) Lösung.
2.3 Die Dirichletsche RWA für die Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 38
Im einfachsten Fall Γ = ΓD erhalten wir die sog. Dirichletsche Randwert-aufgabe (RWA) für die Poisson-Gleichung
−∆u = f auf Ω, (2.7a)u = g auf Γ = ∂Ω. (2.7b)
Neben (2.7) betrachten wir die sog. verallgemeinerte Randwertaufgabe
Bestimme u ∈ C2(Ω) mit u = g längs Γ und∫Ω
∇u · ∇v dx =∫
Ω
fv dx ∀v ∈ C∞0 (Ω).(2.8)
Schließlich betrachten wir noch die Minimierungsaufgabe
Minimiere unter allen Funktionen u ∈ C2(Ω) mit u = g längs Γ
das Funktional J(u) :=1
2
∫Ω
|∇u|2 dx−∫
Ω
fu dx.(2.9)
2.3 Die Dirichletsche RWA für die Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 39
Satz 2.4 Seien g ∈ C(Γ) sowie f ∈ C(Ω) gegebene Funktionen. Sei ferneru ∈ C2(Ω). Dann gilt
(a) Löst u die Minimierungsaufgabe (2.9), so löst u auch die Randwertauf-gabe (2.7).
(b) Die Funktion u ist genau dann Lösung der RWA (2.7), wenn u Lösungder verallgemeinerten RWA (2.8) ist.
Bemerkung 2.5 Nach Satz 2.4 löst jede hinreichend glatte Lösung derVariationsaufgaben (2.8) oder (2.9) auch die RWA (2.7).Schwierigkeit: in vielen Anwendungen tritt der Fall auf, dass keine derartglatte Lösung existiert.
2.3 Die Dirichletsche RWA für die Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 40
Beispiel: Die Situation entspricht der bei den Minimierungsaufgaben
f(x) −→min x ∈ [a, b] (2.10a)und f(x) −→min x ∈ [a, b] ∩Q, (2.10b)
wobei −∞ < a < b
-
Finite Elemente I 41
2.4 Verallgemeinerte Ableitungen
Sei wieder Ω ⊂ Rd, d ∈ N offen und nicht leer.
Für u ∈ C1(Ω) und eine ”Testfunktion“ φ ∈ C∞0 (Ω) gilt die partielle Integra-
tionsformel ∫Ω
u ∂jφdx = −∫
Ω
(∂ju)φdx, 1 ≤ j ≤ d.
Mit anderen Worten: für v := ∂ju gilt die Beziehung∫Ω
u∂jφdx = −∫
Ω
vφ dx ∀φ ∈ C∞0 (Ω). (2.11)
Idee: definiere durch (2.11) die partielle Ableitung auch für den Fall, dassu nicht (im klassischen Sinn) nach xj differenzierbar.
2.4 Verallgemeinerte Ableitungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 42
Definition 2.6 Sei Ω ⊂ Rd, d ∈ N, offen und nicht leer. Gilt für zweiFunktionen u, v ∈ L2(Ω) die Beziehung (2.11), so nennen wir v die verall-gemeinerte Ableitung von u nach xj und schreiben, wie beim klassischenAbleitungsbegriff, v = ∂ju.
Satz 2.7 Besitzt u ∈ L2(Ω) die verallgemeinerte Ableitung v = ∂ju, so istdiese bis auf Funktionswerte auf einer Menge vom Maß Null definiert.
Beispiel 2.8 Sei u(x) := |x|, x ∈ (−1, 1). Dann ist
v(x) :=
−1 falls − 1 < x < 0,α falls x = 0,
1 falls 0 < x < 1
für beliebige reelle Werte von α die verallgemeinerte Ableitung von u.
2.4 Verallgemeinerte Ableitungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 43
2.5 Die Sobolev-Räume H1(Ω) und H10 (Ω)
Definition 2.9 Sei Ω ⊂ Rd, d ∈ N, offen und nicht leer. Der Sobolev-RaumH1(Ω) ist die Menge aller Funktionen u ∈ L2(Ω) mit verallgemeinertenpartiellen Ableitungen erster Ordnung
∂ju ∈ L2(Ω) ∀j = 1, . . . , d.
Für u, v ∈ H1(Ω) definieren wir
(u, v)1 :=
∫Ω
(uv +∇u · ∇v) dx (2.12)
sowie ‖u‖1 := (u, u)1/21 .
Satz 2.10 Durch (2.12) wird auf dem Sobolev-Raum H1(Ω) ein Innenpro-dukt definiert. Mit diesem versehen ist H1(Ω) ein Hilbert-Raum, sofern manwie bei L2(Ω) Funktionen identifiziert, welche fast überall übereinstimmen.
2.5 Die Sobolev-Räume H1(Ω) und H10(Ω) TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 44
Definition 2.11 Sei Ω ⊂ Rd, d ∈ N, offen und nicht leer. Mit H10 (Ω) be-zeichnen wir den Abschluss der Menge C∞0 (Ω) bezüglich ‖ · ‖1.
Satz 2.12 H10 (Ω) ist ein reeller Hilbert-Raum.
Beispiel 2.13 Sei −∞ < a < b < ∞. Ist u ∈ H10 (a, b), so existiert eineeindeutig bestimmte Funktion v ∈ C[a, b] mit u(x) = v(x) fast überall in(a, b) sowie v(a) = v(b) = 0. Ferner gilt die Abschätzung
‖v‖∞ := maxa≤x≤b
|v(x)| ≤ (b− a)1/2(∫ b
a
(u′(x)2 dx
)1/2≤ (b− a)1/2‖u‖1.
2.5 Die Sobolev-Räume H1(Ω) und H10(Ω) TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 45
Beispiel 2.14 Sei −∞ < a < b < ∞ und die Funktion u ∈ C[a, b]stückweise stetig differenzierbar. Es bezeichne K die Menge der Punk-te, in denen u (im klassischen Sinne) differenzierbar ist sowie w : [a, b]→ Rmit
w(x) =
{u′(x) falls x ∈ K,beliebig sonst.
Genauer gesagt: wir nehmen an, dass
(i) u ∈ C[a, b](ii) Es existieren endlich viele Punkte xj , a = x0 < x1 < · · · < xn = b
sodass u auf allen Intervallen (xj , xj+1) stetig differenzierbar ist undu′ stetig auf jeweils [xj , xj+1] fortgesetzt werden kann.
Dann gilt:
(a) w ist die verallgemeinerte Ableitung von u.(b) u ∈ H1(a, b).(c) u ∈ H10 (a, b)⇔ u(a) = u(b) = 0.
2.5 Die Sobolev-Räume H1(Ω) und H10(Ω) TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 46
Verallgemeinerte Randwerte
Definition 2.15 Sei Ω ⊂ Rd, d ∈ N, offen, nicht leer und beschränkt. BeiFunktionen u ∈ H10 (Ω) sagen wir, sie erfüllen die Randbedingung
u = 0 längs ∂Ω (2.13)
im verallgemeinerten Sinne.
Bemerkungen 2.16(a) Beachte: C∞0 (Ω) liegt dicht in H10 (Ω).
(b) Für d = 1 gilt (2.13) sogar im klassischen Sinne (vgl. Beispiel 2.13).
(c) Für d > 1 gilt bei hinreichend glattem Rand Γ = ∂Ω
‖u‖L2(Γ) ≤ C‖u‖1 ∀u ∈ H1(Ω). (2.14)
2.5 Die Sobolev-Räume H1(Ω) und H10(Ω) TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 47
2.6 Die Poincaré-Friedrichs-Ungleichung
Mit Blick auf den abstrakten Lösungssatz 2.3 stellt die Poincaré-Friedrichs-Ungleichung die Koerzivität der Bilinearform a(·, ·) sicher.
Satz 2.17 Sei Ω ⊂ Rd, d ∈ N, offen, nicht leer und beschränkt. Dann exi-stiert eine Konstante C > 0, sodass die Poincaré-Friedrichs-Ungleichung∫
Ω
|u(x)|2 dx ≤ C∫
Ω
|∇u(x)|2 dx ∀u ∈ H10 (Ω) (2.15)
erfüllt ist.
2.6 Die Poincaré-Friedrichs-Ungleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 48
2.7 Ein Existenzsatz für das Dirichlet-Problem
Wir betrachten nach den vorangehenden
-
Finite Elemente I 49
Satz 2.18 (Dirichlet-Prinzip) Sei Ω ⊂ Rd, d ∈ N, offen, nicht leer sowie be-schränkt. Seien ferner f ∈ L2(Ω) sowie g ∈ H1(Ω) gegebene Funktionen.Dann gelten
(a) Die verallgemeinerte Minimierungsaufgabe (2.17) besitzt eine eindeu-tige Lösung u ∈ H1(Ω).
(b) Diese ist auch die eindeutige Lösung der verallgemeinerten Randwert-aufgabe (2.16).
2.7 Ein Existenzsatz für das Dirichlet-Problem TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 50
2.8 Weitere Existenzsätze
Folgender Satz sichert Existenz und Eindeutigkeit der Lösung linearerVariationsaufgaben unter der entscheidenden Annahme der Koerzivität.
Satz 2.19 (Lax-Milgram Lemma, 1954) Sei H ein Hilbert-Raum mit Norm‖ · ‖H , a : H ×H → R eine Bilinearform auf H sowie ` : H → R einlineares Funktional auf H für die es Konstanten C,α > 0 und L gibt mit
|a(u, v)| ≤ C‖u‖H ‖v‖H ∀u, v ∈H , (” a ist stetig “)a(v, v) ≥ α‖v‖2H ∀v ∈H , (” a ist koerziv “)|`(v)| ≤ L‖v‖H ∀v ∈H , (” ` ist stetig “).
Dann besitzt das Variationsproblem
Bestimme u ∈H sodass a(u, v) = `(v) ∀v ∈H
genau eine Lösung.
2.8 Weitere Existenzsätze TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 51
Bemerkungen 2.20(a) Für die Lösung u aus Satz 2.19 gilt die Stabilitätsabschätzung
‖u‖ ≤ 1α‖`‖H ′ .
(b) Ist a(·, ·) zusätzlich noch symmetrisch [Hermitesch], so folgen Existenzund Eindeutigkeit der Lösung aus dem Rieszschen DarstellungssatzA.32.
2.8 Weitere Existenzsätze TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 52
Satz 2.21 (Nečas, 1968; Babuška, 1972) Seien X , Y Banach-Räume und a :X × Y → R eine Bilinearform mit den drei Eigenschaften
(a) a(·, ·) ist stetig, d.h. es existiert eine Konstante C > 0 sodass
|a(x, y)| ≤ C ‖x‖X ‖y‖Y (2.18)
(b) Es existiert eine Konstante γ > 0 sodass
inf0 6=x∈X
sup0 6=y∈Y
|a(x, y)|‖x‖X ‖y‖Y
≥ γ > 0. (2.19)
(c) Es giltsupx∈X
|a(x, y)| > 0 ∀0 6= y ∈ Y . (2.20)
Dann existiert für jedes ` ∈ Y ′ eine eindeutige Lösung der Variationsaufgabe
Bestimme u ∈ X mit a(u, v) = `(v) ∀v ∈ Y . (2.21)
Die Lösung hängt stetig von ` ab, es gilt
‖u‖X ≤1
γ‖`‖Y ′ . (2.22)
2.8 Weitere Existenzsätze TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 53
Bemerkung 2.22 Das Lax-Milgram Lemma ist ein Korollar von Satz 2.21
Folgender Zusatz zu Satz 2.21 ist manchmal nützlich:
Satz 2.23 Werden in Satz 2.21 lediglich (2.18) und (2.19) vorausgesetzt,so ist die Abbildung
A : X → {y ∈ Y : a(x, y) = 0 ∀x ∈X }⊥ ⊂ Y ′
ein Isomorphismus. Ferner ist (2.19) äquivalent zu
‖Ax‖Y ′ ≥ γ‖x‖X ∀x ∈X .
2.8 Weitere Existenzsätze TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 54
2.9 Euler-Lagrange-Gleichungen
Auf Euler und Lagrange geht eine Methode zurück, Funktionen, welche eingegebenes Funktional stationär machen, als Lösung einer Differentialglei-chung zu charakterisieren.
Um zumindest das formale Vorgehen zu beschreiben betrachten wir Funk-tionale der Bauart
J(u) =
∫Ω
F (x, u,∇u) dx
mit einem Gebiet Ω ⊂ Rd und einer hinreichend glatten Funktion
F : Ω× R× Rd → R, (x, u,p) 7→ F (x, u,p).
Eine Funktion u ∈ C2(Ω) ist stationärer ”Punkt“ von J gegenüber glattenStörungen φ ∈ C∞0 (Ω) mit kompaktem Träger, wenn
d
dtJ(u+ tφ)
∣∣∣∣t=0
= 0 ∀φ ∈ C∞0 (Ω).
2.9 Euler-Lagrange-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 55
Differentiation im Integral∫Ω
d
dtF (x, u+ tφ,∇u+ t∇φ)
∣∣t=0
dx
unter Anwendung der Kettenregel führt auf∫Ω
(Fu(x, u,∇u)φ+
d∑j=1
Fpj (x, u,∇u)φxj (x)
)dx = 0.
Partielle Integration im Summenterm ergibt schließlich unter Beachtung von φ ≡ 0auf ∂Ω∫
Ω
(Fu(x, u,∇u)−
d∑j=1
∂xj [Fpj (x, u,∇u)])φdx = 0, ∀φ ∈ C∞0 (Ω),
was nach dem Variationslemma auf die sog. Euler-Lagrange-Gleichung
Fu(x, u,∇u)−d∑
j=1
∂xj [Fpj (x, u,∇u)] = 0, x ∈ Ω, (2.23)
des Funktionals J führt.
2.9 Euler-Lagrange-Gleichungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 56
3 Finite-Elemente Approximationen zur Lösungder Poisson-Gleichung in 2D
3 Finite-Elemente Approximationen zur Lösung der Poisson-Gleichung in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 57
3.1 Aufgabenstellung
Sei Ω ⊂ R2 ein beschränktes polygonales Gebiet mit Rand Γ bestehend aus denTeilmengen ΓD und ΓN sowie g und h hinreichend glatte Randfunktionen. ZurApproximation der Lösung der RWA
−∆u = f in Ω, (3.1a)
u = g längs ΓD, (3.1b)
∂u
∂n= h längs ΓN (3.1c)
betrachten wir die zugehörige schwache Formulierung
Bestimme u ∈ Vg sodass a(u, v) = `(v) ∀v ∈ V0. (3.2a)
Hierbei sind Vg := {u ∈ H1(Ω) : u|ΓD = g} sowie
a(u, v) =
∫Ω
∇u · ∇v dx, `(v) =∫
Ω
fv dx+
∫ΓN
hv ds. (3.2b)
3.1 Aufgabenstellung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 58
Bemerkung 3.1 Oft wird die Variationsaufgabe (3.2) ”homogenisiert“, umVg = V0 zu erreichen.
Dies geschieht mit Hilfe einer beliebigen Funktion ug ∈ H1(Ω) mit derEigenschaft ug|ΓD = g. Dann ist u0 + ug ∈ Vg für alle u0 ∈ V0.
Die homogenisierte Variante von (3.2) lautet somit
Bestimme u0 ∈ V0 sodass a(u0, v) = `(v)− a(ug, v) ∀v ∈ V0, (3.3)
d.h. auf der rechten Seite geht man über zur Linearform ˜̀ := `− a(ug, ·).
3.1 Aufgabenstellung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 59
3.2 Referenzaufgaben
Auf folgende Beispielaufgaben kommen wir im Laufe dieses Abschnittswieder zurück.
3.2 Referenzaufgaben TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 60
Beispiel 1: Ω = (−1, 1)2, ΓD = Γ, f ≡ 1, g ≡ 0.
Dies stellt ein Modell für Wärmeausbreitung auf einer quadratischen Platte dar.Durch Trennung der Veränderlichen bestimmt man die Reihenlösung
u(x, y) =1− x2
2− 16π3
∑k∈N,k ungerade
[sin(kπ(1 + x)/2)
k3 sinh(kπ)(sinh
kπ(1 + y)
2+ sinh
kπ(1− y)2
)].
−1
−0.5
0
0.5
1
−1
−0.5
0
0.5
1−0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
3.2 Referenzaufgaben TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 61
Beispiel 2: Ω = (−1, 1)2 \ [−1, 0]2, ΓD = Γ, f ≡ 1, g ≡ 0.
−1
−0.5
0
0.5
1
−1
−0.5
0
0.5
10
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
3.2 Referenzaufgaben TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 62
Beispiel 3: Ω = (−1, 1)2, ΓD = Γ, f ≡ 0, g = u|Γmit exakter Lösung
u(x, y) =2(1 + y)
(3 + x)2 + (1 + y)2.
−1
−0.5
0
0.5
1
−1
−0.5
0
0.5
10
0.1
0.2
0.3
0.4
0.5
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
3.2 Referenzaufgaben TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 63
Beispiel 4: Ω = (−1, 1)2 \ [−1, 0]2, ΓD = Γ, f ≡ 1, g = u|Γ mit exakterLösung
u(r, θ) = r2/3 sin2θ + π
3, x = r cos θ, y = r sin θ.
−1
−0.5
0
0.5
1
−1
−0.5
0
0.5
10
0.2
0.4
0.6
0.8
1
1.2
1.4
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
3.2 Referenzaufgaben TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 64
3.3 Galerkin-Diskretisierung
Das Galerkin-Verfahren zur Approximation der Lösung von (3.2) bestehtdarin, Ansatz- und Testraum durch endlichdimensionale Räume V hg bzw.V h0 (gleicher Dimension) zu ersetzen.
Hierbei ist h > 0 ein Diskretisierungsparameter.
Gelten V hg ⊂ Vg und V h0 ⊂ V0, so spricht man von einer konformenDiskretisierung.
Werden verschiedene Ansatz- und Testräume verwendet, spricht man vonPetrov-Galerkin-Verfahren, andernfalls von Bubnov-Galerkin-Verfahren.
Wir betrachten zunächst die homogenisierte Variationsaufgabe (3.3) undeinen n-dimensionalen Unterraum V h0 ⊂ V0. (Oder, äquivalent, den Fallg ≡ 0.)Die diskrete Variationsaufgabe lautet somit
Bestimme uh ∈ V h0 sodass a(uh, v) = `(v) ∀v ∈ V h0 . (3.4)
3.3 Galerkin-Diskretisierung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 65
Ist {φ1, φ2, . . . , φn} eine Basis von V h0 sowie uh =∑nj=1 ujφj , so ist (3.4)
äquivalent mit
n∑j=1
uj a(φj , φi) = `(φi), i = 1, 2, . . . , n,
oder, mit A ∈ Rn×n gegeben durch [A]i,j = a(φj , φi), b ∈ Rn durch[b]i = `(φi) sowie u ∈ Rn durch [u ]i = ui, zu dem Galerkin-System
Au = b. (3.5)
Beachte:
• Die diskrete Variationsaufgabe (3.4) bzw. (3.5) besitzt ein eindeutigeLösung.
• Die Galerkin-Matrix A aus (3.5) ist symmetrisch und positiv-definit.
3.3 Galerkin-Diskretisierung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 66
3.4 Finite-Element-Räume in 2D
Galerkin-Verfahren lassen zunächst beliebige Unterräume zu (Spektralver-fahren, Integralgleichungen, etc.).
Die finite-Element-Methode (FEM) zeichnet sich durch die Verwendung vonBasisfunktionen mit kleinem Träger aus. Fast ausschließlich werden hierfürstückweise Polynome eingesetzt.
Deren Konstruktion basiert auf einer Zerlegung von Ω in (disjunkte) Teilge-biete, Elemente genannt, die wir mit K bezeichnen.
Wir betrachten in diesem Abschnitt (Ω ist ein Polygon) Zerlegungen ausDreiecken bzw. Vierecken, welche (in beiden Fällen, auch in 3D) Triangu-lierungen heißen.
Weitere gebräuchliche Bezeichnungen sind Vernetzung, Netz oder Gitter.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 67
Mit V h sei ein zunächst beliebiger endlichdimensionaler Raum auf Ω defi-nierter Funktionen bezeichnet. Entsprechend sei
PK := {v|K : v ∈ V h}
der durch sämtliche Einschränkungen von Funktionen aus V h auf K auf-gespannte Raum.
Bei einer konformen FE-Diskretisierung ist V h ⊂ V erforderlich. EineCharakterisierung von Konformität liefert folgender Satz.
Satz 3.2 Sei T h eine Zerlegung des Gebietes Ω und V h ein endlichdi-mensionaler Funktionenraum. Gilt V h ⊂ C0(Ω) sowie PK ⊂ H1(K) für alleK ∈ T h, so gelten
V h ⊂ H1(Ω), sowie {v ∈ V h : v = 0 auf ∂Ω} ⊂ H10 (Ω).
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 68
Beweis: Aufgrund unserer Annahmen gilt bereits V h ⊂ L2(Ω). Damit auch V h ⊂H1(Ω) müssen wir zeigen, dass jedes v ∈ V h schwache Ableitungen ∂iv, i =1, . . . , d besitzt, d.h. Funktionen vi ∈ L2(Ω) mit∫
Ω
vi φdx = −∫
Ω
v ∂iφ dx ∀φ, φ differenzierbar, φ|∂Ω = 0.
Elementweise gilt∫K
φ∂i(v|K) dx = −∫K
v|K ∂iφdx +
∫∂K
v|K φnK,i ds,
(nK,i die i-te Komponente der äußeren Einheitsnormalen längs ∂K). Summationüber alle K ergibt (mit vi := ∂iv|K∀K)∫
Ω
φ vi dx = −∫
Ω
v ∂iφdx +∑
K∈Th
∫∂K
v|K φnK,i ds.
Die Summe verschwindet jedoch, da φ längs ∂Ω verschwindet und die (orientierten!)Randintegrale längs innerer Ränder je zweimal mit entgegengesetztem Umlaufsinnauftreten. �
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 69
3.4.1 Dreieckelemente
Wir betrachten eine Zerlegung T h von Ω in (abgeschlossene) Dreiecke Kmit folgenden Eigenschaften
(a) Ω = ∪K∈T hK.
(b) Für zwei beliebige Dreiecke K1,K2 ∈ T h ist K1 ∩ K2 entweder leer,ein gemeinsamer Punkt oder eine gemeinsame Kante.
Als Diskretisierungsparameter wird üblicherweise
h := maxK∈T h
diamK,
also ein Maß für die Feinheit der Zerlegung, genommen.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 70
Beispiel 1: Dreieckszerlegung des Äußeren eines Tragflächenquerschnitts.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 71
Besipiel 2: Triangulierungen verschiedener Gebiete(mittels der distmesh-Software von Strang und Persson)
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 72
Beispiel 3: Zerlegung aus Drei- und Vierecken.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 73
Um aus stückweisen Funktionen auf den Teilmengen K von T h einenUnterraum von H1(Ω) zu erhalten, müssen die hieraus konstruierten Funk-tionen nach Satz 3.2 stetig sein. Für stückweise Polynome bezüglich T h
ist dies gewährleistet, sofern diese nur an allen Dreieckskanten stetig sind.
Der einfachste solche Raum ist der aller stückweise linearen Funktionenauf Ω bezüglich T h, d.h.
V h := {v ∈ C(Ω) : v|K ∈P1 ∀K ∈ T h},
wobei mit P1 die Menge alle Polynome (hier in zwei Variablen) vom Gradhöchstens eins bezeichnet sei.
Ein Unterraum von V0 ist gegeben durch
V h0 = {v ∈ V h : v|ΓD = 0}.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 74
Nodale Basis
Funktionen aus V h sind eindeutig durch ihre Funktionswerte an den Knotender Triangulierung (Ecken der Dreiecke) bestimmt. Einen solchen Satz vonParametern, welche eine FE-Funktion vollständig charakterisieren, nenntman Freiheitsgrade (degrees of freedom).
Bei V h0 sind dies alle Knoten bis auf die, die nicht auf ΓD liegen. DerenAnzahl sei n.
Eine besonders nützliche Basis {φ1, . . . , φn}, die sog. nodale Basis, istgegeben durch
φj(xi) = δi,j i, j = 1, . . . , n.
Ist N h = {x1, . . . , xn} die Menge der Knoten mit xj 6∈ ΓD, so gilt
suppφj = {K ∈ T h : xj ∈ K}.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 75
Triangulierung des L-förmigen Gebiets aus Abschnitt 2.1 mit dem Träger dreiernodaler Basisfunktionen.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 76
Konsequenz für das Galerkin-System (3.5):
[b]i = `(φi) =
∫suppφi
fφi dx+
∫suppφi∩ΓN
hφi ds
[A]i,j = a(φj , φi) =
∫suppφi∩suppφj
∇φi · ∇φj dx
Insbesondere: die Galerkin-Matrix A ist dünn besetzt.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 77
Aufbau des Galerkin-Systems
Großer Vorteil der FEM:
• Aufbau des FE Gleichungssystems – man nennt diesen Vorgang As-semblieren – verläuft auch bei komplizierteren Problemen stets nachdem gleichen Grundschema.
• Parametrisierung durch spezielle Eigenschaften der jeweiligen Aufga-be.
• Grundbausteine der Assemblierung stets dieselben, kann bei der Soft-wareumsetzung ausgenutzt werden (z.B. Objektorientierung).
Wir stellen nun die Grundschritte der Assemblierung für unser Modellpro-blem zusammen. Das Vorgehen besteht darin, die Rechnung auf Operatio-nen auf den einzelnen Elementen zurückzuführen, deren Ergebnisse dannzusammengefügt (assembliert) werden.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 78
-
Finite Elemente I 79
Sei K ∈ T h: dann gilt für i, j = 1, 2 . . . , ñ:
a(φj , φi) =
∫Ω
∇φj · ∇φi dx =∑
K∈T h
∫K
∇φj · ∇φi dx =:∑
K∈T haK(φj , φi),
`(φi) =∑
K∈T h
∫K
fφi dx+
∫K∩ΓN
hφi ds =: `K(φi).
Mit
[ÃK ]i,j := aK(φj , φi) i, j = 1, 2, . . . , ñ,
[b̃K ]i := `K(φi, i = 1, 2, . . . , ñ,
folgt alsoà =
∑K∈T h
ÃK , b̃ =∑
K∈T hb̃K .
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 80
Da jedes Element nur zum Träger dreier Basisfunktionen gehört, sind nur(höchstens) neun bzw. drei Einträge von ÃK bzw. b̃K von Null verschieden.
Welche Einträge dies sind kann durch Nachschlagen in einer Elementta-belle ermittelt werden:
[ET (i, j)]i=1,2,3;j=1,...,nK :
Element K1 K2 . . . KnKerster Knoten i(1)1 i
(2)1 . . . i
(nK)1
zweiter Knoten i(1)2 i(2)2 . . . i
(nK)2
dritter Knoten i(1)3 i(2)3 . . . i
(nK)3
Hierbei sei nK die Anzahl der Elemente in T h.
Diese führt neben der bisherigen globalen Knotennumerierung x1, x2, . . . , xñin jedem Element zusätzlich eine lokale Nummerierung x(K)1 , x
(K)2 , x
(K)3 der
zu K gehörenden Knoten ein.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 81
Globale Nummerierung der Knoten (rot) und der Elemente (schwarz)einer Triangulierung des L-Gebiets.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 82
Damit sind die von Null verschiedenen Teilmatrizen bzw. -vektoren von AKund bK gegeben durch
AK :=
aK(φ
(K)1 , φ
(K)1 ) aK(φ
(K)2 , φ
(K)1 ) aK(φ
(K)3 , φ
(K)1 )
aK(φ(K)1 , φ
(K)2 ) aK(φ
(K)2 , φ
(K)2 ) aK(φ
(K)3 , φ
(K)2 )
aK(φ(K)1 , φ
(K)3 ) aK(φ
(K)2 , φ
(K)3 ) aK(φ
(K)3 , φ
(K)3 )
, bK :=`K(φ
(K)1 )
`K(φ(K)2 )
`K(φ(K)3 )
.
Trägt K in der Nummerierung der Elemente den Index k, so ist die Zu-ordnung der lokalen Nummerierung {φ(K)i }i=1,2,3 der zu K gehörendenBasisfunktionen zur globalen Nummerierung {φj}ñj=1 gegeben durch
φ(K)i = φj , j = ET (i, k), i = 1, 2, 3.
Man bezeichnet AK und bK auch als Elementsteifigkeitsmatrix bzw. Elem-entlastvektor.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 83
Nach diesen Überlegungen erhalten wir folgenden Algorithmusa zur As-semblierung:
Initialisiere à := O, b̃ := 0foreach K ∈ Th
berechne AK und bKk := [Index des Elementes K]i1 := ET (1, k), i2 := ET (2, k), i3 := ET (3, k)Ã([i1i2i3], [i1i2i3]) := Ã([i1i2i3], [i1i2i3]) + AKb̃([i1i2i3]) := b̃([i1i2i3]) + bK
end
aHier wird folgende an MATLAB angelehnte Notation verwendet:
A([i1i2i3], [i1i2i3]) =
ai1,i1 ai1,i2 ai1,i3ai2,i1 ai2,i2 ai2,i3ai3,i1 ai3,i2 ai3,i3
, b([i1i2i3]) =bi1bi2bi3
.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 84
Referenzelement
Hilfreich für Implementierung und Analyse: Bezug auf ein ReferenzelementK̂ ⊂ R2. Für alle K ∈ T h gilt dann K = FK(K̂) mit
FK : K̂ → K, K̂ 3 ξ 7→ x ∈ K, x = FK(ξ) = BKξ + bK .
Bei Dreieckelementen üblich: Einheitssimplex
K̂ = {(ξ, η) ∈ R2 : 0 ≤ ξ ≤ 1, 0 ≤ η ≤ 1− ξ}
Für jedes Dreieck K ∈ T h ist die affine Abbildung FK bestimmt durch dieAbbildungsvorschriften
(1, 0) 7→ (x1, y1),(0, 1) 7→ (x2, y2),(0, 0) 7→ (x3, y3), d.h.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 85
K̂
ξ
η
(0, 0) (1, 0)
(0, 1)
x
y
FK
K
(x1, y1)
(x2, y2)
(x3, y3)
[x
y
]=
[x1 − x3 x2 − x3y1 − y3 y2 − y3
]︸ ︷︷ ︸
BK
[ξ
η
]+
[x3
y3
]︸ ︷︷ ︸bK
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 86
Lokale (nodale) Basisfunktionen in K̂:
φ̂1(ξ, η) = ξ, φ̂2(ξ, η) = η, φ̂3(ξ, η) = 1− ξ − η, (ξ, η) ∈ K̂.
Durch die Zuordnung
φ̂ 7→ φ := φ̂ ◦ F−1K , d.h. φ(x ) := φ̂(ξ(x )) = φ̂(F−1K (x ))
wird jeder Funktion φ̂ auf K̂ eine Funktion φ auf K zugeordnet.
Lokale Basisfunktionen auf K:
φj = φ̂j ◦ F−1K : K → R, j = 1, 2, 3.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 87
Rückführung der Integration auf das Referenzelement
Die bei der Assemblierung der Galerkin-Matrix anfallenden Integrale wer-den ebenfalls auf das Referenzelement zurückgeführt. Dies ist auch für dieVerwendung von Quadraturformeln hilfreich.
Aus φ(x ) = φ̂(ξ(x )) folgt (Kettenregela)
∇φ =
[φx
φy
]=
[φ̂ξξx + φ̂ηηx
φ̂ξξy + φ̂ηηy
]=
[ξx ηx
ξy ηy
][φ̂ξ
φ̂η
]= (DF−1K )
>∇̂φ̂.
Wegen x = FK(ξ) = BKξ + bK , d.h. DFK = BK ,
ξ = F−1K (x ) = B−1K (x − bK), d.h. DF
−1K = B
−1K
folgt schließlich∇φ = B−>K ∇̂φ̂.
a∇̂ bedeutet, dass nach den Variablen ξ und η differenziert wird.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 88
Damit ergibt sich (φi = φ(K)i , i = 1, 2, 3)
aK(φj , φi) =
∫K
∇φj · ∇φi dx
=
∫K̂
(B−>K ∇̂φ̂j
)·(B−>K ∇̂φ̂i
)|detBK | dξ.
(3.6)
Hierbei ist (lineares Dreieckelement)
|detBK | = 2|K|,
B−>K =1
2|K|
[y2 − y3 x3 − x2y3 − y1 x1 − x3
],
[∇̂φ̂1 ∇̂φ̂2 ∇̂φ̂3
]=
[1 0 −10 1 −1
].
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 89
Dreieckelemente vom Grad zwei
Nimmt man zu den Ecken des Referenzdreiecks K̂ noch die Kantenmittel-punkte als Knoten hinzu, so lautet die zugehörige nodale Basis von P2 aufK̂:
φ̂1(ξ, η) = ξ(2ξ − 1),
φ̂4(ξ, η) = 4η(1− ξ − η),
φ̂2(ξ, η) = η(2η − 1),
φ̂5(ξ, η) = 4ξ(1− ξ − η),
φ̂3(ξ, η) = (1− (ξ + η))(1− 2(ξ + η)),
φ̂6(ξ, η) = 4ξη
x1
x2
x3
x4
x5
x6
-
Finite Elemente I 90
Die vier nodalen Basisfunktionen im Dreieckelement vom Grad zwei
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 91
3.4.2 Viereckelemente
Auch wenn Dreieckzerlegungen im Allgemeinen flexibler sind, so treten inder Praxis hinreichend oft Gebiete auf, welche einfach in Rechtecke oder(konvexe) Vierecke zerlegbar sind.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 92
Beim einfachsten H1-konformen Viereckelement wird als Polynomraum
Q1 := span{1, ξ, η, ξη},
der Raum der bilinearen Funktionen verwendet.
Im Referenzelement, üblicherweise K̂ = [−1, 1]2, lautet die nodale Basis
φ1(ξ, η) =14 (1− ξ)(1− η),
φ2(ξ, η) =14 (1 + ξ)(1− η),
φ3(ξ, η) =14 (1 + ξ)(1 + η),
φ4(ξ, η) =14 (1− ξ)(1 + η).
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 93
Beim quadratischen Viereckelement ist der Polynomraum
Q2 := span{1, ξ, η, ξη, ξ2, ξ2η, ξη2, η2},
der der biquadratischen Funktionen. Als neue Freiheitsgrade gegenüber dem bi-linearen Viereck kommen hier die Funktionswerte an den Seitenmitten sowie amMittelpunkt hinzu. Die neun Basisfunktionen auf dem Referenzelement lauten hier
φ1(ξ, η) =14ξ(1− ξ)η(1− η),
φ2(ξ, η) =14ξ(1 + ξ)η(1− η),
φ3(ξ, η) =14ξ(1 + ξ)η(1 + η),
φ4(ξ, η) =14ξ(1− ξ)η(1 + η),
φ5(ξ, η) =12(1 + ξ)(1− ξ)η(1− η),
φ6(ξ, η) =12ξ(1 + ξ)(1− η)(1 + η),
φ7(ξ, η) =12(1− ξ)(1 + ξ)η(1 + η),
φ8(ξ, η) =12ξ(1− ξ)(1− η)(1 + η),
φ9(ξ, η) = (1− ξ)(1 + ξ)(1− η)(1 + η).
x1 x2
x3x4
x5
x6
x7
x8 x9
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 94
Biquadratische nodale Basisfunktionen zu einer Ecke, einer Seitenmitte und demMittelpunkt.
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 95
Bei einem achsenparallelen RechteckK mit Kantenlängen hx und hy sowielinkem unteren Eckpunkt (x1, y1) und den weiteren Ecken im Gegenuhrzei-gersinn durchnummeriert lautet die Gebietsabbildung
FK(x ) = BKξ + bK , BK =1
2
[hx 0
0 hy
], b =
1
2
[x1 + x3
y1 + y3
].
In diesem Fall ist (vgl. (3.6))
|detBK | = 14hxhy,
B−>K = 2
[1/hx 0
0 1/hy
].
3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 96
3.5 Approximationsfehler
Wir betrachten zunächst allgemeine konforme Galerkin-Diskretisierungeneiner Variationsaufgabe
Bestimme u ∈ V sodass a(u, v) = `(v) ∀v ∈ V , (3.7)
wobei wir annehmen, dass die Voraussetzungen des Lax-Milgram-Lemmas(Satz 2.19) erfüllt sind. Die diskrete Variationsaufgabe lautet entsprechend
Bestimme uh ∈ V h sodass a(uh, v) = `(v) ∀v ∈ V h, (3.8)
mit einem Unterraum V h ⊂ V .
Frage: Was kann man aussagen über ‖u− uh‖, insbesondere für h→ 0?
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 97
3.5.1 Das Céa-Lemma
Satz 3.3 (Lemma von Céa) Gelten für die Variationsaufgabe (3.7) die Vor-aussetzungen des Lax-Milgram-Lemmas, so gilt für den Fehler u − uh derGalerkin-Approximation
‖u− uh‖ ≤ Cα
infv∈V h
‖u− v‖. (3.9)
Hierbei bezeichnen C und α die Stetigkeits- bzw. Koerzivitätskonstante ausdem Lax-Milgram-Lemma.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 98
Bemerkungen 3.4(a) Ist die Bilinearform zusätzlich noch symmetrisch (Hermitesch), so läßt
sich (3.9) verbessern zu
‖u− uh‖ ≤√C
αinfv∈V h
‖u− v‖. (3.10)
(b) Für die Konvergenz einer Folge von Galerkin-Approximationen (imGrenzwert h→ 0) ist somit hinreichend, dass für die zugehörige Familie{V h}h>0 von Unterräumen gilt
limh→0
infv∈V h
‖u− v‖ = 0 ∀u ∈ V .
(c) In diesem Sinne wird durch das Céa-Lemma die Abschätzung desGalerkin-Fehlers zu einem Approximationsproblem.
(d) Die Tatsache, dass a(u− uh, v) = 0 für alle v ∈ V h wird auch Galerkin-Orthogonalität genannt.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 99
Das Céa-Lemma läßt sich auch auf den Fall verallgemeinern, dass kei-ne Koerzivität, sondern lediglich eine inf-sup-Bedingung vorliegt (vgl.Satz 2.21).
Satz 3.5 Seien X , Y Banach-Räume und a : X × Y → R eine Bilinear-form, welche die Voraussetzungen von Satz 2.21 erfüllt. Ferner seien fürUnterräume X h ⊂ X und Y h ⊂ Y gleicher Dimension die Bedingung(2.19) mit einer Konstanten γh sowie Bedingung (2.20) erfüllt. Dann gilt fürdie Lösung uh der diskreten Variationsaufgabe
Bestimme uh ∈X h sodass a(uh, v) = `(v) ∀v ∈ Y h, (3.11)
und die Lösung u ∈X von (2.21) die Abschätzung
‖u− uh‖X ≤(
1 +C
γh
)inf
v∈X h‖u− v‖X .
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 100
Bemerkung 3.6 Im Gegensatz zum Lax-Milgram Lemma folgen Bedingun-gen (2.19) und (2.20) nicht aus deren Gültigkeit für die Räume X und Y .Insbesondere kann die Konstante γ in (2.19) von h abhängen, d.h. γ = γh.Entscheidend ist, ob γh gleichmäßig nach unten beschränkt ist.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 101
3.5.2 A priori Fehlerabschätzungen
Die einfachsten a priori Fehlerabschätzungen beruhen auf zusätzlicherRegularität:
Definition 3.7 Mit H2(Ω) wird der Raum aller Funktionen in u ∈ L2(Ω)bezeichnet, welche schwache Ableitungen bis einschließlich Ordnung zweibesitzen.
H2(Ω) ist ein Hilbert-Raum bezüglich des Innenproduktes
(u, v)2 := (u, v)1 +2∑
i,k=1
(∂jku, ∂jkv).
Hierzu gehört die Norm
‖u‖2 =(‖u‖21 + |u|22
)1/2, |u|2 :=
( 2∑i,k=1
‖∂jku‖2)1/2
.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 102
Definition 3.8 Die Variationsaufgabe (3.2) heißt H2-regulär falls für hinrei-chend glatte Daten f , g und h die Lösung sogar in Vg ∩H2(Ω) liegt.
Aufgrund des Sobolevschen Einbettungssatzes liegt in der Äquivalenzklassejeder Funktion u ∈ H2(Ω) eine aus C(Ω). Insbesondere kann man bei H2-Funktionen vom Wert u(x ) dieser Funktion an einer Stelle x ∈ Ω sprechen.
Definition 3.9 Sei T h eine zulässige Triangulierung des Gebietes Ω undV h der FE-Raum aus stückweise linearen Funktionen bezüglich T h. Füralle K ∈ T h mit Ecken {xj}3j=1 und u ∈ H2(K) wird die durch
IKu ∈P1(K), (IKu)(xj) = u(xj), j = 1, 2, 3
als lokale Interpolierende zu u bezeichnet. Die globale InterpolierendeIhu ∈ V h zu u ∈ H2(Ω) ist definiert durch
Ihu|K = IKu ∀K ∈ T h.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 103
Lemma 3.10 Für alle u ∈ H2(Ω) und alle K ∈ T h gilt
‖∇(u− IKu)‖2K ≤2h2K|K|‖∇̂(û− IK̂ û)‖
2K̂.
Folgende Aussage ist ein Spezialfall des Bramble-Hilbert-Lemmas, wel-ches wir in allgemeiner Form später beweisen werden:
Lemma 3.11 Es existiert eine Konstante C sodass für alle û ∈ H2(K̂) gilt
‖∇̂(û− IK̂ û)‖K̂ ≤ C|û− IK̂ û|2,K̂ = C|û|2,K̂ .
Lemma 3.12 Für alle u ∈ H2(K) und û = u ◦ FK gilt
|û|22,K̂≤ 12h2K
h2K|K||u|22,K .
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 104
Lemma 3.13 Bezeichnet θK ∈ (0, π/3] den kleinsten Innenwinkel eines(nichtentarteten) Dreiecks K, so gilt
h2K4
sin θK ≤ |K| ≤h2K2
sin θK .
Korollar 3.14 Für die globale Interpolierende Ih eines stückweise linearenFE-Raumes V h ⊂ H1(Ω) bezüglich einer Triangulierung T h gilt
‖∇(u− Ihu)‖2 ≤ C∑
K∈T h
1
sin2 θKh2K |u|22,K ∀u ∈ H2(Ω).
Definition 3.15 Eine Familie von Triangulierungen {T h}h>0 heißt qua-siuniform, falls eine gleichmäßige untere Schranke θ für den minimalenInnenwinkel aller Dreiecke K ∈ T h für alle h existiert.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 105
Satz 3.16 Ist die Variationsaufgabe (3.2)H2-regulär, so gilt für die Galerkin-Approximation uh bezüglich eines Unterraumes V h ⊂ H1(Ω) aus stückweiselinearen Funktionen bezüglich einer quasiuniformen Triangulierung T h mith := maxK∈T h diamK
‖u− uh‖1 ≤ Ch|u|2mit einer von u und h unabhängigen Konstanten C.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 106
3.5.3 A posteriori Fehlerabschätzungen
Gegeben: RWA (3.2), Lösung u, V h basierend auf T h, u ≈ uh ∈ V h
Gesucht: (lokaler) Fehlerschätzer η : T h → R, K 7→ ηK mit
• ηK ≈ ‖∇(u− uh)‖K ,• Aufwand zur Berechnung von {ηK}K∈T h möglichst linear in |T h|,• obere Schranke der Form∑
K∈T h‖∇(u− uh)‖2K ≤ C(θ)
∑K∈T h
η2K ,
(θ untere Schranke für minimalen Innenwinkel),• untere Schranke der Form
ηK ≤ C(θωK )‖∇(u− uh)‖ωK ,
ωK lokale Elementumgebung von K, K ⊂ ωK ⊂ T h.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 107
Ausgangspunkt: für alle Testfunktionen v ∈ V0 gilt
a(u− uh, v) = `(v)− a(uh, v). (3.12)
Die rechte Seite von (3.12) stellt ein Residuumsterm dar, ein Element desDuals von V0.
Man kann zeigen: Norm von u−uh beidseitig durch Dualnorm des Residu-ums beschränkt.
Idee: Schätze u − uh durch Schätzung von ‖` − a(uh, ·)‖ bzw. durchnäherungsweise Lösung von (3.12).
Annahme: homogene Neumann-Daten, d.h. `(v) = (f, v).
Mit (u, v)K :=∫Kuv dx bzw. aK(u, v) :=
∫K∇u · ∇v dx wird (3.12) zu
a(u− uh, v) =∑
K∈T haK(u− uh, v) =
∑K∈T h
(f, v)K −∑
K∈T haK(u
h, v).
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 108
Partielle Integration:
− aK(uh, v) =∫K
v∆uh dx −∑
E∈E (K)
〈n · ∇uh, v
〉E, (3.13)
mit den Bezeichnungen
E (K) Menge der Kanten des Elements K ∈ T h.n = nE,K äußerer Normalenvektor längs E bez. K,
〈·, ·〉E Innenprodukt in L2(E), (Dualform in H−1/2(E)×H1/2(E))
n · ∇uh (diskreter) Fluss über E, i.A. unstetig an Elementgrenzen.
Definiere daher für zwei Elemente K1 und K2 mit gemeinsamer Kante Eden Sprung des Flusses von v über E durch
s∂v
∂n
{:= (∇v|K1 −∇v|K2) · nE,K1 = (∇v|K2 −∇v|K1) · nE,K2 .
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 109
Mit diesen Bezeichnungen sowie (3.13) wird aus (3.12)
∑K∈T h
aK(u− uh, v) =∑
K∈T h
(f + ∆uh, v)K − 12 ∑E∈E (K)
〈s∂uh
∂n
{, v
〉E
,wobei wir den Fluss zu gleichen Teilen auf die angrenzenden Elementeverteilt haben.
2 Anteile:
Element-Residuum RK := (f + ∆uh)|K
Fluss-Sprung RE :=s∂uh
∂n
{
Beachte: RK ≡ RE ≡ 0 für exakte Lösung u.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 110
Hier: betrachte lineare Dreiecke bzw. bilineare Rechtecke
• RE ≡ konstant bei linearen Dreiecken,
• in beiden Fällen RK = f |K .
Weitere Vereinfachung: approximiere RK ≈ R0K , wobei f orthogonal aufden Raum der konstanten Funktionen projiziert wird.
Weitere Notation:
E h :=⋃
K∈T hE (K) = E hΩ ∪ E hD ∪ E hN ,
E hΩ := Eh ∩ Ω, E hD := E h ∩ ΓD, E hN := E h ∩ ΓN .
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 111
Mit der Bezeichnung
R∗E :=
12
r∂uh
∂n
z, E ∈ E hΩ ,
−nE,K · ∇uh, E ∈ E hN ,0, E ∈ E hD,
erhalten wir schließlich∑K∈T h
aK(u− uh, v) =∑
K∈T h
[(RK , v)K −
∑E∈E (K)
〈R∗E , v〉E
], (3.14)
Exakter Fehler durch (3.14) für alle v ∈ V0 bestimmt.
Idee: bestimme eK ≈ u− uh auf K durch Lösen von
aK(eK , v) = (R0K , v)K −
∑E∈E (K)
〈R∗E , v〉E ∀v ∈ ṼhK . (3.15)
und setze ηK := ‖∇eK‖K . Wahl von Ṽ hK ?
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 112
Eine Möglichkeit [Bank & Weiser, 1985]
Ṽ hK := P̃K ⊕ B̃K , (3.16)
P̃K := span{φE : E ∈ E (K) \ E hD},
B̃K := span{φK}
Hierbei sei φE die quadratische bzw. biquadratische Kanten-Bläschenfunktion(BF) mit φE 6≡ 0 auf E sowie φK die kubische bzw. biquadratische BF aufK mit
0 ≤ φK ≤ 1, φK ≡ 0 auf ∂K, φK = 1 am Schwerpunkt von K.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 113
Bei dieser Wahl:
• Berechnung von eK gemäß (3.15) erfordert Lösen eines LGS derDimension 4 bzw. 5.
• Wegen (∇v,∇v)K > 0 ∀v ∈ Ṽ hK sind diese eindeutig lösbar (Konstantekann nicht durch BF dargestellt werden).
• Wichtig, da dies Verträglichkeitsbedingung (reine Neumann-RWA) um-geht.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 114
Beispiel: Wir betrachten die RWA aus Beispiel 3, diskretisiert mit bilinearenRechteckelementen auf einem uniformen achsenparallelen Gitter mit 2`
Elementen in jeder Koordinatenrichtung.
Wir vergleichen für verschiedene Gitterweiten
• den exakten Fehler ‖∇(u− uh)‖,• den geschätzten (globalen) Fehler η := (
∑K∈T h η
2K)
1/2,• den Quotienten Xη := η/‖∇(u− uh)‖.
h ‖∇(u− uh)‖ η Xη5.000× 10−1 4.9542× 10−2 5.0314× 10−2 0.984662.500× 10−1 2.5163× 10−2 2.5114× 10−2 0.998061.250× 10−1 1.2581× 10−2 1.2574× 10−2 0.999376.250× 10−2 6.2909× 10−3 6.2881× 10−3 0.999563.125× 10−2 3.1455× 10−3 3.1446× 10−3 0.99972
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 115
Lemma 3.17 Zu u ∈ V0 existiert eine Quasi-Interpolierende Ĩhu ∈ V h0 mit
‖u− Ĩhu‖K ≤ C1(βωK )hK ‖∇u‖ωK ∀K ∈ T h, (3.17a)
‖u− Ĩhu‖E ≤ C2(βωK )h1/2E ‖∇u‖ωK ∀E ∈ E
h. (3.17b)
Hierbei sei ωK die Vereinigung aller Elemen-te mit mindestens einer gemeinsamen Eckemit K (siehe rechts) und βωK eine obereSchranke für das Elementenseitenverhältnisin ωK .
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 116
Lemma 3.18 Sei φ ein auf einem Drei- oder Rechteckelement K definier-tes Polynom. Dann existiert eine (nur vom Seitenverhältnis abhängige)Konstante C mit
‖∇φ‖K ≤C
hK‖φ‖K . (3.18)
Hierbei bezeichnet hK die Länge der längsten Kante in K.
Satz 3.19 Wird die Lösung der Variationsaufgabe (3.2) approximiert durchbilineare Rechteckelemente auf einem Gitter T h mit oberer Schranke β∗ fürdas Seitenverhältnis, so gilt für den durch (3.15) definierten Fehlerschätzermit lokalem Funktionenraum definiert durch (3.16) die Abschätzung
‖∇(u− uh)‖ ≤ C(β∗)
∑K∈T h
η2K + h2∑
K∈T h‖RK −R0K‖2K
1/2 . (3.19)Hierbei bezeichnet h die Länge der längsten Kante in T h.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 117
Satz 3.20 Wird die Lösung der Variationsaufgabe (3.2) approximiert durchbilineare Rechteckelemente auf einem Gitter T h mit oberer Schranke β∗ fürdas Seitenverhältnis, so gilt für den durch (3.15) definierten Fehlerschätzermit lokalem Funktionenraum definiert durch (3.16) die Abschätzung
ηK ≤ C(βωK )‖∇(u− uh)‖. (3.20)
Hierbei bezeichnet ωK die Vereinigung der (5) Elemente, die mit K ge-meinsame Kanten haben.
3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 118
3.6 Matrixeigenschaften
Vektordarstellung von FE-Funktionen: sei {φ1, . . . , φn} Basis von V h0 .
Identifiziere
v =
n∑j=1
vjφj ∈ V h0 ←→ v =
v1...
vn
∈ Rn.L2(Ω)-Norm:
‖v‖2 = (v, v) = v>Mv , [M ]i,j = (φj , φi) Massenmatrix
Energienorm:
‖v‖2a = a(v, v) = v>Av , [A]i,j = a(φj , φi) Steifigkeitsmatrix
3.6 Matrixeigenschaften TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 119
Lemma 3.21 Für FE-Räume aus P1- oder Q1-Elementen basierend aufeiner quasiuniformen Familie von Zerlegungen {T h}h>0 von Ω ⊂ R2existieren von h unabhängige Konstanten c und C mit
ch2min v>v ≤ v>Mv ≤ Ch2 v>v .
Hierbei bezeichne hmin := minK∈T h hK für jede Zerlegung T h von Ω.
Definition 3.22 Eine Familie von Zerlegungen {T h}h>0 von Ω heißt uni-form, falls für alle T h gilt hmin ≥ ch.
Korollar 3.23 Für FE-Räume aus P1- oder Q1-Elementen basierend aufeiner uniformen Familie von Zerlegungen {T h}h>0 von Ω ⊂ R2 existierenvon h unabhängige Konstanten c und C mit
ch2 v>v ≤ v>Mv ≤ Ch2 v>v .
3.6 Matrixeigenschaften TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 120
Bemerkungen 3.24(a) Für Zerlegungen von Rd ist h2 durch hd zu ersetzen.(b) Die Abschätzung gilt ebenso für die Räume Pk, Qk, k ≥ 2 mit von k
abhängigen Konstanten c und C.(c) Bei der Galerkin-Diskretisierung der RWA (3.2) geht die rechte Seite f
der Differentialgleichung ein durch den Vektor
f =
(f, φ1)
...
(f, φn)
∈ Rn.M.a.W.: Hier wird die Funktion f durch deren L2-Projektion nach V h
repräsentiert. Im Gegensatz zu Ansatz- und Testfunktionen ist f nichtdie Koordinatendarstellung der L2-Projektion von f bezüglich der Basis{φ1, . . . , φn}.
3.6 Matrixeigenschaften TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 121
Satz 3.25 Für P1 bzw. Q1-Elemente bezüglich einer Familie uniformerTriangulierungen {T h}h>0 gilt für die Galerkin-Matrix A in (3.5)
ch2 ≤ v>Av
v>v≤ C ∀0 6= v ∈ Rn (3.21)
mit zwei von h unabhängigen Konstanten c und C. Hierbei bezeichnet h dieLänge der längsten Kante in T h und n die Dimension des FE-Raums V h0 .
Bemerkungen 3.26(a) Die Ungleichungen (3.21) beschränken den Wertebereich, und, da A
symmetrisch, somit das Spectrum der Matrix A.(b) Insbesondere gilt für die Konditionszahl κ(A) von A
κ(A) = ‖A‖ ‖A−1‖ = λmax(A)λmin(A)
≤ Cch−2.
3.6 Matrixeigenschaften TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 122
4 Die Lösung der diskreten Poisson-Gleichung
Die Galerkin-Matrix der FE-Diskretisierung der Poisson-RWA (3.1) istsymmetrisch und positiv-definit (gilt nicht notwendig für andere Dgln.)
Allen FE-Diskretisierungen gemeinsam: die Galerkin-Matrix ist dünn be-setzt (sparse). Wir betrachten die Matrix aus Beispiel 2 (L-Gebiet).
bilinear biquadratisch
N h # Einträge 6= 0 Anteil h # Einträge 6= 0 Anteil65 0.25 251 5.94% 0.5 313 7.4%
225 0.125 1339 2.64% 0.25 2073 4.1%
883 0.0625 6107 0.78% 0.125 10141 1.5%
3201 0.0312 26011 0.25% 0.0625 44767 0.44%
12545 0.0156 107291 0.068% 0.0312 187742 0.12%
4 Die Lösung der diskreten Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 123
0 10 20 30 40 50 60
0
10
20
30
40
50
60
nz = 251
h = 0.25, bilineare Elemente
0 10 20 30 40 50 60
0
10
20
30
40
50
60
nz = 313
h = 0.5, biquadratische Elemente
4 Die Lösung der diskreten Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 124
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1Q1 finite element subdivision
h = 0.25, bilineare Elemente−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1Q2 finite element subdivision
h = 0.5, biquadratische Elemente
4 Die Lösung der diskreten Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 125
Direkte oder iterative Verfahren ?
• In der FE-Praxis wurde lange Zeit mit direkten Verfahren, also Variantender Gauß-Elimination, gearbeitet.
• Stark verbreitet: frontal solvers. Hier wird die Faktorisierung der Matrixso früh wie möglich während der Assemblierung vorgenommen.
• Weiterentwicklungen auch heute noch relevant, auch als eigenständigeEliminationsverfahren, z.B. UMFPACK.
• Faustregel: in 2D reichen ausgereifte direkte, auf dünn-besetzte Ma-trizen spezialisierte Löser aus. Für 3D-Diskretisierungen muss aufIterationsverahren zurückgegriffen werden.
• Vorteil direkter Verfahren: typischerweise als ”black box“ einsetzbar.Bei Iterationsverfahren ist deutlich mehr Expertise des Anwenderserforderlich.
4 Die Lösung der diskreten Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 126
4.1 Erste Iterationsverfahren
Zunächst sei allgemein gegeben
Ax = b, A ∈ Cn×n invertierbar, b ∈ Cn. (4.1)
Grundidee:
Startnäherung x0 ≈ A−1b,Erzeuge rekursiv Folge {xm}m∈N mit lim
m→∞xm = A
−1b.
Bezeichnungen:
rm := b −Axm Residuum, Defekt,x ? := A−1b (exakte) Lösung,
em := x? − xm Fehler.
4.1 Erste Iterationsverfahren TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 127
Die ersten Iterationsverfahren (Handrechnung!) beruhen auf Umstellungeinzelner Gleichungen, Auflösen nach den entsprechenden Unbekanntenund Wiedereinsetzen (Relaxation).
• Carl-Friedrich Gauß (∼ 1832), Normalgleichungen aus Astronomie,Landesvermessung, 20–40 Unbekannte
”Ich empfehle Ihnen diesen Modus zur Nachahmung. Schwerlich wer-den Sie je wieder direkt eliminieren, wenigstens nicht, wenn Sie mehrals 2 Unbekannte haben. Das indirekte Verfahren lässt sich halb imSchlafe ausführen, oder man kann während desselben an andere Dingedenken.“ Brief von Gauß an Gerling, 1823
• Carl Gustav Jacobi (∼ 1846)
• Ludwig Seidel (∼ 1862)
• Christian August Nagel (Sächsische Triangulation, 1867-1886) 159Unbekannte
4.1 Erste Iterationsverfahren TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 128
Formale Beschreibung von Relaxationsverfahren durch Zerlegungen (split-tings):
A = M −N , M invertierbar, (4.2a)
überführt (4.1) in Fixpunktform Mx = Nx +b, zugehörige Fixpunktiteration
xm+1 = Txm + c, T = M−1N , c = M−1b. (4.2b)
Typische algorithmische Umsetzung: (beachte: Txm + c = xm + M−1rm)
Algorithmus 1: Relaxationsverfahren.Gegeben : A invertierbar, b, Startnäherung x0
1 for m = 0, 1, 2 . . . bis Abbruchkriterium erfüllt do2 rm ← b −Axm3 Löse Mh = rm4 xm+1 ← xm + h
4.1 Erste Iterationsverfahren TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 129
Klassische Zerlegungen
• M = IRichardson-Verfahren
• M = D := diag(A)Jacobi-Verfahren, Einzelschrittverfahren
• M = D − L, L = − tril(A)aGauß-Seidel-Verfahren, Einzelschrittverfahren, Liebmann-Verfahren
• M = D − ωL, ω ∈ RSOR-Verfahren (successive-overrelaxation)
• M = ω2−ω (1ωD − L)D
−1( 1ωD −U ), U := − triu(A), ω ∈ RSSOR-Verfahren (symmetric SOR)
aMit tril(A) bzw. triu(A) sei, abweichend von der MATLAB-Semantik, der echte unterebzw. obere Dreieckanteil von A bezeichnet.
4.1 Erste Iterationsverfahren TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 130
Für den Fehler der m-ten Iterierten aus (4.2b) gilt
em+1 = x? − xm+1 = (I −M−1A)(x ? − xm) = Tem = · · · = Tme0.
Lemma 4.1 Für eine beliebige Matrix T ∈ Cn×n gilt
limm→∞
Tm = O ⇔ ρ(T ) < 1.
Satz 4.2 Das durch die Zerlegung (4.2a) definierte Relaxationsverfahren(4.2b) konvergiert genau dann für alle Startnäherungen x0, wenn für dieIterationsmatrix T = M−1N gilt ρ(T ) < 1.
4.1 Erste Iterationsverfahren TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 131
4.2 M -Matrizen und reguläre Zerlegungen
Definition 4.3 Eine Matrix A ∈ Cn×n heißt reduzibel (zerlegbar), falls eseine Permutationsmatrix P gibt sodass
PAP> =
[A1,1 A1,2
O A2,2
]
mit quadratischen Matrizen A1,1, A2,2. Eine Matrix heißt irreduzibel, fallssie nicht reduzibel ist.
Lemma 4.4 Eine Matrix A ∈ Cn×n ist genau dann irreduzibel, wenn esfür jede Zerlegung I = I1 ∪ I2 der Indexmenge I = {1, 2, . . . , n} mitI1 6= ∅, I2 6= ∅, I1 ∩I2 = ∅ stets ein Element ai,j von A gibt mit i ∈ I1,j ∈ I2 und ai,j 6= 0.
4.2 M -Matrizen und reguläre Zerlegungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 132
Definition 4.5 Eine Matrix A ∈ Rm×n heißt nichtnegativ (positiv), falls
ai,j ≥ 0 (> 0) 1 ≤ i ≤ m, 1 ≤ j ≤ n.
Entsprechend bedeutet A ≥ B (A > B), dass A−B ≥ O (> O).
Satz 4.6 (Satz von Perron und Frobenius) Sei A ≥ 0 irreduzibel.Dann gelten:
(a) Der Spektralradius ρ(A) ist einfacher Eigenwert von A (der sog.Perron-Eigenwert von A) und positiv.
(b) Zum Perron-Eigenwert gehört ein positiver Perron-Eigenvektor x > 0 .
(c) Jeder nichtnegative Eigenvektor von A ist Eigenvektor zum Perron-Eigenwert ρ(A) (und damit positiv).
(d) Ist O ≤ B ≤ A, so folgt ρ(B) ≤ ρ(A). Ist zusätzlich A 6= B , so giltρ(B) < ρ(A).
4.2 M -Matrizen und reguläre Zerlegungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 133
Definition 4.7 Eine nichtsinguläre Matrix A ∈ Cn×n heißt M-Matrixa, falls
ai,j ≤ 0 für i 6= j und A−1 ≥ O .
Eine Hermitesche M-Matrix heißt Stieltjes-Matrix.
Satz 4.8 Für B ∈ Cn×n mit ρ(B) < 1 ist I −B nichtsingulär. Es gilt
(I −B)−1 =∞∑j=0
Bj .
Diese Reihe wird Neumannsche Reihe genannt. Konvergiert umgekehrtdie Neumannsche Reihe, so ist ρ(B) < 1.
Satz 4.9 Es sei B ∈ Rn×n und B ≥ O . Die Relation ρ(B) < 1 gilt genaudann, wenn die Inverse (I −B)−1 existiert und (I −B)−1 ≥ O .
aDie Bezeichnung wurde 1937 durch Alexander Ostrowski eingeführt, zu Ehren von Her-mann Minkowski.
4.2 M -Matrizen und reguläre Zerlegungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 134
Satz 4.10 Eine Matrix A ∈ Cn×n mit ai,j ≤ 0 für i 6= j ist genau dann eineM-Matrix, wenn D := diag(A) > O (insbesondere sind M -Matrizen stetsreell) und wenn für die Matrix B := I −D−1A gilt ρ(B) < 1.
Definition 4.11 Die Zerlegung A = M − N der Matrix A ∈ Rn×nheißtreguläre Zerlegung von A, wenn M−1 ≥ O und N ≥ O gelten.
Satz 4.12 Ist A = M−N eine reguläre Zerlegung von A und ist A−1 ≥ O ,so gilt
ρ(M−1N ) =ρ(A−1N )
1 + ρ(A−1N )< 1,
d.h. jede reguläre Zerlegung einer Matrix mit nichtnegativer Inversen führtzu einem konvergenten Iterationsverfahren. Diese Aussage gilt insbeson-dere für M-Matrizen.
4.2 M -Matrizen und reguläre Zerlegungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 135
Satz 4.13 Für A ∈ Rn×n mit A−1 > O seien A = M1 −N1 = M2 −N2zwei reguläre Zerlegungen. Falls N2 ≥ N1 ≥ O und N2 6= N1 ist, so gilt
0 < ρ(M−11 N1) < ρ(M−12 N2) < 1,
d.h. das aus der ersten Zerlegung resultierende Verfahren konvergiertschneller als das aus der zweiten.
4.2 M -Matrizen und reguläre Zerlegungen TU Bergakademie Freiberg, WS 2010/111
-
Finite Elemente I 136
Satz 4.14 (Starkes Zeilensummenkriterium) Gilt für A = D − L −U ∈Cn×n
‖D−1(L + U )‖∞ < 1,
so konvergieren Gauß-Seidel und Jacobi-Verfahren.
Satz 4.15 (Schwaches Zeilensummenkriterium) Ist A ∈ Cn×n irreduzi-bel und gilt ∑
j 6=i
|aij | ≤ |aii|, i = 1, . . . , n,
wobei für mindestens einen Index i echte Ungleichheit gilt, so konvergierendas Jacobi- und sowie das Gauß-Seidel-Verfahren.
Satz 4.16 Ist A eine M-Matrix, dann konvergiert sowohl das Jacobi- alsauch das Gauß-Seidel-Verfahren.Ist A zusätzlich irreduzibel, so konvergiert das Gauß-Se