commit 4ec62c0901f8c9b5479a8a2e196a7420fff5514b
parent 3f3537af772ff2fb6fcb07a2b88f0c757dc1ffb5
Author: miksa234 <milutin@popovic.xyz>
Date: Tue, 5 Apr 2022 15:19:02 +0200
programing p5
Diffstat:
11 files changed, 699 insertions(+), 29 deletions(-)
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/Untitled-checkpoint.ipynb b/num_ana/prog/.ipynb_checkpoints/Untitled-checkpoint.ipynb
@@ -1,6 +0,0 @@
-{
- "cells": [],
- "metadata": {},
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/num_ana/prog/.ipynb_checkpoints/prb5_p-checkpoint.ipynb b/num_ana/prog/.ipynb_checkpoints/prb5_p-checkpoint.ipynb
@@ -0,0 +1,77 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "a365a539",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "0416380f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def damped_richardson(A, b, n, ω):\n",
+ " x_k = np.zeros(b.shape[0])\n",
+ " for i in range(n):\n",
+ " x_k = (np.identity(b.shape[0]) - ω*A)@x_k + ω*b\n",
+ " return x_k"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "301ce785",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "A = np.array([[3, 1, 1],\n",
+ " [1, 3, 1],\n",
+ " [1, 1, 4]])\n",
+ "b = np.array([2, -4, 3])\n",
+ "\n",
+ "# find suitable ω\n",
+ "λ = np.linalg.eig(A)[0]\n",
+ "ω = 2/(max(λ) + min(λ))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "077214f2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "x = damped_richardson(A, b, 20, ω)\n",
+ "np.testing.assert_allclose(A@x, b, rtol=1e-5)"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/num_ana/prog/prb2_p.ipynb b/num_ana/prog/prb2_p.ipynb
@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": 2,
"id": "517aef32",
"metadata": {},
"outputs": [],
@@ -21,7 +21,11 @@
},
{
"cell_type": "code",
+<<<<<<< HEAD
"execution_count": 2,
+=======
+ "execution_count": 3,
+>>>>>>> origin/master
"id": "e14c738f",
"metadata": {},
"outputs": [],
@@ -58,7 +62,11 @@
},
{
"cell_type": "code",
+<<<<<<< HEAD
"execution_count": 6,
+=======
+ "execution_count": 4,
+>>>>>>> origin/master
"id": "caef181e",
"metadata": {},
"outputs": [
@@ -66,6 +74,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
+<<<<<<< HEAD
"1\t44.76606865271526\n",
"2\t59.35975010638207\n",
"3\t22.760834328149066\n",
@@ -96,11 +105,36 @@
"28\t6.194275175956779\n",
"29\t3.495733758960119\n",
"Max. and Min. Singular value are far apart from each other for uneaven p\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"
+>>>>>>> origin/master
]
},
{
"data": {
+<<<<<<< HEAD
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAD4CAYAAABWiRm9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAASiElEQVR4nO3df2xdd3nH8fczNxVWYXO7elGSsqUblRGC0SCrAhWhrqWEAVo9hCrQNgWpUvYHm4qYsib8M5i2EZaNH9IkpoyyBQloqxKSCiRC1R9ik6aCgwuBhqylakWdtDGjFnSyIA3P/vAxczI7vjf2Offrc94vKfK933vt+xwd5X7uec73fG9kJpIklehXBl2AJEnLMaQkScUypCRJxTKkJEnFMqQkScW6pMkXu/LKK3Pr1q1NvqQkqWBHjx79UWaOLvd4oyG1detWJicnm3xJSVLBIuLpCz1uu0+SVCxDSpJULENKklSsnkIqIkYi4t6I+H5EHI+IN0TEFRFxf0Q8Xv28vO5iJUnd0uuR1CeBr2bmK4HXAseB3cADmXkN8EB1X5KkNbPi7L6I+DXgTcB7ATLz58DPI+IW4IbqaQeAh4E76ijyYh2ammbfkROcnJ1j88gwu7aPMbFty6DLkiT1qJcjqauBGeBfI2IqIj4dEZcBGzPzVPWcZ4GNS/1yROyMiMmImJyZmVmbqntwaGqaPQePMT07RwLTs3PsOXiMQ1PTjdUgSVqdXkLqEuB1wKcycxvwP5zX2sv57/tY8js/MnN/Zo5n5vjo6LLXa625fUdOMHfm7Dljc2fOsu/IicZqkCStTi8h9QzwTGY+Ut2/l/nQei4iNgFUP0/XU+LFOTk719e4JKk8K4ZUZj4L/DAixqqhm4DHgPuAHdXYDuBwLRVepM0jw32NS5LK0+vsvj8HPhcR3wGuBf4O2AvcHBGPA2+u7hdj1/YxhjcMnTM2vGGIXdvHlvkNSVJpelq7LzMfBcaXeOimNa1mDS3M4nN2nyStX40uMNu0iW1bDCVJWsdcFkmSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVKxLenlSRDwF/BQ4C7yYmeMRcQVwN7AVeAq4NTOfr6dMSVIX9XMk9XuZeW1mjlf3dwMPZOY1wAPVfUmS1sxq2n23AAeq2weAiVVXI0nSIr2GVAJfi4ijEbGzGtuYmaeq288CG5f6xYjYGRGTETE5MzOzynIlSV3S0zkp4I2ZOR0RvwHcHxHfX/xgZmZE5FK/mJn7gf0A4+PjSz5HkqSl9HQklZnT1c/TwJeA64DnImITQPXzdF1FSpK6acWQiojLIuJlC7eBtwDfBe4DdlRP2wEcrqtISVI39dLu2wh8KSIWnv/5zPxqRHwTuCcibgOeBm6tr0xJUhetGFKZ+STw2iXG/xu4qY6iJEkCV5yQJBXMkJIkFcuQkiQVy5CSJBXLkJIkFcuQkiQVy5CSJBXLkJIkFcuQkiQVy5CSJBWr16/qKMKhqWn2HTnBydk5No8Ms2v7GBPbtgy6LElSTdZNSB2ammbPwWPMnTkLwPTsHHsOHgMwqCSppdZNu2/fkRO/DKgFc2fOsu/IiQFVJEmq27oJqZOzc32NS5LWv3UTUptHhvsalyStf+smpHZtH2N4w9A5Y8Mbhti1fWxAFUmS6rZuJk4sTI5wdp8kdce6CSmYDypDSZK6Y920+yRJ3WNISZKKZUhJkoplSEmSimVISZKKZUhJkoplSEmSimVISZKKZUhJkorVc0hFxFBETEXEl6v7V0fEIxHxRETcHRGX1lemJKmL+jmSuh04vuj+R4GPZ+YrgOeB29ayMEmSegqpiLgKeDvw6ep+ADcC91ZPOQBM1FCfJKnDej2S+gTwl8Avqvu/Dsxm5ovV/WeAJVd+jYidETEZEZMzMzOrqVWS1DErhlREvAM4nZlHL+YFMnN/Zo5n5vjo6OjF/AlJUkf18lUd1wN/EBFvA14C/CrwSWAkIi6pjqauAqbrK1OS1EUrHkll5p7MvCoztwLvBh7MzD8CHgLeVT1tB3C4tiolSZ20muuk7gA+EBFPMH+O6s61KUmSpHl9fTNvZj4MPFzdfhK4bu1LkiRpnitOSJKKZUhJkoplSEmSimVISZKKZUhJkoplSEmSitXXFHQt7dDUNPuOnODk7BybR4bZtX2MiW1LLmUoSeqDIbVKh6am2XPwGHNnzgIwPTvHnoPHAAwqSVol232rtO/IiV8G1IK5M2fZd+TEgCqSpPYwpFbp5OxcX+OSpN4ZUqu0eWS4r3FJUu8MqVXatX2M4Q1D54wNbxhi1/axAVUkSe3hxIlVWpgc4ew+SVp7htQamNi2xVCSpBrY7pMkFcuQkiQVy5CSJBXLkJIkFcuQkiQVy5CSJBXLkJIkFcuQkiQVy5CSJBXLkJIkFcuQkiQVy5CSJBVrxZCKiJdExDci4tsR8b2I+HA1fnVEPBIRT0TE3RFxaf3lSpK6pJcjqZ8BN2bma4FrgbdGxOuBjwIfz8xXAM8Dt9VWpSSpk1YMqZz3QnV3Q/UvgRuBe6vxA8BEHQVKkrqrp3NSETEUEY8Cp4H7gR8As5n5YvWUZ4Alv1ApInZGxGRETM7MzKxByZKkrugppDLzbGZeC1wFXAe8stcXyMz9mTmemeOjo6MXV6UkqZP6mt2XmbPAQ8AbgJGIWPhm36uA6bUtTZLUdb3M7huNiJHq9jBwM3Cc+bB6V/W0HcDhmmqUJHXUJSs/hU3AgYgYYj7U7snML0fEY8BdEfE3wBRwZ411SpI6aMWQyszvANuWGH+S+fNTkiTVwhUnJEnFMqQkScUypCRJxTKkJEnFMqQkScUypCRJxTKkJEnF6uViXhXg0NQ0+46c4OTsHJtHhtm1fYyJbUuu6StJrWFIrQOHpqbZc/AYc2fOAjA9O8eeg8cADCpJrWa7bx3Yd+TELwNqwdyZs+w7cmJAFUlSMwypdeDk7Fxf45LUFobUOrB5ZLivcUlqC0NqHdi1fYzhDUPnjA1vGGLX9rEBVSRJzXDixDqwMDnC2X2SusaQWicmtm0xlCR1ju0+SVKxDClJUrFs953HlR0kqRyG1CKu7CBJZbHdt4grO0hSWQypRVzZQZLKYkgt4soOklQWQ2oRV3aQpLI4cWIRV3aQpLIYUudxZQdJKoftPklSsVYMqYh4eUQ8FBGPRcT3IuL2avyKiLg/Ih6vfl5ef7mSpC7p5UjqReAvMvNVwOuB90XEq4DdwAOZeQ3wQHVf69ihqWmu3/sgV+/+CtfvfZBDU9ODLklSx60YUpl5KjO/Vd3+KXAc2ALcAhyonnYAmKipRjVgYbWN6dk5kv9bbcOgkjRIfZ2TioitwDbgEWBjZp6qHnoW2Li2palJrrYhqUQ9h1REvBT4IvD+zPzJ4scyM4Fc5vd2RsRkREzOzMysqljVx9U2JJWop5CKiA3MB9TnMvNgNfxcRGyqHt8EnF7qdzNzf2aOZ+b46OjoWtSsGrjahqQS9TK7L4A7geOZ+bFFD90H7Khu7wAOr315aoqrbUgqUS8X814P/AlwLCIercY+COwF7omI24CngVtrqVCNcLUNSSWK+dNJzRgfH8/JycnGXk+SVLaIOJqZ48s97ooTkqRiGVKSpGIZUpKkYhlSkqRiGVKSpGL5fVJq3KGpaae6S+qJIaVGLSxku7BO4MJCtoBBJen/sd2nRrmQraR+eCQ1AF1ud7mQraR+eCTVsK5/b5ML2UrqhyHVsK63u1zIVlI/bPc1rOvtLheyldQPQ6phm0eGmV4ikLrU7prYtsVQktQT230Ns90lSb3zSKphtrua0eUZlFKbGFIDYLurXl4wLLWH7T61TtdnUEptYkipdbo+g1JqE0NKreMFw1J7GFJqHWdQSu3hxAm1TlMzKJ1BKNXPkGqxLr+J1j2D0hmEUjNs97VU1xeyrZszCKVmGFIt5ZtovZxBKDXDkGop30Tr5QxCqRmGVEv5JlqvJmcQHpqa5vq9D3L17q9w/d4HbdmqUwyplnIadr0mtm3hI+98DVtGhglgy8gwH3nna2qZQei5RXXZirP7IuIzwDuA05n56mrsCuBuYCvwFHBrZj5fX5nqlwvZ1q+JNRgvdG7Rfaku6GUK+r8B/wR8dtHYbuCBzNwbEbur+3esfXlaDReyXf88t6iuWzGkMvPrEbH1vOFbgBuq2weAhzGkOqnL12I1oakvyXQ/qlQXe05qY2aeqm4/C2xc7okRsTMiJiNicmZm5iJfTiXyfEn9mji36H5UyVY9cSIzE8gLPL4/M8czc3x0dHS1L6eCeC1W/ZqYoOF+VMkudlmk5yJiU2aeiohNwOm1LErrg+dLmlH3ucUm96NtRfXrYo+k7gN2VLd3AIfXphytJ16L1Q5N7UfbiroYK4ZURHwB+E9gLCKeiYjbgL3AzRHxOPDm6r46pqlrsbyYtV5N7UfbiroYvczue88yD920xrVonWniWixXG69fU9fUNdVWtKXYLn5Vh1al7vMlXszajCauqWtiOr0fatrHZZFUNCdntEcTbcWmWoq2oJvjkZSK1tTFrGCbqG5NtBWb+FDj0VqzDCkVbdf2sXPeEKC+yRm+8dSv7rZiEx9qmmxB+8HJdp8K19Rq4848a4cmWopNTgBpYsp+6a1Lj6RUvCZO6jvzrB2aaCk21YJu4ohtPXQQDCkJZ561Sd0fappqQTfxwWk9zJ613SfhzDP1rqkWdBMrgayH2bMeSUk480z9aaIF3cQRW5OzZy+WISVVnHnWH8+v1auJD05NtS5Xw5CSGtLEG0LTM8/qPmLrehDW/cGpqSWxVsOQkhrizLP+NNm67HIYNtG6XA1DSmqQM89611Tr0qPCshlSUos01b5p4oitqdZlm44K2xiEhpTUMs48609bjgrb2h71OilJfWviWqGmvoyxLdcjNXkdXpPfsOyRlKSL0paZZ205KmxTe3QxQ0pSsZpoXbbleqQ2tUcXM6QkdV4bjgqbmtnZ9CoVhpQkNaANQQjNr1JhSElSS7SlPbqYISVJ6kuTq1Q4BV2SVCxDSpJULENKklQsQ0qSVCxDSpJUrMjM5l4sYgZ4+rzhK4EfNVZEWbq87dDt7e/ytkO3t99tP9dvZebocr/QaEgtWUDEZGaOD7SIAenytkO3t7/L2w7d3n63vb9tt90nSSqWISVJKlYJIbV/0AUMUJe3Hbq9/V3eduj29rvtfRj4OSlJkpZTwpGUJElLMqQkScUaWEhFxFsj4kREPBERuwdVx6BExFMRcSwiHo2IyUHXU7eI+ExEnI6I7y4auyIi7o+Ix6uflw+yxross+0fiojpav8/GhFvG2SNdYmIl0fEQxHxWER8LyJur8Zbv+8vsO1d2fcviYhvRMS3q+3/cDV+dUQ8Ur333x0Rl17w7wzinFREDAH/BdwMPAN8E3hPZj7WeDEDEhFPAeOZ2YmL+iLiTcALwGcz89XV2N8DP87MvdUHlcsz845B1lmHZbb9Q8ALmfkPg6ytbhGxCdiUmd+KiJcBR4EJ4L20fN9fYNtvpRv7PoDLMvOFiNgA/AdwO/AB4GBm3hUR/wx8OzM/tdzfGdSR1HXAE5n5ZGb+HLgLuGVAtagBmfl14MfnDd8CHKhuH2D+P3DrLLPtnZCZpzLzW9XtnwLHgS10YN9fYNs7Iee9UN3dUP1L4Ebg3mp8xX0/qJDaAvxw0f1n6NDOqyTwtYg4GhE7B13MgGzMzFPV7WeBjYMsZgD+LCK+U7UDW9fuOl9EbAW2AY/QsX1/3rZDR/Z9RAxFxKPAaeB+4AfAbGa+WD1lxfd+J04Mzhsz83XA7wPvq1pCnZXzfecuXQ/xKeB3gGuBU8A/DrSamkXES4EvAu/PzJ8sfqzt+36Jbe/Mvs/Ms5l5LXAV8x20V/b7NwYVUtPAyxfdv6oa64zMnK5+nga+xPwO7Jrnqr79Qv/+9IDraUxmPlf9B/4F8C+0eP9X5yO+CHwuMw9Ww53Y90tte5f2/YLMnAUeAt4AjETEJdVDK773DyqkvglcU83yuBR4N3DfgGppXERcVp1IJSIuA94CfPfCv9VK9wE7qts7gMMDrKVRC2/QlT+kpfu/Onl+J3A8Mz+26KHW7/vltr1D+340Ikaq28PMT5Q7znxYvat62or7fmArTlTTLj8BDAGfycy/HUghAxARv8380RPAJcDn2779EfEF4Abml+p/Dvgr4BBwD/CbzH+Fy62Z2boJBsts+w3Mt3sSeAr400XnaFojIt4I/DtwDPhFNfxB5s/NtHrfX2Db30M39v3vMj8xYoj5A6J7MvOvq/e/u4ArgCngjzPzZ8v+HZdFkiSVyokTkqRiGVKSpGIZUpKkYhlSkqRiGVKSpGIZUpKkYhlSkqRi/S9bhJKASc+kxQAAAABJRU5ErkJggg==\n",
+=======
+ "text/plain": [
+ "<matplotlib.collections.PathCollection at 0x7fde841e3a00>"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "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",
+>>>>>>> origin/master
"text/plain": [
"<Figure size 504x288 with 1 Axes>"
]
@@ -115,7 +149,11 @@
"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",
+<<<<<<< HEAD
" N = D - Q\n",
+=======
+ " N = D-Q\n",
+>>>>>>> origin/master
" 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",
@@ -124,38 +162,26 @@
" \n",
"n = 20\n",
"Q = poisson_mat(n)\n",
+<<<<<<< HEAD
"P = np.arange(1, 30)\n",
+=======
+ "P = np.arange(1, 50)\n",
+>>>>>>> origin/master
"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\")"
+ "plt.scatter(P, cond_2)"
]
},
{
"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",
+ "id": "2469e597",
"metadata": {},
"outputs": [],
"source": []
@@ -177,7 +203,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.10.2"
+ "version": "3.10.3"
}
},
"nbformat": 4,
diff --git a/num_ana/prog/prb4_p.ipynb b/num_ana/prog/prb4_p.ipynb
@@ -0,0 +1,181 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "84bcf23b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from numpy import testing as testing\n",
+ "from matplotlib import pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f8656109",
+ "metadata": {},
+ "source": [
+ "# Sheet 4, Exercise 1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "0de86192",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def conjugate_gradient(A:np.ndarray, b:np.ndarray, steps: int) -> np.ndarray:\n",
+ " x_k = np.zeros(b.shape) \n",
+ " r_k = b - (A @ x_k)\n",
+ " p_k = r_k\n",
+ " for k in range(steps):\n",
+ " α_k = (p_k @ r_k) / (p_k @ A @ p_k)\n",
+ " x_k = x_k + (α_k * p_k) # x_k+1\n",
+ " r_k = r_k - α_k * (A @ p_k) # r_k+1\n",
+ " if not np.any(r_k): # stop if r_k+1 = 0 \n",
+ " break\n",
+ " β_k = ((A @ p_k) @ r_k) / ((A @ p_k) @ p_k)\n",
+ " p_k = r_k - (β_k * p_k)\n",
+ " return x_k\n",
+ "\n",
+ "def poisson_mat(n:int, m : int =None) -> np.ndarray:\n",
+ " return 2 * np.eye(n, m) + (-1) * np.eye(n, m, k=1) + (-1) * np.eye(n, m, k=-1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "60a4d70d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for n in range(3, 21):\n",
+ " Q = poisson_mat(n)\n",
+ " b = np.ones(n)\n",
+ " x = np.linalg.inv(Q) @ b\n",
+ " x_k = conjugate_gradient(Q, b, 10)\n",
+ " testing.assert_allclose(x_k, x, rtol=1e-5, err_msg=f'CG failed at dim - {n}')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "67accdb6",
+ "metadata": {},
+ "source": [
+ "What happends to a matrix that is not positive definite? Consider the system\n",
+ "\\begin{align}\n",
+ " A = \n",
+ " \\begin{pmatrix}\n",
+ " 1 & 0 & 0\\\\\n",
+ " 0 & 1 & 0\\\\\n",
+ " 0 & 0 & -2\\\\\n",
+ " \\end{pmatrix}, \\qquad\n",
+ " b = \n",
+ " \\begin{pmatrix}\n",
+ " 1 \\\\\n",
+ " 1 \\\\\n",
+ " 1\n",
+ " \\end{pmatrix}.\n",
+ "\\end{align}\n",
+ "Because $A$ is indefinite, the coefficients $\\alpha_k$ are undefined, \n",
+ "i.e. there is an $r^k \\neq 0$ such that $(r^k)^TAr^k = 0$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "12b54df7",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/tmp/ipykernel_219747/1775264609.py:6: RuntimeWarning: divide by zero encountered in double_scalars\n",
+ " α_k = (p_k @ r_k) / (p_k @ A @ p_k)\n",
+ "/tmp/ipykernel_219747/1775264609.py:12: RuntimeWarning: invalid value encountered in subtract\n",
+ " p_k = r_k - (β_k * p_k)\n",
+ "/tmp/ipykernel_219747/1775264609.py:6: RuntimeWarning: invalid value encountered in matmul\n",
+ " α_k = (p_k @ r_k) / (p_k @ A @ p_k)\n",
+ "/tmp/ipykernel_219747/1775264609.py:8: RuntimeWarning: invalid value encountered in matmul\n",
+ " r_k = r_k - α_k * (A @ p_k) # r_k+1\n",
+ "/tmp/ipykernel_219747/1775264609.py:11: RuntimeWarning: invalid value encountered in matmul\n",
+ " β_k = ((A @ p_k) @ r_k) / ((A @ p_k) @ p_k)\n"
+ ]
+ },
+ {
+ "ename": "AssertionError",
+ "evalue": "\nNot equal to tolerance rtol=1e-05, atol=0\nCG failed for non-definite matrix\nx and y nan location mismatch:\n x: array([nan, nan, nan])\n y: array([1., 1., 1.])",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)",
+ "Input \u001b[0;32mIn [4]\u001b[0m, in \u001b[0;36m<cell line: 4>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 2\u001b[0m b \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mones(\u001b[38;5;241m3\u001b[39m)\n\u001b[1;32m 3\u001b[0m x_k \u001b[38;5;241m=\u001b[39m conjugate_gradient(A, b , \u001b[38;5;241m10\u001b[39m)\n\u001b[0;32m----> 4\u001b[0m \u001b[43mtesting\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43massert_allclose\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx_k\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrtol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1e-5\u001b[39;49m\u001b[43m \u001b[49m\u001b[43m,\u001b[49m\u001b[43merr_msg\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43mf\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mCG failed for non-definite matrix\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n",
+ " \u001b[0;31m[... skipping hidden 2 frame]\u001b[0m\n",
+ "File \u001b[0;32m~/.local/lib/python3.10/site-packages/numpy/testing/_private/utils.py:745\u001b[0m, in \u001b[0;36massert_array_compare.<locals>.func_assert_same_pos\u001b[0;34m(x, y, func, hasval)\u001b[0m\n\u001b[1;32m 740\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m bool_(x_id \u001b[38;5;241m==\u001b[39m y_id)\u001b[38;5;241m.\u001b[39mall() \u001b[38;5;241m!=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m:\n\u001b[1;32m 741\u001b[0m msg \u001b[38;5;241m=\u001b[39m build_err_msg([x, y],\n\u001b[1;32m 742\u001b[0m err_msg \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mx and y \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m location mismatch:\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 743\u001b[0m \u001b[38;5;241m%\u001b[39m (hasval), verbose\u001b[38;5;241m=\u001b[39mverbose, header\u001b[38;5;241m=\u001b[39mheader,\n\u001b[1;32m 744\u001b[0m names\u001b[38;5;241m=\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mx\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124my\u001b[39m\u001b[38;5;124m'\u001b[39m), precision\u001b[38;5;241m=\u001b[39mprecision)\n\u001b[0;32m--> 745\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mAssertionError\u001b[39;00m(msg)\n\u001b[1;32m 746\u001b[0m \u001b[38;5;66;03m# If there is a scalar, then here we know the array has the same\u001b[39;00m\n\u001b[1;32m 747\u001b[0m \u001b[38;5;66;03m# flag as it everywhere, so we should return the scalar flag.\u001b[39;00m\n\u001b[1;32m 748\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(x_id, \u001b[38;5;28mbool\u001b[39m) \u001b[38;5;129;01mor\u001b[39;00m x_id\u001b[38;5;241m.\u001b[39mndim \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n",
+ "\u001b[0;31mAssertionError\u001b[0m: \nNot equal to tolerance rtol=1e-05, atol=0\nCG failed for non-definite matrix\nx and y nan location mismatch:\n x: array([nan, nan, nan])\n y: array([1., 1., 1.])"
+ ]
+ }
+ ],
+ "source": [
+ "A = np.array([[1, 0 , 0], [0, 1, 0], [0, 0, -2]])\n",
+ "b = np.ones(3)\n",
+ "x_k = conjugate_gradient(A, b , 10)\n",
+ "testing.assert_allclose(x_k, b, rtol=1e-5 ,err_msg=f'CG failed for non-definite matrix')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5c5e809e",
+ "metadata": {},
+ "source": [
+ "# Sheet 4, Exercise 4 1) check"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "20f80867",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[3. 2. 1.]\n"
+ ]
+ }
+ ],
+ "source": [
+ "A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])\n",
+ "b = np.array([4, 0, 0])\n",
+ "x = conjugate_gradient(A, b, 10)\n",
+ "print(x)"
+ ]
+ }
+ ],
+ "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.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/num_ana/prog/prb4_p.pdf b/num_ana/prog/prb4_p.pdf
Binary files differ.
diff --git a/num_ana/prog/prb5_p.ipynb b/num_ana/prog/prb5_p.ipynb
@@ -0,0 +1,77 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "a365a539",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "0416380f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def damped_richardson(A, b, n, ω):\n",
+ " x_k = np.zeros(b.shape[0])\n",
+ " for i in range(n):\n",
+ " x_k = (np.identity(b.shape[0]) - ω*A)@x_k + ω*b\n",
+ " return x_k"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "301ce785",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "A = np.array([[3, 1, 1],\n",
+ " [1, 3, 1],\n",
+ " [1, 1, 4]])\n",
+ "b = np.array([2, -4, 3])\n",
+ "\n",
+ "# find suitable ω\n",
+ "λ = np.linalg.eig(A)[0]\n",
+ "ω = 2/(max(λ) + min(λ))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "077214f2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "x = damped_richardson(A, b, 20, ω)\n",
+ "np.testing.assert_allclose(A@x, b, rtol=1e-5)"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/num_ana/sheets/sheet_3.pdf b/num_ana/sheets/sheet_3.pdf
Binary files differ.
diff --git a/num_ana/sheets/sheet_4.pdf b/num_ana/sheets/sheet_4.pdf
Binary files differ.
