Blatt 02, 23.10.2019¶

Aktivieren Sie wie gewohnt ihre Arbeitsumgebung und starten Sie den Jupyter Notebook server, siehe zB Blatt 1, Aufgabe 0).
Laden Sie dieses Blatt von der Homepage herunter und speichern Sie es unter ~/NPDGL1/notebooks/blatt_02.ipynb.
Importieren Sie numpy und pymor.basic und machen Sie matplotlib für das Notebook nutzbar.

%matplotlib notebook
import numpy as np
from pymor.basic import *

äquidistante Gitter in 1d¶

Als Partitionierung eines zusammenhängenden beschränkten Gebietes $\Omega \subset \mathbb{R}$ (also eines Intervals $(a, b)$ für $a < b$) betrachten wir ein äquidistantes Gitter $\mathcal{T}_h := \{ T_0, \dots, T_{I - 1}\}$

für $I \in \mathbb{N}$, mit
Gitterweite $h := \frac{|b - a|}{I}$,
$I + 1$ Knotenpunkten $x_i := a + ih$ für $0 \leq i \leq I$, sodass $x_0 = a$ und $x_I = b$ und
$I$ Elementen $T_i := (x_i, x_{i + 1}) \in \mathcal{T}_h$ für $0 \leq i \leq I - 1$.

Das Gitter ist also

eine Überdeckung von $\Omega$, d.h. $\bigcup_{i = 0}^{I - 1}\overline{T_i} = \overline{\Omega}$, und zwar

in disjunkte Elemente, d.h. $\overline{T_i} \cap \overline{T_j} = \begin{cases}\text{Knotenpunkt }x_i,&\text{falls } j = i - 1\\\overline{T_i},&\text{falls } j = i\\\text{Knotenpunkt } x_{i + i},&\text{falls } j = i + 1\\\emptyset,&\text{sonst}\end{cases}.$

stückweise lineare Funktionen in 1d¶

Für ein zusammenhängendes beschränktes Gebiet $T \subset \mathbb{R}$ sei $\mathbb{P}_k(T)$ der Raum der Polynome über $T$ vom Grad höchstens $k \in \mathbb{N}_0$.

Zu einem gegebenen Gitter $\mathcal{T}_h$ betrachten wir den Raum der stetigen stückweise linearen Funktionen bezüglich des Gitters:

$$S_h^1 := \big\{ u \in C^0(\Omega) \;\big|\; u|_T \in \mathbb{P}_1(T) \big\}$$

Jede Funktion $u_h \in S_h^1$ ist dabei eindeutig definiert durch ihre Werte in den Knotenpunkten des Gitters. Wir können also der Funktion $u_h$ ihren Vektor an Freiheitsgeraden (DoF Vektor, von degree of freedom) $\underline{u_h} \in \mathbb{R}^{I + 1}$ zuordnen, sodass $(\underline{u_h})_i = u_h(x_i)$ für $0 \leq i \leq I$. Es besteht also die folgende Isomorphie:

$$\begin{align*}S_h^1 &\cong \mathbb{R}^{I + 1}\\ u_h&\leftrightarrow \underline{u_h}\end{align*}$$

Aufgabe 1: Interpolation¶

Da wir beliebige Funktionen $f \in C^0(\Omega)$ nicht fehlerfrei handhaben können (zum Beispiel Integrieren oder visualisieren) betrachten wir geeignete Approximationen $f_h \in S_h^1$. Dazu definieren wir den zu $S_h^1$ passenden Interpolationsoperator

$$\begin{align*}\Pi_{S_h^1}; C^0(\Omega) &\to S_h^1,\\f&\mapsto \Pi_{S_h^1}[f] =: f_h,\end{align*}$$

wobei $f_h$ eindeutig bestimmt ist durch die Auswertung von $f$ an den Knotenpunkten des Gitters:

$$\begin{align*}f_h(x_i) = f(x_i)&&\text{für } 0 \leq i \leq I\end{align*}$$

Implementieren Sie den Interpolationsoperator $\Pi_{S_h^1}$ als Python Funktion und testen Sie Ihr Program für $\Omega = (0, 2\pi)$, $I = 4$ und $f(x) = sin(x)$!

Hinweis: Sie können der Einfachheit halber for-Schleifen verwenden!

Schreiben Sie eine Python Funktion f_h = interpolate(grid, f), die für ein gegebenes pyMOR Gitter in einer Raumdimension, grid$ = \mathcal{T}_h$, und eine gegebene pyMOR Funktion, f$ = f$, den DoF-Vektor der Interpolierten, f_h$ = \underline{f_h}$, als Numpy-Array entsprechender länge zurück gibt.
- Stellen Sie dazu in der Funktion sicher, dass das gegebene Gitter die Partition eines Gebietes in einer Raumdimension ist.
- Speichern Sie in der Funktion alle Knotenpunkte des Gitters in vertices ab.
- Werten Sie die Funktion f and allen Knotenpunkten aus (siehe Dokumentation von FunctionInterface) und speichern Sie die Werte in values.
- Stellen Sie sicher, dass Sie genauso viele Funktionswerte wie Knotenpunkte berechnet haben.
- Geben Sie values zurück.

Ein äquidistantes Gitter in einer Raumdimension wird in pyMOR durch die Klasse OnedGrid modelliert.
- Legen Sie ein entsprechendes Objekt mit Namen grid an,
- geben Sie die Dimension des Gitters,
- die Anzahl der Knotenpunkte sowie
- die Position jedes Knotenpunktes aus (Wiederholung von Blatt 00).

Analytische Funktionen werden in pyMOR durch ExpressionFunction modelliert. Legen Sie ein entsprechendes Objekt mit Namen f an.

Berechnen Sie die Interpolation von $f$ auf dem Gitter $\mathcal{T}_h$ mittels Ihrer interpolate Funktion.
- Rufen Sie Ihre Funktion auf und speichern Sie das Resultat in f_h ab.
- Geben Sie den Python-Typ von f_h aus.
- Geben Sie die Länge von f_h aus.
- Geben Sie mit Hilfe einer for-Schleife die Koordinaten aller Knotenpunkte des Gitters zusammen mit den Werten von $f$ aus (Hinweis: Nutzen Sie zip).

Aufgabe 2: Visualisierung¶

Schreiben Sie eine Python Funktion visualize(grid, f_h, name=''), welche die Interpolation einer Funktion über einem Gitter mit Hilfe von matplotlib visualisiert und testen Sie Ihr Program für $\Omega = (0, 2\pi)$, $f(x) = sin(x)$ und Gitter mit $I = 2^2, 2^3, 2^4, 2^5$ Elementen!

Machen Sie das benötigte Modul mit dem Befehl from matplotlib import pyplot as plt unter dem Namen plt verfügbar und machen Sie sich mit plt.plot und plt.legend vertraut.

Schreiben Sie die Funktion visualize, welche eine Funktion visualisiert.
- Erstellen Sie in der Funktion mit plt.figure() eine neue (leere) Grafik.
- Stellen Sie in der Funktion sicher, dass das gegebene Gitter die Partition eines Gebietes in einer Raumdimension ist.
- Stellen Sie in der Funktion sicher, dass die Länge von f_h mit der Anzahl der Knotenpunkte des Gitters überein stimmt.
- Nutzen Sie die plot Funktion aus dem plt Modul, um die Werte von f_h an den Knotenpunktes des Gitters zu visualisieren und linear zu interpolieren, und den Namen der Funktion zu übergeben.
- Nutzen Sie die legend Funktion aus dem plt Modul, um die Legende anzeigen zu lassen.

Testen Sie Ihre Funktion für das Gitter und die Interpolation aus Aufgabe 1.

Erweitern Sie Ihre Funktion, sodass Sie mehrere Interpolierte gleichzeitig visualisieren können.
- Erstellen Sie vor Ihrer Funktion folgendermaßen eine Farbskala:
```
colors = list(plt.get_cmap('tab20c').colors)
colors.reverse()
```
- Erstellen Sie in Ihrer Funktion mit plt.figure() einmalig eine neue (leere) Grafik.
- Testen Sie in Ihrer Funktion jedes der Argumente grid, f_h und name, ob es sich um eine Python Liste handelt (Hinweis: nutzen Sie dafür isinstance), und wandeln Sie es gegebenenfalls in eine Liste der Länge der 1 um.
- Stellen Sie in Ihrer Funktion sicher, dass alle drei Listen dieselbe Länge haben.
- Iterieren Sie in Ihrer Funktion über diese Listen (Hinweis: nutzen Sie zip) und gehen Sie wie oben vor, um jede Interpolierte bezüglich ihres Gitters zu visualisieren. Geben Sie in Ihrer Funktion zusätzlich beim Aufruf von plot mit Hilfe von color=colors.pop() jeder visualisierung eine andere Farbe.
- Rufen Sie danchach in Ihrer Funktion einmal legend() auf, um die Legend anzeigen zu lassen.

Testen Sie Ihre Funktion für die oben angegeben Gitter. Erstellen Sie dazu drei leere Listen und füllen Sie diese mit Hilfe einer for-Schleife mit den entsprechenden Gittern, Interpolierten und Namen (welche die Anzahl der Gitterelemente enthalten).

Aufgabe 3: Finite Differenzen Verfahren in einer Raumdimension¶

Wir betrachten das vereinfachte Diffusionsproblems aus Blatt 00 in einer Raumdimension mit $A = 1$ (Poisson Problem), gesucht ist also $u \in C^2(\Omega)$, sodass

$$\begin{align*}\partial_{xx} u(x) &= f(x) &&\text{für } x \in \Omega\\ u(x) &= g_\text{D}(x) &&\text{für } x \in \partial\Omega.\end{align*}$$

Zur Diskretisierung mit dem Verfahren der Finiten Differenzen für ein gegebenes äquidistantes Gitter $\mathcal{T}_h$ gehen wir folgendermaßen vor:

Wir betrachten den Differentialoperator $\partial_{xx}:C^2(\Omega) \to C^0(\Omega)$, $u \mapsto \partial_{xx}u$. Gesucht ist eine Approximation, welche für Funktionen in $S_h^1$ wohldefiniert ist: $\partial_{h,h}: S_h^1 \to S_h^1$.
- Betrachten Sie dazu für eine Funktion $u \in C^\infty(\Omega)$ die Taylor-Entwicklung von $u$ um einen Gitterpunkt $x_i$, ausgewertet an einem beliebigen Punkt $y \in \Omega$:
  
  $$u_(y) = u(x_i) + u'(x_i) \cdot (y - x_i) + \tfrac{1}{2} u''(x_i)\cdot(y - x_i)^2 + \mathcal{O}\big((y - x_i)^3\big)$$
- Nutzen Sie diese Entwicklung, ausgewertet an den benachbarten Gitterpunkten $x_{i - 1}$ und $x_{i + 1}$, um einen Differenzenquotionten als Approximation von $u''(x_i)$ der Güte $\mathcal{O}(h)$ herzuleiten (Hinweis: nuzten Sie die Gitterweite $h$, um die Darstellung zu vereinfachen).
- Definieren Sie $\partial_{h, h}$ mit Hilfe dieses Differenzenquotient.

Leiten Sie eine Approximationen obiger Differentialgleichung für die approximative Lösung $u_h \in S_h^1$ her.
- Formulieren Sie dazu die PDGL um, indem Sie $\partial_{h, h}$ verwenden, und Gleichheit nur noch an allen Knotenpunkten des Gitters fordern.
- Leiten Sie daraus $I - 1$ Gleichungen für den Zusammenhang zwischen der gesuchten Funktion $u_h$ und der gegebenen Funktion $f$ an allen inneren Gitterpunkten $x_i$, für $1 \leq i \leq I$, her.
- Leiten Sie mit Hilfe der Randbedingung zwei Bedingungen für $u_h$ an den Randpunkten des Gitters, $x_0$ und $x_{I + 1}$, her.

Leiten Sie mit Hilfe dieser $I + 1$ Gleichungen ein Gleichungssystem der Form

$$\underline{b_h} \cdot \underline{u_h} = \underline{f_h}$$

für den DoF-Vektor der approximativen Lösung $u_h \in S_h^1$ her, indem Sie
- die Einträge der Matrix $\underline{b_h} \in \mathbb{R}^{I + 1 \times I + 1}$ und
- die Einträge des Vektors $\underline{f_h} \in \mathbb{R}^{I + 1}$
bestimmen.

Assemblieren Sie die Matrix $\underline{b_h}$ und den Vektor $\underline{f_h}$ in Numpy Arrays entsprechender Größe mit Hilfe von for-Schleifen und/oder Ihrer interpolate Funktion.

Lösen Sie das resultierende Gleichungssystem mit Hilfe von numpy.linalg.solve (Hinweis: numpy ist als np verfügbar).

Fassen Sie Ihr Programm in einer Funktion u_h = discretize_elliptic_1d_fd(grid, f, g_D) zusammen und testen Sie Ihr Programm für $\Omega = (0, 1)$, $f = 1$ und $g_\text{D} = 0$ und für Gitter mit $2^2, 2^3, 2^4, 2^5$ Elementen (dabei soll g_D ein Python Tuple mit zwei Werten sein, also in diesem Beispiel g_D = (0, 0)). Visualisieren Sie die Lösungen mit Hilfe der visualize Funktion aus Aufgabe 2.

Testen Sie Ihr Programm für $\Omega = (0, 1)$ und ein Gitter mit $2^7$ Elementen und die folgenden Fälle, visualisieren Sie die Lösungen.
- $f = 1$, $g_D(x) = \begin{cases}0,&x = 0\\1, &x = 1\end{cases}$
- $g_D = 0$, $f(x) = \begin{cases}1, &x \in [0.2, 0.4]\\0, &\text{sonst}\end{cases}$