Wichtig: Damit alle benötigten Pakete richtig eingebunden werden, führen Sie die nächste Zelle einmal aus, sobald Sie das Notebook neu öffnen.
# Bibliotheken einbinden
%matplotlib inline
import numpy as np
import numpy.linalg as linalg
import matplotlib.pyplot as plt
Gegeben sei für $u \in C^2([0,1])$ das Randwertproblem mit homogenen Randdaten \begin{equation} -u''(x) = f(x)\, , \quad x \in (0,1)\, , \quad u(0) = u(1) = 0\, , \end{equation} wobei $f \in C([0,1])$ gelte. Zur näherungsweisen Lösung sei für einen Parameter $n \in N$ das Gitter $x_i = ih$, $i = 0, \ldots, n$ mit $h = 1/n$ gegeben. Verwendet man zur Diskretisierung Finite Differenzen, d.h. $-u''(x_i) \approx (2u(x_i) - u(x_{i-1}) - u(x_{i+1}))/h^2$, so erhält man das Gleichungssystem \begin{equation} \frac{1}{h^2} \begin{pmatrix} 2 & -1 & & & & \\ -1 & 2 & -1 & & & \\ & \ddots & \ddots & \ddots & & \\ & & -1 & 2 & -1 & \\ & & & -1 & 2 \end{pmatrix} \cdot \begin{pmatrix} u_1 \\ u_2 \\ \vdots \\ u_{n-2} \\ u_{n-1} \end{pmatrix} = \begin{pmatrix} f_1 \\ f_2 \\ \vdots \\ f_{n-2} \\ f_{n-1} \end{pmatrix}\, , \qquad \qquad \qquad (1) \end{equation} wobei $f_i = f(x_i)$ gesetzt wird. Dabei ist jede Komponente $u_i$, $i = 1, \ldots, n-1$, der Lösung des Gleichungssystems eine Näherung der exakten Lösung des kontinuierlichen Problems and der Stelle $x_i$, d.h. $u_i \approx u(x_i)$.
Schreiben Sie jeweils Programme zur iterativen Lösung des Gleichungssystems (1) durch das Gauß-Seidel-, das Jacobi- und das Richardson-Verfahren. Abkürzend werde das System (1) als $Ay = b$ geschrieben. Der Aufruf des Programms soll mit \begin{equation} \texttt{y = loeser}(\texttt{b, y0, tol, kmax}); \end{equation} erfolgen, wobei $\texttt{loeser}$ entweder $\texttt{jacobi}$, $\texttt{gaussseidel}$ oder $\texttt{richardson}$ ist. Eingabeparameter sind \begin{align} \texttt{b} \qquad &\text{Vektor der rechten Seite} \\ \texttt{y0} \qquad &\text{Startvektor } y^{(0)} \text{ der Iteration} \\ \texttt{tol} \qquad &\text{Fehlertoleranz} \\ \texttt{kmax} \qquad &\text{maximale Anzahl der Iterationen} \end{align} und zusätzlich $\texttt{omega}$ für die Schrittweite $\omega > 0$ im Richardson-Verfahren. Ausgabeparameter ist $\texttt{y}$, also die Approximation der Lösung $A^{-1}b$.
Dabei soll die Iteration abgebrochen werden, falls die maximale Iterationszahl $\texttt{kmax}$ überschritten wurde oder \begin{equation} | Ay^{(k)} - b|_2 + | y^{(k-1)} - y^{(k)}|_2 < \texttt{tol} \end{equation} gilt.
Verwenden Sie zum Testen z.B. die konstante Funktion $f \equiv 1$. Plotten Sie die Lösung $(u_i)_i$ als Funktion über dem Interval $[0,1]$ und beobachten Sie die Konvergenz für grösser werdendes $n$. Die exakte Lösung $u$ erhält man durch zweimaliges Aufintegrieren von $-f$.
Testen Sie für das Richardson-Verfahren verschiedene Parameter $\omega > 0$, z.B. $\omega = 1$. Laut Vorlesung sollte $\omega < 2/\lambda_\max(A)$ sein, wobei $\lambda_\max$ der grösste Eigenwert von $A$ ist. Der optimale Wert ist laut Vorlesung gegeben durch $$\omega = \frac{2}{\lambda_\max(A) + \lambda_\min(A)}\, .$$ Die Eigenwerte $0 <\lambda_\min(A) \leq \ldots \leq \lambda_\max(A)$ von $A$ wurden auf dem Übungsblatt 11 bestimmt.
Berechnen Sie für festes $n$ die exakte Lösung $y = A^{-1}b$, z.B. mit Hilfe der $\texttt{linalg}$-Bibliothek und dem Befehl $\texttt{linalg.solve(A,b)}$, und plotten Sie den Fehler $e(k) = \| y - y^{(k)} \|_2$ für $k =1,2,\ldots$ in einer semilog-Skala. Interpretieren Sie die Ergebnisse!
$\textit{Hinweis}$: Assemblieren Sie nicht explizit die Matrix $A$, sondern nur die Matrix-Vektor-Multiplikation $Ay$. Nutzen Sie die dünne Bandstruktur aus!
# Fuegen Sie hier Ihre Loesung ein!