Institut für Numerische und Angewandte Mathematik
FB Mathematik und Informatik der Universität Münster
Prof. Dr. Mario Ohlberger, Dr. Stephan Rave und Tobias Wedemeier
%matplotlib notebook
Aufgabe 1 a)¶
Zunächst definieren wir erneut die Diskretisierungen von Blatt 1, Aufgabe 2:
from pymor.basic import *
set_log_levels({'pymor': 'WARN'})
# Aufgabe 2 a)
domain = RectDomain(bottom='neumann')
neumann_data = ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type= {'neum': ()})
diffusion = ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
2, (),
values={'K': 10}
)
problem = StationaryProblem(
domain=domain,
diffusion=diffusion,
neumann_data=neumann_data
)
d_a, _ = discretize_stationary_cg(problem, diameter=1/100)
# Aufgabe 2 b)
diffusion = ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu)',
2, (),
values={'K': 10},
parameter_type= {'diffu': ()}
)
problem = StationaryProblem(
domain=domain,
diffusion=diffusion,
neumann_data=neumann_data
)
d_b, _ = discretize_stationary_cg(problem, diameter=1/100)
# Aufgabe 2 c)
diffusion_R = BitmapFunction('R.png', range=[1, 0])
diffusion_B = BitmapFunction('B.png', range=[1, 0])
diffusion = LincombFunction(
[ConstantFunction(1., 2), diffusion_R, diffusion_B],
[1., ExpressionParameterFunctional('-(1 - R)', {'R': ()}), ExpressionParameterFunctional('-(1 - B)', {'B': ()})]
)
problem = StationaryProblem(
domain=domain,
diffusion=diffusion,
neumann_data=ConstantFunction(-1, 2)
)
d_c, _ = discretize_stationary_cg(problem, diameter=1/100)
Parameterräume geben in pyMOR den zulässigen Wertebereich für die Parameter an, von denen ein gegebenes parametrisches Objekt, hier die Diskretisierung, abhängt. Wenn wir also z.B. $\mathcal{P} \subset \mathbb{R}^p$ schreiben, so würde der parameter_type des Objekts $\mathbb{R}^p$ entsprechen, wohingegen parameter_space die Menge $\mathcal{P}$ definiert.
In pyMOR ist gegenwärtig genau ein Parameterraum implementiert, CubicParameterSpace, der Parameterräume $\mathcal{P}$ beschreibt, die kartesische Produkte von abgeschlossenen Intervallen sind. Beim Anlegen muss der zugehörige parameter_type sowie die gemeinsamen unteren und oberen Schranken der Teilintervalle angegeben werden. Alternativ
kann mit dem ranges-Parameter ein Dictionary von Intervallen für die einzelnen Parameter-Komponenten angegeben werden:
d_a = d_a.with_(parameter_space=CubicParameterSpace(d_a.parameter_type, -10, 10))
d_b = d_b.with_(parameter_space=CubicParameterSpace(d_b.parameter_type,
ranges={'diffu': [0.0001, 1.], 'neum': [-10, 10]}))
d_c = d_c.with_(parameter_space=CubicParameterSpace(d_c.parameter_type, 0.0001, 1.))
Hier wollen wir das parameter_space-Attribut der Diskretisierungen d_a, d_b, d_c auf einen geeigneten Parameterraum setzten. Da jedoch Diskretisierungen wie auch Operatoren in pyMOR unveränderlich sind ('immutable'), können wir nur ein neues Diskretisierungs-Objekt erzeugen, in dem das parameter_space-Attribut passend gesetzt ist. Um dies zu erleichtern, haben solche unveränderliche Objekte in pyMOR eine with_-Methode. Der Aufruf
obj.with_(attrib1=value1, attrib2=value2)gibt dabei eine Kopie von obj zurück, in dem die Attribute attrib1, attrib2 auf value1 bzw. value2 gesetzt sind. In der Anweisung
d_a = d_a.with_(...)wird dabei die modifizierte Kopie von d_a wieder dem Namen d_a zugewiesen und das alte Diskretisierungs-Objekt verworfen.
Die Unveränderlichkeit von Diskretisierungen und Operatoren in pyMOR ist eine Zusicherung an den Progammierer, dass diese Objekte, z.B. nachdem sie einer Funktion als Argument übergeben worden sind, ihr Verhalten nicht (überraschend) verändern können. - Der angehängte Unterstrich in with_ ist vonnöten, da with in Python ein reserviertes Wort ist und nicht an eine Methode vergeben werden kann. Das Anhängen eines Unterstrichs ist eine Standard-Konvention in Python.
Im Folgenden legen wir uns auf eine der drei oben definierten Diskretisierungen fest:
d = d_c
Aufgabe 1 b)¶
Ein wesentlicher Nutzen von Parameterräumen in pyMOR ist, dass mit ihren sample_...-Methoden leicht Mengen von in ihnen liegenden Parametern generiert werden können. CubicParameterSpace verfügt über die Methoden sample_uniformly und sample_randomly. Wir verwenden hier die letztere Methode:
U = d.solution_space.empty()
for mu in d.parameter_space.sample_randomly(30):
U.append(d.solve(mu))
Alle Lösungssnapshots werden hier in dem VectorArray U aneinander gehängt. Neue VectorArrays werden in pyMOR über ihren zugehörigen VectorSpace erzeugt. Alle Lösungen einer Diskretisierung d gehören zu d.solution_space, dessen empty-Methode ein neues leeres Array erzeugt. (Alle Vektoren in einem VectorArray gehören in pyMOR zu einem festen VectorSpace.) Mit der append-Methode des Arrays U werden die in dem von d.solve zurückgegebenen Array enthaltenen Vektoren (in diesem Fall genau ein Vektor) and U angehängt.
Die Länge eines VectorArrays ist die Anzahl der in ihm enthaltenen Vektoren:
len(U)
Ein mehrere Vektoren enthaltendes VectorArray kann als Zeitreihe visualisiert werden:
d.visualize(U)
Schließlich geben wir noch die $H^1$- und $H^1_0$-Normen der Vektoren in U aus:
print('H^1-Norm:', d.h1_0_norm(U))
print('H^1_0-Norm:', d.h1_0_semi_norm(U))
Achtung: d besitzt außerdem die Normen d.h1_norm und d.h1_semi_norm, die sich von obigen nur um die Randwertbehandlung der zugrunde liegenden Systemmatrizen unterscheiden. Bei der Berechnung von Normen von $H^1_0$-Funktionen liefern beide Varianten das selbe Ergebnis. Für die Berechnung von Riesz-Repräsentanten muss jedoch unbedingt das für den jeweiligen Raum ($H^1$ bzw. $H^1_0$) korrekt assemblierte Skalarprodukt verwendet werden.
Aufgabe 2 a)¶
Sei $v_{proj}$ die orthogonale Projektion von $v$ auf den von den Vektoren in $U$ aufgespannten Vektorraum. Dann ist $v - v_{proj}$ orthogonal zu diesem Raum und insbesondere zu jedem Vektor aus $U$. Sind $u_i$, $i = 1, \ldots, N$ die Vektoren in $U$, so gilt also:
$$(v - v_{proj}, u_i) = 0, \qquad i = 1, \ldots, N.$$Sind nun $\lambda_j$, $j = 1, \ldots, N$ die Koeffizienten von $v_{proj}$ bezüglich der $u_j$ basis, d.h. $\sum_{j=1}^N \lambda_j u_j = v_{proj}$, so folgt:
$$ \sum_{j=1}^N \lambda_j (v_i, v_j) = (v, u_i), \qquad i = 1, \ldots, N.$$Setzten wir $M_{i,j} := (v_i, v_j)$, $R := [(v, u_1), \ldots, (v, u_N)]^T$ und $\Lambda := [\lambda_1, \ldots, \lambda_N]^T$, so erhalten wir das lineare Gleichungssystem
$$ M \cdot \Lambda = R,$$mit dem wir die Koeffizienten $\lambda_i$ und somit $v_{proj}$ bestimmen können.
Aufgabe 2 b)¶
Jedes VectorArray besitzt eine gramian-Methode, welche gerade die oben definierte Gram-Matrix $M$ bezüglich des Euklidischen Skalarprodukts berechnet. gramian ist ein Spezialfall der allgemeineren dot-Methode, die im Falle von U.dot(V) die Matrix $(u_i, v_j)$ aller Euklidischen Skalarprodukte zwischen den Vektoren in U und V berechnet.
Die Linearkombination $\sum_{j=1}^N \lambda_j u_j$ kann in pyMOR effizient mit der lincomb-Methode von U bestimmt werden. Diese erwartet ein 1D- oder 2D-Array von Koeffizienten $\lambda_j$. (Im Falle eines 2D-Arrays entspricht dabei jede Zeile einer Folge von Linearkoeffizienten und es werden die entsprechenden Linearkombinationen für jede Zeile des Arrays zurückgeliefert.)
Damit können wir $v_{proj}$ wie folgt berechnen:
V = d.solve(d.parameter_space.sample_randomly(1, seed=99)[0])
M = U.gramian()
rhs = U.dot(V) # beachte die Reihenfolge von U und V um einen Spaltenvektor zu erhalten!
v = np.linalg.solve(M, rhs)
V_proj = U.lincomb(v.T) # hier wird ein Zeilenvektor von Koeffizienten erwartet!
Wie Diskretisierungen und Operatoren, sind auch Parameterräume in pyMOR unveränderlich. Daher liefert auch die sample_randomly-Methode immer wieder die gleichen Parameter. Um eine neue Zufallsfolge von Parametern zu erhalten, kann ein neuer seed für den Zufallsgenerator übergeben werden.
Wird visualize ein Tupel von VectorArrays übergeben, wird ein Plot für jedes Array generiert. Der separate_colorbars Parameter legt fest, ob alle Plots mit der gleichen Colorbar skaliert werden sollen:
d.visualize((V, V_proj, V - V_proj), separate_colorbars=True)
Aufgabe 2 c)¶
Um die Matrix $M$ mit dem korrekten Skalarprodukt zu berechnen verwenden wir die apply2-Methode des entsprechenden Skalarprodukt-Operators. Dabei ist op.apply2(U, V) gerade äquivalent zu U.dot(op.apply(V)). Ist also op durch eine Matrix P gegeben und fassen wir U, V als Matrizen von Spaltenvektoren auf, so erhalten wir gerade die Matrix
$ U^T \cdot P \cdot V$.
Achtung: In NumpyVectorArray sind die einzelnen Vektoren gerade zeilenweise in einem 2D-Array gespeichert. U.data entspricht also immer einer Matrix von Zeilenvektoren!
M = d.h1_0_semi_product.apply2(U, U)
rhs = d.h1_0_semi_product.apply2(U, V)
v = np.linalg.solve(M, rhs)
V_h1_proj = U.lincomb(v.T)
d.visualize((V, V_h1_proj, V - V_h1_proj), separate_colorbars=True)
Schließlich berechnen wir noch die Projektionsfehler:
for norm_name, norm in [('eukld.', None), ('H^1_0', d.h1_0_semi_norm), ('L^2', d.l2_0_norm)]:
print('{:<16}'.format(norm_name + '-Norm:'), end='')
for proj_name, VV in [('eukld.', V_proj), ('H^1_0', V_h1_proj)]:
if norm is None:
err = ((V - VV).l2_norm() / V.l2_norm())[0]
else:
err = (norm(V - VV) / norm(V))[0]
print('{}: {:.2e} '.format(proj_name, err), end='')
print()
Aufgabe 2 d)¶
Zunächst berechnen wir die 100 Test-Snapshots:
V = d.solution_space.empty()
for mu in d.parameter_space.sample_randomly(100, seed=99):
V.append(d.solve(mu))
print('.', end='', flush=True)
Die Berechnung der Projektions-Matrix und rechten Seite für verschiedene Basisgrößen können wir optimieren, indem wir diese zunächst für die volle Basis berechnen und dann entsprechend ausschneiden. Wie bei Numpy-Arrays können wir auch bei VectorArrays slicing verwenden. Dabei erhalten wir ebenfalls einen View in das Original-Array. Es wird kein Kopie der Daten erzeugt.
M = d.h1_0_semi_product.apply2(U, U)
rhs = d.h1_0_semi_product.apply2(U, V)
errors = []
for N in range(len(U)):
try:
v = np.linalg.solve(M[:N, :N], rhs[:N, :])
except np.linalg.LinAlgError:
v = np.zeros((N, len(V)))
V_proj = U[:N].lincomb(v.T)
errors.append(np.max(d.h1_0_semi_norm(V - V_proj)))
Je nach Problem kann es passieren, dass die Vektoren in U (numerisch) linear abhängig werden. In diesem Fall schlägt der Aufruf von np.linalg.solve mit einer np.linalg.LinalgError-Exception fehl. Dies fangen wir in dem try-except-Block ab und setzen in diesem Fall V_proj auf 0.
In der Modellreduktion haben wir es oft mit exponentiellem Fehlerabfall zu tun. Dafür eignet sich ein semi-logarithmischer Plot, den wir mit matplotlib.pyplot.semilogx bzw. matplotlib.pyplot.semilogy erhalten:
from matplotlib import pyplot as plt
plt.figure()
plt.semilogy(errors)
plt.show()
Aufgabe 2 e)¶
Die Konditionen berechnen wir analog mit numpy.linalg.cond. Dabei ist zu beachten, dass wir keine $0 \times 0$-Matrix übergeben dürfen:
conds = []
for N in range(1, len(U)):
conds.append(np.linalg.cond(M[:N, :N]))
plt.figure()
plt.semilogy(range(1, len(U)), conds)
plt.show()
Aufgabe 3 a)¶
Wollen wir mit Hilfe eines Skalarprodukt-Operators Normen berechnen, können wir statt apply2 die pairwise_apply2-Methode verwenden, die gerade die Diagonale der von apply2 berechneten Matrix zurückliefert:
def greedy(M, product, N):
basis = M.space.empty()
max_errors = []
for n in range(N):
lhs = product.apply2(basis, basis)
rhs = product.apply2(basis, M)
try:
m = np.linalg.solve(lhs, rhs)
except np.linalg.LinAlgError:
break
M_proj = basis.lincomb(m.T)
M_err = M - M_proj
errors = np.sqrt(product.pairwise_apply2(M_err, M_err))
max_errors.append(np.max(errors))
basis.append(M[np.argmax(errors)])
return basis
Aufgabe 3 b)¶
Zunächst berechnen wir eine reduzierte Basis mit Hilfe des implementierten Greedy-Algorithmus:
greedy_basis = greedy(U, d.h1_0_product, 30)
Anschließend berechnen wir wie bisher die Approximationsfehler:
M_greedy = d.h1_0_semi_product.apply2(greedy_basis, greedy_basis)
rhs = d.h1_0_semi_product.apply2(greedy_basis, V)
greedy_proj_errors = []
for N in range(len(greedy_basis)):
try:
v = np.linalg.solve(M_greedy[:N, :N], rhs[:N, :])
except np.linalg.LinAlgError:
v = np.zeros((N, len(V)))
V_proj = greedy_basis[:N].lincomb(v.T)
greedy_proj_errors.append(np.max(d.h1_0_semi_norm(V - V_proj)))
plt.figure()
plt.semilogy(errors, label='without greedy')
plt.semilogy(greedy_proj_errors, label='with greedy')
plt.legend()
plt.show()
Schließlich noch die Konditionen der Projektions-Matrizen:
greedy_conds = []
for N in range(1, len(U)):
greedy_conds.append(np.linalg.cond(M_greedy[:N, :N]))
plt.figure()
plt.semilogy(range(1, len(U)), conds, label='without greedy')
plt.semilogy(range(1, len(U)), greedy_conds, label='with greedy')
plt.legend()
plt.show()
