Praktikum zur Vorlesung Numerik partieller Differentialgleichungen I

Mario Ohlberger, Felix Schindler

Blatt 04, 04.12.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_04.ipynb.
• Importieren Sie numpy und pymor.basic und machen Sie matplotlib für das Notebook nutzbar.

%autosave 0
%matplotlib notebook
import numpy as np
from pymor.basic import *
from matplotlib import pyplot as plt

Autosave disabled

Aufgabe 1: Programmteile wiederverwendbar machen

1. Machen Sie Teile von Blatt 2 wiederverwedbar, indem Sie

   • eine Datei interpolations.py erstellen und darin
     • eine Funktion interpolate_lagrange_p1(grid, f) erstellen, welche die Interpolation im Sinne von Blatt 2, Aufgabe 1 berechnet und zurück gibt, und
   • eine Datei visualizations.py erstellen und darin
     • eine Funktion visualize_lagrange_p1(grid, f_h, name, title) erstellen, welche eine beliebige Menge von diskreten Funktionen im Sinne von Blatt 2, Aufgabe 2 visualisiert.

cat interpolations.py

def interpolate_lagrange_p1(grid, f):
    vertices = grid.centers(1)
    num_vertices = len(vertices)
    values = f.evaluate(vertices)
    assert len(values) == len(vertices)
    return values

cat visualizations.py

from matplotlib import pyplot as plt

def visualize_lagrange_p1(grid, f_h, name='f_h', title='title'):
    colors = list(plt.get_cmap('tab20c').colors)
    colors.reverse()

    plt.figure()
    if not isinstance(grid, list):
        grid = [grid,]
    if not isinstance(f_h, list):
        f_h = [f_h,]
    if not isinstance(name, list):
        name = [name,]
    assert len(grid) == len(f_h) == len(name)
    for gg, ff, nn in zip(grid, f_h, name):
        assert gg.dim == 1
        assert len(ff) == gg.size(1)
        plt.plot(gg.centers(1), ff, color=colors.pop(), label=nn)
    plt.legend()
    plt.title(title)

1. Testen Sie diese Funktionen für ein Gitter des Gebiets $\Omega = (0, 1)$ mit 4 Elementen, indem Sie die Funktionen aus den beiden Dateien importieren und aufrufen.

Hinweis: Nennen Sie das Gebiet omega, das Gitter coarse_grid, die Funktion f und ihre Interpolierte f_H.

omega = (0., 2.*np.pi)
coarse_grid = OnedGrid(omega, num_intervals=4)

from interpolations import interpolate_lagrange_p1
from visualizations import visualize_lagrange_p1

f = ExpressionFunction('sin(x)')
f_H = interpolate_lagrange_p1(coarse_grid, f)
visualize_lagrange_p1(coarse_grid, f_H, 'f_H', 'Interpolation on coarse grid')

1. Machen Sie Teile von Blatt 3 wiederverwedbar, indem Sie

   • eine Datei grid_tools.py erstellen und darin
     • eine Funktion F(grid, T, x_hat) erstellen, welche die Auswertung der Referenzabbildung im Sinne von Blatt 3, Aufgabe 1 berechnet und zurück gibt, und
     • eine Funktion F_inv(grid, T, x) erstellen, welche die Auswertung der inversen Referenzabbildung im Sinne von Blatt 3, Aufgabe 1 berechnet und zurück gibt und
   • eine Datei shape_functions.py erstellen und darin
     • eine Funktion lagrangian_p1_basis(T_hat, x_hat) erstellen, welche die Auswertung der $P^1$-Basis im Sinne von Blatt 3, Aufgabe 2 berechnet und zurück gibt.

cat grid_tools.py

import numpy as np

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))

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))

cat shape_functions.py

# +
import numpy as np

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 + 1e-7)) and np.all((x_hat - 1e-7) <= a_hat[1]), 'Points {} not in reference element!'.format(x_hat)
    return np.array([1 - x_hat, x_hat]).reshape(2, len(x_hat))

def lagrangian_p1_basis_jacobians(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 + 1e-7)) and np.all((x_hat - 1e-7) <= a_hat[1]), 'Points {} not in reference element!'.format(x_hat)
    return np.array([-1*np.ones(x_hat.shape), np.ones(x_hat.shape)]).reshape(2, len(x_hat))

Aufgabe 2: Integration

Schreiben Sie eine Funktion L2_norm_lagrange_p1(grid, f_h), welche $\|f_h\|_{L^2(\Omega)}$ berechnet, wobei f_h der Vektor der Freiheitsgrade zum Gitter grid sei und testen Sie Ihre Funktion mit unten stehendem Code.

Leiten Sie dazu zuerst eine entsprechende Integrationsformel her:

• Lokalisieren Sie das Integral bezüglich der Elemente des Gitters
• Nutzen Sie die Referenzabbildung
• Nutzen Sie eine Quadratur

Hinweis:

• Sie erhalten die Quadraturpunkte und Gewichte mit Hilfe von quadrature von einem Referenzelement.
• Sie erhalten die aus dem Transformationssatz resultierende Skalierung mit Hilfe von integration_elements vom Gitter

from shape_functions import lagrangian_p1_basis

def L2_norm_lagrange_p1(grid, f_h):
    assert len(f_h) == grid.size(1)
    l2_norm_sqrd = 0.
    T_hat = grid.reference_element(0)
    points, weights = T_hat.quadrature(order=2)
    for T in range(grid.size(0)):
        ie = grid.integration_elements(0)[T]
        for x_q, w_q in zip(points, weights):
            phi_x_q = lagrangian_p1_basis(T_hat, x_q)
            global_vertex_ids = grid.subentities(0, 1)[T]
            f_h_x_q = phi_x_q.ravel().dot(f_h[global_vertex_ids])
            l2_norm_sqrd += ie * w_q * f_h_x_q**2
    return np.sqrt(l2_norm_sqrd)

unit_interval = OnedGrid()
one = np.ones(unit_interval.size(1))

assert L2_norm_lagrange_p1(unit_interval, one) == 1

Mit der Definition der Lagrange-Basis von Blatt 5 ergibt sich aus der Basisdarstellung von $f_h$, dass

$$\begin{align*} f_h|_T = \sum_{I=1}^{\dim S_h^1}(\underline{f_h})_I\; \varphi_I|_T = \sum_{i = 1}^{d + 1} (\underline{f_h})_I\; \hat{\varphi}_i \circ F_T&&\text{mit } I = \sigma_T(i). \end{align*}$$

Damit ergibt sich für die $L^2$-Norm

$$\begin{align*} \|f_h\|_{L^2}^2 &= \int_\Omega f_h(x)^2 \text{d}x\\ &= \sum_{T \in \tau_h} \int_T f_h(x)^2 \text{d}x\\ &= \sum_{T \in \tau_h} \int_T \Big[\sum_{i = 1}^{d + 1} (\underline{f_h})_I\; \big(\hat{\varphi}_i \circ F_T\big)(x)\Big]^2 \text{d}x\\ &= \sum_{T \in \tau_h} \int_{\hat{T}} |\det\nabla F|\; \Big[\sum_{i = 1}^{d + 1} (\underline{f_h})_I\; \big(\hat{\varphi}_i \circ F_T \circ F_T^{-1}\big)(\hat{x})\Big]^2 \text{d}\hat{x}\\ &= \sum_{T \in \tau_h} \int_{\hat{T}} |\det\nabla F|\; \big(\sum_{i = 1}^{d + 1}(\underline{f_h})_I\; \hat{\varphi}_i(\hat{x})\big)^2\text{d}\hat{x}\\ &\approx \sum_{T \in \tau_h} \sum_{q = 1}^Q \omega_q\, |\det\nabla F|\; \big(\sum_{i = 1}^{d + 1}(\underline{f_h})_I \;\hat{\varphi}_i(\hat{x}_q)\big)^2 \end{align*}$$

Aufgabe 3: Prolongation diskreter Funktionen

Schreiben Sie eine Funktion prolong_lagrange_p1(coarse_grid, coarse_DoFs, fine_grid), welche eine diskrete Funtion (gegeben durch den Vektor an Freiheitsgeraden coarse_DoFs bezüglich des Gitters coarse_grid) auf einem feinerem Gitter fehlerfrei interpoliert.

• Schreiben Sie dazu analog zu Ihrer interpolate_lagrange_p1 Funktion eine Interpolation bezüglich des feinen Gitters.
• Suchen sie an jedem Interpolationspunkt des feinen Gitters das Element des groben Gitters, welches das feine Element enthält.
• Nutzen Sie die Basisdarstellung der diskreten Funktion auf dem groben Gitter, um diese an diesem Punkt auszuwerten.

from grid_tools import F, F_inv

def prolong_lagrange_p1(coarse_grid, coarse_DoFs, fine_grid):
    assert len(coarse_DoFs) == coarse_grid.size(1)
    fine_vertices = fine_grid.centers(1)
    f_h = np.zeros(len(fine_vertices))
    for i, x in enumerate(fine_vertices):
        found_element = False
        for T in range(coarse_grid.size(0)):
            global_vertex_ids = coarse_grid.subentities(0, 1)[T]
            a = coarse_grid.centers(1)[global_vertex_ids]
            if a[0] <= x <= a[1]:
                found_element = True
                break # this is the correct element
        assert found_element, 'Point {} is not contained in grid!'.format(x)
        x_hat = F_inv(coarse_grid, T, x)
        phi_x = lagrangian_p1_basis(coarse_grid.reference_element(0), x_hat)
        f_x = phi_x.ravel().dot(coarse_DoFs[global_vertex_ids])
        f_h[i] = f_x
    return f_h

fine_grid = OnedGrid(omega, num_intervals=8)

f_h = prolong_lagrange_p1(coarse_grid, f_H, fine_grid)

visualize_lagrange_p1([coarse_grid, fine_grid], [f_H, f_h], ['f_H', 'f_h'])

Aufgabe 4: Interpolationsfehler

Bestimmen Sie den Interpolationsfehler, indem Sie zuerst

• ein feines Referenzgitter mit 256 Elementen erstellen,
• die Interpolation von $f$ auf diesem Gitter als Referenzlösung berechnen,

und dann für Gitter der Feinheit 4, 8, 16, 32 jeweils

• die Interpolation von $f$ auf das Gitter berechnen,
• die Interpolation auf das Referenzgitter prolongieren
• die $L^2$-Norm der Differenz zur Referenzlösung berechnen.

Visualisieren Sie die Differenzen in einem Plot.

I = 4
num_refs = 4

reference_grid = OnedGrid(omega, num_intervals=I * 2**(num_refs + 2))
f_ref = interpolate_lagrange_p1(reference_grid, f)

differences = []
labels = []
grid_sizes = []
grid_widths = []

for _ in range(num_refs):
    print(f'prolonging from grid of size {I}...')
    grid = OnedGrid(omega, num_intervals=I)
    grid_sizes.append(grid.size(0))
    grid_widths.append(np.max(grid.diameters(0)))
    f_h = interpolate_lagrange_p1(grid, f)
    f_h = prolong_lagrange_p1(grid, f_h, reference_grid)
    differences.append(f_ref - f_h)
    labels.append(f'f - f_h ({I} elements)')
    I *= 2

visualize_lagrange_p1([reference_grid] * num_refs, differences, labels, title='interpolation error')

prolonging from grid of size 4...
prolonging from grid of size 8...
prolonging from grid of size 16...
prolonging from grid of size 32...

Aufgabe 5: Interpolationsgüte

Wir wissen nach Satz 2.34, dass $$ \|f - f_h\|_{L^2(\Omega)} \leq C h^k $$ für eine geeignete Konstante $C > 0$ (die nicht vom Gitter abhängt). Bestimmen Sie die Konvergenzordnung $k$ experimentell, indem Sie den EOC (experimental order of convergence) berechnen, der folgendermaßen definiert ist: $$ EOC^{(\nu)} := \frac{\ln\Big(\frac{\big\|f - f_h^{(\nu)}\big\|_{L^2(\Omega)}}{\big\|f - f_h^{(\nu - 1)}\big\|_{L^2(\Omega)}}\Big)}{\ln\Big(\frac{h^{(\nu)}}{h^{(\nu - 1)}}\Big)} $$ Dabei bezeichnet $\nu = 0, 1, 2, \dots$ das Level der Verfeinerung und $f_h^{(\nu)}$ die Interpolierte auf das jeweilige Gitter mit Gitterweite $h^{(\nu)}$.

Visualisieren Sie die Fehler geeigneter Weise.

errors = []

for diff in differences:
    l2_norm_sqrd = 0.
    T_hat = reference_grid.reference_element(0)
    points, weights = T_hat.quadrature(order=9)
    #print(f'points = {points}')
    #print(f'weights = {weights}')
    for T in range(reference_grid.size(0)):
        ie = reference_grid.integration_elements(0)[T]
        for x_q, w_q in zip(points, weights):
            # evaluate difference in x_q
            phi_x_q = lagrangian_p1_basis(T_hat, x_q)
            global_vertex_ids = reference_grid.subentities(0, 1)[T]
            diff_x_q = phi_x_q.ravel().dot(diff[global_vertex_ids])
            l2_norm_sqrd += ie * w_q * diff_x_q**2
    errors.append(np.sqrt(l2_norm_sqrd))

plt.figure()
plt.loglog(grid_sizes, errors, label='L^2 error')
plt.xlabel('grid sizes')
plt.ylabel('error')
plt.legend()

<matplotlib.legend.Legend at 0x7f7f72a52f10>

for i in range(len(errors) - 1):
    e_coarse = errors[i]
    e_fine = errors[i + 1]
    h_coarse = grid_widths[i]
    h_fine = grid_widths[i + 1]
    print(np.log(e_fine / e_coarse) / np.log(h_fine / h_coarse))

1.9422269563568915
1.9889383083363068
2.0105654872253864