Praktikum zur Vorlesung Numerik partieller Differentialgleichungen I
Mario Ohlberger, Felix Schindler
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
numpyundpymor.basicund machen Siematplotlibfü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
verticesab. - Werten Sie die Funktion
fand allen Knotenpunkten aus (siehe Dokumentation von FunctionInterface) und speichern Sie die Werte invalues. - Stellen Sie sicher, dass Sie genauso viele Funktionswerte wie Knotenpunkte berechnet haben.
- Geben Sie
valueszurück.
def interpolate(grid, f):
assert grid.dim == 1
vertices = grid.centers(1)
num_vertices = len(vertices)
# loop-based variant
#values = np.zeros(num_vertices)
#for i in range(num_vertices):
# values[i] = f.evaluate(vertices[i])
# vectorized variant
values = f.evaluate(vertices)
assert len(values) == len(vertices)
return values
Ein äquidistantes Gitter in einer Raumdimension wird in pyMOR durch die Klasse
OnedGridmodelliert.- Legen Sie ein entsprechendes Objekt mit Namen
gridan, - geben Sie die Dimension des Gitters,
- die Anzahl der Knotenpunkte sowie
- die Position jedes Knotenpunktes aus (Wiederholung von Blatt 00).
- Legen Sie ein entsprechendes Objekt mit Namen
grid = OnedGrid(domain=(0, 2*np.pi), num_intervals=4)
dim = grid.dim
print('grid lives in a {}-dimension space'.format(dim))
I = grid.size(0)
h = 2*np.pi / I
print('grid has {} elements, h = {}'.format(I, h))
vertices = grid.centers(dim)
print('grid has {} vertices:'.format(len(vertices)))
for i in range(len(vertices)):
x_i = vertices[i]
print(' x_{} = {}'.format(i, x_i))
- Analytische Funktionen werden in pyMOR durch
ExpressionFunctionmodelliert. Legen Sie ein entsprechendes Objekt mit Namenfan.
f = ExpressionFunction('sin(x)', dim_domain=1, shape_range=())
Berechnen Sie die Interpolation von $f$ auf dem Gitter $\mathcal{T}_h$ mittels Ihrer
interpolateFunktion.- Rufen Sie Ihre Funktion auf und speichern Sie das Resultat in
f_hab. - Geben Sie den Python-Typ von
f_haus. - Geben Sie die Länge von
f_haus. - Geben Sie mit Hilfe einer
for-Schleife die Koordinaten aller Knotenpunkte des Gitters zusammen mit den Werten von $f$ aus (Hinweis: Nutzen Siezip).
- Rufen Sie Ihre Funktion auf und speichern Sie das Resultat in
f_h = interpolate(grid, f)
print(type(f_h))
print(f_h.shape)
print(len(f_h))
print()
for x, value in zip(vertices, f_h):
print('f_h({}) = {}'.format(x, value))
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 pltunter dem Namenpltverfügbar und machen Sie sich mitplt.plotundplt.legendvertraut.
from matplotlib import pyplot as plt
# Dokumentation von plot anzeigen
#plt.plot?
# Eines der Beispiele aus der Dokumentation, x = [0, 1, 2, 3], y = x^2
#plt.plot([0, 1, 2, 3], [0, 1, 4, 9])
# Dokumentation von legend anzeigen
#plt.legend?
# Eine Möglichkeit besteht darin, bei plot direkt ein label mit zu geben, anstelle des Aufrufs von oben nutzt man also zB
plt.plot([0, 1, 2, 3], [0, 1, 4, 9], label='f')
plt.legend()
# Andere nützliche Funktionen
plt.title('Plot von f(x) = x^2')
plt.xlabel('x')
plt.ylabel('f(x)')
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_hmit der Anzahl der Knotenpunkte des Gitters überein stimmt. - Nutzen Sie die
plotFunktion aus dempltModul, um die Werte vonf_han den Knotenpunktes des Gitters zu visualisieren und linear zu interpolieren, und den Namen der Funktion zu übergeben. - Nutzen Sie die
legendFunktion aus dempltModul, um die Legende anzeigen zu lassen.
- Erstellen Sie in der Funktion mit
def visualize(grid, f_h, name=''):
plt.figure()
assert grid.dim == 1
assert len(f_h) == grid.size(1)
plt.plot(grid.centers(1), f_h, label=name)
plt.legend()
- Testen Sie Ihre Funktion für das Gitter und die Interpolation aus Aufgabe 1.
visualize(grid, f_h, 'f_h')
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_hundname, ob es sich um eine Python Liste handelt (Hinweis: nutzen Sie dafürisinstance), 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 vonplotmit Hilfe voncolor=colors.pop()jeder visualisierung eine andere Farbe. - Rufen Sie danchach in Ihrer Funktion einmal
legend()auf, um die Legend anzeigen zu lassen.
colors = list(plt.get_cmap('tab20c').colors)
colors.reverse()
def visualize(grid, f_h, name=''):
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()
- 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).
grids = []
funcs = []
names = []
for I in (2**2, 2**3, 2**4, 2**5):
gg = OnedGrid(domain=(0, 2*np.pi), num_intervals=I)
grids.append(gg)
funcs.append(interpolate(gg, f))
names.append('f_h ({} grid elements)'.format(I))
visualize(grids, funcs, names)
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.
Lösung: Obige Taylor-Entwicklung ausgewertet an $x_{i - 1}$ und $x_{i + 1}$ ergibt $$\begin{align*} u(x_{i - 1}) &= u(x_i) + u'(x_i)\cdot -h + \tfrac{1}{2}u''(x_i)\cdot h^2 + \mathcal{O}(h^3)\\ u(x_{i + 1}) &= u(x_i) + u'(x_i)\cdot h + \tfrac{1}{2}u''(x_i)\cdot h^2 + \mathcal{O}(h^3). \end{align*}$$ Nach Addition bleiben nur noch Terme $0$-ter und $2$-ter Ordnung, $$\begin{align*} u(x_{i - 1}) + u(x_{i + 1}) = 2 u(x_i) + u''(x_i) + \mathcal{O}(h^3), \end{align*}$$ sodass wir nach Umordnen und Division durch $h^2$ den Differenzenquotienten $$\begin{align*} u' '(x_{i}) = \frac{u(x_{i - 1}) - 2u(x_i) + u(x_{i + 1})}{h^2} + \mathcal{O}(h) \end{align*}$$ für die zweite Ableitung von $u$ in einem Gitterpunkt mit Fehler der Ordnung $h$ erhalten.
Lösung: Dadurch motiviert definieren wir als Approximation des Differentialoperators $\partial_{xx}$ den diskreten Operator auf dem Unterraum der diskreten Funktionen mit Null-Randwerten (obiger Differenzenquotient ist nur für innere Gitterpunkte sinnvoll) $$\begin{align*} \partial_{h, h}: S_h^1 \cap H^1_0(\Omega) &\to S_h^1 \cap H^1_0(\Omega),\\ u_h &\mapsto v_h := \partial_{h, h}\; u_h, \end{align*}$$ der Aufgrund oben genannter Isomorphie zwischen $S_h^1$ und $\mathbb{R}^{I 1}$ eindeutig durch die Angabe des DoF-Vektors von $v_h$ definiert ist (zur Erinnerung: $\underline{u_h}$ ist der DoF-Vektor von $u_h$): $$\begin{align*} (\underline{v_h})_i := \begin{cases} h^{-2} \big((\underline{u_h})_{i - 1} - 2\, (\underline{u_h})_i + (\underline{u_h})_{i + 1}\big)&\text{, für } 1 \leq i < I,\\ 0&\text{, für } i \in \{0, I\} \end{cases}. \end{align*}$$
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.
Lösung: Gesucht ist $u_h \in S_h^1$, sodass $$\begin{align*} -\partial_{h, h}\, u_h(x_i) &= f(x_i) &&\text{für } 1 \leq i < I,\\ u_h(x_i) &= g_\text{D}(x_i) &&\text{für } i \in \{0, I\}. \end{align*}$$
Daraus ergeben sich die $I - 1$ Gleichungen für die inneren Knoten, $$\begin{align*} -\frac{u_h(x_{i - 1})}{h^2} + 2\, \frac{u_h(x_i)}{h^2} - \frac{u_h(x_{i + 1}}{h^2} = f(x_i)&&\text{für } 1 \leq i < I, \end{align*}$$ und die zwei Gleichungen für die Randknoten, $$\begin{align*} u_h(x_0) = g_\text{D}(x_0) && u_h(x_{I + 1}) = g_\text{D}(x_{I + 1}). \end{align*}$$
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.
Lösung: Gesucht ist $\underline{u_h} \in \mathbb{R}^{I + 1}$, sodass $$\begin{align*} \underline{b_h} \cdot \underline{u_h} = \underline{f_h}, \end{align*}$$ mit $$\begin{align*} (\underline{b_h})_{ij} := \begin{cases} 1,&\text{für } j = i, i \in \{0, I\},\\ -\tfrac{1}{h^2},&\text{für } j = i - 1, 1 \leq i < I,\\ \tfrac{2}{h^2},&\text{für } j = i, 1 \leq i < I,\\ -\tfrac{1}{h^2},&\text{für } j = i + 1, 1 \leq i < I,\\ 0,&\text{sonst} \end{cases}&&\text{und}&& (\underline{f_h})_j := \begin{cases} g_\text{D}(x_i),&\text{für } i \in \{0, I\},\\ f(x_i),&\text{für } 1 \leq i < I\\ \end{cases}. \end{align*}$$
- 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 IhrerinterpolateFunktion.
I = grid.size(0)
h = 2*np.pi / I
b_h = np.zeros((I + 1, I + 1))
f_h = np.zeros(I + 1)
for i in range(1, I):
b_h[i][i - 1] = -1/h**2
b_h[i][i] = 2/h**2
b_h[i][i + 1] = -1/h**2
f_h[i] = f.evaluate(vertices[i])
b_h[0][0] = 1
f_h[0] = 0.
b_h[I][I] = 1
f_h[I] = 0.
print(b_h)
print()
print(f_h)
- Lösen Sie das resultierende Gleichungssystem mit Hilfe von numpy.linalg.solve (Hinweis: numpy ist als
npverfügbar).
u_h = np.linalg.solve(b_h, f_h)
visualize(grid, u_h, 'u_h')
- 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 sollg_Dein Python Tuple mit zwei Werten sein, also in diesem Beispielg_D = (0, 0)). Visualisieren Sie die Lösungen mit Hilfe dervisualizeFunktion aus Aufgabe 2.
def discretize_elliptic_1d_fd(grid, f, g_D):
vertices = grid.centers(1)
h = (grid._domain[1] - grid._domain[0]) / I
b_h = np.zeros((I + 1, I + 1))
f_h = np.zeros(I + 1)
for i in range(1, I):
b_h[i][i - 1] = -1/h**2
b_h[i][i] = 2/h**2
b_h[i][i + 1] = -1/h**2
f_h[i] = f.evaluate(vertices[i])
b_h[0][0] = 1
f_h[0] = g_D[0]
b_h[I][I] = 1
f_h[I] = g_D[1]
u_h = np.linalg.solve(b_h, f_h)
return u_h
omega = (0, 1)
f = ConstantFunction(1.)
g_D = (0., 0.)
grids = []
solutions = []
names = []
for I in (2**2, 2**3, 2**4, 2**5):
grid = OnedGrid(omega, num_intervals=I)
u_h = discretize_elliptic_1d_fd(grid, f, g_D)
grids.append(grid)
solutions.append(u_h)
names.append('u_h ({} elements)'.format(I))
visualize(grids, solutions, names)
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}$
omega = (0, 1)
I = 2**7
for f, g_D, name in ((ConstantFunction(1), (0., 1.), 'u_h (case 1)'),
(ExpressionFunction('(0.2 <= x[..., 0]) * (x[..., 0] <= 0.4) * 1.'), (0., 0.), 'u_h (case 2)')):
grid = OnedGrid(omega, num_intervals=I)
u_h = discretize_elliptic_1d_fd(grid, f, g_D)
visualize(grid, u_h, name)
- Messen Sie die benötigte Zeit sowohl für die Erstellung des Gitters (inklusive einem Aufruf von
grid.centers(1)), die Assemblierung als auch für die Lösung des Problems für verschieden feine Gitter und visualisieren Sie die Zeiten in Abhängigkeit von $I$. Experimentieren Sie mit verschiedenen Algorithmen zur Assemblierung (for-Schleifen, Vektorisiert sparse).
import time
assemble_data = {}
solve_data = {}
def measure_timings(discretizer, name, max_refinements=13):
omega = (0, 1)
f = ConstantFunction(1)
g_D = (0, 0)
data = {'I': [], 'grid': [], 'assemble': [], 'solve': []}
assemble_data[name] = ([], [])
solve_data[name] = ([], [])
for ref in range(2, max_refinements + 1):
I = 2**ref
data['I'].append(I)
timings = discretizer(omega, I, f, g_D)
for kk, vv in timings.items():
data[kk].append(vv)
assemble_data[name][0].append(I)
assemble_data[name][1].append(timings['assemble'])
solve_data[name][0].append(I)
solve_data[name][1].append(timings['solve'])
print()
plt.figure()
for id in ('grid', 'assemble', 'solve'):
plt.semilogy(data['I'], data[id], label=id)
plt.legend()
plt.xlabel('I')
plt.ylabel('time (s)')
def discretize_double_for(omega, I, f, g_D):
timings = {}
print('creating grid with {} elements ... '.format(I), end='')
t = time.time()
grid = OnedGrid(omega, I)
vertices = grid.centers(1) # call centers for time measurements, is created on demand
timings['grid'] = time.time() - t
print('done (took {}s)'.format(timings['grid']))
print('assembling ... ', end='')
t = time.time()
h = (omega[1] - omega[0]) / I
b_h = np.zeros((I + 1, I + 1))
f_h = np.zeros(I + 1)
for i in range(0, I + 1):
for j in range(0, I + 1):
if i == 0 and j == 0:
b_h[i][j] = 1.
elif i == I and j == I:
b_h[i][j] = 1.
elif j == i - 1:
b_h[i][j] = -1/h**2
elif j == i:
b_h[i][j] = 2/h**2
elif j == i + 1:
b_h[i][j] = -1/h**2
else:
b_h[i][j] = 0.
if i == 0:
f_h[i] = g_D[0]
elif i == I:
f_h[I] = g_D[1]
else:
f_h[i] = f.evaluate(vertices[i])
timings['assemble'] = time.time() - t
print('done (took {}s)'.format(timings['assemble']))
print('solving ... ', end='')
t = time.time()
u_h = np.linalg.solve(b_h, f_h)
timings['solve'] = time.time() - t
print('done (took {}s)'.format(timings['solve']))
return timings
measure_timings(discretize_double_for, 'double_for')
def discretize_single_for(omega, I, f, g_D):
timings = {}
print('creating grid with {} elements ... '.format(I), end='')
t = time.time()
grid = OnedGrid(omega, I)
vertices = grid.centers(1) # call centers for time measurements, is created on demand
timings['grid'] = time.time() - t
print('done (took {}s)'.format(timings['grid']))
print('assembling ... ', end='')
t = time.time()
h = (omega[1] - omega[0]) / I
b_h = np.zeros((I + 1, I + 1))
f_h = np.zeros(I + 1)
for i in range(1, I):
b_h[i][i - 1] = -1/h**2
b_h[i][i] = 2/h**2
b_h[i][i + 1] = -1/h**2
f_h[i] = f.evaluate(vertices[i])
b_h[0][0] = 1
f_h[0] = g_D[0]
b_h[I][I] = 1
f_h[I] = g_D[1]
timings['assemble'] = time.time() - t
print('done (took {}s)'.format(timings['assemble']))
print('solving ... ', end='')
t = time.time()
u_h = np.linalg.solve(b_h, f_h)
timings['solve'] = time.time() - t
print('done (took {}s)'.format(timings['solve']))
return timings
measure_timings(discretize_single_for, 'single_for')
def discretize_vectorized_dense(omega, I, f, g_D):
timings = {}
print('creating grid with {} elements ... '.format(I), end='')
t = time.time()
grid = OnedGrid(omega, I)
vertices = grid.centers(1) # call centers for time measurements, is created on demand
timings['grid'] = time.time() - t
print('done (took {}s)'.format(timings['grid']))
print('assembling ... ', end='')
t = time.time()
h = (omega[1] - omega[0]) / I
b_h = (
np.diag(np.ones(I)*(-1/h**2), k=-1)
+ np.diag(np.ones(I + 1)*(2/h**2))
+ np.diag(np.ones(I)*(-1/h**2), k=1)
)
b_h[0][0] = 1
b_h[0][1] = 0
b_h[I - 1][I] = 0
b_h[I][I] = 1
f_h = interpolate(grid, f)
f_h[0] = g_D[0]
f_h[I] = g_D[1]
timings['assemble'] = time.time() - t
print('done (took {}s)'.format(timings['assemble']))
print('solving ... ', end='')
t = time.time()
u_h = np.linalg.solve(b_h, f_h)
timings['solve'] = time.time() - t
print('done (took {}s)'.format(timings['solve']))
return timings
measure_timings(discretize_vectorized_dense, 'vectorized_dense')
import scipy as sp
def discretize_vectorized_sparse(omega, I, f, g_D):
timings = {}
print('creating grid with {} elements ... '.format(I), end='')
t = time.time()
grid = OnedGrid(omega, I)
vertices = grid.centers(1) # call centers for time measurements, is created on demand
timings['grid'] = time.time() - t
print('done (took {}s)'.format(timings['grid']))
print('assembling ... ', end='')
t = time.time()
h = (omega[1] - omega[0]) / I
lower_diag = np.ones(I)*(-1/h**2)
main_diag = np.ones(I + 1)*(2/h**2)
upper_diag = np.ones(I)*(-1/h**2)
main_diag[0] = 1
upper_diag[0] = 0
main_diag[I] = 1
lower_diag[I - 1]
b_h = sp.sparse.diags(
[lower_diag, main_diag, upper_diag],
[-1, 0, 1],
shape=(I + 1, I + 1))
f_h = interpolate(grid, f)
f_h[0] = g_D[0]
f_h[I] = g_D[1]
timings['assemble'] = time.time() - t
print('done (took {}s)'.format(timings['assemble']))
print('solving ... ', end='')
t = time.time()
u_h = sp.sparse.linalg.spsolve(b_h, f_h)
timings['solve'] = time.time() - t
print('done (took {}s)'.format(timings['solve']))
return timings
measure_timings(discretize_vectorized_sparse, 'vectorized_sparse')
plt.figure()
for name, data in assemble_data.items():
plt.semilogy(data[0], data[1], label=name)
plt.legend()
plt.title('assemble')
plt.xlabel('I')
plt.ylabel('time (s)')
plt.figure()
for name, data in solve_data.items():
plt.semilogy(data[0], data[1], label=name)
plt.legend()
plt.title('solve')
plt.xlabel('I')
plt.ylabel('time (s)')
