Praktikum zur Vorlesung Numerik partieller Differentialgleichungen I

Mario Ohlberger, Felix Schindler

Blatt 03, 13.11.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_03.ipynb.
  • Importieren Sie numpy und pymor.basic und machen Sie matplotlib für das Notebook nutzbar.
In [1]:
%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!

In [2]:
omega = (0., 1.)
grid = OnedGrid(omega, num_intervals=4)

Aufgabe 1: Die Referenzabbildung [Lemma 2.11, Übungszettel 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$.

  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

      • sowohl einzelne Zahlen, also F(grid, T, 0.5), als auch
      • eine Python-Liste von Zahlen, also F(grid, T, [0., 0.5, 1.]), als auch
      • ein numpy-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.

In [3]:
def F(grid, T, x_hat):
    if not isinstance(x_hat, np.ndarray):
        x_hat = np.array(x_hat)
    T_hat = grid.reference_element(0)
    a_hat = T_hat.subentities(0, 1)
    assert np.all(a_hat[0] <= x_hat) and np.all(x_hat <= a_hat[1]), 'Points {} not in reference element!'.format(x_hat)
    A, b = grid.embeddings(0)
    A, b = A[T], b[T]
    return (A*x_hat + b).reshape(len(x_hat))
In [4]:
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)))

T_hat = [0 1]
  F_0(T_hat) = [0.   0.25]
  F_1(T_hat) = [0.25 0.5 ]
  F_2(T_hat) = [0.5  0.75]
  F_3(T_hat) = [0.75 1.  ]

  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:

    • 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.

In [5]:
#F(grid, 0, -1)
In [6]:
#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

  • Knotenpunkte $\hat{a}_i$ des Referenzelementes $\hat{T}$ und
  • Knotenpunkte $a_i$ des Elementes $T$,

für $1 \leq i \leq 2$.

  1. Ü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.

In [7]:
for T in range(grid.size(0)):
    global_vertex_ids = grid.subentities(0, 1)[T]
    a = grid.centers(1)[global_vertex_ids]
    print(np.all(np.isclose(a, F(grid, T, a_hat).reshape(a.shape))))

True
True
True
True

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.

  1. Ü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.

In [8]:
rho_hat = 1
for T in range(grid.size(0)):
    h = grid.diameters(0)[T]
    A, _ = grid.embeddings(0)
    A = A[T]
    print(np.linalg.norm(A) <= h/rho_hat)

True
True
True
True

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*}$$

  1. 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.

In [9]:
def F_inv(grid, T, x):
    if not isinstance(x, np.ndarray):
        x = np.array(x)
    global_vertex_ids = grid.subentities(0, 1)[T]
    a = grid.centers(1)[global_vertex_ids]
    assert np.all(a[0] <= x) and np.all(x <= a[1]), 'Points {} not inside element {}!'.format(x, T)
    A, b = grid.embeddings(0)
    A, b = A[T], b[T]
    x_hat = np.zeros(len(x))
    for i in range(len(x)):
        x_hat[i] = np.linalg.solve(A, x[i] - b)
    return x_hat.reshape(len(x))
In [10]:
print('T_hat = {}'.format(a_hat))
for T in range(grid.size(0)):
    global_vertex_ids = grid.subentities(0, 1)[T]
    a = grid.centers(1)[global_vertex_ids]
    print('  F_inv_{}([{}, {}]) = {}'.format(T, a[0][0], a[1][0], F_inv(grid, T, a)))
    assert np.all(np.isclose(a_hat, F_inv(grid, T, a).reshape(a_hat.shape)))

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.]

  1. 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:

    • Es gelten dieselben hinweise wie in Teilaufgabe 2.
In [11]:
#F_inv(grid, 0, -1)
In [12]:
#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.

  1. Ü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.

In [13]:
h_hat = 1
for T in range(grid.size(0)):
    rho = grid.diameters(0)[T]
    A, _ = grid.embeddings(0)
    A = A[T]
    A_inv = np.linalg.inv(A)
    print(np.linalg.norm(A_inv) <= h_hat/rho)

True
True
True
True

Aufgabe 2: Die $P^1$-Lagrange Basis [Definition 2.13]

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.

  1. 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.

In [14]:
def lagrangian_p1_basis(T_hat, x_hat):
    if not isinstance(x_hat, np.ndarray):
        x_hat = np.array(x_hat)
    a_hat = T_hat.subentities(0, 1)
    assert np.all(a_hat[0] <= x_hat) and np.all(x_hat <= a_hat[1]), 'Points {} not in reference element!'.format(x_hat)
    return np.array([1 - x_hat, x_hat]).reshape(2, len(x_hat))
In [15]:
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]))

phi_0([0, 0.5, 1]) = [1.  0.5 0. ]
phi_1([0, 0.5, 1]) = [0.  0.5 1. ]

  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:

    • Es gelten dieselben hinweise wie in Teilaufgabe 1.2.
In [16]:
#lagrangian_p1_basis(T_hat, -1)
In [17]:
#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 $$

  1. Überprüfen Sie diese Bedingung programmatisch (d.h. nicht allein durch print) mit Hilfe von np.eye.
In [18]:
np.allclose(lagrangian_p1_basis(T_hat, a_hat), np.eye(2, 2))
Out[18]:
True
  1. Visualisieren Sie $\hat{\varphi}_0$ und $\hat{\varphi}_1$ mit Hilfe von matplotlib zusammen in einem Plot, indem Sie die Funktionen an genügend vielen Punkten auswerten. Geben Sie geeignete Legende und Titel an.
In [19]:
x_hat = [0, 1] # two points are enough for linear functions
phi = lagrangian_p1_basis(grid.reference_element(0), x_hat)
plt.plot(x_hat, phi[0], label='phi_0')
plt.plot(x_hat, phi[1], label='phi_1')
plt.legend()
plt.title('P1 Lagrangian basis on reference element')
Out[19]:
Text(0.5, 1.0, 'P1 Lagrangian basis on reference element')