Praktikum zur Vorlesung Numerik partieller Differentialgleichungen I
Mario Ohlberger, Felix Schindler
~/NPDGL1/notebooks/blatt_03.ipynb.numpy und pymor.basic und machen Sie matplotlib für das Notebook nutzbar.%matplotlib notebook
import numpy as np
from pymor.basic import *
from matplotlib import pyplot as plt
Nutzen Sie für diesen Zettel das folgende Gitter mit 4 Elementen!
omega = (0., 1.)
grid = OnedGrid(omega, num_intervals=4)
Ein Merkmal affiner Gitter $\mathcal{T}_h$ mit Referenzelement $\hat{T} \subset \mathbb{R}^d$ ist es, dass zu jedem Gitterelement $T \in \mathcal{T}_h$ eine affine Abbildung, die Referenzabbildung, $$\begin{align*} F_T: \hat{T}&\to T\\ \hat{x}&\mapsto x := F_T(\hat{x}) := A_T\,\hat{x} + b_T, \end{align*}$$ gibt, mit einer positiv definiten Matrix $A_T \in \mathbb{R}^{d \times d}$ und einem Vektor $b_T \in \mathbb{R}^d$, sodass $$ T = F_T(\hat{T}). $$
In einer Raumdimension gilt $\hat{T} := [\hat{a}_0, \hat{a}_1]$ mit $\hat{a}_0 = 0$ und $\hat{a}_1 = 1$.
Schreiben Sie eine Funktion F(grid, T, x_hat), welche die Auswertung von $F_T(\hat{x})$ für alle $\hat{x} \in $x_hat berechnet und testen Sie Ihre Funktion mit dem Code in der nächsten Code-Zelle im Notebook!
Hinweise:
Sie können davon ausgehen, dass Sie als erstes Argument grid ein affines Gitter im pyMOR Sinne übergeben wird, und als zweites Argument T$\in \mathbb{N}$ der Index eines Gitterelementes $T \in \mathcal{T}_h$.
Als drittes Argument soll es möglich sein
F(grid, T, 0.5), als auchF(grid, T, [0., 0.5, 1.]), als auchnumpy-array, also F(grid, T, np.linspace(0., 1., 3)),übergeben zu können. Sie können dies mit Hilfe von isinstance innerhalb der Funktion erreichen:
if not isinstance(x_hat, np.ndarray):
x_hat = np.array(x_hat)Nach diesem Aufruf können Sie davon ausgehen, dass x_hat immer ein antsprechendes numpy-array ist
Nutzen Sie die embeddings Funktion des affinen Gitters (siehe Dokumentation oben), um $A_T$ und $b_T$ zu erhalten.
Geben Sie ein numpy-array der Größe len(x_hat) zurück.
T_hat = grid.reference_element(0)
a_hat = T_hat.subentities(0, 1)
print('T_hat = {}'.format(a_hat))
for T in range(grid.size(0)):
print(' F_{}(T_hat) = {}'.format(T, F(grid, T, a_hat)))
Stellen Sie in Ihrer Funktion außerdem sicher, dass alle Punkte in x_hat im Referenzelement $\hat{T}$ liegen und testen Sie, dass der folgende Aufruf fehlschlägt! Geben Sie eine informative Fehlermeldung, falls die Punkte nicht im Referenzelement liegen.
Hinweis:
Nutzen Sie dazu assert mit zwei Argumenten: assert condition, 'error message if condition is False'
Kommentieren Sie nach erfolgreichem Test die nächsten Code-Zellen mit Hilfe von # aus.
F(grid, 0, -1)
F(grid, 0, 2)
Eine Aussage von Lemma 2.11 ist dass $F_T(\hat{T}) = T$, es gilt also insbesondere
$$ F_T(\hat{a}_i) = a_i $$für alle
für $1 \leq i \leq 2$.
Überprüfen Sie diese Aussage für alle Element des Gitters.
Hinweis:
Nutzen Sie die subentities Methode des affinen Gitters (siehe Dokumentation oben), um die globalen Indices der Knotenpunkte eines Elementes zu erhalten.
Nutzen Sie diese globalen Indices, um die Koordinaten der Knotenpunkte eines Elementes zu erhalten.
Vergleichen Sie die so erhaltenen Koordinaten der Knotenpunkte eines Elementes mit den Ergebnissen der Anwendung der Referenzabbildung auf a_hat.
Eine weitere Aussage von Lemma 2.11 und Übungszettel 4, Aufgabe 2 ist die folgende Abschätzung für den Gradienten der Referenzabbildung:
$$ \|\nabla F_T\| = \|A_T\| \leq \frac{h(T)}{\rho(\hat{T})}\quad,\forall T \in \mathcal{T}_h, $$wobei $h(T)$ der Durchmesser eines Gitterelementes und $\rho(\hat{T})$ der Durchmesser der Größten in $\hat{T}$ enthaltenen Kugel ist.
Überprüfen Sie diese Aussage für alle Elemente des Gitters.
Hinweis:
Nutzen Sie np.linalg.norm, um die Norm einer Matrix zu berechnen.
Den Durchmesser eines Elementes erhalten Sie vom Gitter (siehe Dokumentation oben), $\rho(\hat{T})$ können Sie analytisch herleiten.
Außerdem existiert zu jedem Gitterelement $T \in \mathcal{T}_h$ die inverse Referenzabbildung $$\begin{align*} F^{-1}_T: T &\to \hat{T}\\ x&\mapsto \hat{x} := F^{-1}_T(x). \end{align*}$$
Schreiben Sie eine Funktion F_inv(grid, T, x), welche die Auswertung von $F^{-1}_T(x)$ berechnet und testen Sie Ihre Funktion analog zum Test von der Funktion F oben!
Die Ausgabe des Testaufrufs soll folgendermaßen aussehen:
T_hat = [0 1]
F_inv_0([0.0, 0.25]) = [0. 1.]
F_inv_1([0.25, 0.5]) = [0. 1.]
F_inv_2([0.5, 0.75]) = [0. 1.]
F_inv_3([0.75, 1.0]) = [0. 1.]Hinweise:
Es gelten dieselbe Hinweise wie für Teilaufgabe 1.
Nutzen Sie np.linalg.solve.
Testen Sie Ihre Funktion, indem Sie für jedes Element die Knotenpunkte in a abspeichern und das Resultat von F_inv(grid, T, a) betrachten.
Stellen Sie in Ihrer Funktion außerdem sicher, dass alle Punkte in x im Element $T$ liegen und testen Sie, dass der folgende Aufruf fehlschlägt! Geben Sie eine informative Fehlermeldung, falls die Punkte nicht im Element liegen.
Hinweis:
F_inv(grid, 0, -1)
F_inv(grid, 0, 2)
Eine weitere Aussage von Lemma 2.11 und Übungszettel 4, Aufgabe 2 ist die folgende Abschätzung für den Gradienten der inversen Referenzabbildung:
$$ \|\nabla F_T^{-1}\| = \|A_T^{-1}\| \leq \frac{h(\hat{T})}{\rho(T)}\quad,\forall T \in \mathcal{T}_h, $$wobei $h(\hat{T})$ der Durchmesser des Referenzelementes und $\rho(T)$ der Durchmesser der Größten in $T$ enthaltenen Kugel ist.
Überprüfen Sie diese Aussage für alle Elemente des Gitters.
Hinweis:
Nutzen Sie np.linalg.inv um die Inverse einer Matrix zu bestimmen und np.linalg.norm, um die Norm einer Matrix zu berechnen.
In 1d können Sie $\rho(T) = h(T)$ verwenden, $h(\hat{T})$ können Sie analytisch herleiten.
Die Lagrange-Basis erster Ordnung auf dem 1d Referenzelement $\hat{T}$ ist durch die Funktionen $\hat{\varphi}_0, \hat{\varphi}_1: \hat{T} \to \mathbb{R}$ mit $$\begin{align*} \hat{\varphi}_0(\hat{x}) &:= 1 - \hat{x}\\ \hat{\varphi}_1(\hat{x}) &:= \hat{x} \end{align*}$$ für gegeben.
Schreiben Sie eine Funktion lagrangian_p1_basis(T_hat, x_hat), welche die Auswertung von $\hat{\varphi}_0(\hat{x})$ und $\hat{\varphi}_1(\hat{x})$ für alle Punkte $\hat{x} \in $x_hat berechnet und testen Sie ihre Funktion mit Hilfe des Codes in der nächsten Code-Zelle.
Hinweise:
Es gelten dieselben Hinweise wie in Aufgabe 1.1, was die Argumente der Funktion betrifft.
Geben Sie ein numpy.array der Größe (2, len(x_hat)) zurück, sodass die erste Zeile die Auswertungen von $\hat{\varphi}_0$, und die zweite Zeile die Auswertungen von $\hat{\varphi}_1$ enthält.
x_hat = [0, 0.5, 1]
phi = lagrangian_p1_basis(T_hat, x_hat)
print('phi_0({}) = {}'.format(x_hat, phi[0]))
print('phi_1({}) = {}'.format(x_hat, phi[1]))
Stellen Sie in Ihrer Funktion außerdem sicher, dass alle Punkte in x_hat im Referenzelement $\hat{T}$ liegen und testen Sie, dass der folgende Aufruf fehlschlägt! Geben Sie eine informative Fehlermeldung, falls die Punkte nicht im Referenzelement liegen.
Hinweis:
lagrangian_p1_basis(T_hat, -1)
lagrangian_p1_basis(T_hat, 2)
Nach Definition ist die $P^1$-Lagrange Basis eine Nodale (knotenbasierte) Basis, es gilt also: $$ \hat{\varphi}_i(\hat{a}_j) = \delta_{ij}\quad\forall 1\leq i, j \leq 2 $$
print) mit Hilfe von np.eye.matplotlib zusammen in einem Plot, indem Sie die Funktionen an genügend vielen Punkten auswerten. Geben Sie geeignete Legende und Titel an.