commit fb582e38d41ae807e1907b2352f66fe2582a36c3
parent 9ae6dd03cde73e87788f86a561b7c25935e43517
Author: miksa <milutin@popovic.xyz>
Date: Mon, 21 Mar 2022 07:00:31 +0100
done sheet 3
Diffstat:
8 files changed, 474 insertions(+), 40 deletions(-)
diff --git a/num_ana/build/prb2.pdf b/num_ana/build/prb2.pdf
Binary files differ.
diff --git a/num_ana/build/prb3.pdf b/num_ana/build/prb3.pdf
Binary files differ.
diff --git a/num_ana/prb2.tex b/num_ana/prb2.tex
@@ -161,6 +161,56 @@ 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.
+\subsubsection{}
+Now let $r \in \mathbb{R}$ and
+\begin{align}
+ A_r =
+ \begin{pmatrix}
+ 1 & r & r\\
+ r & 1 & r\\
+ r & r & 1\\
+ \end{pmatrix}
+\end{align}
+We can show that the Gauss-Siedel method for $A_r$ converges, provided that
+$t \in (-\frac{1}{2}, 1)$ and the Jacobi methods for $A_r$ does not converge
+for $r \in (\frac{1}{2}, 2)$. We start of with calculating the iteration
+matrices, then calculating their eigenvalues. Then $r$ is determined by
+$\rho(A) < 1$. For the Jacobi method we have
+\begin{align}
+ B_J = D^{-1}(E+F) =
+ \begin{pmatrix}
+ 0 & -r & -r\\
+ -r & 0 & -r\\
+ -r & -r & 0\\
+ \end{pmatrix}.
+\end{align}
+Then the eigenvalues of the matrix are
+\begin{align}
+ \det(B_J - \lambda I) &= -\lambda^3 - 2r^4 +3r^3\lambda = 0\\
+ \lambda_1 &= r \qquad \lambda_2 = -2r
+\end{align}
+Then $\rho(A) < 1$ determines the range of $r$
+\begin{align}
+ |\lambda_max| = 2|r| < 1 \Rightarrow |r| < \frac{1}{2}
+\end{align}
+We conclude that for $r\in (\frac{1}{2} , 1)$ does not converge.
+Now for the Gauss-Siedel method
+\begin{align}
+ B_G = (D-E)^{-1}F =
+ \begin{pmatrix}
+ 0 & -r & -r\\
+ 0 & r^2 & r^2-r\\
+ 0 & -r^3+r^2 & -r^3+2r^2\\
+ \end{pmatrix}.
+\end{align}
+then
+\begin{align}
+ \det\left(B_G - \lambda I \right) &= -\lambda^3-r^3\lambda^2
+ 3r^2\lambda^2 - r^3\lambda = 0\\
+ \Rightarrow \lambda_{1 / 2}&=
+ \frac{1}{2}\left(\pm\sqrt{r-4}(r-1)r^{\frac{3}{2}}-r^{3}+3r^{2}\right),
+\end{align}
+$\lambda < 1$ for $r \in \left( -\frac{1}{2}, 1\right)$.
\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,
@@ -175,6 +225,9 @@ the finite difference method on a $n \times n$ grid,
\end{pmatrix}
\end{align}
\subsubsection{}
+(Can also be done with Gershorin disks)
+\newline
+
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}
@@ -201,7 +254,6 @@ thereby we can conclude that $\lambda_k \in [0, 4] \;\; \forall k \in \{1,
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}
@@ -219,6 +271,9 @@ Now let $P_n \in \mathbb{R}^{n\times n}$ be the matrix
\end{pmatrix}
\end{align}
\subsubsection{}
+(Can also be done with Gershorin disks)
+\newline
+
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
diff --git a/num_ana/prb3.tex b/num_ana/prb3.tex
@@ -0,0 +1,260 @@
+\include{preamble.tex}
+
+\usepackage[final]{pdfpages}
+
+\begin{document}
+\maketitle
+\tableofcontents
+\section{Sheet 3}
+\subsection{Problem 1}
+Take a linear system $Ax = b$, where $A \in \mathbb{R}^{n\times n}$ and $b
+\in \mathbb{R}^{n}$. We want to solve it using the Gradient decent method, an
+iteration
+\begin{align}
+ x^{k+1} = x^{k} + \alpha_k r^{k},
+\end{align}
+where $r^{k} = b - Ax^{k}$ and the residual $\alpha_k= \frac{(r^{k})^T
+r^k}{(r^k)^T Ar^k}$.
+\subsubsection{}
+We compute $x^1$ for
+\begin{align}
+ A =
+ \begin{pmatrix}
+ 2 & -1 \\
+ -1 & 2
+ \end{pmatrix},
+ \qquad
+ b = \begin{pmatrix}1 \\ 1 \end{pmatrix},
+\end{align}
+with an initial guess $x^0 = 0$.
+\begin{align}
+ r^0 = \begin{pmatrix} 1 & 1 \end{pmatrix} \qquad
+ Ar^0 = \begin{pmatrix} 1 & 1 \end{pmatrix} \quad \Rightarrow \quad
+ \alpha_0 = 1.
+\end{align}
+Then for $x^1$ we have
+\begin{align}
+ x^1 = \begin{pmatrix} 0 \\ 0 \end{pmatrix} + \begin{pmatrix}1 \\
+1 \end{pmatrix} = b
+\end{align}
+\subsubsection{}
+Suppose the $k$-th error $e^k = x - x^k$ is an eigenvector of $A$ to the
+eigenvalue $\lambda$, then
+\begin{align}
+ A e^{k}=\lambda e^{k} = \lambda (x^{k}-x) = \lambda x^{k} - \lambda x.
+\end{align}
+For the next iteration step we need $r^{k}$ which is
+\begin{align}
+ r^{k} &= b - Ax^{k} + Ax - Ax = (b - Ax) - A(x^{k} - x) = -\lambda e^{k}\\
+ (r^k)^{T} r^k &= \lambda^2 (e^k)^{T}e^{k}\\
+ (r^k)^{T} A r^k &= \lambda^3 (e^k)^{T}e^{k}\\
+ \Rightarrow \alpha_k &= \frac{1}{\lambda}
+\end{align}
+Then the next step $x^{k+1}$ is
+\begin{align}
+ x^{k+1} = x^{k}+\alpha_k r^{k} = x^{k}-\frac{\lambda}{\lambda}e^{k} =
+ x^{k} - e^{k} = x^{k}-x^{k}+ x = x,
+\end{align}
+i.e. $x^{k+1}$ is then the solution.
+\subsection{Exercise 2}
+We show the norm equivalence of the vector norms $ \|\cdot\|_\infty$ and
+$\|\cdot\|_2$ and the optimality of the constants $C, C' > 0$and the
+optimality of the constants $C, C' > 0$
+\subsubsection{}
+We show
+\begin{align}
+ \|v\|_\infty = \max_i |v_i| \leq \sqrt{\sum_{i}^{} v_i^2} \quad
+ \Rightarrow \quad C = 1.
+\end{align}
+Then
+\begin{align}
+ \|v\|_2 = \sqrt{\sum_i v_i^2} \le \sqrt{\sum_i \max_i |v_i|^2} =
+ \sqrt{\sum_i \|v\|^2_\infty} = \sqrt{n} \|v\|_\infty.
+\end{align}
+We get
+\begin{align}
+ \|v\|_\infty \le \|v\|_2 \le \sqrt{n} \|v\|_\infty.
+\end{align}
+For $u := \frac{v}{\|v\|_\infty} \; \Rightarrow \; \|v\|_\infty = 1$ we get
+\begin{align}
+ 1 \le \|u\|_2 \le \sqrt{n}
+\end{align}
+which states optimality of the constants.
+\subsection{Problem 3}
+Now we do the same as in Problem 2 but with matrix norms induced by vector
+norms, specifically for $\|\cdot\|_2$ and $\|\cdot\|_\infty$ matrix norms
+induced by vector norms.
+\subsubsection{}
+We know that
+\begin{align}
+ \|A\|_2 = \max_{v\neq 0} \frac{\|Ax\|_2}{\|x\|_2} \qquad \|A\|_2 &=
+ \rho(A^*A)\\
+ \|A\|_\infty = \max_{i,j} |a_{ij}|
+\end{align}
+Together with the norm inequalities
+\begin{align}
+ \|Av\|_\infty &\le \|A\|_\infty \|v\|_\infty\\
+ \|A\|_2^2 = \rho(A^TA) \le \|A^TA\|_\infty\\
+ \Rightarrow \|\rho(A^TA)\|_\infty = \|A^TA x\|_\infty \le \|A^T
+ A\|_\infty \|x\|_\infty,
+\end{align}
+thereby
+\begin{align}
+ \|A^T A\|_\infty &= \max_{i,j} \left|(A^T A)_{ij})\right| = \max_{i,j}
+ \left| \sum_l A_{il} A_{lj} \right| \\
+ &=\max_{i,j} \sum_l \left| A_{il} A_{lj} \right| \leq n
+ \|A\|_\infty^2\\
+ \nonumber\\
+ \Rightarrow \|A\|_2 &\leq \sqrt{n} \|A\|_\infty
+\end{align}
+\subsubsection{}
+Let $A \in \mathbb{C}^{n\times n}$, $b \in \mathbb{C}^n$ defined as
+\begin{align}
+ A =
+ \begin{pmatrix}
+ 1 & \cdots & 1 \\
+ 0 & \cdots & 0 \\
+ \vdots & \vdots &\vdots \\
+ 0 & \cdots & 0
+ \end{pmatrix} \qquad b = \begin{pmatrix} 1 \\ \vdots \\ 1
+ \end{pmatrix}
+\end{align}
+Then $Ab = \begin{pmatrix} n & 0 & \cdots & 0 \end{pmatrix}^T$ and the
+$\|\cdot\|_2$ of $A$ is
+\begin{align}
+ \|A\|_2 = \max_{v\neq 0} \frac{\|Av\|_2}{\|v\|} = \frac{\|Ab\|_2}{\|b\|}
+ = \frac{n}{\sqrt{n}} = \sqrt{n}
+\end{align}
+\subsubsection{}
+Let $A$ be defined as above, we show that $\|A\|_\infty = \sqrt{n} \|A\|_2$,
+where $C = \sqrt{n} $ is optimal
+\begin{align}
+ \|A\|_\infty = \max_i \sum_j |A_{ij}| = n = \sqrt{n} \sqrt{n} = \sqrt{n}
+ \|A\|_2,
+\end{align}
+thereby $C$ is optimal.
+\subsubsection{}
+We show that $\|A\|_2 \le \sqrt{n} \|A\|_\infty$ for all $A \in
+\mathbb{C}^{n\times n}$, and equality holds for $A$ whose columns are all
+zero, accept the first one filled with ones. We know the norms
+$\|\cdot\|_2$ and $\|\cdot\|_\infty$ are induced by vector norms which are
+consistent. Then for all $v \in \mathbb{C}^{n}$ we have
+\begin{align}
+ \|v\|_\infty &\le\sqrt{n} \|v\|_2\\
+ \Leftrightarrow \|A\|_\infty &\le \sqrt{n} \|A\|_2 \qquad\ \forall \ A
+ \in \mathbb{C}^{n\times n}.
+\end{align}
+\subsection{Problem 4}
+Let $\langle x, y \rangle$ be the standard Euclidean scalar product on
+$\mathbb{R}^n$ and $\|\cdot\|_2$ the Euclidean norm. Let $B\in
+\mathbb{R}^{n\times n}$ be an antisymmetric matrix, i.e. $B^T = -B$. And let
+$A:= I - B$
+\subsubsection{}
+We show that for all $x \in \mathbb{R}^{n}$ we have $\langle Bx, x\rangle =0$
+and $\langle Ax, x\rangle = \|x\|_2^2$.
+\begin{align}
+ \langle Bx, x\rangle = (Bx)^T x = x^T B^T x = -x^T B x\\
+ (-x^TBx)^T &= x^TBx\; \Rightarrow\; x^T Bx = jx^TBx\\
+ \Rightarrow 2x^T B x &= 2 \langle Bx, x\rangle = 0.
+\end{align}
+And the second statement follows from the previous equation
+\begin{align}
+ \langle Ax, x\rangle &= \rangle\left(I-B\right)x,x\rangle \\
+ &=\langle x, x\rangle - \underbrace{\langle Bx,
+ x\rangle}_{=0} = \langle x, x\rangle =
+ \|x\|_2^2
+\end{align}
+\subsubsection{}
+We show that $\|Ax\|_2^2 = \|x\|_2^2 + \|Bx\|_2^2$ and that $\|A\|_2 =
+\sqrt{1+\|B\|_2^2}$. We start with the vector norm
+\begin{align}
+ \|Ax\|^2_2
+ &= \|x - Bx\|_2^2 \\
+ &= \langle x -Bx, x-Bx\rangle \\
+ &= (x^T - x^T B^T)(x - Bx)\\
+ &= x^T x - x^T B^T x - x^T Bx + x^T B^T Bx\\
+ &= x^T x + x^TB^TBx = \|x\|_2^2 + \|Bx\|_2^2
+\end{align}
+where $i, j$ run from $1$ to $n$. Then the vector induced norm follows
+\begin{align}
+ \|A\|_2^2
+ &= \|I-B\|_2^2 = \sup_{x\neq 0} \frac{\|Ax\|_2^2}{\|x\|^2_2}\\
+ &= \sup_{x\neq 0} \frac{\|x\|^2_2 +\|Bx\|_2^2 }{\|x\|_2^2} = 1+
+ \sup_{x\neq 0}\frac{\|Bx\|_2^2}{\|x\|^2_2}= 1+\|B\|_2^2,
+\end{align}
+pulling the square root on both sides gives us the answer.
+\subsubsection{}
+We show that $A$ is invertible and the inverse matrix norm is given
+\begin{align}
+ \|A^{-1}\|^2 = \max_{x\neq}\frac{\|x\|^2}{\|Ax\|^2}
+\end{align}
+To show that $A$ is invertible we use the Rank-Nullity Theorem
+\begin{align}
+ n &=\text{rg}(A) + \dim\text{Im}(A)\\
+ \Rightarrow n - \dim\text{Im}(A) = rg(A).
+\end{align}
+If $A$ is invertible it has maximal rank that means $\dim \text{Im}(A) = 0$.
+Let's look what the image of A is
+\begin{align}
+ Ax = 0 \Rightarrow \|Ax\|_2^2 = \|x\|_2^2 + \|Bx\|_2^2 = 0
+\end{align}
+holds only for $x = 0$, thereby $\dim \text{Im}(A) = 0$ and the rank is
+maximal. Now let $A^{-1}y =x$ then
+\begin{align}
+ \frac{\|A^{-1}y\|_2}{\|y\|_2} &= \frac{\|x\|_2}{\|Ax\|_2}\\
+ \Rightarrow \|A^{-1}\|_2 &= \max_{x\neq 0} \frac{\|x\|_2}{\|Ax\|_2}.
+\end{align}
+\subsubsection{}
+Next we show that $\|A^{-1}\|_2 \le 1$.
+\begin{align}
+ \|A^{-1}\|_2 = \max_{x\neq 0} \frac{\|x\|_2}{\|Ax\|_2} =
+ \frac{1}{\sqrt{1+\|B\|_2^2}} \le 1
+\end{align}
+\subsection{Problem 5}
+We take $A$ and $B$ from last exercise.
+\subsubsection{}
+We show that for $k \in \{1,\ldots,m\}$ and $\mathcal{W} \subset
+\mathbb{R}^n$ with $\dim \mathcal{W} = k$ spanned by $w_1,\ldots,w_k \in
+\mathbb{R}^n$. We show that if $x \in \mathcal{W}$ is such that
+\begin{align}
+ \langle Ax, w\rangle = \langle b , w \rangle \qquad \forall w\in
+ \mathcal{W}.
+\end{align}
+then $\|x\|_2 \le \|b\|_2$. We can choose $b = x$ then
+\begin{align}
+ \langle Ax, x\rangle &= \|x\|_2^2 = \langle b, x\rangle \leq
+ \|b\|_2\|x\|_2.\\
+ \Rightarrow \|x\|_2 \le \|b\|_2
+\end{align}
+\subsubsection{}
+Next we show that for $x :=\sum_j c_j w_j$ and
+\begin{align}
+ \sum_j c_j \langle Aw_j, w_i\rangle = \langle b , w_i\rangle
+\end{align}
+for $i = 1, \ldots , k$, we have a unique solution for $x \in \mathcal{W}$.
+We do this by showing that every solution is the 0 solution $b=0$
+\begin{align}
+ \langle 0 , w\rangle = 0 = \langle Ax ,x \rangle = \|x\|_2^2.
+\end{align}
+\subsubsection{}
+Lastly for $x^* := A^{-1}b$ we show and inequality relation
+\begin{align}
+ \|x^* - x\|_2 \le \|A\|_2 \min_{x \in \mathcal{W}}\|x^* -w\|_2.
+\end{align}
+(We take the minimum $w = x$ and calculate)
+\begin{align}
+ \|x^* - x\|_2^2
+ &= \langle A(x^* -x), x^* - x\rangle\\
+ &= \langle A(x^* -x), x^* - x +w -w\rangle\\
+ &= \langle A(x^* -x), x^* - w\rangle + \langle A(x^* -x), w - x\rangle\\
+ &= \langle A(x^* -x), x^* - w\rangle +\underbrace{\langle A(x^* -x), w -
+ x\rangle}_{=0}\\
+ \le \|A\|_2 \|x^* - x\|_2\|x^* - w\|.
+\end{align}
+When we take the minimum and divide by $\|x^* - x\|_2$ we get
+\begin{align}
+ \|x^* - x\| \le \|A\|_2\min_{x\in\mathcal{W}} \|x^* - w\|_2.
+\end{align}
+
+%\printbibliography
+\end{document}
diff --git a/num_ana/prog/.ipynb_checkpoints/prb2_p-checkpoint.ipynb b/num_ana/prog/.ipynb_checkpoints/prb2_p-checkpoint.ipynb
@@ -0,0 +1,142 @@
+{
+ "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": 2,
+ "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": 5,
+ "id": "caef181e",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "1\t44.76606865271519\n",
+ "2\t59.35975010638131\n",
+ "3\t22.760834328149002\n",
+ "4\t35.62184004487854\n",
+ "5\t15.344146462132624\n",
+ "6\t25.450588757787248\n",
+ "7\t11.636387156050118\n",
+ "8\t19.801556558635184\n",
+ "9\t9.41247817754299\n",
+ "Max. and Min. Singular value are far apart from each other for uneaven p\n"
+ ]
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAD4CAYAAAC5S3KDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUAUlEQVR4nO3db4xcV3nH8e9TxxEroN2ELJa9JnUQkVHUNHG1ioISoeAQnAIiVoQiKJVcKZLf0ApUarB5U1GVJtQSkBdVJYvQGok/iYLjRBRhIicRNKoCazbggHET0kRk7cRLyQqoVuCYpy/mrrN29s/s3J2dM3O/H8mauWfu3Tk+szO/PeeeeyYyE0mSSvMHva6AJEnzMaAkSUUyoCRJRTKgJElFMqAkSUW6YDWf7JJLLslNmzat5lNKkgp25MiRX2TmyHyPrWpAbdq0ifHx8dV8SklSwSLiuYUec4hPklQkA0qSVCQDSpJUpLYCKiKGI+K+iPhpRByLiLdFxMUR8VBEPFXdXtTtykqSmqPdHtRdwLcy863AVcAxYDdwODMvBw5X25IkrYglZ/FFxB8Bbwf+CiAzfwf8LiJuAW6odtsPPAp8ohuVbMfBiUn2HjrOiekZNgwPsWvbZrZvGe1VdSRJNbXTg7oMmAL+LSImIuILEfFaYF1mnqz2eQFYN9/BEbEzIsYjYnxqamplan2egxOT7DlwlMnpGRKYnJ5hz4GjHJyY7MrzSZK6r52AugD4M+BfM3ML8H+cN5yXre/smPd7OzJzX2aOZebYyMi812LVtvfQcWZOnzmnbOb0GfYeOt6V55MkdV87AfU88HxmPl5t30crsF6MiPUA1e2p7lRxaSemZ5ZVLkkq35IBlZkvAD+PiM1V0Y3AT4AHgR1V2Q7gga7UsA0bhoeWVS5JKl+7s/j+BvhyRPwIuBr4J+BO4KaIeAp4Z7XdE7u2bWZo7ZpzyobWrmHXts0LHCFJKl1ba/Fl5hPA2DwP3biitenQ7Gw9Z/FJ0uBY1cViu2n7llEDSZIGiEsdSZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKdEE7O0XEs8CvgTPAy5k5FhEXA/cAm4Bngdsy86XuVFOS1DTL6UG9IzOvzsyxans3cDgzLwcOV9uSJK2IOkN8twD7q/v7ge21ayNJUqXdgErg2xFxJCJ2VmXrMvNkdf8FYN18B0bEzogYj4jxqampmtWVJDVFW+eggOszczIi3gg8FBE/nftgZmZE5HwHZuY+YB/A2NjYvPtIknS+tnpQmTlZ3Z4C7geuAV6MiPUA1e2pblVSktQ8SwZURLw2Il4/ex94F/Ak8CCwo9ptB/BAtyopSWqedob41gH3R8Ts/l/JzG9FxPeBeyPiduA54LbuVVOS1DRLBlRmPgNcNU/5/wI3dqNSkiS5koQkqUgGlCSpSAaUJKlIBpQkqUgGlCSpSAaUJKlIBpQkqUgGlCSpSAaUJKlI7a5mXoSDE5PsPXScE9MzbBgeYte2zWzfMtrrakmSuqBvAurgxCR7Dhxl5vQZACanZ9hz4CiAISVJA6hvhvj2Hjp+NpxmzZw+w95Dx3tUI0lSN/VNQJ2YnllWuSSpv/VNQG0YHlpWuSSpv/VNQO3atpmhtWvOKRtau4Zd2zb3qEaSpG7qm0kSsxMhnMUnSc3QNwEFrZAykCSpGfpmiE+S1CwGlCSpSAaUJKlIBpQkqUgGlCSpSAaUJKlIBpQkqUgGlCSpSG0HVESsiYiJiPhGtX1ZRDweEU9HxD0RcWH3qilJaprl9KA+Ahybs/0Z4HOZ+RbgJeD2layYJKnZ2gqoiNgIvAf4QrUdwFbgvmqX/cD2LtRPktRQ7fagPg98HPh9tf0GYDozX662nwfmXSQvInZGxHhEjE9NTdWpqySpQZYMqIh4L3AqM4908gSZuS8zxzJzbGRkpJMfIUlqoHZWM78OeF9EvBt4DfCHwF3AcERcUPWiNgKT3aumJKlpluxBZeaezNyYmZuADwAPZ+aHgEeA91e77QAe6FotJUmNU+c6qE8AfxsRT9M6J3X3ylRJkqRlfmFhZj4KPFrdfwa4ZuWrJEmSK0lIkgplQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkoq0rOugBtXBiUn2HjrOiekZNgwPsWvbZrZvmXftW0nSKml8QB2cmGTPgaPMnD4DwOT0DHsOHAUwpCSphxo/xLf30PGz4TRr5vQZ9h463qMaSZLAgOLE9MyyyiVJq6PxAbVheGhZ5ZKk1dH4gNq1bTNDa9ecUza0dg27tm3uUY0kSeAkibMTIZzFJ0llaXxAQSukDCRJKkvjh/gkSWUyoCRJRTKgJElFMqAkSUUyoCRJRTKgJElFMqAkSUUyoCRJRTKgJElFMqAkSUUyoCRJRVoyoCLiNRHxvYj4YUT8OCI+VZVfFhGPR8TTEXFPRFzY/epKkpqinR7Ub4GtmXkVcDVwc0RcC3wG+FxmvgV4Cbi9a7WUJDXOkgGVLb+pNtdW/xLYCtxXle8HtnejgpKkZmrrHFRErImIJ4BTwEPAz4DpzHy52uV5YN7vq4iInRExHhHjU1NTK1BlSVITtBVQmXkmM68GNgLXAG9t9wkyc19mjmXm2MjISGe1lCQ1zrJm8WXmNPAI8DZgOCJmv/BwIzC5slWTJDVZO7P4RiJiuLo/BNwEHKMVVO+vdtsBPNClOkqSGqidr3xfD+yPiDW0Au3ezPxGRPwE+FpE/CMwAdzdxXpKkhpmyYDKzB8BW+Ypf4bW+ShJklacK0lIkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkorUzoW6WsTBiUn2HjrOiekZNgwPsWvbZrZvmXfdXEnSMhhQNRycmGTPgaPMnD4DwOT0DHsOHAUwpCSpJof4ath76PjZcJo1c/oMew8d71GNJGlwGFA1nJieWVa5JKl9BlQNG4aHllUuSWqfAVXDrm2bGVq75pyyobVr2LVtc49qJEmDw0kSNcxOhHAWnyStPAOqpu1bRg0kSeoCh/gkSUVqRA/Ki2klqf8MfEB5Ma0k9aeBH+LzYlpJ6k8DH1BeTCtJ/WngA8qLaSWpPw18QHkxrST1p4GfJOHFtJLUnwY+oMCLaSWpHw38EJ8kqT8ZUJKkIi05xBcRbwK+BKwDEtiXmXdFxMXAPcAm4Fngtsx8qXtVHSyubiFJi2unB/Uy8LHMvAK4FvhwRFwB7AYOZ+blwOFqW22YXd1icnqG5JXVLQ5OTPa6apJUjCUDKjNPZuYPqvu/Bo4Bo8AtwP5qt/3A9i7VceC4uoUkLW1Z56AiYhOwBXgcWJeZJ6uHXqA1BDjfMTsjYjwixqempurUdWC4uoUkLa3tgIqI1wFfBz6amb+a+1hmJq3zU6+Smfsycywzx0ZGRmpVdlC4uoUkLa2tgIqItbTC6cuZeaAqfjEi1lePrwdOdaeKg8fVLSRpaUsGVEQEcDdwLDM/O+ehB4Ed1f0dwAMrX73BtH3LKHfceiWjw0MEMDo8xB23XuksPkmaI1qjc4vsEHE98F3gKPD7qviTtM5D3QtcCjxHa5r5Lxf7WWNjYzk+Pl63zpKkARERRzJzbL7HlrwOKjP/E4gFHr6xTsUkSVqIK0lIkopkQEmSimRASZKKZEBJkorUiO+DGjQuNCupCQyoPjO70OzsWn6zC80ChpSkgeIQX59xoVlJTWEPagmlDae50KykprAHtYgSv7fJhWYlNYUBtYgSh9NcaFZSUzjEt4gSh9NmhxdLGnaUpG4woBaxYXiIyXnCqNfDadu3jBpIkgaeQ3yLGMThtIMTk1x358Nctvs/uO7Oh3t6Pk2SFmMPahGDNpzmNVSS+okBtYRBGk5bbNLHoPwfJQ0Oh/gapMRJH5K0EAOqQbyGSlI/MaAaZBAnfUgaXJ6DapBBm/QhabAZUA1TZ9JHaesSShpsBlQXDdIHulPUJa02z0F1SYkLzdZR4rqEkgabAdUlg/aB7hR1SavNgOqSQftAd4q6pNVmQHXJoH2g152i7hqAkpZryYCKiC9GxKmIeHJO2cUR8VBEPFXdXtTdavafQbvmaPuWUe649UpGh4cIYHR4iDtuvbKtCRKDdj5O0uqIzFx8h4i3A78BvpSZf1KV/TPwy8y8MyJ2Axdl5ieWerKxsbEcHx9fgWr3h0GaxVfHdXc+PO/XlowOD/HY7q09qJGkUkTEkcwcm++xJaeZZ+Z3ImLTecW3ADdU9/cDjwJLBlTTeM1RS53zcYPUDpKWp9ProNZl5snq/gvAuoV2jIidwE6ASy+9tMOna5ZBu+ao0y9+HLR2kLQ8tSdJZGuMcMFxwszcl5ljmTk2MjJS9+kaYdCmqHd6Pm7Q2kHS8nTag3oxItZn5smIWA+cWslKNd2gTVHvdA3Auu3g8KDU3zoNqAeBHcCd1e0DK1YjdTwkVrJOzsfVaQeHB6X+1840868C/wVsjojnI+J2WsF0U0Q8Bbyz2tYK8Zqjljrt4PCg1P/amcX3wQUeunGF66JKna/FGKSeQ512cHhQ6n+uZl6oTqeoL9Zz6McP2E7bweFBqf+51NGAGbQJFp3q5fDgoAyxSr1mD2rA1O05DMqwVq+GB+19SSvHgBowu7ZtPucDEtrrOQziB2svhgfrDrEO0h8JUl0O8Q2YThd1ddbbK+oMD65E78tFdaUWe1ADqJOeg7PeXlFneLBXva9Ban9plgElwFlv5+t0eLDTIVbo/I+Euu1vuKlUDvEJcNbbSqnzvVmdfsllnfavO6w4SK+dymMPSoCz3lbSave+6rR/3WFFe27qJgNKZznrrbc6/SOhTvsbbiqZAaXaenHeBQbzQ66TPxLqtH/Twq3E11wL8xyUauvFeRfo3bmX0s671Gn/Ouce67x23Qq3xfTyfFtpvzP9wh6UVkQ/zXqDzv+CL7XX1mn71zn32G89t1722kr8nekHBpR6qlfXHPXbh9zs8YZb5697r4Yk+3E4s5RQNKDUc73offXbh5zh9opOX/denW/rt55+STNrDSj1rSZ9yPVjuC113GqHW6+GJPutp1/SV/YYUOprTfmQ67dw63YoLva6L3RsO6/5QsfWed37radf0lf2GFBqrE7Crd/Ou0Bvwq3UCQlLBdtSz7vQ675YoNY5dqnfmYWO7VWPb6UZUNIy9dN5F+hNuPXjhISljl3odW8n2Do9drHfmcWObef3pU5vcbUmURhQ0ipqSrj144SEUs/3LPQ7s9ixj+3eenafhXptdXqLqzWJwoCS+kQ/hVs/Tkjot/M9Sx272O9Lp73Fdo5dSa4kITXA9i2jPLZ7K/9z53t4bPfWtj9IOl2lolerW/Ti2DoravTq2H6ZRGEPStKi6vTcVrvH14tjezXDrx97qcsVmbniP3QhY2NjOT4+vmrPJ0mroVerNtS5Rm2+cGunl1vn2PlExJHMHJv3MQNKkpqnlKWQDChJUpEWC6hakyQi4uaIOB4RT0fE7jo/S5KkuToOqIhYA/wL8OfAFcAHI+KKlaqYJKnZ6vSgrgGezsxnMvN3wNeAW1amWpKkpqsTUKPAz+dsP1+VnSMidkbEeESMT01N1Xg6SVKTdP1C3czcl5ljmTk2MjLS7aeTJA2IOhfqTgJvmrO9sSpb0JEjR34REc+18bMvAX5Ro25NYTu1x3Zqj+3UHtupPe220x8v9EDH08wj4gLgv4EbaQXT94G/yMwfd/QDz/3Z4wtNO9QrbKf22E7tsZ3aYzu1ZyXaqeMeVGa+HBF/DRwC1gBfXIlwkiQJaq7Fl5nfBL65QnWRJOmsUlcz39frCvQJ26k9tlN7bKf22E7tqd1Oq7rUkSRJ7Sq1ByVJajgDSpJUpOICygVo5xcRX4yIUxHx5JyyiyPioYh4qrq9qJd1LEFEvCkiHomIn0TEjyPiI1W5bTVHRLwmIr4XET+s2ulTVfllEfF49f67JyIu7HVdey0i1kTERER8o9q2jeYREc9GxNGIeCIixquyWu+7ogLKBWgX9e/AzeeV7QYOZ+blwOFqu+leBj6WmVcA1wIfrn6HbKtz/RbYmplXAVcDN0fEtcBngM9l5luAl4Dbe1fFYnwEODZn2zZa2Dsy8+o51z/Vet8VFVC4AO2CMvM7wC/PK74F2F/d3w9sX806lSgzT2bmD6r7v6b1wTKKbXWObPlNtbm2+pfAVuC+qrzx7RQRG4H3AF+otgPbaDlqve9KC6i2FqDVWesy82R1/wVgXS8rU5qI2ARsAR7HtnqVaujqCeAU8BDwM2A6M1+udvH9B58HPg78vtp+A7bRQhL4dkQciYidVVmt912tC3VVjszMiPCagUpEvA74OvDRzPxV6w/fFtuqJTPPAFdHxDBwP/DW3taoLBHxXuBUZh6JiBt6XJ1+cH1mTkbEG4GHIuKncx/s5H1XWg9q2QvQNtyLEbEeoLo91eP6FCEi1tIKpy9n5oGq2LZaQGZOA48AbwOGq3U2wfffdcD7IuJZWqcbtgJ3YRvNKzMnq9tTtP7guYaa77vSAur7wOXVLJkLgQ8AD/a4TiV7ENhR3d8BPNDDuhShOkdwN3AsMz875yHbao6IGKl6TkTEEHATrfN1jwDvr3ZrdDtl5p7M3JiZm2h9Fj2cmR/CNnqViHhtRLx+9j7wLuBJar7viltJIiLeTWvcd3YB2k/3tkZliIivAjfQWsL+ReDvgYPAvcClwHPAbZl5/kSKRomI64HvAkd55bzBJ2mdh7KtKhHxp7ROWq+h9YfqvZn5DxHxZlq9hYuBCeAvM/O3vatpGaohvr/LzPfaRq9Wtcn91eYFwFcy89MR8QZqvO+KCyhJkqC8IT5JkgADSpJUKANKklQkA0qSVCQDSpJUJANKklQkA0qSVKT/B002BwKIALUUAAAAAElFTkSuQmCC\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 = D-Q\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, 50)\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",
+ " if p in np.arange(1, 10):\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\")"
+ ]
+ }
+ ],
+ "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.ipynb b/num_ana/prog/prb2_p.ipynb
@@ -21,7 +21,7 @@
},
{
"cell_type": "code",
- "execution_count": 53,
+ "execution_count": 2,
"id": "e14c738f",
"metadata": {},
"outputs": [],
@@ -58,7 +58,7 @@
},
{
"cell_type": "code",
- "execution_count": 60,
+ "execution_count": 5,
"id": "caef181e",
"metadata": {},
"outputs": [
@@ -66,21 +66,21 @@
"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",
+ "1\t44.76606865271519\n",
+ "2\t59.35975010638131\n",
+ "3\t22.760834328149002\n",
+ "4\t35.62184004487854\n",
+ "5\t15.344146462132624\n",
+ "6\t25.450588757787248\n",
+ "7\t11.636387156050118\n",
+ "8\t19.801556558635184\n",
+ "9\t9.41247817754299\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",
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAD4CAYAAAC5S3KDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUAUlEQVR4nO3db4xcV3nH8e9TxxEroN2ELJa9JnUQkVHUNHG1ioISoeAQnAIiVoQiKJVcKZLf0ApUarB5U1GVJtQSkBdVJYvQGok/iYLjRBRhIicRNKoCazbggHET0kRk7cRLyQqoVuCYpy/mrrN29s/s3J2dM3O/H8mauWfu3Tk+szO/PeeeeyYyE0mSSvMHva6AJEnzMaAkSUUyoCRJRTKgJElFMqAkSUW6YDWf7JJLLslNmzat5lNKkgp25MiRX2TmyHyPrWpAbdq0ifHx8dV8SklSwSLiuYUec4hPklQkA0qSVCQDSpJUpLYCKiKGI+K+iPhpRByLiLdFxMUR8VBEPFXdXtTtykqSmqPdHtRdwLcy863AVcAxYDdwODMvBw5X25IkrYglZ/FFxB8Bbwf+CiAzfwf8LiJuAW6odtsPPAp8ohuVbMfBiUn2HjrOiekZNgwPsWvbZrZvGe1VdSRJNbXTg7oMmAL+LSImIuILEfFaYF1mnqz2eQFYN9/BEbEzIsYjYnxqamplan2egxOT7DlwlMnpGRKYnJ5hz4GjHJyY7MrzSZK6r52AugD4M+BfM3ML8H+cN5yXre/smPd7OzJzX2aOZebYyMi812LVtvfQcWZOnzmnbOb0GfYeOt6V55MkdV87AfU88HxmPl5t30crsF6MiPUA1e2p7lRxaSemZ5ZVLkkq35IBlZkvAD+PiM1V0Y3AT4AHgR1V2Q7gga7UsA0bhoeWVS5JKl+7s/j+BvhyRPwIuBr4J+BO4KaIeAp4Z7XdE7u2bWZo7ZpzyobWrmHXts0LHCFJKl1ba/Fl5hPA2DwP3biitenQ7Gw9Z/FJ0uBY1cViu2n7llEDSZIGiEsdSZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKdEE7O0XEs8CvgTPAy5k5FhEXA/cAm4Bngdsy86XuVFOS1DTL6UG9IzOvzsyxans3cDgzLwcOV9uSJK2IOkN8twD7q/v7ge21ayNJUqXdgErg2xFxJCJ2VmXrMvNkdf8FYN18B0bEzogYj4jxqampmtWVJDVFW+eggOszczIi3gg8FBE/nftgZmZE5HwHZuY+YB/A2NjYvPtIknS+tnpQmTlZ3Z4C7geuAV6MiPUA1e2pblVSktQ8SwZURLw2Il4/ex94F/Ak8CCwo9ptB/BAtyopSWqedob41gH3R8Ts/l/JzG9FxPeBeyPiduA54LbuVVOS1DRLBlRmPgNcNU/5/wI3dqNSkiS5koQkqUgGlCSpSAaUJKlIBpQkqUgGlCSpSAaUJKlIBpQkqUgGlCSpSAaUJKlI7a5mXoSDE5PsPXScE9MzbBgeYte2zWzfMtrrakmSuqBvAurgxCR7Dhxl5vQZACanZ9hz4CiAISVJA6hvhvj2Hjp+NpxmzZw+w95Dx3tUI0lSN/VNQJ2YnllWuSSpv/VNQG0YHlpWuSSpv/VNQO3atpmhtWvOKRtau4Zd2zb3qEaSpG7qm0kSsxMhnMUnSc3QNwEFrZAykCSpGfpmiE+S1CwGlCSpSAaUJKlIBpQkqUgGlCSpSAaUJKlIBpQkqUgGlCSpSG0HVESsiYiJiPhGtX1ZRDweEU9HxD0RcWH3qilJaprl9KA+Ahybs/0Z4HOZ+RbgJeD2layYJKnZ2gqoiNgIvAf4QrUdwFbgvmqX/cD2LtRPktRQ7fagPg98HPh9tf0GYDozX662nwfmXSQvInZGxHhEjE9NTdWpqySpQZYMqIh4L3AqM4908gSZuS8zxzJzbGRkpJMfIUlqoHZWM78OeF9EvBt4DfCHwF3AcERcUPWiNgKT3aumJKlpluxBZeaezNyYmZuADwAPZ+aHgEeA91e77QAe6FotJUmNU+c6qE8AfxsRT9M6J3X3ylRJkqRlfmFhZj4KPFrdfwa4ZuWrJEmSK0lIkgplQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkoq0rOugBtXBiUn2HjrOiekZNgwPsWvbZrZvmXftW0nSKml8QB2cmGTPgaPMnD4DwOT0DHsOHAUwpCSphxo/xLf30PGz4TRr5vQZ9h463qMaSZLAgOLE9MyyyiVJq6PxAbVheGhZ5ZKk1dH4gNq1bTNDa9ecUza0dg27tm3uUY0kSeAkibMTIZzFJ0llaXxAQSukDCRJKkvjh/gkSWUyoCRJRTKgJElFMqAkSUUyoCRJRTKgJElFMqAkSUUyoCRJRTKgJElFMqAkSUUyoCRJRVoyoCLiNRHxvYj4YUT8OCI+VZVfFhGPR8TTEXFPRFzY/epKkpqinR7Ub4GtmXkVcDVwc0RcC3wG+FxmvgV4Cbi9a7WUJDXOkgGVLb+pNtdW/xLYCtxXle8HtnejgpKkZmrrHFRErImIJ4BTwEPAz4DpzHy52uV5YN7vq4iInRExHhHjU1NTK1BlSVITtBVQmXkmM68GNgLXAG9t9wkyc19mjmXm2MjISGe1lCQ1zrJm8WXmNPAI8DZgOCJmv/BwIzC5slWTJDVZO7P4RiJiuLo/BNwEHKMVVO+vdtsBPNClOkqSGqidr3xfD+yPiDW0Au3ezPxGRPwE+FpE/CMwAdzdxXpKkhpmyYDKzB8BW+Ypf4bW+ShJklacK0lIkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkorUzoW6WsTBiUn2HjrOiekZNgwPsWvbZrZvmXfdXEnSMhhQNRycmGTPgaPMnD4DwOT0DHsOHAUwpCSpJof4ath76PjZcJo1c/oMew8d71GNJGlwGFA1nJieWVa5JKl9BlQNG4aHllUuSWqfAVXDrm2bGVq75pyyobVr2LVtc49qJEmDw0kSNcxOhHAWnyStPAOqpu1bRg0kSeoCh/gkSUVqRA/Ki2klqf8MfEB5Ma0k9aeBH+LzYlpJ6k8DH1BeTCtJ/WngA8qLaSWpPw18QHkxrST1p4GfJOHFtJLUnwY+oMCLaSWpHw38EJ8kqT8ZUJKkIi05xBcRbwK+BKwDEtiXmXdFxMXAPcAm4Fngtsx8qXtVHSyubiFJi2unB/Uy8LHMvAK4FvhwRFwB7AYOZ+blwOFqW22YXd1icnqG5JXVLQ5OTPa6apJUjCUDKjNPZuYPqvu/Bo4Bo8AtwP5qt/3A9i7VceC4uoUkLW1Z56AiYhOwBXgcWJeZJ6uHXqA1BDjfMTsjYjwixqempurUdWC4uoUkLa3tgIqI1wFfBz6amb+a+1hmJq3zU6+Smfsycywzx0ZGRmpVdlC4uoUkLa2tgIqItbTC6cuZeaAqfjEi1lePrwdOdaeKg8fVLSRpaUsGVEQEcDdwLDM/O+ehB4Ed1f0dwAMrX73BtH3LKHfceiWjw0MEMDo8xB23XuksPkmaI1qjc4vsEHE98F3gKPD7qviTtM5D3QtcCjxHa5r5Lxf7WWNjYzk+Pl63zpKkARERRzJzbL7HlrwOKjP/E4gFHr6xTsUkSVqIK0lIkopkQEmSimRASZKKZEBJkorUiO+DGjQuNCupCQyoPjO70OzsWn6zC80ChpSkgeIQX59xoVlJTWEPagmlDae50KykprAHtYgSv7fJhWYlNYUBtYgSh9NcaFZSUzjEt4gSh9NmhxdLGnaUpG4woBaxYXiIyXnCqNfDadu3jBpIkgaeQ3yLGMThtIMTk1x358Nctvs/uO7Oh3t6Pk2SFmMPahGDNpzmNVSS+okBtYRBGk5bbNLHoPwfJQ0Oh/gapMRJH5K0EAOqQbyGSlI/MaAaZBAnfUgaXJ6DapBBm/QhabAZUA1TZ9JHaesSShpsBlQXDdIHulPUJa02z0F1SYkLzdZR4rqEkgabAdUlg/aB7hR1SavNgOqSQftAd4q6pNVmQHXJoH2g152i7hqAkpZryYCKiC9GxKmIeHJO2cUR8VBEPFXdXtTdavafQbvmaPuWUe649UpGh4cIYHR4iDtuvbKtCRKDdj5O0uqIzFx8h4i3A78BvpSZf1KV/TPwy8y8MyJ2Axdl5ieWerKxsbEcHx9fgWr3h0GaxVfHdXc+PO/XlowOD/HY7q09qJGkUkTEkcwcm++xJaeZZ+Z3ImLTecW3ADdU9/cDjwJLBlTTeM1RS53zcYPUDpKWp9ProNZl5snq/gvAuoV2jIidwE6ASy+9tMOna5ZBu+ao0y9+HLR2kLQ8tSdJZGuMcMFxwszcl5ljmTk2MjJS9+kaYdCmqHd6Pm7Q2kHS8nTag3oxItZn5smIWA+cWslKNd2gTVHvdA3Auu3g8KDU3zoNqAeBHcCd1e0DK1YjdTwkVrJOzsfVaQeHB6X+1840868C/wVsjojnI+J2WsF0U0Q8Bbyz2tYK8Zqjljrt4PCg1P/amcX3wQUeunGF66JKna/FGKSeQ512cHhQ6n+uZl6oTqeoL9Zz6McP2E7bweFBqf+51NGAGbQJFp3q5fDgoAyxSr1mD2rA1O05DMqwVq+GB+19SSvHgBowu7ZtPucDEtrrOQziB2svhgfrDrEO0h8JUl0O8Q2YThd1ddbbK+oMD65E78tFdaUWe1ADqJOeg7PeXlFneLBXva9Ban9plgElwFlv5+t0eLDTIVbo/I+Euu1vuKlUDvEJcNbbSqnzvVmdfsllnfavO6w4SK+dymMPSoCz3lbSave+6rR/3WFFe27qJgNKZznrrbc6/SOhTvsbbiqZAaXaenHeBQbzQ66TPxLqtH/Twq3E11wL8xyUauvFeRfo3bmX0s671Gn/Ouce67x23Qq3xfTyfFtpvzP9wh6UVkQ/zXqDzv+CL7XX1mn71zn32G89t1722kr8nekHBpR6qlfXHPXbh9zs8YZb5697r4Yk+3E4s5RQNKDUc73offXbh5zh9opOX/denW/rt55+STNrDSj1rSZ9yPVjuC113GqHW6+GJPutp1/SV/YYUOprTfmQ67dw63YoLva6L3RsO6/5QsfWed37radf0lf2GFBqrE7Crd/Ou0Bvwq3UCQlLBdtSz7vQ675YoNY5dqnfmYWO7VWPb6UZUNIy9dN5F+hNuPXjhISljl3odW8n2Do9drHfmcWObef3pU5vcbUmURhQ0ipqSrj144SEUs/3LPQ7s9ixj+3eenafhXptdXqLqzWJwoCS+kQ/hVs/Tkjot/M9Sx272O9Lp73Fdo5dSa4kITXA9i2jPLZ7K/9z53t4bPfWtj9IOl2lolerW/Ti2DoravTq2H6ZRGEPStKi6vTcVrvH14tjezXDrx97qcsVmbniP3QhY2NjOT4+vmrPJ0mroVerNtS5Rm2+cGunl1vn2PlExJHMHJv3MQNKkpqnlKWQDChJUpEWC6hakyQi4uaIOB4RT0fE7jo/S5KkuToOqIhYA/wL8OfAFcAHI+KKlaqYJKnZ6vSgrgGezsxnMvN3wNeAW1amWpKkpqsTUKPAz+dsP1+VnSMidkbEeESMT01N1Xg6SVKTdP1C3czcl5ljmTk2MjLS7aeTJA2IOhfqTgJvmrO9sSpb0JEjR34REc+18bMvAX5Ro25NYTu1x3Zqj+3UHtupPe220x8v9EDH08wj4gLgv4EbaQXT94G/yMwfd/QDz/3Z4wtNO9QrbKf22E7tsZ3aYzu1ZyXaqeMeVGa+HBF/DRwC1gBfXIlwkiQJaq7Fl5nfBL65QnWRJOmsUlcz39frCvQJ26k9tlN7bKf22E7tqd1Oq7rUkSRJ7Sq1ByVJajgDSpJUpOICygVo5xcRX4yIUxHx5JyyiyPioYh4qrq9qJd1LEFEvCkiHomIn0TEjyPiI1W5bTVHRLwmIr4XET+s2ulTVfllEfF49f67JyIu7HVdey0i1kTERER8o9q2jeYREc9GxNGIeCIixquyWu+7ogLKBWgX9e/AzeeV7QYOZ+blwOFqu+leBj6WmVcA1wIfrn6HbKtz/RbYmplXAVcDN0fEtcBngM9l5luAl4Dbe1fFYnwEODZn2zZa2Dsy8+o51z/Vet8VFVC4AO2CMvM7wC/PK74F2F/d3w9sX806lSgzT2bmD6r7v6b1wTKKbXWObPlNtbm2+pfAVuC+qrzx7RQRG4H3AF+otgPbaDlqve9KC6i2FqDVWesy82R1/wVgXS8rU5qI2ARsAR7HtnqVaujqCeAU8BDwM2A6M1+udvH9B58HPg78vtp+A7bRQhL4dkQciYidVVmt912tC3VVjszMiPCagUpEvA74OvDRzPxV6w/fFtuqJTPPAFdHxDBwP/DW3taoLBHxXuBUZh6JiBt6XJ1+cH1mTkbEG4GHIuKncx/s5H1XWg9q2QvQNtyLEbEeoLo91eP6FCEi1tIKpy9n5oGq2LZaQGZOA48AbwOGq3U2wfffdcD7IuJZWqcbtgJ3YRvNKzMnq9tTtP7guYaa77vSAur7wOXVLJkLgQ8AD/a4TiV7ENhR3d8BPNDDuhShOkdwN3AsMz875yHbao6IGKl6TkTEEHATrfN1jwDvr3ZrdDtl5p7M3JiZm2h9Fj2cmR/CNnqViHhtRLx+9j7wLuBJar7viltJIiLeTWvcd3YB2k/3tkZliIivAjfQWsL+ReDvgYPAvcClwHPAbZl5/kSKRomI64HvAkd55bzBJ2mdh7KtKhHxp7ROWq+h9YfqvZn5DxHxZlq9hYuBCeAvM/O3vatpGaohvr/LzPfaRq9Wtcn91eYFwFcy89MR8QZqvO+KCyhJkqC8IT5JkgADSpJUKANKklQkA0qSVCQDSpJUJANKklQkA0qSVKT/B002BwKIALUUAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 504x288 with 1 Axes>"
]
@@ -95,7 +95,7 @@
"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",
+ " N = D-Q\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",
@@ -104,41 +104,18 @@
" \n",
"n = 20\n",
"Q = poisson_mat(n)\n",
- "P = np.arange(1, 10)\n",
+ "P = np.arange(1, 50)\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",
+ " if p in np.arange(1, 10):\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": {
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_3.pdf b/num_ana/sheets/sheet_3.pdf
Binary files differ.
