20. Finite Elemente in der Grundvorlesung

20.1. Eine kleine Einführung in Partielle Differentialgleichungen und die Finite Elemente Methode

Wir wollen die Poisson-Gleichung lösen:

Gegeben ist eine offene Menge \(D \subset {\mathbb R}^2\) und eine Funktion \(f : D \rightarrow {\mathbb R}\).

Gesucht ist die Funktion \(u : D \rightarrow {\mathbb R}\) sodass

\[ \Delta u (x,y) = f(x,y) \quad \forall \, (x,y) \in D, \]

wobei der Laplace-Operator definiert ist als

\[ \Delta u (x,y) = \frac{\partial^2 u}{\partial x^2}(x,y) + \frac{\partial^2 u}{\partial y^2}(x,y) \]

Die Lösung ist eindeutig wenn noch Randbedingungen gefordert werden:

\[ u(x,y) = 0 \quad \forall \, (x,y) \in \partial D \]

Mit solchen Partiellen Differentialgleichungen können viele Aufgabenstellungen in den Naturwissenschaften, in technischen Anwendungen, als auch in der Finanzmathematik modelliert werden. Eine exakte Lösung zu finden ist oft nicht möglich, und man berechnet mit Computereinsatz näherungsweise Lösungen (Numerik).

20.2. Approximation von Funktionen auf Triangulierungen

Wir verwenden dazu das Python - Paket NGSolve (www.ngsolve.org).

Wir wollen zuerst ein paar einfache numerische Operationen (differenzieren, integrieren) mit Funktionen durchführen:

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))

Wir verwenden die Funktion (den Gauß-peak):

\[ f(x,y) = \exp (-100 \, ( (x-\tfrac{1}{2})^2 + (y-\tfrac{1}{2})^2 ) \]

func = exp(-100* ((x-0.5)**2+(y-0.5)**2))
Draw (func, mesh, deformation=True);

Diese Funktion wird jetzt auf dem Finite-Elemente Netz interpoliert: Wir berechnen dazu die Funktionswerte in den Knoten, und ersetzen auf jedem einzelnen Dreieck die komplizierte Funktion durch eine affin-lineare Funktion der Bauart \(a + b x + cy\), mit \(a,b,c \in {\mathbb R}\), und mit den gleichen Funktionswerten in den Knoten.

fes = H1(mesh, order=1)
gfu = GridFunction(fes)
gfu.Interpolate (func)
Draw (gfu, deformation=True);

20.3. Berechnung von Integralen:

Wir wollen nun das Integral

\[ \int_D f(x,y) \, d(x,y) \]

berechnen. Das ist das Volumen unter dem Graphen der Funktion.

Wir nutzen die interpolierte Funktion um das Integral näherungsweise (numerisch) zu berechnen: Wir ersetzen die Funktion durch das arithmetische Mittel der Werte in den Eckpunkten, berechnen das Volumen dieses Prismas nach Grundfläche \(\times\) Höhe, und summieren über alle Dreiecke:

\[ \sum_{\text{Dreiecke}\, T} \text{Fläche(T)} \sum_{\text{Eckpunkte}\, x_i} \frac{f(x_i)}{3} \]

Integrate (gfu, mesh)

0.031263650037777674

Um die Genauigkeit zu erhöhen können wir ein feineres Netz verwenden. Eine effizientere Alternative ist es nicht mit affin-linearen, sondern mit Polynomen höherer Ordnung zu interpolieren:

fes = H1(mesh, order=3, dirichlet=".*")
gfu = GridFunction(fes)
gfu.Interpolate (func)
Draw (gfu, deformation=True);

Integrate (gfu, mesh)

0.03141592651374485

Der Fehler zum exakten Wert \(\pi/100\) beträgt:

Integrate (gfu, mesh) - pi/100

-2.2153084922038602e-11

20.4. Berechnung von partiellen Ableitungen:

Wir wollen die partiellen Ableitungen

\[ \frac{\partial f}{\partial x} (x,y) = \lim_{\varepsilon \rightarrow 0} \frac{f(x+\varepsilon,y) - f(x,y) }{\varepsilon} \]

und

\[ \frac{\partial f}{\partial y} (x,y) = \lim_{\varepsilon \rightarrow 0} \frac{f(x, y+\varepsilon) - f(x,y) }{\varepsilon} \]

berechnen. Das sind die Steigungen der Funktion in \(x\)- bzw. in \(y\)-Richtung. Beide partiellen Ableitungen als Vektor geschrieben wird als Gradient bezeichnet:

\[ \nabla f := \operatorname{grad} \, f := \left( \frac{ \partial f}{\partial x}, \frac{\partial f}{\partial y} \right) \]

print ("df/dx:")
Draw (grad(gfu)[0], mesh)
print ("df/dy:")
Draw (grad(gfu)[1], mesh);

df/dx:

df/dy:

Und gleich der Gradient als Vektor gezeichnet, dargestellt durch Kegel: Die Spitze zeigt in in Richtung des steilsten Anstieges:

Draw (grad(gfu), mesh, vectors={"grid_size":50});

20.5. Projektion auf Teilräume

Wir wollen jetzt die gegebene Funktion nicht durch Interpolation annähern, sondern durch orthogonale Projektion bzgl. eines Skalarproduktes.

Sei \((V, \left<.,.\right>)\) ein Euklidscher Vektorraum, und \(V_h \subset V\) ein linearer Teilraum. Für \(u \in V\) ist die orthogonale Projektion \(Pu \in V_h\) auf den Teilraum definiert sodass:

\[ \left< u - Pu, v \right> \quad \forall \, v_h \in V_h \]

Der Projektionsfehler \(u - Pu\) steht orthogonal auf den Teilraum. Diese Projektionsaufgabe können wir auch schreiben als:

\[ \text{Ges.:} \; u_h \in V_h : \left< u_h, v_h \right> = \left< u, v_h \right> \quad \forall \, v_h \in V_h \]

Das Skalarprodukt auf der linken Seite ist bi-linear, d.h. linear in beiden Argumenten. Darum wird es auch Bilinearform genannt:

\[ A : V \times V \rightarrow {\mathbb R} : \quad A(u,v) = \left< u, v \right> \]

Für ein gegebenes, festes \(u\) ist der Term auf der rechten Seite linear in \(v\), also eine Linearform:

\[ f : V \rightarrow {\mathbb R} : \quad f(v) = \left< u, v \right> \]

Das Problem lautet nun:

\[ \text{Ges.:} \; u_h \in V_h : A(u_h, v_h) = f(v_h) \quad \forall \, v_h \in V_h \]

Wir verwenden für die Projektion das \(L_2\) Skalarprodukt:

\[ \left< u, v \right>_{L_2(D)} := \int_D u(x,y) \, v(x,y) \, d(x,y) \]

Dieses Skalarprodukt kann in NGSolve als u*v*dx angegeben werden. Die Projektion wird durch Aufstellen (Assemblierung) und Lösen eines linearen Gleichungssystems berechnet:

u,v = fes.TnT()

A = BilinearForm(u*v*dx).Assemble().mat
f = LinearForm(func*v*dx).Assemble().vec

gfu.vec.data = A.Inverse() * f

Draw (gfu);

Die berechnete Funktion minimiert den Fehler in der \(L_2\)-Norm:

\[ \| u - u_h \|_{L_2} = \min_{v_h \in V_h} \| u - v_h \|_{L_2} \]

sqrt(Integrate ( (gfu-func)**2, mesh))

2.4802985545274948e-05

Die Berechnung der Projektion funktioniert so:

Wir wählen eine Basis \(p_i\) für den Vektorraum \(V_h\). Da \(V_h\) ein Vektorraum von Funktionen ist, sind diese Basiselemente ebenfalls Funktionen \(p_i(x,y)\). Im einfachsten Fall wählt man Funktionen die in einem Knoten 1 sind, in allen anderen 0, und dazwischen linear interpoliert:

scene = Draw (gfu)
from time import sleep
sleep(5)

for i in range(200,220):
    gfu.vec[:] = 0
    gfu.vec[i] = 1
    scene.Redraw()
    sleep(1)

Wir stellen die gesucht Lösung in dieser Basis dar (mit \(N = \dim V_h\)):

\[ u_h(x,y) = \sum_{i=1}^{N} u_i {\bf p}_i(x,y) \]

Der Koeffizientenvektor \(u = (u_1, \ldots u_n) \in {{\mathbb R}^N}\) wird durch die Projektionseigenschaft bestimmt:

\[ A( \sum u_i {\bf p}_i, {\bf p}_j ) = f( {\bf p}_j ) \quad \forall \, j = 1, \ldots N \]

Wir nutzen die Linearität von \(A(.,.)\)

\[ \sum_{i=1}^N A({\bf p}_i, {\bf p}_j ) \; u_i = f( {\bf p}_j ) \quad \forall \, j = 1, \ldots N, \]

und erhalten das lineare Gleichungssystem für den Koeffizientenvektor \(u\)

\[ A u = f \]

mit den Einträgen der Matrix A und des rechte-Seite Vektors f:

\begin{eqnarray*} A_{j,i} & = & A(p_i, p_j) = \int p_i p_j , d(x,y) \ f_j & = & f(p_j) = \int u , p_j , d(x,y) \end{eqnarray*}

20.6. Projektion mit dem \(H^1\)-Skalarprodukt

Wir kehren zurück zur Poisson-Gleichung:

Gesucht ist die Funktion \(u(x,y)\) sodass:

\[ \Delta u(x,y) = f(x,y) \quad \forall \, (x,y) \in D \]

und

\[ u(x,y) = 0 \qquad \forall \, (x,y) \in \partial D \]

Wir wollen die Projektion \(u_h = P u\) der unbekannten Funktion \(u\) auf den endlich-dimensionalen Raum \(V_h\) berechnen. Als Skalarprodukt wählen wir jetzt (der Raum \(H^1\) heißt Sobolev-Raum und ist ein Hilbert-Raum):

\[ \left< u, v \right>_{H^1} = \int_D \nabla u \cdot \nabla v \; d(x,y), \]

wobei die Notation \(\nabla u \cdot \nabla v = \left< \nabla u, \nabla v\right>_{{\mathbb R}^2}\) verwendet wird. Analog zu oben bestimmen wir die Matrix \(A\) mit Einträgen

\[ A_{ij} = \int_D \nabla p_i \cdot \nabla p_j \, d(x,y) \]

und den Vektor \(f\) mit Einträgen

\[ f_j = \int_D \nabla u \cdot \nabla p_j \, d(x,y) \]

Die Pointe ist dass dieser Vektor \(f\) OHNE die unbekannte Funktion berechnet werden kann. Dazu nutzen wir die Green’sche Formel, ähnlich zur partiellen Integration:

\[ \int_D \nabla u \cdot \nabla p_j \, d(x,y) = -\int_D \underbrace{\Delta u}_{f} \, p_j \, d(x,y) + \int_{\partial D} \nabla u \cdot \nu \, p_j \]

Hier steht nun die gegebene Funktion \(f = \Delta u\), die Terme am Rand \(\partial D\) verschwinden wegen der Randbedingung.

A = BilinearForm(grad(u)*grad(v)*dx).Assemble().mat
f = LinearForm(100*func*v*dx).Assemble().vec

gfu.vec.data = A.Inverse(fes.FreeDofs()) * f

Draw(gfu);

Zur Kontrolle berechnen wir \(-\Delta u_h\):

Draw (-Trace(gfu.Operator("hesse")), mesh);

20.7. Ein paar weitere Differentialgleichungen:

• Elastizität beschreibt die Deformation elastischer Obejekte

• Navier-Stokes Gleichung beschreiben die Strömung inkompressibler Flüssigkeiten

• Maxwell’schen Gleichungen beschreiben die Interaktion von elektrischen und magnetischen Feldern