commit f548f0cca8d28d3cd03067fd168425e74e1b34a4
parent 17517f503a61c4cddf911958f88baa83c9ed35b3
Author: miksa <milutin@popovic.xyz>
Date: Sat, 12 Mar 2022 14:51:44 +0100
sheet 2 done
Diffstat:
8 files changed, 470 insertions(+), 1 deletion(-)
diff --git a/num_ana/build/prb1.pdf b/num_ana/build/prb1.pdf
Binary files differ.
diff --git a/num_ana/build/prb2.pdf b/num_ana/build/prb2.pdf
Binary files differ.
diff --git a/num_ana/prb2.tex b/num_ana/prb2.tex
@@ -0,0 +1,299 @@
+\include{preamble.tex}
+
+\usepackage[final]{pdfpages}
+
+\begin{document}
+\maketitle
+\tableofcontents
+\section{Sheet 2}
+\subsection{Problem 1}
+\subsubsection{}
+We let $\rho(A)$ be the spectral radius of a matrix $A \in
+\mathbb{R}^{n\times n}$. A matrix norm $ \|\cdot\|_{1}$ is consistent with the
+vector norm $\|\cdot\|_{2}$ if
+\begin{align}
+ \|Ax\|_{2} \le \|A\|_1\|x\|_2
+\end{align}
+for all $x \in \mathbb{R}^n$ and all $A \in \mathbb{R}^{n\times n}$. Indeed
+every matrix norm induced by a vector norm is consistent. To show this let
+$\|\cdot\|_M$ be a matrix norm and $\|\cdot\|_v$ be a vector norm, defined as
+\begin{align}
+ \|x\|_v = \|xv^T\|_M
+\end{align}
+for all $x \in \mathbb{R}^n$ and some $v \neq 0$, of $\mathbb{R}^{n \times
+n}$ and $\mathbb{R}^{n}$ respectively. Then we have
+\begin{align}
+ \|Ax\|_v = \|Axv^T\| \le \|A\|_M \|xv^T\|_M = \|A\|_M \|v\|_v \qed
+\end{align}
+Note that for $\mathbb{C}^{n \times n}$ and $\mathbb{C}^{n}$ use $v^* \neq 0$ the
+conjugate transpose.
+\subsubsection{}
+Now we consider a splitting of $A = D - (L + U)$, were $D, L$ and $U$ are
+defined as
+\begin{align}
+ D &= \text{diag}\left(a_{11},\ldots,a_{nn}\right), \\
+ \left(L \right)_{ij} &=
+ \begin{cases}
+ -\left( A \right)_{ij}\quad i > j\\
+ 0 \quad i\leq 0
+ \end{cases},
+ \qquad
+ \left(U \right)_{ij} =
+ \begin{cases}
+ -\left( A \right)_{ij}\quad i < j\\
+ 0 \quad i \ge 0
+ \end{cases},
+\end{align}
+Then the matrix of a single Jacobi iteration method is
+\begin{align}
+ B_J = D^{-1}(L+U)
+\end{align}
+We can show that if $A$ is strictly diagonally dominant then
+\begin{align}
+ \rho(B_J) \le \|B_J\|_\infty < 1.
+\end{align}
+If A is strictly diagonally dominant, this means that
+\begin{align}
+ &\mid A_{ii}\mid > \sum_{j\neq i}^{n} \mid A_{ij} \mid \qquad \forall\
+ i\in\{1, \ldots , n\} \\
+ \Leftrightarrow &\sum_{j\neq i} \frac{ \mid A_{ij} \mid}{ \mid A_{ii} \mid}
+ < 1.
+\end{align}
+Now let $(\lambda, v)$ be an eigen-pair of $B_J$, then
+\begin{align}
+ B_J v &= D^{-1}(L+U)v = \lambda v\\
+ (L+U)v &= \lambda D v.
+\end{align}
+For a chosen $i$ this means
+\begin{align}
+ \left| \lambda \right| \left| A_{ii} \right| \left| v_i \right|
+ &=
+ \left|
+ -\sum_{j>i} A_{ij}v_i - \sum_{j<i} A_{ij}v_i
+ \right|\\
+ &\le \sum_{j>i}
+ \left|A_{ij} \right| \left| v_i \right| + \sum_{j<i}
+ \left|A_{ij} \right| \left| v_i \right|\\
+ &= \sum_{j\neq i} \left| A_{ij} \right| \left|v_i\right|.
+\end{align}
+We can choose and $i$ such that $\left| v_i \right| \le \|v\|_\infty$, then
+\begin{align}
+ \left| \lambda \right| \left| A_{ii} \right| \left| v_i \right|
+ &\le \sum_{j\neq i} \left| A_{ij} \right| \|v\|_\infty\\
+ \Rightarrow
+ \left| \lambda \right| \left| A_{ii} \right|
+ &\le \sum_{j\neq i} \left| A_{ij} \right|\\
+ \Rightarrow
+ \left| \lambda \right|
+ &\le \sum_{j\neq i}\frac{\left| A_{ij} \right|}{\left| A_{ii} \right|} <
+ 1. \qed
+\end{align}
+\subsubsection{}
+Next we show that the Jacobi method converges for every initial guess
+$x^0$ to the solution of the equation $Ax = b$, given that $A$ is
+strictly diagonally dominant. So with any initial guess $x^0$ at the $k$-th
+iteration we have
+\begin{align}
+ Dx^{(k)} &= (L+U)x^{(k-1)} + b\\
+ \Leftrightarrow x^{(k)} &= D^{-1}(L+U)x^{(k-1)}+D^{-1}B \\
+ &= B_J x^{(k-1} + D^{-1}b
+\end{align}
+Now let $x$ be the exact solution, then the error at the $k$-th iteration is
+\begin{align}
+ e^{(k)} &= x - x^{(k)} = B_J\left( x - x^{(k-1)} \right) = \ldots\\
+ &= B_J^k e^{(0)}
+\end{align}
+Assume that $e^{(0)} \neq 0$, then we need $\lim_{n \to \infty}B_J^k = 0$. If
+$A$ is \textbf{diagonally dominant} we have $\rho(B_J) < 1$ this means for an
+eigen-pair $\left(\lambda, v \right) $ of $B_J$ we have
+\begin{align}
+ \lim_{n \to \infty}B_J^k v = \lim_{n \to \infty}\lambda^k v \\
+ \Rightarrow \lim_{n \to \infty}\lambda^k = 0,
+\end{align}
+for all $\lambda$ because $\rho(B_J) < 1$. \qed
+\subsection{Problem 2}
+Now consider a $A \in \mathbb{R}^{n \times n}$,
+\begin{align}
+ A = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix}.
+\end{align}
+Let $A = D - (L + U)$ like in the above problem.
+The Gauss-Siedel method has the iteration matrix $B_G = (D - L)^{-1} U$ and
+the Jacobi method has the iteration matrix $B_J = D^{-1}(L + U)$.
+\subsubsection{}
+We show that the spectral radii of $B_J$ and $B_G$ satisfy
+\begin{align}
+ \rho (B_J) = \sqrt{\left| \rho(B_G) \right| },
+\end{align}
+by directly calculating the eigenvalues of $B_J$ and $B_G$ respectively.
+We start with $B_J$,
+\begin{align}
+ \det(B_J - \lambda I)
+ &= \det
+ \begin{pmatrix} -\lambda &
+ -\frac{a_{12}}{a_{11}} \\ -\frac{a_{21}}{a_{22}} & -\lambda
+ \end{pmatrix}\\
+ &= \lambda^2 - \frac{a_{12}a_{21}}{a_{11}a_{22}} = 0\\
+ \Rightarrow \lambda^2 &= \frac{a_{12}a_{21}}{a_{11}a_{22}}
+\end{align}
+For $B_G$ we have
+\begin{align}
+ \det(B_G - \lambda I)
+ &= \det
+ \begin{pmatrix}
+ -\lambda & -\frac{a_{12}}{a_{11}} \\ 0 & -\lambda -
+ \frac{a_{21}a_{12}}{a_{11}a_{22}}\\
+ \end{pmatrix}\\
+ &= -\lambda \left(-\lambda - \frac{a_{21}a_{12}}{a_{11}a_{22}}\right) - 0
+ \\
+ \Rightarrow \lambda_1 &= 0, \qquad \lambda_2 = -
+ \frac{a_{21}a_{12}}{a_{11}a_{22}}.
+\end{align}
+Which satisfies the above condition, remember
+\begin{align}
+ \rho(A) = \max\{ \left| \lambda \right| : \text{$\lambda$ is eigenvalue
+ $A$}\},
+\end{align}
+especially note the absolute value.
+
+\qed
+\subsubsection{}
+Indeed the error of the Gauss-Siedel iteration converges to $0$ if $\rho(B_G)
+< 1$. Which is the case, because exactly then $\rho(B_G) = \rho(B_J)^2 < 1$.
+Where we can also conclude that the Gauss-Siedel method converges twines as
+fast as the Jacobi method for $2 \times 2$ matrices.
+\subsection{Problem 3}
+Let $Q \in \mathbb{R}^{n \times n}$ be the Poisson matrix, the matrix of
+the finite difference method on a $n \times n$ grid,
+\begin{align}
+ Q =
+ \begin{pmatrix}
+ 2 & -1 & & & \\
+ -1& 2 & -1& & \\
+ & \ddots & \ddots & \ddots & \\
+ & & -1 & 2 & -1\\
+ & & & -1 & 2
+ \end{pmatrix}
+\end{align}
+\subsubsection{}
+The eigenvalues of $Q$ lie in the interval $[0 ,4]$. To show this let
+$(\lambda , v)$ be an eigen-pair of $Q$, they satisfy the equation
+\begin{align}
+ Qv = \lambda v
+\end{align}
+At the $k$-th ($k \in \{1,\ldots,n\}$) step we have $v^{(0)} = 0$, $v^{(1)} =
+1$ and $v^{(n+1)} = 0$
+\begin{align}
+ -v^{(k+1} + 2v^{(k)} - v^{(k+1)} &= \lambda v^{(k)}\\
+ \Rightarrow v^{(k+1)} &= (2-\lambda) v^{(k)} - v^{(k-1)},
+\end{align}
+which are the Chebyshev polynomials of the second kind, where $(2-\lambda_k)$
+satisfies
+\begin{align}
+ (2-\lambda_k) &= 2\cdot \cos\left(\frac{k\pi}{n+1}\right)\\
+ \Rightarrow \lambda_k &= 2\left( 1 - \cos\left( \frac{k\pi}{n+1} \right)
+ \right) \\
+ &=4\cdot\sin^2\left( \frac{k\pi}{n+1} \right),
+\end{align}
+thereby we can conclude that $\lambda_k \in [0, 4] \;\; \forall k \in \{1,
+\dots n\}$. \qed
+\subsubsection{}
+\subsubsection{}
+In this section we write a Python script that returns the matrix $Q$ given
+$n$ and for $b = (1, \ldots ,1)^T$ $n=20$ implements the Gauss-Siedel method
+for to solve $Qx = b$ for $x$.
+%\includepdf[pages=1, pagecommand={}, width=\textwidth]{./prog/prb2_p.pdf}
+\begin{figure}[htpb]
+ \centering
+ \includegraphics[width=0.8\textwidth, clip, trim=0cm 5cm 0cm 10cm]{./prog/prb2_p.pdf}
+\end{figure}
+\subsection{Problem 4}
+Now let $P_n \in \mathbb{R}^{n\times n}$ be the matrix
+\begin{align}
+ P_n =
+ \begin{pmatrix}
+ 2 & -1 & & & -1 \\
+ -1& 2 & -1& & \\
+ & \ddots & \ddots & \ddots & \\
+ & & -1 & 2 & -1\\
+ -1 & & & -1 & 2
+ \end{pmatrix}
+\end{align}
+\subsubsection{}
+We can show that all the eigenvalues of $P_n$ are in $[0, 4]$ because $P_n$
+is the finite differenece/laplacian matrix for periodic boundary conditions
+and can be diagonalized by a DFT. So let $\left( \lambda, v \right)$ be the
+eigen-pair of $P_n$, which satisfy the equation
+\begin{align}
+ P_n v = \lambda v
+\end{align}
+This she standard Possion with periodic boundary conditions, with the
+eigenvector $v_j = \omega^{jk} = e^{2\pi i \frac{jk}{n}}$ for a $j \in
+\{1,\ldots,n\}$.
+\begin{align}
+ (P_n v)_j
+ &= 2\omega^{jk}- \omega^{(j-1)k} - \omega^{(j+1)k}\\
+ &= \omega^{jk}(2 - \omega^{-k} - \omega^{k})\\
+ &= \omega^{jk}(2 - 2\cos\left( \frac{2\pi k}{n} \right)\\
+ &= 4\sin^2\left( \frac{2\pi k}{n} \right) \omega^{jk} = \lambda_k v^k_j,
+\end{align}
+thereby we can conclude that $\lambda_k \in [0, 4]$ forall $k$.
+\subsubsection{}
+Because $P_n$ is a real, symmetric, cercular matrix the orthogonal components
+of the eigenvalues are also eigenvectors, i.e. $\text{Re}\left(v \right)$ and
+$\text{Im}\left(v \right) $. We may conclude this by pure calculation. In the $j-th$
+component we have $k$ eigenvalues
+\begin{align}
+ \left(P_n\text{Re}\left( v^k \right)\right) _j
+ &= 2\text{Re}\left( \omega^{jk} \right) -\text{Re}\left( \omega^{(j-1)k}
+ \right) -\text{Re}\left( \omega^{(j+1)k} \right)\\
+ &= \omega^{jk} +\omega^{-jk} - \frac{1}{2}\omega^{(j-1)k}
+ -\frac{1}{2}\omega^{-(j-1)k}
+ -\frac{1}{2}\omega^{(j+1)k} -\frac{1}{2}\omega^{-(j+1)k}\\
+ &= \left( \omega^{jk} + \omega^{-jk} \right)
+ -\frac{1}{2}\left( \omega^{jk}\left(\omega^{k} + \omega^{-k} \right)
+ +\omega^{-jk}\left(\omega^{k} + \omega^{-k} \right) \right)\\
+ &= \frac{1}{2}\left( \omega^{jk} + \omega^{-jk} \right)\left(
+ 2-(\omega^{k} + \omega^{-k}\right) \\
+ &= \text{Re}\left(\omega^{jk} \right) \left(2-2\cos\left( \frac{2\pi
+ k}{n}\right)\right)\\
+ &= 4\sin^2\left( 2\pi \frac{k}{n} \right) \text{Re}\left( \omega^{jk}
+ \right) = \lambda_k\text{Re}(v^k_j).
+\end{align}
+\subsubsection{}
+We define the quantity $m(n) = \min\{|\lambda|:\lambda\; \text{eigenvalue of
+$P_n$}\}$. We can show that the quantity converges to 0 as $n$ goes to
+infinity by calculating for a $k$ that minimises $\lambda_k$ which is for a
+$k\neq$.
+\begin{align}
+ \lim_{n \to \infty}m(n)
+ &= \lim_{n \to \infty}\min_{k\in\{1,\ldots,n\}}\{|\lambda_k(n)|\} = \lim_{n
+ \to \infty}4\sin^2\left(2\pi \frac{k}{n}\right)\\
+ &= \lim_{n \to \infty}4\cdot\left(\frac{x^2}{n^{2}} + \frac{x^{4}}{n^{4}}
+ + O(\frac{1}{n^{6}})\right) = 0,
+\end{align}
+where $x = 2\pi k$.
+\subsection{Problem 5}
+Let $Q$ be like in Problem 3. And split $Q$ as
+\begin{align}
+ Q = D-N,
+\end{align}
+where $D$ consists of diagonal entries of $Q$. For $p \in \mathbb{N}$ let $C_p$ be the
+Neumann polynomial preconditioner , defined as
+\begin{align}
+ C_p = D^{-1}\sum_{k=0}^{p} \left( ND^{-1} \right)^{k}
+\end{align}
+\subsubsection{}
+\subsubsection{}
+The following is a Phython script, that takes $n, p$ as an Input and returns
+$C_p$ furthermore calculates the spectral condition number of the matrix
+$C_pQ$
+\begin{figure}[htpb]
+ \centering
+ \includegraphics[page=2 ,width=\textwidth, clip=true, trim=0cm 7cm 0cm
+ 3cm]{./prog/prb2_p.pdf}
+ \includegraphics[page=3 ,width=0.8\textwidth, clip=true, trim=0cm 17cm 0cm
+ 2.5cm]{./prog/prb2_p.pdf}
+\end{figure}
+%\printbibliography
+\end{document}
+
diff --git a/num_ana/preamble.tex b/num_ana/preamble.tex
@@ -39,7 +39,6 @@
\usepackage[parfill]{parskip}
\usepackage{lipsum}
-
\usepackage{tcolorbox}
\tcbuselibrary{skins,breakable}
diff --git a/num_ana/prog/.ipynb_checkpoints/Untitled-checkpoint.ipynb b/num_ana/prog/.ipynb_checkpoints/Untitled-checkpoint.ipynb
@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/num_ana/prog/prb2_p.ipynb b/num_ana/prog/prb2_p.ipynb
@@ -0,0 +1,165 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "517aef32",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "85da4812",
+ "metadata": {},
+ "source": [
+ "# Sheet 2, Exercise 3"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 53,
+ "id": "e14c738f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def gauss_siedel(A, b, k):\n",
+ " n, m = Q.shape\n",
+ " D = np.reshape([Q[i][j] if i==j else 0 for i in range(n) for j in range(n)], (n ,n))\n",
+ " L = np.reshape([Q[i][j] if i>j else 0 for i in range(n) for j in range(n)], (n ,n))\n",
+ " U = np.reshape([Q[i][j] if i<j else 0 for i in range(n) for j in range(n)], (n ,n))\n",
+ "\n",
+ " x = np.random.rand(n)\n",
+ " for i in range(k):\n",
+ " x = np.linalg.inv(D)@(b - (L + U)@x)\n",
+ " return x\n",
+ "\n",
+ "def poisson_mat(n, m=None):\n",
+ " return 2 * np.eye(n, m) + (-1) * np.eye(n, m, k=1) + (-1) * np.eye(n, m, k=-1)\n",
+ "\n",
+ "# test\n",
+ "for n in range(5, 20):\n",
+ " Q = poisson_mat(n)\n",
+ " b = np.ones(n)\n",
+ " x = gauss_siedel(Q, b, k=1000)\n",
+ " np.testing.assert_allclose(Q@x, b, rtol=1e-5, err_msg=f'GS failed at dim - {n}')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "51428a2b",
+ "metadata": {},
+ "source": [
+ "# Sheet 2, Exercise 5"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 60,
+ "id": "caef181e",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "1\t31706.88589269506\n",
+ "2\t534.147900485917\n",
+ "3\t31706.885892666996\n",
+ "4\t890.0968016459952\n",
+ "5\t31706.885892699214\n",
+ "6\t1245.8213126008322\n",
+ "7\t31706.885892666185\n",
+ "8\t1601.2319940002324\n",
+ "9\t31706.88589270954\n",
+ "Max. and Min. Singular value are far apart from each other for uneaven p\n"
+ ]
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbsAAAD4CAYAAAB10khoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVYUlEQVR4nO3df6zddZ3n8efLFrTqOEW5Q6DFLRkbZupMFvQEmWXWuLjSwk6mHWMMJquNIVM3AxvcnbAD/sOMmoyGHd2QKAkjrGVXrSwiNC5jbZCsu3+A3FLGUpDlDsLQK9I7loKujfzY9/5xPnUPzG3vbe89Pfd+7/ORnNzveX9/nPf5pu3rfr7fT89JVSFJUpe9ZtQNSJI0bIadJKnzDDtJUucZdpKkzjPsJEmdt3zUDRyvU089tdasWTPqNiRJC8SuXbv+oarGplu3aMNuzZo1jI+Pj7oNSdICkeTJI63zMqYkqfMMO0lS5xl2kqTOM+wkSZ1n2EmSOm/Rzsacizt2T3Ldjkf58cFDnLFyBVetP5tN564adVud4jkePs/x8HmOh+9EneMlF3Z37J7kmtv3cOjFlwGYPHiIa27fA+Af4nniOR4+z/HweY6H70Se4yV3GfO6HY/+6sQedujFl7lux6Mj6qh7PMfD5zkePs/x8J3Ic7zkwu7HBw8dU13HznM8fJ7j4fMcD9+JPMdLLuzOWLnimOo6dp7j4fMcD5/nePhO5DlecmF31fqzWXHSslfUVpy0jKvWnz2ijrrHczx8nuPh8xwP34k8x0tugsrhm57OsBoez/HweY6Hz3M8fCfyHKeq5v2gJ0Kv1ys/CFqSdFiSXVXVm27dkruMKUlaegw7SVLnGXaSpM6bMeySvC7J95P8bZK9Sf6i1c9Kcl+SiSRfT3Jyq7+2PZ9o69cMHOuaVn80yfqB+oZWm0hy9RDepyRpCZvNyO6XwIVV9U+Bc4ANSc4HPgt8vqreBjwLXNa2vwx4ttU/37YjyTrgUuDtwAbgi0mWJVkGfAG4GFgHfKhtK0nSvJgx7Krv5+3pSe1RwIXAba2+FdjUlje257T1702SVt9WVb+sqh8BE8B57TFRVY9X1QvAtratJEnzYlb37NoI7EFgP7AT+DvgYFW91DbZBxz+jxGrgKcA2vrngLcM1l+1z5Hq0/WxJcl4kvGpqanZtC5J0uzCrqperqpzgNX0R2K/NcymjtLHjVXVq6re2NjYKFqQJC1CxzQbs6oOAvcAvwesTHL4E1hWA5NteRI4E6Ct/3Xgp4P1V+1zpLokSfNiNrMxx5KsbMsrgPcBj9APvQ+0zTYDd7bl7e05bf13q/8xLduBS9tszbOAtcD3gfuBtW1258n0J7Fsn4f3JkkSMLvPxjwd2NpmTb4GuLWqvpXkYWBbkk8Du4Gb2vY3Af8lyQRwgH54UVV7k9wKPAy8BFxeVS8DJLkC2AEsA26uqr3z9g4lSUuen40pSeoEPxtTkrSkGXaSpM4z7CRJnWfYSZI6z7CTJHWeYSdJ6jzDTpLUeYadJKnzDDtJUucZdpKkzjPsJEmdZ9hJkjrPsJMkdZ5hJ0nqPMNOktR5hp0kqfMMO0lS5xl2kqTOM+wkSZ1n2EmSOs+wkyR1nmEnSeo8w06S1Hkzhl2SM5Pck+ThJHuTXNnqf55kMsmD7XHJwD7XJJlI8miS9QP1Da02keTqgfpZSe5r9a8nOXm+36gkaemazcjuJeBPq2odcD5weZJ1bd3nq+qc9rgLoK27FHg7sAH4YpJlSZYBXwAuBtYBHxo4zmfbsd4GPAtcNk/vT5KkmcOuqp6uqgfa8s+AR4BVR9llI7Ctqn5ZVT8CJoDz2mOiqh6vqheAbcDGJAEuBG5r+28FNh3n+5Ek6R85pnt2SdYA5wL3tdIVSX6Q5OYkp7TaKuCpgd32tdqR6m8BDlbVS6+qT/f6W5KMJxmfmpo6ltYlSUvYrMMuyRuBbwAfr6rngRuA3wTOAZ4G/moYDQ6qqhurqldVvbGxsWG/nCSpI5bPZqMkJ9EPuq9U1e0AVfXMwPq/Br7Vnk4CZw7svrrVOEL9p8DKJMvb6G5we0mS5mw2szED3AQ8UlWfG6ifPrDZHwEPteXtwKVJXpvkLGAt8H3gfmBtm3l5Mv1JLNurqoB7gA+0/TcDd87tbUmS9P/NZmR3AfBhYE+SB1vtE/RnU54DFPAE8DGAqtqb5FbgYfozOS+vqpcBklwB7ACWATdX1d52vD8DtiX5NLCbfrhKkjQv0h9YLT69Xq/Gx8dH3YYkaYFIsquqetOt8xNUJEmdZ9hJkjrPsJMkdZ5hJ0nqPMNOktR5hp0kqfMMO0lS5xl2kqTOM+wkSZ1n2EmSOs+wkyR1nmEnSeo8w06S1HmGnSSp8ww7SVLnGXaSpM4z7CRJnWfYSZI6z7CTJHWeYSdJ6jzDTpLUeYadJKnzZgy7JGcmuSfJw0n2Jrmy1d+cZGeSx9rPU1o9Sa5PMpHkB0neMXCszW37x5JsHqi/M8mets/1STKMNytJWppmM7J7CfjTqloHnA9cnmQdcDVwd1WtBe5uzwEuBta2xxbgBuiHI3At8C7gPODawwHZtvnjgf02zP2tSZLUN2PYVdXTVfVAW/4Z8AiwCtgIbG2bbQU2teWNwC3Vdy+wMsnpwHpgZ1UdqKpngZ3AhrbuTVV1b1UVcMvAsSRJmrNjumeXZA1wLnAfcFpVPd1W/QQ4rS2vAp4a2G1fqx2tvm+a+nSvvyXJeJLxqampY2ldkrSEzTrskrwR+Abw8ap6fnBdG5HVPPf2j1TVjVXVq6re2NjYsF9OktQRswq7JCfRD7qvVNXtrfxMuwRJ+7m/1SeBMwd2X91qR6uvnqYuSdK8mM1szAA3AY9U1ecGVm0HDs+o3AzcOVD/SJuVeT7wXLvcuQO4KMkpbWLKRcCOtu75JOe31/rIwLEkSZqz5bPY5gLgw8CeJA+22ieAzwC3JrkMeBL4YFt3F3AJMAH8AvgoQFUdSPIp4P623Ser6kBb/hPgy8AK4G/aQ5KkeZH+7bbFp9fr1fj4+KjbkCQtEEl2VVVvunV+gookqfMMO0lS5xl2kqTOM+wkSZ1n2EmSOs+wkyR1nmEnSeo8w06S1HmGnSSp8ww7SVLnGXaSpM4z7CRJnWfYSZI6z7CTJHWeYSdJ6jzDTpLUeYadJKnzDDtJUucZdpKkzjPsJEmdZ9hJkjrPsJMkdZ5hJ0nqvBnDLsnNSfYneWig9udJJpM82B6XDKy7JslEkkeTrB+ob2i1iSRXD9TPSnJfq389ycnz+QYlSZrNyO7LwIZp6p+vqnPa4y6AJOuAS4G3t32+mGRZkmXAF4CLgXXAh9q2AJ9tx3ob8Cxw2VzekCRJrzZj2FXV94ADszzeRmBbVf2yqn4ETADntcdEVT1eVS8A24CNSQJcCNzW9t8KbDq2tyBJ0tHN5Z7dFUl+0C5zntJqq4CnBrbZ12pHqr8FOFhVL72qPq0kW5KMJxmfmpqaQ+uSpKXkeMPuBuA3gXOAp4G/mq+GjqaqbqyqXlX1xsbGTsRLSpI6YPnx7FRVzxxeTvLXwLfa00ngzIFNV7caR6j/FFiZZHkb3Q1uL0nSvDiukV2S0wee/hFweKbmduDSJK9NchawFvg+cD+wts28PJn+JJbtVVXAPcAH2v6bgTuPpydJko5kxpFdkq8B7wFOTbIPuBZ4T5JzgAKeAD4GUFV7k9wKPAy8BFxeVS+341wB7ACWATdX1d72En8GbEvyaWA3cNN8vTlJkgDSH1wtPr1er8bHx0fdhiRpgUiyq6p6063zE1QkSZ1n2EmSOs+wkyR1nmEnSeo8w06S1HmGnSSp8ww7SVLnGXaSpM4z7CRJnWfYSZI6z7CTJHWeYSdJ6jzDTpLUeYadJKnzDDtJUucZdpKkzjPsJEmdZ9hJkjrPsJMkdZ5hJ0nqPMNOktR5hp0kqfMMO0lS580YdkluTrI/yUMDtTcn2ZnksfbzlFZPkuuTTCT5QZJ3DOyzuW3/WJLNA/V3JtnT9rk+Seb7TUqSlrbZjOy+DGx4Ve1q4O6qWgvc3Z4DXAysbY8twA3QD0fgWuBdwHnAtYcDsm3zxwP7vfq1JEmakxnDrqq+Bxx4VXkjsLUtbwU2DdRvqb57gZVJTgfWAzur6kBVPQvsBDa0dW+qqnurqoBbBo4lSdK8ON57dqdV1dNt+SfAaW15FfDUwHb7Wu1o9X3T1KeVZEuS8STjU1NTx9m6JGmpmfMElTYiq3noZTavdWNV9aqqNzY2diJeUpLUAccbds+0S5C0n/tbfRI4c2C71a12tPrqaeqSJM2b4w277cDhGZWbgTsH6h9pszLPB55rlzt3ABclOaVNTLkI2NHWPZ/k/DYL8yMDx5IkaV4sn2mDJF8D3gOcmmQf/VmVnwFuTXIZ8CTwwbb5XcAlwATwC+CjAFV1IMmngPvbdp+sqsOTXv6E/ozPFcDftIckSfMm/Vtui0+v16vx8fFRtyFJWiCS7Kqq3nTr/AQVSVLnGXaSpM4z7CRJnWfYSZI6z7CTJHWeYSdJ6jzDTpLUeYadJKnzDDtJUucZdpKkzjPsJEmdZ9hJkjrPsJMkdZ5hJ0nqPMNOktR5hp0kqfMMO0lS5xl2kqTOM+wkSZ1n2EmSOs+wkyR1nmEnSeq8OYVdkieS7EnyYJLxVntzkp1JHms/T2n1JLk+yUSSHyR5x8BxNrftH0uyeW5vSZKkV5qPkd2/qKpzqqrXnl8N3F1Va4G723OAi4G17bEFuAH64QhcC7wLOA+49nBASpI0H4ZxGXMjsLUtbwU2DdRvqb57gZVJTgfWAzur6kBVPQvsBDYMoS9J0hI117Ar4DtJdiXZ0mqnVdXTbfknwGlteRXw1MC++1rtSHVJkubF8jnu//tVNZnkN4CdSX44uLKqKknN8TV+pQXqFoC3vvWt83VYSVLHzWlkV1WT7ed+4Jv077k90y5P0n7ub5tPAmcO7L661Y5Un+71bqyqXlX1xsbG5tK6JGkJOe6wS/KGJL92eBm4CHgI2A4cnlG5GbizLW8HPtJmZZ4PPNcud+4ALkpySpuYclGrSZI0L+ZyGfM04JtJDh/nq1X17ST3A7cmuQx4Evhg2/4u4BJgAvgF8FGAqjqQ5FPA/W27T1bVgTn0JUnSK6Rq3m6pnVC9Xq/Gx8dH3YYkaYFIsmvgv8G9gp+gIknqPMNOktR5hp0kqfMMO0lS5xl2kqTOM+wkSZ1n2EmSOs+wkyR1nmEnSeo8w06S1HmGnSSp8ww7SVLnGXaSpM6b6zeVS5IWiDt2T3Ldjkf58cFDnLFyBVetP5tN564adVsLgmEnSR1wx+5Jrrl9D4defBmAyYOHuOb2PQAGHoadpBPEUcdwXbfj0V8F3WGHXnyZ63Y86nnGsJN0AjjqGL4fHzx0TPWlxgkqkobuaKMOzY8zVq44pvpSY9hJzR27J7ngM9/lrKv/Oxd85rvcsXty1C11hqOO4btq/dmsOGnZK2orTlrGVevPHlFHC4uXMSW8zDZsZ6xcweQ0weaoY/4c/nPqfdHpGXaLhDf3h8ub+8N11fqzX/HLBDjqGIZN567yz+sRGHaLgKOO4fMy23A56tCoGXaLgKOO4fMy2/A56tAoLZgJKkk2JHk0yUSSq0fdz0LiqGP4vLkvdduCCLsky4AvABcD64APJVk32q4WDqcUD9+mc1fxl+//XVatXEGAVStX8Jfv/11HIlJHLJTLmOcBE1X1OECSbcBG4OGRdrVAeHP/xPAym9RdC2JkB6wCnhp4vq/VXiHJliTjScanpqZOWHOj5qhDkuZmoYzsZqWqbgRuBOj1ejXidk4oRx2SdPwWyshuEjhz4PnqVpMkac4WStjdD6xNclaSk4FLge0j7kmS1BEL4jJmVb2U5ApgB7AMuLmq9o64LUlSRyyIsAOoqruAu0bdhySpexbKZUxJkoYmVYtzUmOSKeDJOR7mVOAf5qGdE2mx9bzY+oXF17P9Dt9i63mx9Qvz0/M/qaqx6VYs2rCbD0nGq6o36j6OxWLrebH1C4uvZ/sdvsXW82LrF4bfs5cxJUmdZ9hJkjpvqYfdjaNu4Dgstp4XW7+w+Hq23+FbbD0vtn5hyD0v6Xt2kqSlYamP7CRJS4BhJ0nqvCUZdkluTrI/yUOj7mU2kpyZ5J4kDyfZm+TKUfc0kySvS/L9JH/bev6LUfc0G0mWJdmd5Fuj7mU2kjyRZE+SB5OMj7qfmSRZmeS2JD9M8kiS3xt1T0eS5Ox2Xg8/nk/y8VH3NZMk/679nXsoydeSvG7UPR1Nkitbr3uHeX6X5D27JO8Gfg7cUlW/M+p+ZpLkdOD0qnogya8Bu4BNVbVgv9w2SYA3VNXPk5wE/C/gyqq6d8StHVWSfw/0gDdV1R+Mup+ZJHkC6FXVovgPxEm2Av+zqr7UPvT99VV1cMRtzSjJMvrfxPKuqprrh1kMTZJV9P+urauqQ0luBe6qqi+PtrPpJfkdYBv9L/B+Afg28G+qamK+X2tJjuyq6nvAgVH3MVtV9XRVPdCWfwY8wjRfbruQVN/P29OT2mNB/2aVZDXwr4AvjbqXLkry68C7gZsAquqFxRB0zXuBv1vIQTdgObAiyXLg9cCPR9zP0fw2cF9V/aKqXgL+B/D+YbzQkgy7xSzJGuBc4L4RtzKjdknwQWA/sLOqFnrP/wn4D8D/HXEfx6KA7yTZlWTLqJuZwVnAFPCf26XiLyV5w6ibmqVLga+NuomZVNUk8B+BvweeBp6rqu+Mtqujegj450nekuT1wCW88rtN541ht4gkeSPwDeDjVfX8qPuZSVW9XFXn0P8y3vPaJYsFKckfAPurateoezlGv19V7wAuBi5vl+gXquXAO4Abqupc4P8AV4+2pZm1y61/CPy3UfcykySnABvp/2JxBvCGJP96tF0dWVU9AnwW+A79S5gPAi8P47UMu0Wi3ff6BvCVqrp91P0ci3ap6h5gw4hbOZoLgD9s98C2ARcm+a+jbWlm7Td5qmo/8E369z4Wqn3AvoER/m30w2+huxh4oKqeGXUjs/AvgR9V1VRVvQjcDvyzEfd0VFV1U1W9s6reDTwL/O9hvI5htwi0yR43AY9U1edG3c9sJBlLsrItrwDeB/xwpE0dRVVdU1Wrq2oN/UtW362qBfsbMUCSN7QJS7TLgRfRvyy0IFXVT4CnkpzdSu8FFuwkqwEfYhFcwmz+Hjg/yevbvxvvpX+Pf8FK8hvt51vp36/76jBeZ8F8eeuJlORrwHuAU5PsA66tqptG29VRXQB8GNjT7oEBfKJ94e1CdTqwtc1iew1wa1Utiun8i8hpwDf7/6axHPhqVX17tC3N6N8CX2mXBh8HPjrifo6q/RLxPuBjo+5lNqrqviS3AQ8ALwG7WfgfHfaNJG8BXgQuH9akpSX5Xw8kSUuLlzElSZ1n2EmSOs+wkyR1nmEnSeo8w06S1HmGnSSp8ww7SVLn/T96gU+siVi/EAAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 504x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "def neumann_polynomial_preconditioner(n, p):\n",
+ " Q = poisson_mat(n)\n",
+ " D = np.reshape([Q[i][j] if i==j else 0 for i in range(n) for j in range(n)], (n ,n))\n",
+ " N = Q - D\n",
+ " C_p = np.zeros([n, n])\n",
+ " for k in range(p+1):\n",
+ " C_p += np.linalg.matrix_power(N @ np.linalg.inv(D), k)\n",
+ " return np.linalg.inv(D) @ C_p\n",
+ " \n",
+ " \n",
+ "n = 20\n",
+ "Q = poisson_mat(n)\n",
+ "P = np.arange(1, 10)\n",
+ "cond_2 = []\n",
+ "for p in P:\n",
+ " C_p = neumann_polynomial_preconditioner(n, p)\n",
+ " cond_2.append(np.linalg.cond(C_p @ Q, p=2))\n",
+ " print(p, cond_2[p-1], sep='\\t')\n",
+ " \n",
+ "plt.figure(figsize=[7, 4])\n",
+ "plt.scatter(P, cond_2)\n",
+ "print(\"Max. and Min. Singular value are far apart from each other for uneaven p\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "73ce1ef6",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "53064abd",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "26476d82",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/num_ana/prog/prb2_p.pdf b/num_ana/prog/prb2_p.pdf
Binary files differ.
diff --git a/num_ana/sheets/sheet_2.pdf b/num_ana/sheets/sheet_2.pdf
Binary files differ.
