Praktikum zur Vorlesung Numerik partieller Differentialgleichungen I

Mario Ohlberger, Felix Schindler

Blatt 05, 08.01.2020

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

Aufgabe 1

Sei $\Omega \subset \mathcal{R}^d$ für $d = 1, 2, 3$ ein Gebiet mit Lipschitz-Rand und sei $u \in H^1_0(\Omega)$ die schwache Lösung von

$$\begin{align*} -\nabla\cdot(A \nabla u) &= f &\text{in } \Omega \end{align*}$$

für eine gegebene Diffusion $A \in L^\infty(\Omega)$ und rechte Seite $f \in L^2(\Omega)$.

Sei $\mathcal{T}_h$ ein Simplexgitter von $\Omega$ und $S_h^1$ der zugehörige Lagrange-Raum erster Ordnung. Berechnen Sie $u_h \in S_h^1$ als Lösung von

$$\begin{align*} a_h(u_h, v_h) &= l_h(v_h) &&\text{für alle } v_h \in S_h^1, \end{align*}$$

wobei $a_h: S_h^1 \times S_h^1 \to \mathbb{R}$ und $l_h: S_h^1 \to \mathbb{R}$ die diskreten Varianten von $a: H^1 \times H^1 \to \mathbb{R}$ und $l: H^1 \to \mathbb{R}$, gegeben durch

$$\begin{align*} a(u, v) &:= \int_\Omega A(x) \nabla u(x) \nabla v(x) \text{d}x &&\text{und}\\ l(v) &:= \int_\Omega f(x) v(x)\text{d}x, \end{align*}$$

sind, wenn man alle Integrale durch geeignete Quadraturen ersetzt.

Testen Sie ihr Program für

  • $d = 1$
  • $A = 1$
  • $f = 1$

auf einem Gitter mit 4 Elementen.

In [2]:
A = ConstantFunction(1)
f = ConstantFunction(1)

Bemerkungen zur Basis des Lagrange-Raumes $S_h^1$

Sei $\hat{T}$ das $d$-dimensionale Referenzsimplex ($d=1$: Interval, $d=2$: Dreieck, ...) mit Knotenpunkten $\hat{a}_1, \dots, \hat{a}_{d + 1}$ und sei $T$ ein affin äquivalentes Simplex mit bijektiver Referenzabbildung $F_T: \hat{T} \to T$.

Zu jedem $T \in \mathcal{T}_h$ sei $\sigma_T: \{1, \dots, d+1\} \to \{1, \dots, \dim S_h^1\}$ die Abbildung, welche der lokalen Nummer $i \in \{1, \dots, d + 1\}$ eines Knotenpunktes $\hat{a}_i$ die globale Nummer $I \in \{1, \dots, \dim S_h^1\}$ des Knotens $a_I = F_T(\hat{a}_i)$ von $T$ im Gitter $\mathcal{T}_h$ zuordnet.

Sei $\{\hat{\varphi}_i\}_{i = 1}^{d+1}$ die Lagrange-Basis von $\mathbb{P}_1(\hat{T})$, also $\hat{\varphi}_i(\hat{a}_j) = \delta_{ij}$ für $1 \leq i, j, \leq d + 1$ und sei $\{\varphi_I\}_{I = 1}^{\dim S_h^1}$ die Basis des Lagrange-Raumes erster Ordnung $S_h^1$, also $\varphi_I(a_J) = \delta_{IJ}$ für $1 \leq I,J \leq \dim S_h^1$, wobei $a_J$ für $1 \leq J \leq \dim S_h^1$ die Knoten des Gitters $\mathcal{T}_h$ seien. Dann lässt sich jede globale Basisfunktion, eingeschränkt auf ein Element des Gitters, durch eine lokale Basisfunktion mit Hilfe der Referenzabbildung darstellen:

$$\begin{align*} \varphi_I|_T &= \hat{\varphi_i} \circ F_T^{-1}&&\forall 1 \leq I \leq \dim S_h^1\;\;\forall T \in \mathcal{T}_h, &&\text{mit } I = \sigma_T(i) \end{align*}$$
In [3]:
grid = OnedGrid(domain=(0, 1), num_intervals=4)
d = grid.dim
dim_S_h_1 = grid.size(d)

from grid_tools import F
from shape_functions import lagrangian_p1_basis

T_hat = grid.reference_element(0)
points, weights = T_hat.quadrature(order=1)

Assemblierung der rechten Seite

Mit den Überlegungen zur Basis von $S_h^1$ ergibt sich für die Rechte Seite für $1 \leq I \leq \dim S_h^1$ mit $I = \sigma_T(i)$

$$\begin{align*} l(\varphi_I) &:= \int_\Omega f(x)\, \varphi_I(x)\, \text{d}x\\ &= \sum_{T \in \tau_h} \int_T f(x)\, \big(\hat{\varphi}_i \circ F_T^{-1}\big)(x)\, \text{d}x\\ &= \sum_{T \in \tau_h} \int_{\hat{T}} |\det\nabla F_T|\; \big(f \circ F_T\big)(\hat{x})\, \big(\hat{\varphi}_i \circ F_T^{-1} \circ F_T\big)(\hat{x}) \,\text{d}\hat{x}\\ &= \sum_{T \in \tau_h} \int_{\hat{T}} |\det\nabla F_T|\; \big(f \circ F_T\big)(\hat{x})\, \hat{\varphi}_i(\hat{x}) \,\text{d}\hat{x}. \end{align*}$$

Ersetzen wir das Integral durch eine Quadratur geeigneter Ordnung auf dem Referenzelement mit $Q \in \mathbb{N}$ Quadraturpunkten $\hat{x}_q \in \hat{T}$ und Quadraturgewichten $\omega_q > 0$, für $1 \leq q \leq Q$, so erhalten wir den $I$-ten Eintrag des assemblierten Vektors der rechten Seite, $\underline{l_h} \in \mathbb{R}^{\dim S_h^1}$,

$$\begin{align*} (\underline{l_h})_I := \sum_{T \in \tau_h} \int_{\hat{T}} \sum_{q = 1}^Q |\det\nabla F_T|\; \omega_q\; f\big(F_T(\hat{x}_q)\big)\; \hat{\varphi}_i(\hat{x}_q) &&\text{für } 1 \leq I \leq \dim S_h^1. \end{align*}$$

Es gibt mindestens zwei Möglichkeiten, den Vektor $\underline{l_h}$ algorithmisch zu assemblieren.

Variante 1:

  • für alle $0 \leq I \leq \text{dim}S_h^1$
    • für alle $T \in \mathcal{T}$
      • für alle $1 \leq i \leq d + 1$
        • finde $i$ zu $I$ und $T$
        • füge dem Eintrag $(\underline{f_h})_I$ den entsprechenden Wert hinzu

Komplexität Variante 1: $\mathcal{O}\big(\dim S_h^1 \cdot |\mathcal{T}_h| \cdot (d + 1) \big) = \mathcal{O}(|\mathcal{T}_h|^2)$

Variante 2:

  • für alle $T \in \mathcal{T}_h$:
    • für alle $0 \leq i \leq d + 1$:
      • finde $I$ zu $i$ und $T$
      • füge dem Eintrag $(\underline{l_h})_I$ den entsprechenden Wert hinzu

Komplexität Variante 2: $\mathcal{O}\big(|\mathcal{T}_h| \cdot (d+1)\big) = \mathcal{O}(|\mathcal{T}_h|)$

In [4]:
l_h = np.zeros(dim_S_h_1)

for T in range(grid.size(0)):
    det_grad_F = grid.integration_elements(0)[T]
    for x_q, w_q in zip(points, weights):
        x = F(grid, T, x_q)
        f_x = f.evaluate(x)
        phi_x_q = lagrangian_p1_basis(T_hat, x_q)
        for i in range(d + 1):
            I = grid.subentities(0, 1)[T][i]
            l_h[I] += det_grad_F * w_q * f_x * phi_x_q[i]

print(f'l_h = {l_h}')

l_h = [0.125 0.25  0.25  0.25  0.125]

Bemerkungen zum Gradienten der Lagrange-Basis von $S_h^1$

Mit obiger Definition lässt sich der Gradient jeder globale Basisfunktion, eingeschränkt auf ein Element des Gitters, durch den Gradienten einer lokale Basisfunktion, mit Hilfe der Referenzabbildung und der Kettenregel, darstellen:

$$\begin{align*} \nabla \varphi_I|_T &= \nabla \big( \hat{\varphi}_i \circ F_T^{-1} \big)\\ &= \nabla F_T^{-1}\,\cdot\,\nabla \hat{\varphi}_i \circ F_T^{-1}&&\forall 1 \leq I \leq \dim S_h^1\;\;\forall T \in \mathcal{T}_h, &&\text{mit } I = \sigma_T(i) \end{align*}$$

Assemblierung der linken Seite

Mit den Überlegungen zum Gradienten der Basis von $S_h^1$ ergibt sich für die linke Seite für $1 \leq I,J \leq \dim S_h^1$ mit $I = \sigma_T(i)$ und $J = \sigma_T(j)$ und durch den Umstand dass $\nabla F_T^{-1}$ konstant ist

$$\begin{align*} a(\varphi_I, \varphi_J) &:= \int_\Omega A(x)\,\nabla\varphi_I(x)\,\nabla\varphi_J(x)\,\text{d}x\\ &= \sum_{T \in \mathcal{T}_h} \int_{T}A(x)\,\underbrace{\nabla F_T^{-1}(x)}_{\text{const.}}\,\big(\nabla \hat{\varphi}_i \circ F_T^{-1}\big)(x) \,\underbrace{\nabla F_T^{-1}(x)}_{\text{const.}}\,\big(\nabla \hat{\varphi}_j \circ F_T^{-1}\big)(x)\,\text{d}x\\ &= \sum_{T \in \mathcal{T}_h} \int_{\hat{T}}|\text{det}\nabla F_T|\,\big(A \circ F_T\big)(\hat{x})\,\nabla F_T^{-1}\,\big(\nabla \hat{\varphi}_i \circ F_T^{-1} \circ F_T\big)(\hat{x}) \,\nabla F_T^{-1}\,\big(\nabla \hat{\varphi}_j \circ F_T^{-1} \circ F_T\big)(\hat{x})\,\text{d}\hat{x}\\ &= \sum_{T \in \mathcal{T}_h} \int_{\hat{T}}|\text{det}\nabla F_T|\,\big(A \circ F_T\big)(\hat{x})\,\nabla F_T^{-1}\,\nabla \hat{\varphi}_i(\hat{x}) \,\nabla F_T^{-1}\,\nabla \hat{\varphi}_j(\hat{x})\,\text{d}\hat{x}. \end{align*}$$

Ersetzen wir das Integral durch eine Quadratur geeigneter Ordnung auf dem Referenzelement mit $Q \in \mathbb{N}$ Quadraturpunkten $\hat{x}_q \in \hat{T}$ und Quadraturgewichten $\omega_q > 0$, für $1 \leq q \leq Q$, so erhalten wir den $I,J$-ten Eintrag der assemblierten Matrix der linken Seite $\underline{a_h} \in \mathbb{R}^{\dim S_h^1 \times \dim S_h^1}$,

$$\begin{align*} (\underline{a_h})_{I,J} := \sum_{T \in \tau_h} \sum_{q = 1}^Q |\det\nabla F_T|\; \omega_q\; A\big(F_T(\hat{x}_q)\big)\; \nabla F_T^{-1}\,\nabla \hat{\varphi}_i(\hat{x}_q)\,\nabla F_T^{-1}\,\nabla \hat{\varphi}_j(\hat{x}_q) &&\text{für } 1 \leq I,J \leq \dim S_h^1. \end{align*}$$
In [5]:
from shape_functions import lagrangian_p1_basis_jacobians

a_h = np.zeros((dim_S_h_1, dim_S_h_1))

for T in range(grid.size(0)):
    det_grad_F = grid.integration_elements(0)[T]
    grad_F = grid.embeddings(0)[0][T]
    for x_q, w_q in zip(points, weights):
        x = F(grid, T, x_q)
        A_x = A.evaluate(x)
        # entweder
        #grad_phi_x_q = lagrangian_p1_basis_jacobians(T_hat, x_q)
        # oder
        grad_phi_x_q = np.array([[-1], [1]])
        for i in range(d + 1):
            I = grid.subentities(0, 1)[T][i]
            for j in range(d + 1):
                J = grid.subentities(0, 1)[T][j]
                a_h[I][J] += det_grad_F * w_q * A_x * (grad_phi_x_q[i]/grad_F) * (grad_phi_x_q[j]/grad_F)

print(f'a_h =\n{a_h}')

a_h =
[[ 4. -4.  0.  0.  0.]
 [-4.  8. -4.  0.  0.]
 [ 0. -4.  8. -4.  0.]
 [ 0.  0. -4.  8. -4.]
 [ 0.  0.  0. -4.  4.]]

Randwerte

Die assemblierte linke und rechte Seite wurde bisher ohne die Beachtung der Randwerte assembliert. Um die homogene Dirichlet-Randbedingung auf $\partial\Omega$ zu erzwingen, betrachten wir die Basisdarstellung der Lösung,

$$\begin{align*} u_h = \sum_{I = 1}^{\dim S_h^1} (\underline{u_h})_I \varphi_I &&\text{wobei}&&\underline{u_h} \in \mathbb{R}^{\dim S_h^1} \end{align*}$$

der zu $u_h \in S_h^1$ korrespondierende Vektor an Freiheitsgerade sei. Da es sich bei $S_h^1$ um einen Lagrange-Raum handelt, sind die Freiheitsgerade $(\underline{u_h})_I$ gerade die Werte der Funktion an einem Gitterpunkt, $u_h(a_I)$.

Aus der Randbedingung $u_h|_{\partial\Omega} = 0$ (im schwachen Sinne) folgt also dass

$$\begin{align*} (\underline{u_h})_I = u_h(a_I) = 0&&\text{für alle Inidizes } 1 \leq I \leq \dim S_h^1\text{, sodass }a_I \in \partial\Omega \end{align*}$$

ein Randknoten ist. Zu diesen Indices ersetzten wir also die $I$-te Zeile der Matrix $\underline{a_h}$ durch eine Einheitzeile und setzen den $I$-ten Eintrag von $\underline{l_h}$ auf 0.

In einer Raumdimension sind das gerade der erste und letzte Knoten des Gitters.

In [6]:
if d == 1:
    boundary_dofs = (0, grid.size(d) - 1)

for I in boundary_dofs:
    a_h[I] *= 0
    a_h[I][I] = 1
    l_h[I] = 0

print(f'a_h =\n{a_h}')
print(f'l_h = {l_h}')

a_h =
[[ 1. -0.  0.  0.  0.]
 [-4.  8. -4.  0.  0.]
 [ 0. -4.  8. -4.  0.]
 [ 0.  0. -4.  8. -4.]
 [ 0.  0.  0. -0.  1.]]
l_h = [0.   0.25 0.25 0.25 0.  ]

Lösung des linearen Gleichungssystems

Wir erhalten den Vektor an Freiheitsgeraden der gesuchten Lösung $u_h$ als Lösung des Gleichungssystems

$$\underline{a_h} \cdot \underline{u_h} = \underline{l_h}.$$
In [7]:
u_h = np.linalg.solve(a_h, l_h)
print(f'u_h = {u_h}')

u_h = [-0.       0.09375  0.125    0.09375  0.     ]
In [8]:
from visualizations import visualize_lagrange_p1

visualize_lagrange_p1(grid, u_h, f'{grid.size(0)} grid elements', 'P1 solution')

Aufgabe 2

Bestimmen Sie den EOC von $\|u - u_h\|_{L^2(\Omega)}$ für eine gegebene Referenzlösung $u \in H^1_0(\Omega)$.

Testen Sie ihr Programm für das Beispiel aus Aufgabe 1 und Gitter mit $2^2, 2^3, 2^4, 2^5$ Elementen, indem Sie $u$ analytisch bestimmen, und senden Sie ihr Jupyter Notebook unter Angabe aller bearbeitenden an Felix Schindler.