\documentclass[12pt]{article}

\usepackage[T1]{fontenc} 
\usepackage[ngerman]{babel}
\usepackage{fullpage,amsmath,amsfonts}

\begin{document}

\title{Mathematisches Programmieren}
\author{Sommersemester 2022}
\date{25.05.2022, 01.06.2022, 08.06.2022 und 15.06.2022}
\maketitle

\subsubsection*{Die Poisson-Gleichung}

Sei $\Omega \subset \mathbb R^n$ offen, beschränkt und einfach
zusammenhängend.  Wir betrachten die \emph{Poisson-Gleichung}
\begin{equation}
  - \Delta u = f \quad \text{in } \Omega
\end{equation}
mit homogenen \emph{Dirichlet-Randbedingungen}
\begin{equation}
  u = 0 \quad \text{auf } \partial \Omega \,.
\end{equation}
Hierbei bezeichnet $\Delta$ den \emph{Laplace-Operator}
\begin{equation}
  \Delta u = \frac{\partial^2 u}{\partial x_1^2} + \dots
  + \frac{\partial^2 u}{\partial x_n^2} \,.
\end{equation}
Diese Gleichung beschreibt z.B. die Abhängigkeit des elektrischen
Potentials $u$ in Abhängigkeit von einer gegebenen Ladungsverteilung
$f$.  Die homogene Dirichlet-Randbedingung drückt aus, dass das Gebiet
von einem elektrischen Leiter umschlossen ist, dessen Potential auf
Null liegt.  Die Poisson-Gleichung taucht aber auch in vielfältigen
anderen Zusammenhängen in Mathematik und Physik auf.  Sie ist der
Prototyp einer \emph{elliptischen partiellen Differentialgleichung}.

\subsubsection*{Finite Differenzen im Eindimensionalen ($n=1$)}

Hier ist das Gebiet $\Omega$ einfach ein Intervall.  Mit
entsprechender Skalierung können wir $\Omega = (0, \pi)$ wählen.  Dann
lässt sich das homogene Dirichletproblem schreiben als
\begin{subequations}
  \label{e.poisson}
\begin{gather}
  - u''(x) = f(x) \quad \text{in } (0, \pi) \,, \\
  u(0) = u(\pi) = 0 \,.
\end{gather}
\end{subequations}
Um diese Gleichung zu diskretisieren, unterteilen wir das Intervall $(0,
L)$ in $N$ Teilintervalle der Länge $h = L/N$, wobei hier $L = \pi$,
und definieren die Stützstellen
\begin{gather}
  x_i = i \, h \,.
\end{gather}
Sei nun $u_i \approx u(x_i)$, also eine Näherung zur Lösung $u(x)$ an
der Stützstelle $x_i$, die folgendes lineares Gleichungssystem
erfüllt:
\begin{subequations}
  \label{e.discrete1d}
\begin{gather}
  - \frac{u_{i-1} - 2 \, u_i + u_{i+1}}{h^2} = f(x_i)
  \quad \text{für } i = 1, \dots, N-1 \,, \\
  u_0 = u_N = 0 \,.
\end{gather}
\end{subequations}
Man kann zeigen, dass, falls die exakte Lösung $u$ viermal stetig
differenzierbar ist, die linke Seite dieser Näherung gegen $-u''(x)$
konvergiert, und zwar quadratisch, im Grenzfall $h \to 0$.

\subsubsection*{Finite Differenzen für $n=2$}

In zwei Raumdimensionen kann die Näherung \eqref{e.discrete1d} für
jede der partiellen Ableitungen verwendet werden.  Diese Näherung
funktioniert am Besten, wenn die Ränder des Gebiets achsenparallele
Geradensegmente sind.  Wir beschränken uns daher auf den Fall eines
Rechtecks $\Omega = (0,L_x) \times (0,L_y)$.  In jeder Raumrichtung
nehmen wir $N_x$ bzw.\ $N_y$ Teilintervalle, so dass die jeweiligen
Schrittweiten $h_x = L_x/N_x$ bzw. $h_y = L_y/N_y$ sind.  Die
Stützstellen sind dann
\begin{gather}
  x_{ij} = (i \, h_x, j \, h_y)
\end{gather}
und die Näherung $u_{ij} \approx u(x_{ij})$ löst
\begin{subequations}
  \label{e.discrete2d}
\begin{gather}
  - \frac{u_{i-1,j} - 2 \, u_{ij} + u_{i+1,j}}{h_x^2}
  - \frac{u_{i,j-1} - 2 \, u_{ij} + u_{i,j+1}}{h_y^2}
  = f(x_{ij})
\intertext{für $i = 1, \dots, N_x-1$, $j = 1, \dots, N_y-1$, mit
Randbedingungen} 
  u_{i,0} = u_{i,N_y} = 0
  \quad \text{für } i = 0, \dots, N_x \,, \\
  u_{0,j} = u_{N_x,j} = 0
  \quad \text{für } j = 0, \dots, N_y \,.
\end{gather}
\end{subequations}

\subsubsection*{Aufgaben}



\begin{enumerate}
\item Schreibe das lineare Gleichungssystem \eqref{e.discrete1d} in
Matrixform.  Erzeuge, für gegebenes $N$, diese Matrix in Python
(Hinweis: verwende \texttt{diag}), berechne ihre Eigenwerte (Hinweis:
\texttt{eig} oder \texttt{eigvals} -- weil die Matrix symmetrisch ist,
kann man auch \texttt{eigh} oder \texttt{eigvalsh} verwenden) und
trage die Eigenwerte als Funktion ihres Index in ein Koordinatensystem
zusammen mit den ersten $N-1$ exakten Eigenwerten des
Laplaceoperators, $k^2$ für $k=1, \dots, N-1$, auf.

\item Das lineare Gleichungssystem \eqref{e.discrete1d} hat eine
Tridiagonalstruktur.  Weil das Gauß'sche Eliminationsverfahren diese
Tridiagonalstruktur erhält, kann man das Gleichungssystem leicht in
$O(N)$ Recheneoperationen lösen.  (Das Gauß'sche Verfahren benötigt
$O(N^3)$ Rechenoperationen im allgemeinen Fall und ist damit für sehr
große Matrizen nicht geeignet.)  Dieses vereinfachte Verfahren wird
oft als "`Thomas-Algorithmus"' bezeichnet.

Implementiere den Thomas-Algorithmus in Python (Hinweis: Folge der
Notation im Wikipedia-Artikel über den Thomas-Algorithmus).
Überprüfe die Richtigkeit des Lösers im Vergleich mit dem 
\textsf{Scipy}-Löser
\verb+scipy.linalg.solve_banded+ und vergleiche die benötigte
Rechenzeit für eine sehr große Tridiagonalmatrix.

\item Verwende einen der Löser aus Aufgabe~2, um die Poisson-Gleichung
\eqref{e.poisson} zu lösen.  Teste den Löser in folgenden Schritten:
\begin{enumerate}
\item Wähle eine Funktion $u(x)$ mit $u(0)=u(\pi)=0$, die sich
explizit differenzieren lässt und setze diese in die Poisson-Gleichung
\eqref{e.poisson} ein, um $f(x)$ zu berechnen.
\item Löse das diskrete System \eqref{e.discrete1d} zu dem jetzt
gegebenen $f(x)$ für verschiedene Werte von $N$.  Trage die exakte
Lösung $u(x)$ sowie die berechneten Näherungslösungen $u_i$ in ein
Koordinatensystem auf.
\item Berechne den Diskretisierungsfehler
\[
  E(N) = \max_{i=1, \dots, N-1} \lvert u_i - u(x_i) \rvert
\]
und trage $E(N)$ gegen $N$ in ein doppelt-logarithmisches
Koordinatensystem auf.  Was stellen Sie fest?
\end{enumerate}

\item Bei der Poisson-Gleichung in zwei Raumdimensionen bekommt man
eine Matrix mit $4$ Nebendiagonalen.  Allerdings sind zwei dieser
Nebendiagonalen $O(N)$ weit von der Hauptdiagonalen entfernt.  Damit
funktioniert der Thomas-Algorithmus nicht mehr: Ein normaler
Gauß'sches Eliminationsprozess würde im Allgemeinen alle Einträge, die
zwischen diesen Nebendiagonalen liegen, mit von Null verschiedenen
Werten füllen.  Glücklicherweise gibt es in \textsf{Scipy} effiziente
Löser für lineare Gleichungssysteme mit dünnbesetzten Matrizen.

Im Folgenden ist die Poisson-Gleichung auf dem Rechteck $\Omega = (0,1)
\times (0,1)$ zu lösen.  Der Einfachheit halber kann $N_x=N_y \equiv
N$ gewählt werden.

\begin{enumerate}
\item Erzeuge die Matrix auf der linken Seite der Gleichung
\eqref{e.discrete2d} im CSR-Format (Hinweis:
\verb+scipy.sparse.csr_matrix+).

\item Definiere ein geeignetes Testproblem, wie im Eindimensionalen,
und löse dieses numerisch (Hinweis:
\texttt{scipy.sparse.linalg.spsolve}).

\item Wie skaliert die Rechenzeit von \texttt{spsolve} mit der Größe
des Gleichungssystems, das $(N_x-1)*(N_y-1)$ Variablen enthält?

\end{enumerate}


\end{enumerate}

Geben Sie ihren lauffähigen Code per Email in einem ZIP-Archiv
bis zum 24.06.2022 ab.


\end{document}
