Einführung in das Wissenschaftliche Rechnen

nationales Forschungszentrum in der Helmholtz-Gemeinschaft

Institut für Angewandte und Numerische Mathematik

KIT - Universität des Landes Baden-Württemberg und

nationales Forschungszentrum in der Helmholtz-Gemeinschaft

www.kit.edu

Einführung in das Wissenschaftliche Rechnen

Skript zur Vorlesung im Sommersemester 2019

Christian Wieners

Version vom 6. Februar 2020

Inhaltsverzeichnis

Einführung 2

1 Die Potentialströmung 5

2 Finite-Elemente-Approximation 9

3 Lösen von dünn besetzten LGS 15

4 Gemischte und hybride Finite Elemente 28

5 Finite Volumen für die lineare Transportgleichung 39 5.1 Die lineare Transportgleichung . . . 39 5.2 Das Finite-Volumen-Verfahren . . . 41 5.3 Discontinuous-Galerkin-Verfahren . . . 46

6 Konvektions-Reaktions-Diffusions-Systeme 51

6.1 Die Konvektions-Reaktions-Diffusions-Gleichung . . . 51 6.2 Linear-implizite Verfahren (Rosenbrock-Typ) . . . 60 7 Discontinuous Galerkin-Methode für die Diffusionsgleichung 64 8 UQ (Uncertainty-Quantification) für die Diffusionsgleichung 68

Einführung

• Wie wird Wasserverteilung im Boden beschrieben?

• Wie wird die Ausbreitung von Verschmutzung beschrieben?

• Wie wird ein biologischer Abbauprozess beschrieben?

• Wie werden Unsicherheiten (Materialparameter) quantifiziert?

• Wie können die Parameter gemessen werden?

• Wie kann der Prozess (etwa durch Abpumpen) gesteuert werden? Vereinfachung möglich?

Numerische Methodenentwicklung:

• Wie können die Modelle approximiert werden?

• Ist das approximierte Modell eindeutig lösbar?

• Konvergieren die Approximationen?

• Wie berechnet man die Approximation effizient?

• Wie werden die Algorithmen realisiert?

• Wie werden die Ergebnisse dargestellt?

• Wie kann man die Effizienz/Qualität der Approximation/Algorithmen quantifizie- ren?

Struktur der Lösung einer Praktikumsaufgabe:

1. Problemstellung

2. Mathematisches Modell

3. Konfiguration eines Beispiels (Daten) 4. Numerische Methoden

5. Visualisierung der Ergebnisse 6. Qualitative Auswertung

7. Diskussion der Ergebnisse (Insgesamt 1-2 Seiten) Ergänzung B.Sc.-Arbeit (Literaturrecherche)

• Eigenschaften des Modells

• Eigenschaften der numerischen Verfahren (Insgesamt 40-60 Seiten) Anforderung Promotion: neue Methode / neues Ergebnis.

1 Die Potentialströmung

Problem

Ω ⊂ R³ beschreibe eine poröse Bodenschicht. Von oben sickert Regenwasser langsam nach unten ins Grundwasser. Ziel ist die Bestimmung der transportierten Wassermenge sowie deren Verteilung.

Modell

Sei q: Ω × [0, T] −→ R³ der Fluss mit divq = 0 und %: Ω ×[0, T] −→ R≥0 die Massendichte, d. h. für jedes KontrollvolumenY ⊂Ωist durch

Z

Y

%(x, t) dx

die Wassermenge inY zum Zeitpunkttgegeben. Die zeitliche Änderung der Wassermenge inY ergibt sich durch den Zu- und Abfluss am Rand vonY. Dieser Zusammenhang lässt sich mithilfe der Bilanzgleichung

d dt

Z

Y

%(x, t) dx=− Z

∂Y

%(x, t)q(x, t)·nda .

mathematisch formulieren. Dabei bezeichne ∂Y den Rand von Y und n die nach außen orientierte Einheitsnormale auf∂Y. Der Gaußsche Integralsatz liefert

d dt

Z

Y

%(x, t) dx=− Z

∂Y

div(%q) da .

Damit erhält man die partielle Differentialgleichung

∂_t%+ div(%q) = 0.

Im Fall ∂_t% ≡ 0 spricht man von einer stationären Strömung. Ist außerdem die Dichte

%≡%0 konstant, vereinfacht sich die partielle Differentialgleichung zudivq= 0.

(1.1) Lemma

Seidivq = 0 inΩ. DefiniereΓ_in = {x ∈∂Ω : q(x)·n(x) <0}undΓ_out = {x ∈∂Ω : q(x)·n(x)>0}. Dann gilt

Z

Γin

q·nda=− Z

Γout

q·nda

Beweis.

0 = Z

Ω

divqdx= Z

∂Ω

q·nda= Z

Γin

q·nda+ Z

Γout

q·nda

Darcy-Gesetz Das Darcy-Gesetz

q =−κ(∇p−G)

beschreibt den Fluss eines Fluids in einem porösen Medium, dabei sind p: Ω−→R hydrostatischer Druck

κ: Ω −→R^3×3sym Permeabilitätstensor G= (0,0, %₀g₀)^T Gravitation.

Ist die Permeabilität richtungsunabhängig, d. h.κ =κ₀I, spricht man von isotropem Ver- halten. Wenn zusätzlichκ ≡κ₀ konstant ist, dann ist die Permeabilität ortsunabhängig (=

homogen).

Mit der Hilfsgrößeu(x) =−p(x)+%₀g₀x₃vereinfacht sich das Darcy-Gesetz zuq =κ∇u.

Damit ergibt sich die partielle Differentialgleichung divκ∇u= 0.

Bemerkung

Im Fall ortsunabhängiger Permeabilität entspricht das der Laplace-Gleichung (div∇u=)∆u= 0.

Randwertaufgabe

Sei∂Ω =Γ_N∪Γ_D. Durch die Vorgabe von Funktionen u_D: Γ_D−→R (Dirichlet-Werte)

g_N: Γ_N−→R (Neumann-Werte) ergibt sich die Randwertaufgabe: Bestimmeu: Ω−→Rmit

divκ∇u = 0 inΩ u = u_D aufΓ_D κ∇u·n = g_N aufΓ_N. 2D-Reduktion

Sei Ω2D = {(x1, x3) ∈ R² : (x1,0, x3) ∈ Ω3D} ein Querschnitt von Ω2D. Die Lösung sei lokal bei Ω₀ symmetrisch um x₂ = 0. Dann gilt ∂₂u_3D = 0. Mit u_2D(x₁, x₃) = u_3D(x₁,0, x₃)folgt

divκ_2D∇u_2D = 0,

wobei κ_2D(x₁, x₃) ∈ R^2×2sym. Man kann also allgemein von Ω ⊂ R^D und κ ∈ R^D×Dsym

ausgehen, wobeiDdie Werte2und3annehmen kann.

Beispiel

SeiΩ = (0,1)², Γ_D = [0,1]× {0}, Γ_N=∂Ω\Γ_D, u_D ≡0, κ=Iund g_N(x) =

(−1 fallsx₃ = 1

0 sonst .

Die exakte Lösung der zugehörigen Randwertaufgabe lautetu(x₁, x₃) =−x₃, denn dann gilt

∆u= 0, sowie

∇u·n =







−1 fallsx₂ = 1 1 fallsx₂ = 0 0 sonst

.

Schwache Formulierung (1.2) Lemma

a) Seidivq = 0inΩ, dann gilt für jede Testfunktionφ: Ω−→Rmitφ= 0aufΓ_D Z

Ω

q· ∇φdx= Z

∂Ω

q·nφda

b) Sei

Z

Ω

q· ∇φdx= Z

∂Ω

gφda

für alle Testfunktionen. Dann giltdivq = 0inΩundq·n=g auf∂Ω.

Beweis. a)

div(φq) =∇φ·q+φdivq

⇒ Z

∂Ω

φq·nda= Z

Ω

∇φ·qdx+ Z

Ω

φdivqdx

b)

Z

∂Ω

φ(q·n−g) da= Z

Ω

φdivqdx

⇒divq= 0inΩ

⇒q·n−g = 0auf∂Ω.

Anwendung auf Randwertaufgabe Starke Form:

divq= 0, q=−κ∇uinΩ

−q·n=g_NaufΓ_N,u=u_DaufΓ_D Schwache Form:

Z

Ω

κ∇u· ∇φdx= Z

ΓN

g_Nφda

∀φmitφ = 0aufΓ_D

2 Finite-Elemente-Approximation

Um numerisch eine näherungsweise Lösung der schwachen Formulierung der Potential- strömung zu bestimmen, werden ein endlich-dimensionaler Funktionenraum V_h konstru- iert und einu_h ∈V_h mitu_h ≈u_DaufΓ_Dund

Z

Ω

κ∇uh· ∇φhdx= Z

ΓN

gNφhda ∀φh ∈Vhmitφh = 0aufΓD (1) bestimmt.

Wir betrachten im Folgenden den Spezialfall, dass Ω ⊂ R² ein Polygongebiet ist. Die hier vorgestellte Vorgehensweise lässt sich auf dreidimensionale Problemstellungen über- tragen. Es sei eine ZerlegungΩ = S

K∈KK vonΩ in Dreiecke oder Vierecke gegeben.

Diese sogenannten Zellen werden über

K = conv(V_K) mit

V_K =

({zK,0, zK,1, zK,2} ⊂R² Dreieck {zK,0, zK,1, zK,2, zK,3} ⊂R² Viereck beschrieben. Die Menge aller Eckpunkte der Zerlegung wird mit

V_K = [

K∈K

V_K ={z₁, . . . , z_N} notiert. Im Folgenden sei die ZerlegungKzulässig.

(2.1) Definition

Eine ZerlegungKheißt zulässig, wenn für alleK, K⁰ ∈ K conv(VK∩ VK⁰) = convVK ∩convVK⁰

gilt, d. h. wenn der Durchschnitt von zwei (verschiedenen) Zellen leer, eine Ecke oder eine Kante ist.

Idee

Wähle Knotenwerteu(z)für z ∈ VK und interpoliere die Werte in jedemK ∈ K linear (bei Dreiecken) oder bilinear (bei Vierecken). Um effizienter rechnen zu können, wird zunächst eine Referenzzelle festgelegt.

Referenzdreieck

Bei einer Zerlegung in Dreiecke wird das Referenzdreieck Kˆ = conv

ˆ z₀ =

0 0

,zˆ₁ =

1 0

,zˆ₂ =

0 1

betrachtet. Die affin-lineare Koordinatentransformationϕ_K: ˆK −→Ksei durch ϕ_K(ξ) = z_K,0+B_KξmitB_k = Dϕ_K

gegeben und liefert für jedesK ∈ K den Übergang vom ReferenzdreieckKˆ aufK. Auf dem ReferenzdreieckKˆ wird der Ansatzraum

Vˆ = span{λˆ0,λˆ1,λˆ2}=P¹( ˆK) mit

λˆ₀(ξ) = 1−ξ₁−ξ₂ λˆ₁(ξ) = ξ₁

λˆ₂(ξ) = ξ₂ definiert. Nach Konstruktion gilt dann

λˆ_i(ˆz_j) =

(1 i=j 0 i6=j .

Für den Ansatzraum auf dem DreieckK ∈ Kergibt sich über die Koordinatentransforma- tionϕ_K

V_K = span{λˆ_i◦ϕ⁻¹_K :i= 0,1,2}=P1(K).

Die Kroneckereigenschaft

λˆ_i◦ϕ⁻¹_K (z_j) =

(1, i=j 0, i6=j bleibt so erhalten.

Referenzviereck

Bei einer Zerlegung in Vierecke werden auf die gleiche Weise ein Referenzviereck Kˆ = conv

ˆ z0 =

0 0

,zˆ1 =

1 0

,zˆ2 =

1 1

,zˆ3 =

0 1

sowie die zugehörige KoordinatentransformationϕK: ˆK −→K ϕ_K(ξ) =z_K,0+B_K(ξ)mitB_K(ξ) = Dϕ_K(ξ) definiert. Mit den Ansatzfunktionen

λˆ₀₀(ξ) = (1−ξ₁)(1−ξ₂) λˆ₀₁(ξ) = (1−ξ₁)ξ₂ λˆ10(ξ) = ξ1(1−ξ2) λˆ₁₁(ξ) = ξ₁ξ₂

ergeben sich wie eben die lokalen Ansatzräume Vˆ = span{λˆ₀₀,ˆλ₀₁,λˆ₁₀,ˆλ₁₁} undV_K = span{λˆij ◦ϕ⁻¹_K : i, j = 0,1}.

Legt man zu jedemK ∈ Kmithilfe einer Abbildung

µ_K: IKˆ −→I_h ={1, . . . , N}

und

IKˆ =

({0,1,2} Dreieck {0,1,2,3} Viereck

die Zuordnung zwischen den Ecken von K und denen der Referenzzelle Kˆ fest, ergibt sich inK der Zusammenhang

φ_k(x) = ˆλ_i◦ϕ⁻¹_K (x) mitµ_K(i) = k. Folglich gilt inK

Dφ_k(x) = ˆDˆλ_i(ˆx)Dϕ⁻¹_K (x)mitϕ_K(ˆx) =x bzw.

∇φ_k(x) =B_K^−T(ˆx) ˆ∇λˆ_i(ˆx).

Finite-Elemente-Approximation

Zu einer gegebenen ZerlegungKmit der maximalen Kantenlängeh= maxK∈Kdiam(K) werden die Finite-Elemente-Ansatzräume

V_h ={φ_h ∈C(Ω) :φ_h|_K ∈V_K ∀K ∈ K}

und

V_h(u_D) = {φ_h ∈V_h: φ_h(z) = u_D(z), z ∈ VK∩Γ_D}

definiert. Aus der schwachen Formulierung (1) ergibt sich damit die Finite-Elemente- Formulierung: Bestimmeu_h ∈V_h(u_D)mit

X

K∈K

Z

K

κ∇u_h· ∇φ_hdx= Z

ΓN

g_Nφ_hda ∀φ_h ∈V_h(0).

Unter Verwendung der Knotenbasis{λ₁, . . . , λ_N}vonV_hmit ˆλ_n(ˆz_k) =

(1 n=k 0 n6=k

erhält man die äquivalente algebraische Formulierung: Sei I = {1, . . . , N} und I_D = {h ∈ I: z_h ∈ Γ_D} die Indexmenge der Knoten, die auf dem Dirichletrand Γ_D liegen.

Bestimmeu∈R^N mitu_h =P

u_nλ_n ∈ V_h(u_D)und

N

X

n=1

u_nX

K∈K

Z

K

κ∇λ_n· ∇λ_kdx= Z

ΓN

g_Nλ_kda, ∀k ∈ I \ I_D.

Wir definieren nun die SteifigkeitsmatrixA ∈R^N×N und den Lastvektorb∈ R^N über

A[k, n] =









 P

K∈K

R

K

κ∇λ_n· ∇λ_kdx, k /∈ I_D

1, n =k ∈ I_D

0, n 6=k ∈ I_D

und

b[k] =





 R

ΓN

gNλkda k /∈ ID

u_D(z_k) k ∈ I_D . Damit ergibt sich das zu lösende lineare Gleichungssystem

A u=b.

Nach Konstruktion gilt u_k = u_D(z_k) für k ∈ I_D. Zur Bestimmung von A wird auf den Zellen der Zerlegung K und den Basisfunktionen λ_n mit Einschränkungen λ_K,n auf K gerechnet:

A[n, k] = Z

Ω

κ∇λ_K,k· ∇λ_K,ndx= X

K∈K

Z

K

κ∇λ_K,k· ∇λ_K,ndx.

Den Vorgang des Zusammenfügens vonAbezeichnet man alsAssemblieren. Die Integra- tion auf jeder Zelle wird über

X

K∈K

Z

K

κ∇λK,k· ∇λK,ndx= X

K∈K

Z

Kˆ

det DϕKκ∇λK,k◦ϕK· ∇λK,n◦ϕKdˆx

auf die Referenzzelle transformiert. Dort wird mit Hilfe eines Quadraturverfahrens nume- risch integriert. Fallsκ|_K für alleK ∈ Kkonstant ist, ergibt sich

A[n, k] = X

K∈K

X

ξ∈Ξ

w_ξdetB_K(ξ)κ◦ϕ_KB_K^−T(ξ) ˆ∇ˆλ_K,k(ξ)·B_K^−T(ξ) ˆ∇λˆ_K,n(ξ).

Für das Referenzdreieck werden Ξ =

1 3 1 3

und ω_ξ = 1 2 verwendet, für das Referenzviereck

Ξ = α

α

,

1−α α

,

α 1−α

,

1−α 1−α

mitα=

√2−1 2√

2 undω_ξ = 1 4. Algorithmus

A0) A= 0,b = 0 A1) ∀K ∈ K:

A = A+R^T_KA_KR_K, b = b+R^T_Kb_Kd.h.A[r_k,n, r_k,k] = A[r_k,n, r_k,k] +A_K[n, k], b[rk,n] =b[rk,n] +b_K[n]

A2) ∀k ∈ I_D :b[k] =u_D(z_K), A[k, k] = 1, A[k, n] = 0fürn 6=k A3) ∀k ∈ I_D, n6=k :b[n] =b[n]−A[n, k]u_D(z_k)

(2.2) Lemma

Seiκuniform symmetrisch positiv definit, d. h. κy·y ≥ κ₀|y|² ∀y ∈ R^D mitκ₀ >0und seiΓ_D∩ Vh6=∅. Dann istAregulär.

Beweis. Per Konstruktion gilt:Aist symmetrisch positiv definit, also zeige A v·v = 0⇔v = 0

Sei alsoA v·v = 0. Außerdem gilt A v·v =X

K

X

ξ

w_ξ|detB_K|κ∇v(ϕ_K(ξ))· ∇v(ϕ_K(ξ))

≥κ₀X

K

X

ξ

w_ξdetB_K|∇v(ϕ_K(ξ))|²

=κ₀ Z

Ω

|∇v|²dx

⇒ ∇v = 0inΩ

⇒v ≡const inΩ

Mit z_K ∈Γ_D∩ V_h ⇒ v[k] = 0 ⇒ v(z_k) = 0 ⇒ v ≡0.

Bemerkung

Im Folgenden seiA₀ ∈R^N×N gegeben durch

A₀[n, k] =











A[n, k] n, k /∈ ID

1 n=k∈ ID

0 n6=k∈ ID

0 k6=n∈ ID

.

Nach Satz 2.2 istA₀ positiv definit. Zerlegt manuin u=u₀+u_D ergibt sich das äquivalente Gleichungssystem

A₀u₀ =b−Au_D.

Man kann also ohne Einschränkung annehmen, dassApositiv definit ist.

3 Lösen von dünn besetzten LGS

Ein Vorteil der finiten Elemente Methode ist, dass die auftretenden Matrizen in der Regel dünnbesetzt sind. Deswegen ist es von Vorteil spezielle Lösungsverfahren für dünnbesetz- te Gleichungssysteme zu verwenden.

Betrachte eine Familie von linearen GleichungssystemenAu=binR^N. SeiI ={1, ..., N} die Index- Menge,|I|= N undG= G_A ={(n, k) ∈ I × I : A[n, k]6= 0}der Matrix- graph.

(3.1) Definition

Eine Matrix heißt dünn besetzt , wenn|G_A|=O(N). Sonst heißt sie voll besetzt.

(3.2) Lemma

SeiI_K ={k ∈ I :K ⊂suppφ_k}={k ∈ I :z_k ∈ V_K}.

a) Dann gilt für den Graph der FE- Matrix

G⊂[

I_K× I_K und |G| ≤N(max

k∈I |{K :z_k∈ V_K}|+ 1).

Bemerkung (D=2): Seiγ₀ der kleinste Innenwinkel aller ZellenK. Dann ist|G| ≤ (^2π_γ

0 + 1)N.

b) Das Assemblieren benötigtO(N)Operationen.

c) Die MultiplikationAubenötigtO(N)Operationen.

d) Das compressed row storage format benötigt O(N) Speicher: Sei M = |G| die Anzahl der Matrixeinträge ungleich Null. Sei außerdem

a[m], m= 1..., M, Vektor der Matrixeinträge 6= 0 c[m]in{1, ..., N}Spaltenindex vona[m]

d[n]∈ {1, ..., M + 1}, n= 1, N+ 1mitd[1] = 1, d[N + 1] =M + 1Diagonalpointer.

Dann werden die Einträge der MatrixAfolgendermaßen gespeichert:

A[n, c[m]] = a[m] fürm=d[n], ..., d[n+ 1]−1, n= 1, ..., N

⇒(Au)[n] =

d[n+1]−1

X

m=d[n]

a[m]u[c[m]].

Assemblieren

SeiAK[i, j] = Z

K

K∇(ˆλj ◦ϕ⁻¹_K ) dx

V_K = {z_µK(i) : 1≤i≤V_K} ⊂ V, V_K =|V_K|

A = X

K

R^T_KA_KR_K mitR_K ∈R^VK×V

R_K[i, n] =

(1 µ_K(i) =n

0 sonst .

Nested Dissection SeiΩ= S

p∈P

Ω^p,P ={1, ..., P}, Ω^p∩Ω^q =∅mitp, q ∈ P, p6=qeine nichtüberlappende Gebietszerlegung(Ω^p∩Ω^q6=∅)mit SkeletonΓP =Ω∩(∪∂Ω^p).

SeienI^p ={n:z_n ∈Ω^p\ΓP}Indizes zu inneren Knotenpunkten inΩ^p undI_Γ ={n ∈ I :zn ∈ΓP}Interface- Indizes. Dann gilt:A[n, k] = 0fürn ∈Ω^p, k ∈Ω^q(p6=q).

A=











A₁₁ 0 A_Γ1

. .. ... 0 A_pp A_{Γ p} A_1Γ · · · A_pΓ A_{Γ Γ}











Aufwand in 2D:O((^N_P)²) +O((P(^N_P)¹²)³), fallsp < √ N.

Das folgende Lemma zeigt, dass das Schurkomplement die Invertierbarkeit der ursprüng- lichen Matrix erbt.

(3.3) Lemma

WennAsymmetrisch positiv definit ist, dann ist das Schurkomplement S =A_{Γ Γ} −

P

X

p=1

A_{Γ p}A⁻¹_ppA_pΓ ebenfalls symmetrisch positiv definit, und es gilt füru=A⁻¹b:

u_Γ = S⁻¹(b_Γ −X

A⁻¹_ppA_pΓb_p), u_p = A⁻¹_pp(b_p−A_pΓu_Γ).

Beweis. AusA_ppu_p =b_p−A_pΓu_Γ folgt A_{Γ Γ}u_Γ = b_Γ −X

A_{Γ p}u_p =b_Γ −X

A_{Γ p}A⁻¹_pp(b_p−A_pΓu_Γ)

⇒Su_Γ = b_Γ −X

A_{Γ p}A⁻¹_ppb_p. Seiv_Γ 6= 0undv_p =−A⁻¹_ppA_pΓv_Γ

⇒SvΓ ·vΓ =AΓ ΓvΓ ·vΓ −X

(A⁻¹_ppApΓvΓ)(ApΓ ·vΓ)

=A_{Γ Γ}v_Γ ·v_Γ −X

v_p(A_ppv_p)

=Av·v

⇒Sv_Γ ·v_Γ =Av·v = 0

⇒v = 0 ⇒ v_Γ = 0

⇒S ist positiv definit.

Nun wähle

P = 2^S, Ω^0,1 =Ω Ω^1,1∪Ω^1,2 =Ω rekursiv

Ω^s,q = Ω^s+1,2q−1∪Ω^s+1,2q, q= 1, ...,2^S Ω^p = Ω^S,p, p= 1, ...,2^S =P.

Nested Dissection: Rekursive Anwendung von Lemma 3.3 fürs=S−1, S−2, ...,0.

Beispiele für Partitionierungsverfahren

Sobald die Finite Elemente Methode parallel ausgeführt wird, was in der Praxis der Regel entspricht, muss das zugrunde liegende Gebiet auf die einzelnen Prozesse verteilt werden.

Dies geschieht mit Hilfe eines Partitionierungsverfahrens.

A) Rekursive Koordination - Bisektion (RCB) SeizK = _|K|¹ P

z∈VK

z ∈R²Zellenmittelpunkt Definiere

z <₁ y⇔z[1]< y[1]oderz[1] =y[1]undz[2]< y[2]) z <₂ y⇔z[2]< y[2]oderz[2] =y[2]undz[1]< y[1]) s= 0: I⁰ = 1, ..., N ,VK =z₁, ..., z_N fürN =|K|, z_n=z_Kn

I⁰ =I^1,1∪ I^2,2 mitz <₁ y fürz ∈ I^1,1, y ∈ I^1,2

s= 2: I⁰ =I^2,1∪ I^2,2∪ I^2,3∪ I^2,4 mit I^s−1,t =I^s,2t−1∪ I^s,2tundz <2 y

Fahre nun mit abwechselnder Sortierung rekursiv fort, bis das Gebiet weit ge- nug zerlegt ist. Das Verfahren ist sehr schnell, allerdings nicht invariant gegen- über Gebietstransformationen und in ungünstigen Fällen wird das InterfaceI_Γ sehr groß.

B) Rekursiv Trägheits- Bisektion (RIB)

Definiere diskreten Laplace- Operator zum Graphen

ΓK = {(K, K⁰) :∂K∩∂K⁰ =f gemeinsame Seite} ⊂ K × K d_K = |{K⁰ 6=K : (K, K⁰)∈ΓK}|

L[K, K⁰] =







d_K K =K⁰

−1 (K, K⁰)∈ΓK

0 sonst

L·e= 0 mite= (1, ...,1)^> ∈R^NK ⇒symmetrisch, positiv definit

Wir approximieren den Eigenvektor w ∈ W = v ∈R^NK :v·e= 0 zum kleinsten Eigenwertµ >0vonLw =µw.

w·e= 0 ⇒Definiere

K^1,1 ={K :w[K]>0},K^1,2 =K \ K^1,1. Inverse Iteration zur Eigenwert- Berechnung:

I0) Wählew⁰ ∈R^NK \ {0, e}, k = 0, ε >0

I1) Approximierew^k≈(L|w)⁻¹w^k−1durch einige cg-Schritte.

I2) berechnew^k:= _|w¹k|w^k, µ_k =Lw^k·w^k I3) falls|µ_k−µk−1|< εSTOP

I4) gehe zu I1)

Dieses Verfahren ist invariant gegenüber Gebietstransformationen.

C) Space-Filling Curve

Definiere eine Kurveγ : [a, b]→Ω, die jede Zelle genau einmal durchläuft. Damit wird eine Nummerierung inKdefiniert. Setze

I^q ={K_[(q−1)^M

p]+1, ..., K_[q^M

p]}.

Damit ensteht eine hierarchische Konstruktion, die auf unregelmäßige Gitter er- weiterbar ist. Allerdings ist das Verfahren aufwendig und nicht optimal.

Netzgenerierung

Wenn ein GebietΩgegeben ist, muss dieses erst in Zellen zerlegt werden, bevor die Finite Elemente Methode anwendbar ist. Da nicht jedes Gebiet eine regelmäßige Struktur hat, sind verschiedene Methoden zur Netzgenerierung nötig.

A) Strukturierte Gitter

1) Wähle BoxΩ ⊃[α₁, β₁]×[α₂, β₂].

Wähleh >0und GitterΓh =hZ²∩Ω.

Konstruiere regelmäßiges Vierecksgitter zuΓ_h.

2) Wähle weitere Punkte auf dem Rand ∂Ω und konstruiere (unregelmäßige) Vierecke am Rand.

3) Verschiebung einiger Knoten, um gleichmäßigere Vierecke zu erhalten (bzw.

Einfügen von Dreiecken).

Dieses Verfahren hat den Vorteil, dass es schnell durchführbar ist, allerdings ist es nicht geeignet den Rand gut zu approximieren.

B) Advancing Front

1) Wähle (gleichmäßige) Punkte auf dem Rand.

2) Konstruiere eine Elementschicht entlang dem Rand, so dass ein neuer innerer Rand entsteht. Wiederhole diesen Prozess so lange die inneren Ränder nicht zusammenstoßen.

3) Verbinde die inneren Ränder durch Elemente.

4) Glätten der entstandenen Zellen.

Bei diesem Verfahren entstehen bessere Elemente am Rand, allerdings ist die Kolli- sionskontrolle im Inneren aufwendig.

C) Voronoi- Konstruktion

1) Wähle gleichmäßig verteilte Punkteznauf dem Rand und im Inneren.

2) Bilde Voronoi-Zellen

C_n ={z∈Ω :|z−z_n| ≤ |z−z_k| ∀k 6=n}konvex!

3) Bilde Kanten:

E ={(z_n, z_k) :n 6=k, C_n∩C_kgemeinsame Seite}

4) Bilde zugehöriges Dreiecksnetz (Tetraedernetz).

Der Vorteil dieses Verfahrens sind die guten Netze die entstehen, allerdings ist man auf Dreiecke eingeschränkt.

Iterative Lösungsverfahren

SeiB ∈R^N×N ein Vorkonditionierer zuA=L+D+RmitB ≈A⁻¹, z.B.

A) Jacobi,B_Jac =D⁻¹

B) Gauß-SeidelB_GS = (L+D)⁻¹ C) Sym. Gauß-Seidelc=B_SGSrmit

c₁ = B_GSr r1 = r−Ac1

c = c₁+B_GS^> r₁

⇒B_SGS = B_GS+B_GS^> −B_GS^>AB_GS D) Parallele Vorkonditionierer zuA= P

p∈P

Ap

WähleB^p zuRpAR^>_p. Definierec=BPrmit c_p = A⁻¹_p R_pr

c = X

R^>_pc_p

E) Zweigitter- Verfahren mitA=A₁, A₀ =R₀A₁R^>₀ und

RestriktionenR₀: R^N →R^N0 mitN₀ < N und ProlongationR^>₀ : R^N0 →R^N Definierec=B_ZGrdurch:

Glätten:

c₁ = B_GSr r₁ = r−Ac₁ r₀ = R₀r₁ Grobgitterkorrektur:

c₀ = A⁻¹₀ r₀ c = c₁+R^>₀c₀ Symmetrische Variante: Nachglätten mitB_GS^>.

F) Mehrgitter:A=A_L, AL−1 =RL−1A_LR^>_L−1, ...rekursiv glätten.

G) Least Squares:

1) Wähleε > 0als Abbruchbedingung undΘ > 0 als Reduktionsfaktor. Seiu₀ gegeben. Setze

k = 0, r⁰ =b−Au⁰, ϕ₀ =|r⁰| 2) falls%k<max{Θ%0, ε}STOP

3) Wähle

c^k−1 = Br^k r^k+1 = r^k−Ac^k u^k+1 = u^k+c^k 4) k:=k+ 1, gehe zu 2)

Es gilt

r^k = (I−BA)^kr⁰

⇒ %k≤ kIN −BAk^k%0

⇒ FallskI_N −BAk<1, so konvergiert LS FallsI_N −BAnormal und| · |=| · |₂, so gilt

%k

%₀ →%(I_N −BA) Krylov- Raum- Verfahren

Krylov- Raum- Verfahren sind eine Klasse von iterativen Verfahren zum Lösen von dünn- besetzten Gleichungssystemen.

Berechne eine ONB im Krylov- RaumK_k(BA, Br) = {P(BA)Br : P ∈ Pk−1}, wobei A, B symmetrisch, positiv definit sind. Es gilt in der Energienorm|u|_A=√

Au·u:

1) Wähleε,Θ>0, gegebenx⁰ Setze

k = 0, r⁰ =b−Ax⁰, ϕ₀ = 1, p⁰ =u⁰ 2) Falls|r^k|<max{Θ|r⁰|, ε} STOP

3) Setze

c^k = Br^k

%₁ = c·r p^k := %₁

%0

p^k p^k+1 := p^k+c^k

%₀ = %₁ q^k+1 = Ap^k+1

α = %₀

p^k+1·q^k+1 u^k+1 = u^k+α_kp^k+1

r^k+1 = r^k−αq^k+1 4) k :=k+ 1, gehe zu 2)

Es gilt in der Energienorm|u|_A =√

Au·u:

|u^k−u⁰|_A≤2 √

κ−1

√ κ+ 1

k

|u¹−u⁰|_Amitκ =κ(BA).

Das Folgende wurde im SS’15 in der Vorlesung notiert nach Space-Filling Curve:

Subspace Correction Methode

SeiA ∈R^N×N, N_p < N, p∈1, ..., P , R_p ∈R^Np×N, A_p =R_pAR^>_p. Beispiel:

A ist FE-Matrix, sei weiter N = |V_K|, V_h = span{λ₁, ..., λ_N},V_K = {z₁, ..., z_N} =

∪ VK. Sei zudemK= S

p∈P

K^pnicht überlappend undV^P = S

k∈K^P

V_K überlappend.

RP[n, k] ={zP,1, ..., zP,N_P}

=

(1 z_P,n=z_K 0 sonst Nun wähleB_P ≈A⁻¹_P .

Beispiel:

1) B_P =A⁻¹_P LU-Zerlegung

2) B_P = (diagA_P)⁻¹ Jacobi-Verfahren

3) B_P = (diagA_P + lower A_P)⁻¹ Gauß-Seidel Parallel Subspace Correction

P0) Wähleu⁰ ∈R^N (u⁰[n] =u_D(z_n), z_n∈Γ_D)

Berechner⁰ =b−Au⁰, setze k= 0, und wähle ε >0, θ >0

P1) Falls|r^k|< ε STOP

P2) ∀p: berechne

c_p =B_PR_Pr^k und c^k =X R^>_Pc_P P3)

u^k+1 =u^k+θc^k r^k+1 =r^k−θAc^k

P4) Setzek :=k+ 1, gehe zu P1).

Sucsessive Subspace Correction

P2)* Setzer˜=r^k für p= 1, ..., P : c_P =B_PR_Pr˜ dann ˜r:= ˜r−AR^>_Pc_P. Symmetrisierung

Für ungeradekersetzeB_P durchB_P^>.

⇒fallsA, B_P positiv definit: CG-Verfahren.

Poisson Problem

Sei Ω ⊂ R^N beschränkt. Suche u : Ω → R, sodassu ∈ C²(Ω)∩C(Ω)und∂Ω = Γ_D∪Γ_Nmit

−∆u=f inΩ

u=gaufΓ_D (P)

∇un =haufΓ_N wobeif ∈C(Ω) und g, h∈C(∂Ω) gegeben.

Spezialfall: Fürf = 0 ⇒ Laplace-Gleichung

Schwache Formulierung

u∈L₁(Ω) ist schwache Lösung von (P), falls Z

Ω

−∆uvdx= Z

Ω

∇u· ∇vdx+ Z

∂Ω

∇un·vda

= Z

Ω

∇u· ∇vdx+ Z

ΓN

∇un·vda

= Z

Ω

f·vdx

∀v ∈ V(0) wobeiV(0) =H_0,Γ¹

D :={v ∈V :v = 0 auf Γ_D} mit v =H¹(Ω) :={v ∈ L₂(Ω) :v besitzt schwache Ableitung ∂_iv ∈L₂(Ω)}.

Sucheu∈V(g) := {u∈V :u=g auf Γ_D}, sodass a(u, v) :=

Z

Ω

∇u· ∇vdx

= Z

Ω

f vdx+ Z

ΓN

∇un·vda (VP)

:=hl, vi ∀v ∈V(0).

Galerkin-Verfahren

Wähle endlich dimensionalen Ansatzraum bzw. TestraumV_h ⊂ V mi Knotenbasis{λ_z : z ∈ V}, V ={z₁, ..., z_n}. Dann ist

v_h =X

z∈V

v_zλ_z,

a(uh, vh) = hl, vhi ∀vh ∈Vh(0) mit u_h =g auf Γ_D und u_z =g(z)∀z∈Γ_D∩ V.

Das heißt u^>_NA_N =l^> mit A_N = (a(λ_z, λ_z))_z,z und l = (hl, λ_zi)_z.

Au=

A_N 0 0 I

u_N u_D

=









 l g(z₁^D)

... g(z_n^D)











Ritz-Galerkin-Verfahren

Genau dann istuh ∈Vh Lösung von (VP), wenn sie das Energiefunktional E(u_h) = 1

2a(u_h, u_h)− hl, u_hi

= 1

2u^>Au−l^>u mitu_h =g aufΓ_D minimiert, das heißt suche Nullstelle von

∇E(u_h) =A^>u−l = 0^!

∇²E(u_h) =A^>=A symmetrisch positiv definit.

Dann linearisieren:

∇E(u^k)≈ ∇E(u^k−1) +∇²E(u^k−1)∆w^k = 0^! Im Newton-Schritt:

∆w^k=−(∇²E(u^k−1))⁻¹∇E(u^k−1) Setzeu^k =u^k−1 −∆w^k

STOP:k∇E(u^k)k< ε.

(3.4) Satz

SeienA, B_p symmetrisch positiv definit und seiB ∈R^N×N mitu^k+1 =u^k+θBr^k, d.h.

additiv/parallel:B =P

R^>_pB_pR_p multiplikativ/successive:I_N −BA=

P

Q

p=1

(I_N −R^>_pBR₀a).

Dann existiert

Θ₀ >0 mit %=||I_N −ΘBA||_A<1 für Θ∈(0,Θ₀) und es gilt

|u−u^k|_A≤%^k|u−u⁰|_A mit |y|A = p

y^>Ay, ||C||A = sup

y6=0

|C_y|_A

|y|_A = ||L^>CL^−T||2 und A = LL^> Cholesky- Zerlegung.

Beweis.

|u−u^k|_A=|(I_N −ΘBA)^k(u−u⁰)|_A≤ ||I_N −ΘBA||^k_A|u−u⁰|

||I_N −ΘBA||_A=||I_N −ΘL^>BL||₂ = sup

λ∈σ(L^>BL)

|1−Θλ|

wobei gilt σ(L^>BL) = σ(BA)∈[λ_min, λ_max]⊂[0,∞]

⇒ Θ_opt = 2 λ_min+λ_max

⇒ ||I_N −ΘBA||_A<1 für Θ< 1 λ_max. Krylovraum:

(BA)u=Bb soll berechnet werden

⇒ berechneu^k∈u⁰ und span{P(BA)Br⁰ : P ∈P^k−1}

⇒ berechne optimalesP (3.5) Satz

a) SeienA, B symmetrisch positiv definit undαy^>Ay ≤ y^>ABAy ≤ Cy^>Ay. Dann gilt für das cg-Verfahren:

|u^k−u|_A≤2 √

κ−1

√κ+ 1 k

|u−u⁰|_A

mitκ= C α

b) SeienA, B regulär undαy^>y ≤y^>BAy,||BA||₂ ≤C. Dann gilt für GMRES

|u^k−u|₂ ≤κ₂(BA)(1−κ²)^k2|u−u⁰|₂. Multilevel-Vorkonditionierer

SeiV_H ⊂ V_h ein Finite-Elemente-Raum in einer gröberen Triangulierung. Wir elimi- nieren alle Dirichlet-Freiheitsgrade

{λ₁, ..., λ_N} KnotenbasisV_h(0) N = dimV_h(0) {λ^H₁ , ..., λ^H_N

0} KnotenbasisV_H(0) N₀ = dimV_H(0) < N.

Es existiert einR₀ ∈R^N0×N mitλ^H_k =

H

P

n=1

R₀[k, n]λ_n. Wir habenA∈R^N×N undA₀ ∈R^N0×N mit Einträgen

A[k, n] = Z

κ· ∇λ_n· ∇λ_kdx=:a(λ_n, λ_k) A₀[k, n] =a(λ^H_n, λ^H_k)

(3.6) Lemma

a) Es gilt für die so definierten MatrizenA₀ =R₀AR^>₀ (Galerkin-Produkt) b) P₀ =R^>₀A⁻¹₀ R₀A (Galerkin-Projektion) Orthogonal-Projektion bzgl.

hy, ziA=y^>Az.

Beweis.zu a)

A₀[k, n] =a(λ^H_n, λ^H_k) = a X

n⁰

R₀[n, n⁰]λ_n⁰,X

k⁰

R₀[k, k⁰]λ_k⁰

!

=X

n⁰,k⁰

R₀[n, n⁰]R₀[k, k⁰]a(λ_n⁰, λ_k⁰)

= (R₀AR^>₀)[k, n]

zu b)

P₀² = (R₀^>A⁻¹₀ R₀A)² =R^>₀A⁻¹₀ R₀AR^>₀A⁻¹₀ R₀A

=R^>₀A⁻¹₀ R₀A=P₀ DaP₀² =P: Projektion aufP₀(R^N) =R^>₀(R^N0)

(I_N −P₀)² =I_N −2P₀−P₀² =I_N −P₀

Ist es auch Orthogonalprojektion bzgl- Skalarprodukt?

hP₀v,(I_N −P₀)wi_A=v^>P₀^>A(I_N −P₀)w

=v^>(P₀^>A−P₀^>AP₀)w Nebenrechnung:

P₀^>AP₀ =AR^>₀A⁻¹₀ R₀AR^>₀A⁻¹₀ RA=P₀^>A

Zweilevel-Verfahren

Z0) Wähleu⁰ ∈R^N, berechner⁰ =b−Au⁰, setzek = 0, wähleε >0 undΘ>0 Z1) Falls|r^k|< ε STOP

Z2) c^k₀ =A⁻¹₀ R₀r^k (Grobgitterkorrektur) Z3) c^k =R^>₀c^k₀ +B(r^k−AR₀^>c^k₀) (Glätten)

Z4) u^k+1 =u^k+c^k, r^k+1 =r^k−Ac^k und setzek =k+ 1, gehe zu Z1) (3.7) Satz

a) Es gilt

u−u^k+1 = (I_N−BA)(I_N−P₀)(u−u^k)

b) SeiA, B symmetrisch positiv definit, ||BA||_A ≤ 1 undv^>Av ≤ Cv^>ABAv für P₀v = 0. Dann gilt

||(I_N−BA)(I_N−P₀)||_A ≤ r

1− 1 c. Bemerkung

C >0 ist unabhängig vonN, wennB ≈A⁻¹ aufN(P₀) ={v ∈R^N :P₀v = 0}.

Beweis. A=LL^> Cholesky-Zerlegung, w₁, ..., w_NOrthonormalbasis zu Eigenwertµ_nvon(L^>BL)w_n = µ_nw_n>0.

1≥ ||BA||²_A= sup

v6=0

v^>ABABAv

v^>Av = sup

w6=0

w^>L^∗BLw w^>w

= maxµ_n ∈(0,1]

⇒ |(I_N−BA)v|²_A =v^>(A−2ABA+ABABA)v

=w^>(I_N−2L^>BL+ (L^>BL)²)w

=X

(w_n^>w)²(1−2µn+µ²_n)

≤X

w^>_nw(1−µ_n)

=w^>(I_N−L^>BL)w

=v^>(A−ABA)v Betrachte weiter

||(IN−BA)(IN−P0)||A = sup

w6=0

|(I_N−BA)(I_N−P₀)w|²

|w|²_A

= sup

v=(1−P0)w

|(IN−BA)v|²_A

|v|²_A+|v −w|²_A

≤ sup

P0v=0, v6=0

|v|²_A−v^>ABAv

|v|²_A

≤1− 1 c

Anwendung:

B = Θ diag(A⁻¹)mitΘ ∈ (0,1), M ∈ R^N×N, M[n, k] = R

Ωλ_kλ_ndx. V_h ist ein Finite-Elemente-Raum zu regelmäßigen Gitternh= diamK, H = 2h. Dann gilt:

a) inf

φH∈V_H

||φ_h−φ_H||₀ ≤ch||∇φ_h||₀ ∀φ_h ∈V_h, ||φ₀||= qR

Ωφ²dx b) ||λ_h||₀ ≤c₁h⁻¹||∇λ_h||₀

||h⁻¹||₀ ≤c₂||λ_h||₀

⇒v^>M v ≤c₁h⁻²v^>B⁻¹v und h²v^>B⁻¹v ≤c₂v^>M v mit B⁻¹[n, n] =A[n, n] =

Z

κ∇λ_n· ∇λ_ndx c) φ_h =P

v[n]λ_n ⇒cv^>Av≤ ||∇φ_h||²₀ ≤Cv^>Av, cv^>M v ≤ ||φ_h||² ≤Cv^>M v

d) φ_h =P

v[n]λ_n, P₀v = 0

⇒a(φh, φh) =v^>Av= inf

w0∈R^N0

v^>A(v−R^>₀w0)

= inf

w0

v^>AB¹²B⁻¹²(v−R^>₀w0)

≤√

v^>ABAv r

inf

w0

(v −R^>₀w0)B⁻¹(v−R^>₀w0)

| {z }

≤Ch⁻¹ inf ||φ_h−φ_H||₀

≤C||∇φ_H||₀

≤Cv^>Av Multilevel-Verfahren

V₀(0) ⊂V₁(0)⊂ · · · ⊂V_L(0), h_l=h₀2^−l, l= 0, ..., L

N_l = dimV_l(0), A_l ∈R^Nl×Nl, R_l ∈R^Nl×Nl+1, A_l=R_lA_l+1R_l^>

M0) Wähleu⁰_L, r⁰_L=b_L−A_Lu_L M1) Falls|r^k_L|< ε STOP

M2) Setzel =L, c_L = 0, r_L =r_L^k

M3) Fürν = 1, ..., ν_l: w_l =B_lr_l, c_l :=c_l+w_l, r_l :=r_l−A_lw_l M4) rl−1 =Rl−1rl, setze l:=l−1

M5) Fallsl >0 gehe zu M3) M6) c₀ =A⁻¹₀ r₀

M7) Setzel :=l+ 1, c_l =R^>_l−1c_l−1 M8) Glätten (wie (M3))

M9) Fallsl < L gehe zu (M7)

M10) u^k+1_L =u^k_L+c_L, r^k+1_L =r^k_L−A_Lc_L, k :=k+ 1

→ gehe zu (M1)

4 Gemischte und hybride Finite Elemente

Sei ein Differentialgleichungssystems gegeben, das von mehreren Variablen abhängt. Hier kann es sinnvoll sein, die einzelnen Unbekannten in verschiedenen Räumen zu diskretisie- ren. In diesem Fall spricht man von gemischten Finiten Elementen. Betrachte hierzu das Problem

divq = 0 inΩ q = −κ∇uinΩ

−q·n = g_NaufΓ_N u = u_DaufΓ_D wobeiκsymmetrisch positiv definit ist. Es gilt

Z

Ω

divqφdx= 0 für alle Testfunktionenφ :Ω→R. und

Z

Ω

κ⁻¹(q+κ∇u)·ψdx= 0 für alle Testfunktionenψ :Ω→R² Mit dem Satz von Gauß folgt zudem füruψ:

Z

∂Ω

(uψ)·nda= Z

Ω

div(uψ) dx= Z

Ω

∇u·ψdx+ Z

Ω

u·divψdx.

Dies führt auf das Sattelpunktproblem: Bestimme(q, u)mitq·n=−g_NaufΓ_Nund Z

Ω

κ⁻¹q·ψdx− Z

Ω

udivψdx = − Z

ΓD

u_Dψ·nda Z

Ω

divqφdx = 0

für alle(ψ, φ)mitψ·n= 0aufΓ_N. Diskretisierung

SeiKeine Triangulierung undF die Menge aller Seiten∂K ∩∂K⁰ oder∂K ∩∂Ω.

Wähle FE-RäumeW_h, Q_hund W_h(g) ={ψ_h ∈W_h :

Z

F

ψ_h·nda= Z

F

gda ∀F ⊂Γ_N}.

Bestimme(qh, uh)∈Wh(−g_N)×Qh mit Z

Ω

κ⁻¹q_h·ψ_hdx− Z

Ω

u_hdivψ_hdx=− Z

ΓD

u_Dψ_h ·nda

− Z

Ω

divq_hφ_hdx= 0 ∀(ψ_h, φ_h)∈W_h(0)×Q_h.

ZuKˆ = conv{ ⁰₀ , ¹₀

, ⁰₁

}definiereWˆ = span{ψˆ₀,ψˆ₁,ψˆ₂}und Fˆ₀ = convn

0 0

,

1 0

o

Fˆ₁ = convn 1 0

,

0 1

o

Fˆ₂ = convn 0 1

,

0 0

o

mit

ψˆ₀(ξ) =

ξ₁ ξ₂−1

ψˆ₁(ξ) = ξ₁

ξ₂

ψˆ₂(ξ) =

ξ₁−1 ξ₂

.

ZuKˆ = conv{ ⁰₀ , ¹₀

, ¹₁ , ⁰₁

}definiereWˆ = span{ψˆ0,ψˆ1,ψˆ2,ψˆ3}und Fˆ₀ = convn

0 0

,

1 0

o

Fˆ₁ = convn 1 0

,

1 1

o

Fˆ₂ = convn 1 1

,

0 1

o

Fˆ₃ = convn 0 1

,

0 0

o

mit

ψˆ₀(ξ) =

0 ξ₂−1

ψˆ₁(ξ) = ξ₁

0

ψˆ₂(ξ) = 0

ξ₂

ψˆ₃(ξ) =

ξ₁−1 0

. Es gilt

Z

Fˆk

ψˆ·nˆdˆa=

(1 j =k 0 sonst.

(4.1) Lemma

Seinˆ Normale an∂K. Seiˆ ϕ_K : ˆK →K¯ linear affin mit

ϕ_K(ξ) = z_0,K +B_Kξ, J_K = detB_K >0.