Finite Elemente I - tu-freiberg.deernst/Lehre/FEM_I/Folien12/...Finite Elemente I 2 1 Einleitung...

Institut für Numerische Mathematik und Optimierung

Finite Elemente IWintersemester 2012/13

Erste von zwei Vorlesungen im Modul

”Finite-Element Methoden für Mathematiker“

Hörerkreis: 5. Mm, 7. Mm, 9. Mm, 1. MWM

Oliver ErnstInstitut fü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-Räumen

5. Konvergenztheorie

6. Gleichungslö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-Rä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 Größen erfü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-schränkt, offen und zusammenhä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 Größen zusammen, welcheder Gleichung (1.1) genügen:

Anwendung u k f

Elektrostatik elektr. Potential Permittivität Ladungsdichte

Magnetostatik magn. Potential Permeablilität Ladungsdichte

Wärmetransport Temperatur Leitfähigkeit Wärmequelle

Grundwasserströmung Spiegelhö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 Hilfsgröße und der Flussvektor ∇u die wirklichinteressierende Größe. Letzterer wird oft durch (numerische) Differentiation der (nu-merisch berechneten) Lösung u gewonnen, ein instabiler Vorgang. Sog. gemischteFE-Formlierungen gestatten die direkte numerische Berechnung des Flusses.


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) Lösung.


Finite Elemente I 6

Wir betrachten als Beispiel das elektrische Feld im Außenraum zweier Elek-troden in einem Gehäuse. Die Anordnung sei uniform in einer Richtung: dieElektroden mögen recheckigen Querschnitt besitzen, das Gehäuse sei einKreiszylinder.

Auf den Elektroden Γ1 und Γ2 legen wir jeweils den konstanten Potential-wert φ = 1 bzw. φ = −1 fest, am Gehä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.


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.


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 für den rä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 Rä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 ∂Ω.


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, Lösung



1.3 Minimalflächen

Die durch eine Funktion u : Ω → R, Ω ⊂ R2 ein beschränktes Gebiet,definierte Flä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 Fläche A(u) minimal wird, wenn u folgende Randwertaufgabelöst:

−∇·

(1√

1 + |∇u|2∇u

)= 0, in Ω,

u = g, längs ∂Ω,

wobei g die vorgeschriebenen Randwerte darstellt.

1.3 Minimalflächen TU Bergakademie Freiberg, WS 2010/111


Wir betrachten konkret das Gebiet Ω := {(x, y) ∈ R2 : x2 + y2 < 1} undschreiben auf dessen Rand fü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

Minimalfläche

1.3 Minimalflächen TU Bergakademie Freiberg, WS 2010/111


1.4 Die Navier-Stokes-Gleichungen

Die Strö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


Die Größen u und p sind eindeutig gegeben als Lösung des Systemspartieller Differentialgleichungen

∂tu + (u · ∇)u − ν∆u +∇p = f ,∇·u = 0 ,

(ν bezeichnet die kinematische Viskosität des Fluids) zusammen mit ge-eigneten Anfangsbedingungen

u0 = u(x , 0), p0 = p(x , 0), x ∈ Ω,

sowie Randbedingungen längs ∂Ω.



Als Besipiel betrachten wir das sog. DFG-Benchmarkproblem der Strömungin einem Kanal um einen Kreiszylinder:

Am linken Rand wird ein parabolisches Einströmprofil[u

v

]=

[1.2y(0.41− y)/0.412

0

]vorgeschrieben, am rechten Rand Neumann-Randbedingungen und anden übrigen Rändern u = 0 . Ferner ist ν = 0.001.



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



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

Strömungsgeschwindigkeitsfeld bei t = 7s.



1.5 Die Maxwell-Gleichungen

Makroskopische elektromagnetische Phänomene werden beschriebendurch die Maxwell-Gleichungen

∇×E + Ḃ = 0 , (1.3a)

∇×H − Ḋ = J , (1.3b)∇·D = ρ, (1.3c)∇·B = 0, . (1.3d)

Hierbei bezeichnen E die elektrische Feldstärke, D die Verschiebungs-stromdichte, H die magnetische Feldstärke, B die magnetische Fluß-dichte, J die Leitungsstromdichte sowie ρ die elektrische Ladungsdichte.Gleichung (1.3a) ist das Faradaysche Induktionsgesetz, (1.3b) die Maxwell-AmpË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


In linearen und isotropen Materialien gelten die konstitutiven Gleichungen

B = µH , D = �E , J = σE + J i, (1.4)

mit der elektrischen Permittivität �, magnetischen Permeabilität µ und elek-trischer Leitfähigkeit σ als skalare Proportionalitätsfaktoren.

Wir betrachten ein Beispiel aus der geoelektrischen Erkundung, bei demdurch die Messung elektrischer Felder an der Erdoberfläche Rückschlüsseauf die örtliche Verteilung der Leitfähigkeit σ im Untergrund und damit aufdessen Zusammensetzung gezogen werden.



• In geophysikalischen Anwendungen kann der Verschiebungsstrom Ḋvernachlässigt werden. Ferner kann µ = µ0 als konstant angenommenwerden.

• Um nur mit einem Feld rechnen zu mü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 äußere Einheitsnormalelängs ∂Ω.

• Zur Zeit t = t0 wird eine Anfangsbedingung E(x , 0) = E0(x ), x ∈ Ω,benötigt.



Man erhält so die Anfangs-Randwertaufgabe:

µ0σĖ +∇×∇×E = −µ0J̇ i in Ω, (1.5a)n ×E = 0 längs ∂Ω, (1.5b)E(x , 0) = E0(x ). (1.5c)



2 Variationstheorie

Die FEM wird angewandt auf die sog. Variationsformulierung von Rand-und Anfangswertaufgaben für Differentialgleichungen.

In diesem Abschnitt stellen wir die wichtigsten Existenz- und Eindeutig-keitssätze bereit und beleuchten die Beziehung zwischen klassischer undVariationsformulierung von Differentialgleichungen.

Die Charakterisierung von Funktionen, für die ein Funktional stationär wird,als Lösung von Differentialgleichungen ist das Thema der Variationsrech-nung, einer mathematischen Disziplin, deren Ausgangspunkt üblicherweiseauf das Jahr 1696 datiert wird, in dem Johann Bernoulli das Brachistochro-nenproblem der mathematischen Weltöffentlichkeit stellte.

2 Variationstheorie TU Bergakademie Freiberg, WS 2010/111


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 Körper allein unter demEinfluss der Schwerkraft in kü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


Galileo Galilei

GALILEO zeigte 1638, dass ein Körperauf dem Weg von A nach B schnel-ler über den Umweg ADB rutscht alsauf direktem (geraden) Wege AB. Eben-so zeigt er, dass ACDB, ACDEB undACDEFB jeweils zu kürzeren Laufzei-ten führen, und schliesst, dass ein Kreis-bogen die Lösung sei.

A

B

C

D

EF

Falsche Brachistochrone desGALILEO



Gottfried Wilhelm Leibniz Jacob Bernoulli Leonhard Euler

• 1684, Acta Eruditorum: LEIBNIZens Abhandlung über Differential- undIntegralrechnung erscheint.

• Die Brüder JAKOB UND JOHANN BERNOULLI bauen diese Ideen syste-matisch zu praktischen und universell einsetzbaren Lösungtechnikenaus.

• Später im 18. Jahrhundert entwickelt dann LEONHARD EULER formaleLösungstechniken für Variationsprobleme.



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



Aufsummieren aller Teilzeiten ∆t zur Gesamtzeit T gefolgt vom Grenzüber-gang ∆t→ 0 führt auf den Ausdruck

T =1√2g

∫ xB0

√1 + y′(x)2√−y(x)

dx,

der die Laufzeit T in Abhä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.



Das Brachistochronenproblem besitzt folgende allgemeinere Struktur: ge-sucht ist eine zulässige Funktion, (Glattheit, Randbedingungen)

x : [0, tE ]→ Rn, t 7→ x (t),

welche ein Funktional (”Zielfunktional“)

ψ(x ) :=

∫ tE0

L(t,x (t), ẋ (t)) dt

minimiert, mit einer gegebenen Funktion L : R2n+1 → R.

Beim Brachistochronenproblem:

n = 1, x (t) = x(t), L(t, x, ẋ) =

√1 + ẋ2√−x

.



Lösungsansatz: Auf EULER und LAGRANGE geht folgender Ansatz zurück,der notwendige Bedingungen für das Vorliegen einer Minimallösung liefert:

Zu einer Minimallö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 für das Variationsproblem zulässigsind. Insbesondere entspricht x ? = x0.

Einsetzen in das Zielfunktional ψ führt auf eine Funktion J = J(�) := ψ(x�)die an der Stelle � = 0 ein lokales Minimum besitzt, für die also

0 =dJ(�)

d�

∣∣∣∣�=0

=

∫ tE0

[∂L

∂x(t,x ?(t), ẋ ?(t))− ∂

2L

∂t∂ẋ(t,x ?(t), ẋ ?(t))

]φ(t) dt .



Dieses Integral kann nur dann für eine hinreichend große Klasse (z.B. C∞)Funktionen verschwinden, wenn auch

∂L

∂x(t,x ?(t), ẋ ?(t))− ∂

2L

∂t∂ẋ(t,x ?(t), ẋ ?(t)) = 0 . (2.1)

Gleichung (2.1) heißt Euler-Lagrange-Gleichung des Variationsproblemsund liefert eine notwendige Bedingung für das Vorliegen einer minimieren-den Funktion x ?(t).

Beim Brachistochronenproblem gilt ∂tL(t, x, ẋ) = 0 (autonomes Problem).

Folge: Für die Hamilton-Funktion

H(t, x, ẋ) := L(t, x, ẋ)− ẋ∂ẋL(t, x, ẋ)

gilt

d

dtH(t, x?(t), ẋ?(t)) ≡ 0, d.h. H(t, x?(t), ẋ?(t)) = const ∀t.

Man nennt H daher auch ein erstes Integral von L.



Speziell beim Brachistochronenproblem erhält man:

H(t, x, ẋ) =1

√−x√

1 + ẋ2 ≡ C(> 0),

oder äquivalent−x(1 + ẋ2) =: 2r > 0

mit einer neuen Konstanten r.Auflösen nach ẋ unter Beachtung von ẋ(t) < 0 führt auf die AWA

ẋ(t) = −

√2r + x(t)

−x(t), x(0) = 0,

wofür die Methode der Trennung der Veränderlichen auf

t = t(x) = −∫ x

0

√−ξ

2r + ξdξ

führt.



Die geschickte Substitution

ξ = −2r sin2 ϕ2, ϕ ∈

[0,π

2

], dξ = −2r sin ϕ

2cos

ϕ

2dϕ

ξ = 0⇔ ϕ = 0, ξ = x⇔ ϕ = ϑ

führt schließlich auf die parametrisierte Lösungsdarstellung

t(ϑ) = r(ϑ− sinϑ), x(ϑ) = −r(1− cosϑ)

bzw., nach Rückbenennung auf die alten Variablen x und y, auf

x(ϑ) = r(ϑ− sinϑ), y(ϑ) = −r(1− cosϑ), 0 ≤ ϑ ≤ ϑB .

Der Parameterwert ϑB am Endpunkt lässt sich durch Lö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 .



Die Lö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.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


Bemerkungen 2.2(a) Im Fall eines komplexen Vektorraums V fordert man Antilinearitä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(·, ·) gehörende Energie-Innenprodukt genannt.

(c) Die Eigenschaft der Stetigkeit (2.3a) zieht die Stetigkeit der Bilinearforma(·, ·) als Funktion ihrer beiden Argumente nach sich.



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 Lösung.

(b) Die Minimierungsaufgabe (2.4) ist äquivalent zur Variationsaufgabe

Bestimme u ∈H sodass a(u, v) = `(v) ∀v ∈H . (2.5)



2.3 Die Dirichletsche RWA fü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) Lösung.

2.3 Die Dirichletsche RWA für die Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111


Im einfachsten Fall Γ = ΓD erhalten wir die sog. Dirichletsche Randwert-aufgabe (RWA) fü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 lä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 längs Γ

das Funktional J(u) :=1

2

∫Ω

|∇u|2 dx−∫

Ω

fu dx.(2.9)



Satz 2.4 Seien g ∈ C(Γ) sowie f ∈ C(Ω) gegebene Funktionen. Sei ferneru ∈ C2(Ω). Dann gilt

(a) Löst u die Minimierungsaufgabe (2.9), so löst u auch die Randwertauf-gabe (2.7).

(b) Die Funktion u ist genau dann Lösung der RWA (2.7), wenn u Lösungder verallgemeinerten RWA (2.8) ist.

Bemerkung 2.5 Nach Satz 2.4 löst jede hinreichend glatte Lösung derVariationsaufgaben (2.8) oder (2.9) auch die RWA (2.7).Schwierigkeit: in vielen Anwendungen tritt der Fall auf, dass keine derartglatte Lösung existiert.



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


2.4 Verallgemeinerte Ableitungen

Sei wieder Ω ⊂ Rd, d ∈ N offen und nicht leer.

Für u ∈ C1(Ω) und eine ”Testfunktion“ φ ∈ C∞0 (Ω) gilt die partielle Integra-

tionsformel ∫Ω

u ∂jφdx = −∫

Ω

(∂ju)φdx, 1 ≤ j ≤ d.

Mit anderen Worten: für v := ∂ju gilt die Beziehung∫Ω

u∂jφdx = −∫

Ω

vφ dx ∀φ ∈ C∞0 (Ω). (2.11)

Idee: definiere durch (2.11) die partielle Ableitung auch für den Fall, dassu nicht (im klassischen Sinn) nach xj differenzierbar.

2.4 Verallgemeinerte Ableitungen TU Bergakademie Freiberg, WS 2010/111


Definition 2.6 Sei Ω ⊂ Rd, d ∈ N, offen und nicht leer. Gilt fü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

für beliebige reelle Werte von α die verallgemeinerte Ableitung von u.

2.4 Verallgemeinerte Ableitungen TU Bergakademie Freiberg, WS 2010/111


2.5 Die Sobolev-Rä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.

Fü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 überall übereinstimmen.

2.5 Die Sobolev-Räume H1(Ω) und H10(Ω) TU Bergakademie Freiberg, WS 2010/111


Definition 2.11 Sei Ω ⊂ Rd, d ∈ N, offen und nicht leer. Mit H10 (Ω) be-zeichnen wir den Abschluss der Menge C∞0 (Ω) bezü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 überall in(a, b) sowie v(a) = v(b) = 0. Ferner gilt die Abschä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.



Beispiel 2.14 Sei −∞ < a < b < ∞ und die Funktion u ∈ C[a, b]stü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.



Verallgemeinerte Randwerte

Definition 2.15 Sei Ω ⊂ Rd, d ∈ N, offen, nicht leer und beschränkt. BeiFunktionen u ∈ H10 (Ω) sagen wir, sie erfüllen die Randbedingung

u = 0 längs ∂Ω (2.13)

im verallgemeinerten Sinne.

Bemerkungen 2.16(a) Beachte: C∞0 (Ω) liegt dicht in H10 (Ω).

(b) Für d = 1 gilt (2.13) sogar im klassischen Sinne (vgl. Beispiel 2.13).

(c) Für d > 1 gilt bei hinreichend glattem Rand Γ = ∂Ω

‖u‖L2(Γ) ≤ C‖u‖1 ∀u ∈ H1(Ω). (2.14)



2.6 Die Poincaré-Friedrichs-Ungleichung

Mit Blick auf den abstrakten Lösungssatz 2.3 stellt die Poincaré-Friedrichs-Ungleichung die Koerzivität der Bilinearform a(·, ·) sicher.

Satz 2.17 Sei Ω ⊂ Rd, d ∈ N, offen, nicht leer und beschränkt. Dann exi-stiert eine Konstante C > 0, sodass die Poincaré-Friedrichs-Ungleichung∫

Ω

|u(x)|2 dx ≤ C∫

Ω

|∇u(x)|2 dx ∀u ∈ H10 (Ω) (2.15)

erfüllt ist.

2.6 Die Poincaré-Friedrichs-Ungleichung TU Bergakademie Freiberg, WS 2010/111


2.7 Ein Existenzsatz für das Dirichlet-Problem

Wir betrachten nach den vorangehenden


Satz 2.18 (Dirichlet-Prinzip) Sei Ω ⊂ Rd, d ∈ N, offen, nicht leer sowie be-schränkt. Seien ferner f ∈ L2(Ω) sowie g ∈ H1(Ω) gegebene Funktionen.Dann gelten

(a) Die verallgemeinerte Minimierungsaufgabe (2.17) besitzt eine eindeu-tige Lösung u ∈ H1(Ω).

(b) Diese ist auch die eindeutige Lösung der verallgemeinerten Randwert-aufgabe (2.16).

2.7 Ein Existenzsatz für das Dirichlet-Problem TU Bergakademie Freiberg, WS 2010/111


2.8 Weitere Existenzsätze

Folgender Satz sichert Existenz und Eindeutigkeit der Lösung linearerVariationsaufgaben unter der entscheidenden Annahme der Koerzivitä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 fü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 Lösung.

2.8 Weitere Existenzsätze TU Bergakademie Freiberg, WS 2010/111


Bemerkungen 2.20(a) Für die Lösung u aus Satz 2.19 gilt die Stabilitätsabschätzung

‖u‖ ≤ 1α‖`‖H ′ .

(b) Ist a(·, ·) zusätzlich noch symmetrisch [Hermitesch], so folgen Existenzund Eindeutigkeit der Lösung aus dem Rieszschen DarstellungssatzA.32.



Satz 2.21 (Nečas, 1968; Babuška, 1972) Seien X , Y Banach-Rä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 für jedes ` ∈ Y ′ eine eindeutige Lösung der Variationsaufgabe

Bestimme u ∈ X mit a(u, v) = `(v) ∀v ∈ Y . (2.21)

Die Lösung hängt stetig von ` ab, es gilt

‖u‖X ≤1

γ‖`‖Y ′ . (2.22)



Bemerkung 2.22 Das Lax-Milgram Lemma ist ein Korollar von Satz 2.21

Folgender Zusatz zu Satz 2.21 ist manchmal nü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) äquivalent zu

‖Ax‖Y ′ ≥ γ‖x‖X ∀x ∈X .



2.9 Euler-Lagrange-Gleichungen

Auf Euler und Lagrange geht eine Methode zurück, Funktionen, welche eingegebenes Funktional stationär machen, als Lö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 stationärer ”Punkt“ von J gegenüber glattenStörungen φ ∈ C∞0 (Ω) mit kompaktem Träger, wenn

d

dtJ(u+ tφ)

∣∣∣∣t=0

= 0 ∀φ ∈ C∞0 (Ω).

2.9 Euler-Lagrange-Gleichungen TU Bergakademie Freiberg, WS 2010/111


Differentiation im Integral∫Ω

d

dtF (x, u+ tφ,∇u+ t∇φ)

∣∣t=0

dx

unter Anwendung der Kettenregel fü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 führt.

2.9 Euler-Lagrange-Gleichungen TU Bergakademie Freiberg, WS 2010/111


3 Finite-Elemente Approximationen zur Lösungder Poisson-Gleichung in 2D

3 Finite-Elemente Approximationen zur Lösung der Poisson-Gleichung in 2D TU Bergakademie Freiberg, WS 2010/111


3.1 Aufgabenstellung

Sei Ω ⊂ R2 ein beschränktes polygonales Gebiet mit Rand Γ bestehend aus denTeilmengen ΓD und ΓN sowie g und h hinreichend glatte Randfunktionen. ZurApproximation der Lösung der RWA

−∆u = f in Ω, (3.1a)

u = g längs ΓD, (3.1b)

∂u

∂n= h längs ΓN (3.1c)

betrachten wir die zugehö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


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 fü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 über zur Linearform ˜̀ := `− a(ug, ·).

3.1 Aufgabenstellung TU Bergakademie Freiberg, WS 2010/111


3.2 Referenzaufgaben

Auf folgende Beispielaufgaben kommen wir im Laufe dieses Abschnittswieder zurück.

3.2 Referenzaufgaben TU Bergakademie Freiberg, WS 2010/111


Beispiel 1: Ω = (−1, 1)2, ΓD = Γ, f ≡ 1, g ≡ 0.

Dies stellt ein Modell für Wärmeausbreitung auf einer quadratischen Platte dar.Durch Trennung der Veränderlichen bestimmt man die Reihenlö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



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



Beispiel 3: Ω = (−1, 1)2, ΓD = Γ, f ≡ 0, g = u|Γmit exakter Lö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



Beispiel 4: Ω = (−1, 1)2 \ [−1, 0]2, ΓD = Γ, f ≡ 1, g = u|Γ mit exakterLö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.3 Galerkin-Diskretisierung

Das Galerkin-Verfahren zur Approximation der Lösung von (3.2) bestehtdarin, Ansatz- und Testraum durch endlichdimensionale Rä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 Testräume verwendet, spricht man vonPetrov-Galerkin-Verfahren, andernfalls von Bubnov-Galerkin-Verfahren.

Wir betrachten zunächst die homogenisierte Variationsaufgabe (3.3) undeinen n-dimensionalen Unterraum V h0 ⊂ V0. (Oder, ä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


Ist {φ1, φ2, . . . , φn} eine Basis von V h0 sowie uh =∑nj=1 ujφj , so ist (3.4)

ä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 eindeutigeLösung.

• Die Galerkin-Matrix A aus (3.5) ist symmetrisch und positiv-definit.

3.3 Galerkin-Diskretisierung TU Bergakademie Freiberg, WS 2010/111


3.4 Finite-Element-Räume in 2D

Galerkin-Verfahren lassen zunächst beliebige Unterräume zu (Spektralver-fahren, Integralgleichungen, etc.).

Die finite-Element-Methode (FEM) zeichnet sich durch die Verwendung vonBasisfunktionen mit kleinem Träger aus. Fast ausschließlich werden hierfürstü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 Fällen, auch in 3D) Triangu-lierungen heißen.

Weitere gebräuchliche Bezeichnungen sind Vernetzung, Netz oder Gitter.

3.4 Finite-Element-Räume in 2D TU Bergakademie Freiberg, WS 2010/111


Mit V h sei ein zunächst beliebiger endlichdimensionaler Raum auf Ω defi-nierter Funktionen bezeichnet. Entsprechend sei

PK := {v|K : v ∈ V h}

der durch sämtliche Einschränkungen von Funktionen aus V h auf K auf-gespannte Raum.

Bei einer konformen FE-Diskretisierung ist V h ⊂ V erforderlich. EineCharakterisierung von Konformitä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) für alleK ∈ T h, so gelten

V h ⊂ H1(Ω), sowie {v ∈ V h : v = 0 auf ∂Ω} ⊂ H10 (Ω).



Beweis: Aufgrund unserer Annahmen gilt bereits V h ⊂ L2(Ω). Damit auch V h ⊂H1(Ω) mü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 äußeren Einheitsnormalen längs ∂K). Summationü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 φ längs ∂Ω verschwindet und die (orientierten!)Randintegrale längs innerer Ränder je zweimal mit entgegengesetztem Umlaufsinnauftreten. �



3.4.1 Dreieckelemente

Wir betrachten eine Zerlegung T h von Ω in (abgeschlossene) Dreiecke Kmit folgenden Eigenschaften

(a) Ω = ∪K∈T hK.

(b) Für zwei beliebige Dreiecke K1,K2 ∈ T h ist K1 ∩ K2 entweder leer,ein gemeinsamer Punkt oder eine gemeinsame Kante.

Als Diskretisierungsparameter wird üblicherweise

h := maxK∈T h

diamK,

also ein Maß für die Feinheit der Zerlegung, genommen.



Beispiel 1: Dreieckszerlegung des Äußeren eines Tragflächenquerschnitts.



Besipiel 2: Triangulierungen verschiedener Gebiete(mittels der distmesh-Software von Strang und Persson)



Beispiel 3: Zerlegung aus Drei- und Vierecken.



Um aus stückweisen Funktionen auf den Teilmengen K von T h einenUnterraum von H1(Ω) zu erhalten, müssen die hieraus konstruierten Funk-tionen nach Satz 3.2 stetig sein. Für stückweise Polynome bezüglich T h

ist dies gewährleistet, sofern diese nur an allen Dreieckskanten stetig sind.

Der einfachste solche Raum ist der aller stückweise linearen Funktionenauf Ω bezü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 Gradhöchstens eins bezeichnet sei.

Ein Unterraum von V0 ist gegeben durch

V h0 = {v ∈ V h : v|ΓD = 0}.



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 vollstä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 nü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}.



Triangulierung des L-förmigen Gebiets aus Abschnitt 2.1 mit dem Träger dreiernodaler Basisfunktionen.



Konsequenz fü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 dünn besetzt.



Aufbau des Galerkin-Systems

Großer Vorteil der FEM:

• Aufbau des FE Gleichungssystems – man nennt diesen Vorgang As-semblieren – verlä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 für unser Modellpro-blem zusammen. Das Vorgehen besteht darin, die Rechnung auf Operatio-nen auf den einzelnen Elementen zurückzuführen, deren Ergebnisse dannzusammengefügt (assembliert) werden.



Sei K ∈ T h: dann gilt für i, j = 1, 2 . . . , ñ:

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

[ÃK ]i,j := aK(φj , φi) i, j = 1, 2, . . . , ñ,

[b̃K ]i := `K(φi, i = 1, 2, . . . , ñ,

folgt alsoÃ =

∑K∈T h

ÃK , b̃ =∑

K∈T hb̃K .



Da jedes Element nur zum Träger dreier Basisfunktionen gehört, sind nur(höchstens) neun bzw. drei Einträge von ÃK bzw. b̃K von Null verschieden.

Welche Einträ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 führt neben der bisherigen globalen Knotennumerierung x1, x2, . . . , xñin jedem Element zusätzlich eine lokale Nummerierung x(K)1 , x

(K)2 , x

(K)3 der

zu K gehörenden Knoten ein.



Globale Nummerierung der Knoten (rot) und der Elemente (schwarz)einer Triangulierung des L-Gebiets.



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 )

.

Trä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 gehörendenBasisfunktionen zur globalen Nummerierung {φj}ñ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.



Nach diesen Überlegungen erhalten wir folgenden Algorithmusa zur As-semblierung:

Initialisiere Ã := 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)Ã([i1i2i3], [i1i2i3]) := Ã([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

.



Referenzelement

Hilfreich für Implementierung und Analyse: Bezug auf ein ReferenzelementK̂ ⊂ R2. Fü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 üblich: Einheitssimplex

K̂ = {(ξ, η) ∈ R2 : 0 ≤ ξ ≤ 1, 0 ≤ η ≤ 1− ξ}

Fü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.



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



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.



Rückführung der Integration auf das Referenzelement

Die bei der Assemblierung der Galerkin-Matrix anfallenden Integrale wer-den ebenfalls auf das Referenzelement zurückgeführt. Dies ist auch fü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.



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

].



Dreieckelemente vom Grad zwei

Nimmt man zu den Ecken des Referenzdreiecks K̂ noch die Kantenmittel-punkte als Knoten hinzu, so lautet die zugehö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


Die vier nodalen Basisfunktionen im Dreieckelement vom Grad zwei



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.



Beim einfachsten H1-konformen Viereckelement wird als Polynomraum

Q1 := span{1, ξ, η, ξη},

der Raum der bilinearen Funktionen verwendet.

Im Referenzelement, ü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 + η).



Beim quadratischen Viereckelement ist der Polynomraum

Q2 := span{1, ξ, η, ξη, ξ2, ξ2η, ξη2, η2},

der der biquadratischen Funktionen. Als neue Freiheitsgrade gegenü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



Biquadratische nodale Basisfunktionen zu einer Ecke, einer Seitenmitte und demMittelpunkt.



Bei einem achsenparallelen RechteckK mit Kantenlä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.5 Approximationsfehler

Wir betrachten zunä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) erfü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 über ‖u− uh‖, insbesondere für h→ 0?

3.5 Approximationsfehler TU Bergakademie Freiberg, WS 2010/111


3.5.1 Das Céa-Lemma

Satz 3.3 (Lemma von Céa) Gelten für die Variationsaufgabe (3.7) die Vor-aussetzungen des Lax-Milgram-Lemmas, so gilt für den Fehler u − uh derGalerkin-Approximation

‖u− uh‖ ≤ Cα

infv∈V h

‖u− v‖. (3.9)

Hierbei bezeichnen C und α die Stetigkeits- bzw. Koerzivitätskonstante ausdem Lax-Milgram-Lemma.



Bemerkungen 3.4(a) Ist die Bilinearform zusätzlich noch symmetrisch (Hermitesch), so läßt

sich (3.9) verbessern zu

‖u− uh‖ ≤√C

αinfv∈V h

‖u− v‖. (3.10)

(b) Für die Konvergenz einer Folge von Galerkin-Approximationen (imGrenzwert h→ 0) ist somit hinreichend, dass für die zugehörige Familie{V h}h>0 von Unterräumen gilt

limh→0

infv∈V h

‖u− v‖ = 0 ∀u ∈ V .

(c) In diesem Sinne wird durch das Céa-Lemma die Abschätzung desGalerkin-Fehlers zu einem Approximationsproblem.

(d) Die Tatsache, dass a(u− uh, v) = 0 für alle v ∈ V h wird auch Galerkin-Orthogonalität genannt.



Das Céa-Lemma läßt sich auch auf den Fall verallgemeinern, dass kei-ne Koerzivität, sondern lediglich eine inf-sup-Bedingung vorliegt (vgl.Satz 2.21).

Satz 3.5 Seien X , Y Banach-Räume und a : X × Y → R eine Bilinear-form, welche die Voraussetzungen von Satz 2.21 erfüllt. Ferner seien fürUnterräume X h ⊂ X und Y h ⊂ Y gleicher Dimension die Bedingung(2.19) mit einer Konstanten γh sowie Bedingung (2.20) erfüllt. Dann gilt fürdie Lösung uh der diskreten Variationsaufgabe

Bestimme uh ∈X h sodass a(uh, v) = `(v) ∀v ∈ Y h, (3.11)

und die Lösung u ∈X von (2.21) die Abschätzung

‖u− uh‖X ≤(

1 +C

γh

)inf

v∈X h‖u− v‖X .



Bemerkung 3.6 Im Gegensatz zum Lax-Milgram Lemma folgen Bedingun-gen (2.19) und (2.20) nicht aus deren Gültigkeit für die Räume X und Y .Insbesondere kann die Konstante γ in (2.19) von h abhängen, d.h. γ = γh.Entscheidend ist, ob γh gleichmäßig nach unten beschränkt ist.



3.5.2 A priori Fehlerabschätzungen

Die einfachsten a priori Fehlerabschätzungen beruhen auf zusätzlicherRegularitä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 bezüglich des Innenproduktes

(u, v)2 := (u, v)1 +2∑

i,k=1

(∂jku, ∂jkv).

Hierzu gehört die Norm

‖u‖2 =(‖u‖21 + |u|22

)1/2, |u|2 :=

( 2∑i,k=1

‖∂jku‖2)1/2

.



Definition 3.8 Die Variationsaufgabe (3.2) heißt H2-regulär falls für hinrei-chend glatte Daten f , g und h die Lösung sogar in Vg ∩H2(Ω) liegt.

Aufgrund des Sobolevschen Einbettungssatzes liegt in der Ä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 zulässige Triangulierung des Gebietes Ω undV h der FE-Raum aus stückweise linearen Funktionen bezüglich T h. Fü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.



Lemma 3.10 Für alle u ∈ H2(Ω) und alle K ∈ T h gilt

‖∇(u− IKu)‖2K ≤2h2K|K|‖∇̂(û− IK̂ û)‖

2K̂.

Folgende Aussage ist ein Spezialfall des Bramble-Hilbert-Lemmas, wel-ches wir in allgemeiner Form später beweisen werden:

Lemma 3.11 Es existiert eine Konstante C sodass für alle û ∈ H2(K̂) gilt

‖∇̂(û− IK̂ û)‖K̂ ≤ C|û− IK̂ û|2,K̂ = C|û|2,K̂ .

Lemma 3.12 Für alle u ∈ H2(K) und û = u ◦ FK gilt

|û|22,K̂≤ 12h2K

h2K|K||u|22,K .



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 Für die globale Interpolierende Ih eines stückweise linearenFE-Raumes V h ⊂ H1(Ω) bezü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 gleichmäßige untere Schranke θ für den minimalenInnenwinkel aller Dreiecke K ∈ T h für alle h existiert.



Satz 3.16 Ist die Variationsaufgabe (3.2)H2-regulär, so gilt für die Galerkin-Approximation uh bezüglich eines Unterraumes V h ⊂ H1(Ω) aus stückweiselinearen Funktionen bezüglich einer quasiuniformen Triangulierung T h mith := maxK∈T h diamK

‖u− uh‖1 ≤ Ch|u|2mit einer von u und h unabhängigen Konstanten C.



3.5.3 A posteriori Fehlerabschätzungen

Gegeben: RWA (3.2), Lösung u, V h basierend auf T h, u ≈ uh ∈ V h

Gesucht: (lokaler) Fehlerschätzer η : T h → R, K 7→ ηK mit

• ηK ≈ ‖∇(u− uh)‖K ,• Aufwand zur Berechnung von {ηK}K∈T h möglichst linear in |T h|,• obere Schranke der Form∑

K∈T h‖∇(u− uh)‖2K ≤ C(θ)

∑K∈T h

η2K ,

(θ untere Schranke für minimalen Innenwinkel),• untere Schranke der Form

ηK ≤ C(θωK )‖∇(u− uh)‖ωK ,

ωK lokale Elementumgebung von K, K ⊂ ωK ⊂ T h.



Ausgangspunkt: fü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 beschränkt.

Idee: Schätze u − uh durch Schätzung von ‖` − a(uh, ·)‖ bzw. durchnäherungsweise Lö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).



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 äußerer Normalenvektor längs E bez. K,

〈·, ·〉E Innenprodukt in L2(E), (Dualform in H−1/2(E)×H1/2(E))

n · ∇uh (diskreter) Fluss über E, i.A. unstetig an Elementgrenzen.

Definiere daher für zwei Elemente K1 und K2 mit gemeinsamer Kante Eden Sprung des Flusses von v über E durch

s∂v

∂n

{:= (∇v|K1 −∇v|K2) · nE,K1 = (∇v|K2 −∇v|K1) · nE,K2 .



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 für exakte Lösung u.



Hier: betrachte lineare Dreiecke bzw. bilineare Rechtecke

• RE ≡ konstant bei linearen Dreiecken,

• in beiden Fä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 .



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) für alle v ∈ V0 bestimmt.

Idee: bestimme eK ≈ u− uh auf K durch Lösen von

aK(eK , v) = (R0K , v)K −

∑E∈E (K)

〈R∗E , v〉E ∀v ∈ ṼhK . (3.15)

und setze ηK := ‖∇eK‖K . Wahl von Ṽ hK ?



Eine Möglichkeit [Bank & Weiser, 1985]

Ṽ 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-Blä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.



Bei dieser Wahl:

• Berechnung von eK gemäß (3.15) erfordert Lösen eines LGS derDimension 4 bzw. 5.

• Wegen (∇v,∇v)K > 0 ∀v ∈ Ṽ hK sind diese eindeutig lösbar (Konstantekann nicht durch BF dargestellt werden).

• Wichtig, da dies Verträglichkeitsbedingung (reine Neumann-RWA) um-geht.



Beispiel: Wir betrachten die RWA aus Beispiel 3, diskretisiert mit bilinearenRechteckelementen auf einem uniformen achsenparallelen Gitter mit 2`

Elementen in jeder Koordinatenrichtung.

Wir vergleichen für verschiedene Gitterweiten

• den exakten Fehler ‖∇(u− uh)‖,• den geschä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



Lemma 3.17 Zu u ∈ V0 existiert eine Quasi-Interpolierende Ĩhu ∈ V h0 mit

‖u− Ĩhu‖K ≤ C1(βωK )hK ‖∇u‖ωK ∀K ∈ T h, (3.17a)

‖u− Ĩ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 für das Elementenseitenverhältnisin ωK .



Lemma 3.18 Sei φ ein auf einem Drei- oder Rechteckelement K definier-tes Polynom. Dann existiert eine (nur vom Seitenverhältnis abhängige)Konstante C mit

‖∇φ‖K ≤C

hK‖φ‖K . (3.18)

Hierbei bezeichnet hK die Länge der längsten Kante in K.

Satz 3.19 Wird die Lösung der Variationsaufgabe (3.2) approximiert durchbilineare Rechteckelemente auf einem Gitter T h mit oberer Schranke β∗ fürdas Seitenverhältnis, so gilt für den durch (3.15) definierten Fehlerschätzermit lokalem Funktionenraum definiert durch (3.16) die Abschätzung

‖∇(u− uh)‖ ≤ C(β∗)

∑K∈T h

η2K + h2∑

K∈T h‖RK −R0K‖2K

1/2 . (3.19)Hierbei bezeichnet h die Länge der längsten Kante in T h.



Satz 3.20 Wird die Lösung der Variationsaufgabe (3.2) approximiert durchbilineare Rechteckelemente auf einem Gitter T h mit oberer Schranke β∗ fürdas Seitenverhältnis, so gilt für den durch (3.15) definierten Fehlerschätzermit lokalem Funktionenraum definiert durch (3.16) die Abschätzung

ηK ≤ C(βωK )‖∇(u− uh)‖. (3.20)

Hierbei bezeichnet ωK die Vereinigung der (5) Elemente, die mit K ge-meinsame Kanten haben.



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


Lemma 3.21 Für FE-Räume aus P1- oder Q1-Elementen basierend aufeiner quasiuniformen Familie von Zerlegungen {T h}h>0 von Ω ⊂ R2existieren von h unabhängige Konstanten c und C mit

ch2min v>v ≤ v>Mv ≤ Ch2 v>v .

Hierbei bezeichne hmin := minK∈T h hK für jede Zerlegung T h von Ω.

Definition 3.22 Eine Familie von Zerlegungen {T h}h>0 von Ω heißt uni-form, falls für alle T h gilt hmin ≥ ch.

Korollar 3.23 Für FE-Räume aus P1- oder Q1-Elementen basierend aufeiner uniformen Familie von Zerlegungen {T h}h>0 von Ω ⊂ R2 existierenvon h unabhängige Konstanten c und C mit

ch2 v>v ≤ v>Mv ≤ Ch2 v>v .



Bemerkungen 3.24(a) Für Zerlegungen von Rd ist h2 durch hd zu ersetzen.(b) Die Abschätzung gilt ebenso für die Räume Pk, Qk, k ≥ 2 mit von k

abhä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

repräsentiert. Im Gegensatz zu Ansatz- und Testfunktionen ist f nichtdie Koordinatendarstellung der L2-Projektion von f bezüglich der Basis{φ1, . . . , φn}.



Satz 3.25 Für P1 bzw. Q1-Elemente bezüglich einer Familie uniformerTriangulierungen {T h}h>0 gilt für die Galerkin-Matrix A in (3.5)

ch2 ≤ v>Av

v>v≤ C ∀0 6= v ∈ Rn (3.21)

mit zwei von h unabhängigen Konstanten c und C. Hierbei bezeichnet h dieLänge der längsten Kante in T h und n die Dimension des FE-Raums V h0 .

Bemerkungen 3.26(a) Die Ungleichungen (3.21) beschränken den Wertebereich, und, da A

symmetrisch, somit das Spectrum der Matrix A.(b) Insbesondere gilt für die Konditionszahl κ(A) von A

κ(A) = ‖A‖ ‖A−1‖ = λmax(A)λmin(A)

≤ Cch−2.



4 Die Lösung der diskreten Poisson-Gleichung

Die Galerkin-Matrix der FE-Diskretisierung der Poisson-RWA (3.1) istsymmetrisch und positiv-definit (gilt nicht notwendig für andere Dgln.)

Allen FE-Diskretisierungen gemeinsam: die Galerkin-Matrix ist dünn be-setzt (sparse). Wir betrachten die Matrix aus Beispiel 2 (L-Gebiet).

bilinear biquadratisch

N h # Einträge 6= 0 Anteil h # Einträ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 Lösung der diskreten Poisson-Gleichung TU Bergakademie Freiberg, WS 2010/111


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



−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



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 früh wie möglich während der Assemblierung vorgenommen.

• Weiterentwicklungen auch heute noch relevant, auch als eigenständigeEliminationsverfahren, z.B. UMFPACK.

• Faustregel: in 2D reichen ausgereifte direkte, auf dünn-besetzte Ma-trizen spezialisierte Löser aus. Für 3D-Diskretisierungen muss aufIterationsverahren zurückgegriffen werden.

• Vorteil direkter Verfahren: typischerweise als ”black box“ einsetzbar.Bei Iterationsverfahren ist deutlich mehr Expertise des Anwenderserforderlich.



4.1 Erste Iterationsverfahren

Zunächst sei allgemein gegeben

Ax = b, A ∈ Cn×n invertierbar, b ∈ Cn. (4.1)

Grundidee:

Startnä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) Lösung,

em := x? − xm Fehler.

4.1 Erste Iterationsverfahren TU Bergakademie Freiberg, WS 2010/111


Die ersten Iterationsverfahren (Handrechnung!) beruhen auf Umstellungeinzelner Gleichungen, Auflö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 lässt sich halb imSchlafe ausführen, oder man kann während desselben an andere Dingedenken.“ Brief von Gauß an Gerling, 1823

• Carl Gustav Jacobi (∼ 1846)

• Ludwig Seidel (∼ 1862)

• Christian August Nagel (Sächsische Triangulation, 1867-1886) 159Unbekannte



Formale Beschreibung von Relaxationsverfahren durch Zerlegungen (split-tings):

A = M −N , M invertierbar, (4.2a)

überführt (4.1) in Fixpunktform Mx = Nx +b, zugehö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, Startnäherung x0

1 for m = 0, 1, 2 . . . bis Abbruchkriterium erfüllt do2 rm ← b −Axm3 Löse Mh = rm4 xm+1 ← xm + h



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.



Fü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 Fü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 für alle Startnäherungen x0, wenn für dieIterationsmatrix T = M−1N gilt ρ(T ) < 1.



4.2 M -Matrizen und regulä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 esfü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 reguläre Zerlegungen TU Bergakademie Freiberg, WS 2010/111


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 gehö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 zusätzlich A 6= B , so giltρ(B) < ρ(A).



Definition 4.7 Eine nichtsinguläre Matrix A ∈ Cn×n heißt M-Matrixa, falls

ai,j ≤ 0 für i 6= j und A−1 ≥ O .

Eine Hermitesche M-Matrix heißt Stieltjes-Matrix.

Satz 4.8 Für B ∈ Cn×n mit ρ(B) < 1 ist I −B nichtsingulä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 eingeführt, zu Ehren von Her-mann Minkowski.



Satz 4.10 Eine Matrix A ∈ Cn×n mit ai,j ≤ 0 für i 6= j ist genau dann eineM-Matrix, wenn D := diag(A) > O (insbesondere sind M -Matrizen stetsreell) und wenn für die Matrix B := I −D−1A gilt ρ(B) < 1.

Definition 4.11 Die Zerlegung A = M − N der Matrix A ∈ Rn×nheißtreguläre Zerlegung von A, wenn M−1 ≥ O und N ≥ O gelten.

Satz 4.12 Ist A = M−N eine reguläre Zerlegung von A und ist A−1 ≥ O ,so gilt

ρ(M−1N ) =ρ(A−1N )

1 + ρ(A−1N )< 1,

d.h. jede reguläre Zerlegung einer Matrix mit nichtnegativer Inversen führtzu einem konvergenten Iterationsverfahren. Diese Aussage gilt insbeson-dere für M-Matrizen.



Satz 4.13 Für A ∈ Rn×n mit A−1 > O seien A = M1 −N1 = M2 −N2zwei regulä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.



Satz 4.14 (Starkes Zeilensummenkriterium) Gilt fü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 fü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 zusätzlich irreduzibel, so konvergiert das Gauß-Se

Finite Elemente I - tu-freiberg.deernst/Lehre/FEM_I/Folien12/...Finite Elemente I 2 1 Einleitung...

Documents

Transcript of Finite Elemente I - tu-freiberg.deernst/Lehre/FEM_I/Folien12/...Finite Elemente I 2 1 Einleitung...