From d876321aac6b5c5f7fa8b1988e9e8a2144c46220 Mon Sep 17 00:00:00 2001
From: Frederic Weidling <fweidli@client61.num.math.uni-goettingen.de>
Date: Fri, 23 Dec 2016 08:52:32 +0100
Subject: [PATCH] blatt9
---
09_Iterative.ipynb | 280 ++++++++++++++++++++++++++++++++++++
10_Poisson1d_function.ipynb | 168 ++++++++++++++++++++++
2 files changed, 448 insertions(+)
create mode 100644 09_Iterative.ipynb
create mode 100644 10_Poisson1d_function.ipynb
diff --git a/09_Iterative.ipynb b/09_Iterative.ipynb
new file mode 100644
index 0000000..374da5d
--- /dev/null
+++ b/09_Iterative.ipynb
@@ -0,0 +1,280 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Aufgabe 34\n",
+ "## Teil 1\n",
+ "Implementation der Jacobi Iteration:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "from numpy import array, diag\n",
+ "from numpy.linalg import norm\n",
+ "\n",
+ "def jacobi_iteration(A,b,x0,tol):\n",
+ " d=diag(A)\n",
+ " B=-(A-diag(d))\n",
+ " d=d**(-1)\n",
+ " k=1\n",
+ " x=d*(B@x0+b)\n",
+ " x_old=x0\n",
+ " q=0\n",
+ " \n",
+ " while norm(x-x_old)>tol:\n",
+ " k=k+1\n",
+ " x_old=x.copy()\n",
+ " x=d*(B@x+b)\n",
+ "\n",
+ " return x,k"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Wir testen an einem Beispiel mit starken Zeilensummenkriterium (d.h. wir müssen die $\\infty$-Norm verwenden):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Startwert: [ 2.69261466 -17.14628008 2.2221343 ],\n",
+ "Berechnete Lösung: [-0.07144352 0.57141562 -0.64283772],\n",
+ "Wahre Lösung: [-0.07142857 0.57142857 -0.64285714],\n",
+ "Differenz: 1.692e-04,\n",
+ "Iterationen: 24\n"
+ ]
+ }
+ ],
+ "source": [
+ "from numpy.random import randn\n",
+ "from numpy.linalg import solve\n",
+ "\n",
+ "A=array([[3.,1,-1],[3,5,1],[2,1,-4]])\n",
+ "b=array([1.,2,3])\n",
+ "xtrue=solve(A,b)\n",
+ "x0=randn(3)*10\n",
+ "x,k = jacobi_iteration(A,b,x0,1e-4)\n",
+ "\n",
+ "print('Startwert: {0},\\nBerechnete Lösung: {1},\\nWahre Lösung: {2},\\nDifferenz: {3:1.3e},\\nIterationen: {4}'.format(x0,x,xtrue,norm(A@x-b),k))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Teil 2\n",
+ "Implementation des Gauß-Seidel Verfahren:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "from numpy import tril, triu\n",
+ "\n",
+ "def gauss_seidel_iteration(A,b,x0,tol):\n",
+ " k=1\n",
+ " x=solve(tril(A),-triu(A,1)@x0+b)\n",
+ " x_old=x0\n",
+ " q=0\n",
+ " \n",
+ " while norm(x-x_old)>tol:\n",
+ " k=k+1\n",
+ " x_old=x.copy()\n",
+ " x=solve(tril(A),-triu(A,1)@x+b)\n",
+ " \n",
+ " return x,k"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Wir testen am selben Beispiel wie oben:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Startwert: [ 2.69261466 -17.14628008 2.2221343 ],\n",
+ "Berechnete Lösung: [-0.07139824 0.57140518 -0.64284782],\n",
+ "Wahre Lösung: [-0.07142857 0.57142857 -0.64285714],\n",
+ "Differenz: 6.062e-05,\n",
+ "Iterationen: 13\n"
+ ]
+ }
+ ],
+ "source": [
+ "x,k = gauss_seidel_iteration(A,b,x0,1e-4)\n",
+ "\n",
+ "print('Startwert: {0},\\nBerechnete Lösung: {1},\\nWahre Lösung: {2},\\nDifferenz: {3:1.3e},\\nIterationen: {4}'.format(x0,x,xtrue,norm(A@x-b),k))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Teil 3\n",
+ "Anwendung auf die Poissongleichung. Zunächst Import der benötigten Funktionen."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "from numpy import eye\n",
+ "\n",
+ "def tridiag(a,b,c,n):\n",
+ " A=zeros([n,n]) \n",
+ " for iii in range(n):\n",
+ " if iii==0:\n",
+ " A[iii,iii]=b \n",
+ " A[iii,iii+1]=c\n",
+ " elif iii==n-1:\n",
+ " A[iii,iii]=b\n",
+ " A[iii,iii-1]=a\n",
+ " else:\n",
+ " A[iii,iii]=b\n",
+ " A[iii,iii-1]=a\n",
+ " A[iii,iii+1]=c\n",
+ " return A\n",
+ "\n",
+ "def poisson2d_matrix(n):\n",
+ " A=zeros([(n-1)**2,(n-1)**2]) # Hier speichern wie die Poissonmatrix rein\n",
+ " B=tridiag(-1,4,-1,n-1) # Die Matrix entlang der Diagonalen\n",
+ " C=-eye(n-1) # Die Matrix auf den Matrixnebendiagonalen\n",
+ " # Jetzt bauen wir Matrix nach dem selben Prinzip zusammen, wie in tridiag:\n",
+ " for iii in range(n-1): # für \"blockweise\" durchlaufen der Matrix\n",
+ " ind=iii*(n-1) # Umrechnen von Blockindex zu Matrixindex\n",
+ " A[ind:ind+(n-1),ind:ind+n-1]=B #Setzen des Diagonalblocks\n",
+ " # Setzen der Nebendiagonalblöcke\n",
+ " if iii==0:\n",
+ " A[ind:ind+(n-1),ind+(n-1):ind+2*(n-1)]=C\n",
+ " elif iii==n-2:\n",
+ " A[ind:ind+(n-1),ind-(n-1):ind]=C\n",
+ " else:\n",
+ " A[ind:ind+(n-1),ind+(n-1):ind+2*(n-1)]=C\n",
+ " A[ind:ind+(n-1),ind-(n-1):ind]=C\n",
+ " return n**2*A\n",
+ "\n",
+ "from numpy import array,zeros,ogrid,broadcast_arrays\n",
+ "\n",
+ "def rechte_seite(n):\n",
+ " rechteSeite=zeros([n-1,n-1])\n",
+ " x,y=ogrid[1/n:1-1/n:(n-1)*1j,1/n:1-1/n:(n-1)*1j]\n",
+ " x, y = broadcast_arrays(x, y)\n",
+ " rechteSeite[(1/8<=x) & (x<=1/4) & (1/8<=y) & (y<=5/8)]=2\n",
+ " rechteSeite[(1/8<=y) & (y<=1/4) & (1/8<=x) & (x<=5/8)]=2\n",
+ " z=(x-2/3)**2+(y-2/3)**2\n",
+ " rechteSeite[z<=1/16]=-1\n",
+ " return rechteSeite.ravel()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "N=[5,10,15,20,25,30,35,40]\n",
+ "k_jacobi=[]\n",
+ "k_gauss_seidel=[]\n",
+ "\n",
+ "for n in N:\n",
+ " A=poisson2d_matrix(n)\n",
+ " b=rechte_seite(n)\n",
+ " x0=zeros(b.shape)\n",
+ " tol=0.001/(n**(1/2))\n",
+ " \n",
+ " x,k=jacobi_iteration(A,b,x0,tol)\n",
+ " k_jacobi.append(k)\n",
+ " \n",
+ " x,k=gauss_seidel_iteration(A,b,x0,tol)\n",
+ " k_gauss_seidel.append(k)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg4AAAFkCAYAAABIPLOYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XmcjvX+x/HXh+yFolCyRKFVMxEVibKlg6Ni4shS0pFq\nWn7aS5tSDVpUWiQxHelIdSzJWtbM0CIUWSqMJcbOLN/fH99bjcky+3XPzPv5eNwPM9d13df9ue9z\ndL99V3POISIiIpIRRYIuQERERPIPBQcRERHJMAUHERERyTAFBxEREckwBQcRERHJMAUHERERyTAF\nBxEREckwBQcRERHJMAUHERERyTAFBxEREcmwbAUHM3vAzFLNLCbNsZGhY2kfk9I9r4SZvWZmW81s\nl5mNN7PTslOLiIiI5L4sBwczawD0Ab49wunJQCWgcugRle78UOBaoBPQFDgd+DirtYiIiEjeyFJw\nMLMTgQ+AW4AdR7jkgHNui3Nuc+iRmOa5ZYFeQLRzbrZzbgnQE7jczBpmpR4RERHJG1ltcXgN+Mw5\nN+Mo55uZWYKZrTCz4WZ2SppzkcAJwPRDB5xzK4H1QOMs1iMiIiJ54ITMPsHMugD1gUuOcslkfLfD\nGqAWMAiYZGaNnd/DuzJw0Dm3M93zEkLnjvSaFYBWwFpgf2ZrFhERKcRKAjWAqc65bdm9WaaCg5lV\nxY9PuNo5l3Ska5xz49L8uszMvgdWA82AmVmssxUwJovPFREREegKjM3uTTLb4hAJnArEm5mFjhUF\nmprZHUCJUKvCn5xza8xsK1AbHxw2AcXNrGy6VodKoXNHshbggw8+oF69epksueCJjo5myJAhQZcR\nOH0Onj6Hv+iz8PQ5/EWfBSxfvpxu3bpB6Ls0uzIbHL4ELkh37D1gOfBc+tAAf7ZSVAA2hg7FAclA\nC2BC6Jo6QDVg/lFedz9AvXr1iIiIyGTJBU+5cuX0OaDP4RB9Dn/RZ+Hpc/iLPovD5EhXf6aCg3Nu\nD/Bj2mNmtgfY5pxbbmZlgMfxYxw24VsZngd+AqaG7rHTzN4BYsxsO7ALeBmY65xblM33IyIiIrko\n04MjjyBtK0MKcCHQHSgPbMAHhsfSjYmIDl07HigBTAH65UAtIiIikouyHRycc83T/LwfaJ2B5xwA\n+oceIiIikk9or4p8KCoq/UKchZM+B0+fw1/0WXj6HP6izyLn2RHGM4YdM4sA4uLi4jTIRUREJBPi\n4+OJjIwEiHTOxWf3fmpxEBERkQxTcBAREZEMU3AQERGRDFNwEBERkQxTcBAREZEMU3AQERGRDFNw\nEBERkQxTcBAREZEMU3AQERGRDFNwEBERkQxTcBAREZEMU3AQEREpwHbuzNn7KTiIiIgUUAsXQk5v\nEKrgICIiUsA4BzExcMUVcOqpOXtvBQcREZEC5I8/oH17uPdeiI6Gt9/O2fufkLO3ExERkaDMmwdd\nusCePfD553DttRAfn7OvoRYHERGRfC41FQYPhqZNoVo1WLrUh4bcoOAgIiKSj23dCu3awYABcP/9\nMHMmnHlm7r2euipERETyqa+/9l0TBw7A5MnQunXuv6ZaHERERPKZ1FQYNAiaNYNatXzXRF6EBlBw\nEBERyVc2b4a2beHhh+HBB2H6dDjjjLx7/WwFBzN7wMxSzSwm3fEnzWyDme01s2lmVjvd+RJm9pqZ\nbTWzXWY23sxOy04tIiIiBd3s2VC/PixZAlOnwlNPwQl5POggy8HBzBoAfYBv0x0fANwROtcQ2ANM\nNbPiaS4bClwLdAKaAqcDH2e1FhERkYIsJcWHhObNoW5d3zVxzTXB1JKl4GBmJwIfALcAO9Kdvgt4\nyjn3uXPuB6A7Phh0CD23LNALiHbOzXbOLQF6ApebWcOsvQ0REZGCKSEBWrWCxx+HRx+FadOgSpXg\n6slqi8NrwGfOuRlpD5pZTaAyMP3QMefcTmAh0Dh06BL8bI6016wE1qe5RkREpNCbMQMuugh++AG+\n/BKeeAKKFg22pkwHBzPrAtQHHjzC6cqAAxLSHU8InQOoBBwMBYqjXSMiIlJopaT4Foarr4bzz/dd\nE82bB12Vl6khFWZWFT8+4WrnXFLulCQiIlJ4bdwIXbv6gZADB8JDDwXfypBWZsdiRgKnAvFmZqFj\nRYGmZnYHUBcwfKtC2laHSsCS0M+bgOJmVjZdq0Ol0Lmjio6Oply5cocdi4qKIiqn9wwVEREJwLRp\n0K2bDwozZsCVV2bu+bGxscTGxh52LDExMQcrBHPOZfxiszJA9XSH3wOWA88555ab2QbgBefckNBz\nyuJDRHfn3Eeh37cAXZxzE0LX1Ando5FzbtERXjcCiIuLiyMiIiKz71FERCSsJSf78QvPPutnS4we\nDafl0CIF8fHxREZGAkQ657K95VWmWhycc3uAH9MeM7M9wDbn3PLQoaHAI2a2ClgLPAX8BkwM3WOn\nmb0DxJjZdmAX8DIw90ihQUREpCD7/Xe46SaYOxeeecbvOVEkjJdnzIllIw5rsnDODTaz0sCbQHng\nK6CNc+5gmsuigRRgPFACmAL0y4FaRERE8o0pU+Bf/4ISJWDWLLjiiqArOr5sBwfn3N/GeTrnngCe\nOMZzDgD9Qw8REZFCJTnZr8nw3HPQpg28/z5UrBh0VRmj3TFFRETy0K+/QlQULFgAzz8P990X3l0T\n6Sk4iIiI5JH//Q+6d4cyZWDOHLjssqAryrx8lHFERETyp6Qk+L//g3bt4PLL/SZV+TE0gFocRERE\nctX69dClC3zzDbz0EkRHw58rIeVDCg4iIiK55NNPoUcPKFsWvv4aLr006IqyT10VIiIiOezgQbjn\nHmjf3q/+uGRJwQgNoBYHERGRHLVmje+aWLIEhg2D/v3zd9dEegoOIiIiOWTCBOjZE045xa8E2aBB\n0BXlPHVViIiIZNOBA3DnnfDPf/qtsOPjC2ZoALU4iIiIZMvq1dC5M3z/Pbz6Kvz73wWrayI9BQcR\nEZEsGj8eeveGU0+F+fOhMGzgrK4KERGRTNq/H/r1gxtugNatfddEYQgNoBYHERGRTPn5Z9818eOP\n8PrrcNttBbtrIj21OIiIiGTQhx9CZCTs3u03qerbt3CFBlBwEBEROa59+3xIiIry+03ExUH9+kFX\nFQx1VYiIiBzDypVw443w008wYgTcckvha2VISy0OIiIiRzFmjO+aOHAAFi6EW28t3KEBFBxERET+\nZu9e37LQrZtf1GnxYrjwwqCrCg/qqhAREUlj+XLfNbF6Nbz7rt/dsrC3MqSlFgcREZGQ99+HSy6B\n1FT45hu/74RCw+EUHEREpNDbs8eHhJtv9ms0LFoE550XdFXhSV0VIiJSqC1b5rsm1q6FUaOge/eg\nKwpvanEQEZFCyTk/hqFBAyhSxA+AVGg4PgUHEREpdHbv9iGhd2/o2tVPtaxXL+iq8odMBQcz62tm\n35pZYugxz8xapzk/0sxS0z0mpbtHCTN7zcy2mtkuMxtvZqfl1BsSERE5lu++8wMgP/nEr9Pw1ltQ\nunTQVeUfmW1x+BUYAEQAkcAMYKKZpc1pk4FKQOXQIyrdPYYC1wKdgKbA6cDHma5cREQkE5zzIeHS\nS6FkSb9s9E03BV1V/pOpwZHOuf+lO/SImd0ONAKWh44dcM5tOdLzzaws0Avo4pybHTrWE1huZg2d\nc4syVb2IiEgG7Nrld7GMjfV7TsTEQKlSQVeVP2V5jIOZFTGzLkBpYF6aU83MLMHMVpjZcDM7Jc25\nSHxYmX7ogHNuJbAeaJzVWkRERI5m6VK/bPTnn/vdLV9/XaEhOzIdHMzsfDPbBRwAhgMdQ1/+4Lsp\nugPNgf8DrgQmmf25fEZl4KBzbme62yaEzomIiOQI53xIaNQITjwR4uP9Gg2SPVlZx2EFcBFQDrge\neN/MmjrnVjjnxqW5bpmZfQ+sBpoBM7NbbHR0NOXKlTvsWFRUFFFR6YdRiIhIYZaU5JeKHjsW+vWD\nF1/04xoKutjYWGJjYw87lpiYmKOvYc657N3AbBqwyjl3+1HObwYeds69ZWZXAV8CJ6dtdTCztcAQ\n59ywo9wjAoiLi4sjIiIiW/WKiEjBlpzsBz1OnAijR/vFnQqz+Ph4IiMjASKdc/HZvV9OrONQBChx\npBNmVhWoAGwMHYoDkoEWaa6pA1QD5udALSIiUoilpkKvXjBhAowbp9CQGzLVVWFmz+LHMawHTgK6\n4scxtDSzMsDj+KmVm4DawPPAT8BUAOfcTjN7B4gxs+3ALuBlYK5mVIiISHY452dMjBnjZ0+0bx90\nRQVTZsc4nAaMAqoAicB3QEvn3AwzKwlciB8cWR7YgA8MjznnktLcIxpIAcbjWyqmAP2y8yZERKRw\ncw7uugvefhvee08tDbkps+s43HKMc/uB1kc7n+a6A0D/0ENERCRbnIMHHoBXXoE339R+E7lNe1WI\niEi+9uSTMHgwDB0KffoEXU3Bp+AgIiL51uDB8MQTMGiQ76qQ3KfgICIi+dLLL8OAAfDYY76rQvKG\ngoOIiOQ7I0b4Fob77/ctDpJ3FBxERCRfGT3aT7vs3x+efx7+3NRA8oSCg4iI5BvjxvmlpHv39oMh\nFRrynoKDiIjkCxMnQteufjnpN96AIvoGC4Q+dhERCXtTp/pFnTp0gJEjoWjRoCsqvBQcREQkrM2a\n5QNDq1Z+OekTsrKvs+QYBQcREQlb8+ZBu3bQpIkf31C8eNAViYKDiIiEpcWLoU0buOQS+OQTKFky\n6IoEFBxERCQMffcdtGwJ550Hn30GpUsHXZEcouAgIiJhZflyuPpqqFkTJk2Ck04KuiJJS8FBRETC\nxqpV0KIFVKoEX3wB5csHXZGkp+AgIiJhYd06HxrKloUvv4QKFYKuSI5EwUFERAL3++/QvLmfajl9\num9xkPCk2bAiIhKozZv9mIakJPjqKzjjjKArkmNRcBARkcBs2+ZDQ2IizJkD1asHXZEcj4KDiIgE\nYscOvxrkpk0wezbUrh10RZIRCg4iIpLndu2Ctm3hl19g5kyoVy/oiiSjFBxERCRP7d0L110Hy5b5\n2RMXXRR0RZIZCg4iIpJn9u+Hjh39ctJTp0KDBkFXJJml4CAiInkiKclvjT1njl8R8vLLg65IsiJT\n6ziYWV8z+9bMEkOPeWbWOt01T5rZBjPba2bTzKx2uvMlzOw1M9tqZrvMbLyZnZYTb0ZERMJTcjJ0\n7epbGSZMgKuuCroiyarMLgD1KzAAiAAigRnARDOrB2BmA4A7gD5AQ2APMNXM0m6EOhS4FugENAVO\nBz7OxnsQEZEwlpoKvXr5wDBuHLRuffznSPjKVFeFc+5/6Q49Yma3A42A5cBdwFPOuc8BzKw7kAB0\nAMaZWVmgF9DFOTc7dE1PYLmZNXTOLcrWuxERkbDiHPTtC2PGwNix0L590BVJdmV5yWkzK2JmXYDS\nwDwzqwlUBqYfusY5txNYCDQOHboEH1bSXrMSWJ/mGhERKQCcg7vugrfegnffhc6dg65IckKmB0ea\n2fnAfKAksAvo6JxbaWaNAYdvYUgrAR8oACoBB0OB4mjXiIhIPuccPPAAvPIKvPEG3Hxz0BVJTsnK\nrIoVwEVAOeB64H0za5qjVYmISL42cCAMHgxDh8JttwVdjeSkTAcH51wy8Evo1yVm1hA/tmEwYPhW\nhbStDpWAJaGfNwHFzaxsulaHSqFzxxQdHU25cuUOOxYVFUVUVFRm34aIiOSS55/3wWHQIN9VIXkn\nNjaW2NjYw44lJibm6GuYcy57NzCbDqxzzvUysw3AC865IaFzZfEhortz7qPQ71vwgyMnhK6pgx9Y\n2ehogyPNLAKIi4uLIyIiIlv1iohI7nn5ZR8WHnvMhwcJXnx8PJGRkQCRzrn47N4vUy0OZvYsMBk/\nmPEkoCtwJdAydMlQ/EyLVcBa4CngN2Ai+MGSZvYOEGNm2/FjJF4G5mpGhYhI/jZihA8N998PTzwR\ndDWSWzLbVXEaMAqoAiQC3wEtnXMzAJxzg82sNPAmUB74CmjjnDuY5h7RQAowHigBTAH6ZedNiIhI\nsEaP9tMu77jDd1WYBV2R5JbMruNwSwaueQJ44hjnDwD9Qw8REcnnxo2DHj2gd28YNkyhoaDL8joO\nIiIiEyf6paSjovy0yyL6Vinw9D+xiIhkydSpftOqDh3gvfegaNGgK5K8oOAgIiKZNnOmDwytWvnl\npE/QXsuFhoKDiIhkyrx5cN110KSJH99QvPjxnyMFh4KDiIhk2OLF0KYNREbCJ59AyZJBVyR5TcFB\nREQy5LvvoGVLOPdc+PxzKF066IokCAoOIiJyXMuXw9VXQ82aMHkynHRS0BVJUBQcRETkmFatghYt\noFIl+OILKF8+6IokSAoOIiJyVOvW+dBQtix8+SVUqBB0RRI0BQcRETmi33+H5s39VMvp032Lg4hm\n3oqIyN8kJPiWhqQkmDMHzjgj6IokXCg4iIjIYbZtg2uugZ07YfZsqFEj6IoknCg4iIjIn3bs8KtB\nbtzoQ8PZZwddkYQbBQcREQFg1y5o2xZ++cUvKX3uuUFXJOFIwUFERNi71y8jvWyZnz1x0UVBVyTh\nSsFBRKSQ278fOnaEb77x6zQ0aBB0RRLOFBxERAqxgwf91thz5sD//geXXx50RRLuFBxERAqp5GTo\n2hWmTIFPP/VrNogcj4KDiEghlJoKPXvChAkwfjy0bh10RZJfKDiIiBQyzkHfvjB2rH906BB0RZKf\nKDiIiBQizsFdd8Fbb8F770HnzkFXJPmN9qoQESkknIMHHoBXXoE33oCbbw66IsmPFBxERAqJgQNh\n8GAYMgRuuy3oaiS/UnAQESkEnn/eB4dnn4W77w66GsnPMjXGwcweBDoCdYF9wDxggHPupzTXjATS\nN4BNcc61TXNNCSAG6AyUAKYC/3bObc7KmxARkaMbNsx3UTz6KDz4YNDVSFalulT2J+9nz8E97Ena\nw96kvX/7eW/SXvYk7Tns53XL1+VoHZkdHNkEeAVYHHruIOALM6vnnNuX5rrJQA/AQr8fSHefoUAb\noBOwE3gN+Dh0fxERySEjRvgWhvvu8y0OkntSXSp7k/Ye80v8SF/2f/58nPN7k/ZmqI4TipxAmWJl\nKFO8DKWLlcY22fGflAmZCg5pWw0AzKwHsBmIBL5Oc+qAc27Lke5hZmWBXkAX59zs0LGewHIza+ic\nW5SZmkRE5MhGj/bTLu+4w49tsJz9/sh3UlJT/vpiP96XePov/uTj/+t+X/K+4xcBFC9anNLFSh/2\n5V6mWOjP4mU4uezJxzyf9ucjXVesaLHDXi8+Pp7IQZE59jlmdzpmecABf6Q73szMEoDtwAzgEefc\noWsiQ687/dDFzrmVZrYeaAwoOIiIZNO4cdCjB/Tq5bsqClto2LBrAz0+6cGKrSv+/LLfn7w/Q88t\nUbTEUb+USxcrTcXSFf/68j7KF/qxvtxPKJK/V0LIcvVmZvguh6+dcz+mOTUZ3+2wBqiF786YZGaN\nnXMOqAwcdM7tTHfLhNA5ERHJhokT/VLSXbrAm29CkUI2DH7Z5mW0GdMGh6Nn/Z6Z+td66WKlKVqk\naNBvIaxlJ/YMB84FDtsSxTk3Ls2vy8zse2A10AyYmY3XIzo6mnLlyh12LCoqiqioqOzcVkSkwBg7\n1i8l3b49jBoFRQvZd+CstbPo8GEHqpevzqSbJnFG2TOCLilPxcbGEhsbe9ixxMTEHH0N840AmXyS\n2avAdUAT59z6DFy/GXjYOfeWmV0FfAmcnLbVwczWAkOcc8OO8PwIIC4uLo6IiIhM1ysiUtDt2QN3\n3gnvvgvdusE770Dx4kFXlbc+/OFDbv7kZppUa8LHN35MuZLljv+kQiA+Pp7IyEiASOdcfHbvl+kG\nrFBoaA9clcHQUBWoAGwMHYoDkoEWaa6pA1QD5me2HhGRwu7776FBA/jwQxg5Et5/v3CFBuccL857\nkaiPo+h8XmcmdZ2k0JCLMruOw3AgCvgHsMfMKoVOJTrn9ptZGeBx/BiHTUBt4HngJ/xaDTjndprZ\nO0CMmW0HdgEvA3M1o0JEJOOc82MYoqPhnHMgLg7q1g26qryVkprC3VPu5tVvXuXhJg/z1FVPYYVt\nJGgey+wYh774WRSz0h3vCbwPpAAXAt3xMy424APDY865pDTXR4euHY9fAGoK0C+TtYiIFFo7dsCt\nt/otsW+/HV56CUqVCrqqvLUvaR83/fcmPl35KW+2e5M+kX2CLqlQyOw6Dsfs2nDO7QeOu6u7c+4A\n0D/0EBGRTFi40M+Y2L7dB4dOnYKuKO9t3buV62Kv47uE75jYZSLtzmkXdEmFRiGbpCMikn+lpvqF\nnK64AipXhqVLC2doWP3Hai575zJ+2f4Ls26epdCQxxQcRETygc2boW1bGDAA7r0X5syBGjWCrirv\nLfp9EY3faQzA/N7zaXBGg4ArKnzy9/JVIiKFwPTpfoplaipMnQotWwZdUTA+W/kZncd3pn7l+nwa\n9SkVS1cMuqRCSS0OIiJhKjkZHnkErrkGzjsPvv228IaGNxe/SYf/dKB17dZM7z5doSFACg4iImFo\n/Xpo1gyeew6eftq3NFQuhIvyO+d4aPpD9P1fX/o16MdHN3xEqWKFbPpImFFXhYhImJk40S8bfeKJ\nMHs2XH758Z9TEB1MOUjvT3vzwXcf8OI1L3JP43u0RkMYUIuDiEiYOHDALxvdoQNceaWfNVFYQ0Pi\n/kTajmnLuGXj+LDTh9x72b0KDWFCLQ4iImHgp5/82gzLlsErr0C/foVvK+xDftv5G23HtOXXnb8y\n7V/TaFq9adAlSRoKDiIiAfvgA+jbF844wy/uVL9+0BUF5/uE72kzpg1FixRlbq+5nHvquUGXJOmo\nq0JEJCC7d0OPHvCvf8E//+n3mijMoWHGmhlcMfIKTi1zKvN7z1doCFNqcRARCcC330LnzvDbbzBq\nFHTvHnRFwRrz3Rh6TuxJ85rN+eiGjzipxElBlyRHoRYHEZE85BwMHw6XXgolS/pWhsIcGpxzDPpq\nEN0mdKPbhd34LOozhYYwp+AgIpJHtm+H66/3Ax9vvRUWLIA6dYKuKjjJqcn8+3//5qEZD/H4lY/z\nzj/eoVjRYkGXJcehrgoRkTwwbx5ERcHOnfDf/0LHjkFXFKw9B/cQ9XEUk36exNvXvU3viN5BlyQZ\npBYHEZFclJoKgwZB06ZQtapfm6Gwh4bNezbT/P3mzFgzg8+iPlNoyGfU4iAikksSEvyMiS+/hAcf\nhCeegGKFvCX+520/02ZMG3Yf3M2cnnOIqBIRdEmSSQoOIiK5YNo0HxrA7zNxzTXB1hMOFvy2gOti\nr6NCqQosuGUBNcrXCLokyQJ1VYiI5KCkJN+60KoVXHSRn3ap0AATV0zkqlFXUbdiXeb1nqfQkI8p\nOIiI5JB16/weEy+84Mc1TJ4MlSoFXVXwXlv0Gh3/05F257Rj2r+mcUqpU4IuSbJBXRUiIjngv/+F\n3r2hXDn46ito3DjoioKX6lJ58MsHGTxvMNGNonmx5YsUMf17Nb/T/4IiItmwf79fl6FTJ2jeHJYs\nUWgAOJB8gK7/7coL815gSKshxLSKUWgoINTiICKSRStX+mWjV6zwq0H27Vt4d7RMa8f+HXT4sAML\nflvARzd8RKdzOwVdkuQgBQcRkSwYNcq3NFSt6ne0vOiioCsKD+sT19N2TFs27t7I9O7Tubza5UGX\nJDksU+1GZvagmS0ys51mlmBmE8zsnCNc96SZbTCzvWY2zcxqpztfwsxeM7OtZrbLzMab2WnZfTMi\nIrlt1y6/t0SPHnDDDX6vCYUGb+mmpTR+pzF7kvYwt9dchYYCKrMdTk2AV4BLgauBYsAXZlbq0AVm\nNgC4A+gDNAT2AFPNrHia+wwFrgU6AU2B04GPs/geRETyxJIlEBkJEybA6NEwciSUKRN0VeFh2upp\nNB3ZlConVmF+7/nUrVg36JIkl2QqODjn2jrnRjvnljvnvgd6ANWAyDSX3QU85Zz73Dn3A9AdHww6\nAJhZWaAXEO2cm+2cWwL0BC43s4bZfkciIjnMOXjlFWjUCE480bcydOsWdFXhY9TSUbQd25Yrql3B\nrB6zqHxi5aBLklyU3SGu5QEH/AFgZjWBysD0Qxc453YCC4FD44wvwY+tSHvNSmB9mmtERMLCH3/4\nvSXuvNMPfpw/H875Wwdt4eSc4+k5T9NjYg96XNSDT6M+5cTiJwZdluSyLA+ONDPDdzl87Zz7MXS4\nMj5IJKS7PCF0DqAScDAUKI52jYhI4L7+Gm66CXbvhokT4R//CLqi8HFoS+y34t/iqaue4uEmD2Oa\nUlIoZKfFYThwLtAlh2oREQkLKSnwzDPQrBlUr+6XjVZo+Mvug7tp/2F7Ri4dyXvt3+ORpo8oNBQi\nWWpxMLNXgbZAE+fcxjSnNgGGb1VI2+pQCViS5priZlY2XatDpdC5o4qOjqZcuXKHHYuKiiIqKior\nb0NE5G82bvSbU82YAQ8/DI8/Dido4vqfEnYncO3Ya/lp209MumkS19TSRhzhJDY2ltjY2MOOJSYm\n5uhrmHMuc0/woaE9cKVz7pcjnN8AvOCcGxL6vSw+RHR3zn0U+n0L0MU5NyF0TR1gOdDIObfoCPeM\nAOLi4uKIiNAWrCKSO6ZM8VMtixaFMWP8SpDyl5VbV9JmTBv2J+9nUtdJ1K9cP+iSJAPi4+OJjIwE\niHTOxWf3fpldx2E40BW4CdhjZpVCj5JpLhsKPGJm15nZBcD7wG/ARPhzsOQ7QIyZNTOzSOBdYO6R\nQoOISG5LSoIBA6BNG4iI8F0TCg2Hm7t+Lpe9exmlipViwS0LFBoKscw2wPXFD36cle54T3xAwDk3\n2MxKA2/iZ118BbRxzh1Mc300kAKMB0oAU4B+mS1eRCS71qyBqCg/xXLwYLj3XiiiLRUO8/GPH9P1\nv11pVLUREzpP4ORSJwddkgQoU8HBOZehv07OuSeAJ45x/gDQP/QQEQnE+PFwyy1w8sl+BsWllwZd\nUfgZtmAY0VOj6Xx+Z95r/x4lTigRdEkSMOVqESl09u2D22/3S0Zfc41fEVKh4XCpLpV7p97L3VPv\n5v7L7mf3DcDhAAAgAElEQVTMP8coNAigTa5EpJBZvtzvaPnzz/DGG9Cnj3a0TG9/8n66T+jO+B/H\n80qbV7ij4R1BlyRhRMFBRAoF5/zeEv37+7UZFi2CCy4Iuqrw88e+P2j/YXsWb1jMxzd+TMd6HYMu\nScKMgoOIFHg7d/quibFjoXdvGDZMm1Mdydoda2kzpg1b9mxhRvcZND5TuwDI3yk4iEiBFhcHXbpA\nQoIPDlov7sjiN8Zz7dhrKV2sNPN6z+OcCtqQQ45MgyNFpEByDoYOhcaNoVw5iI9XaDiaKaum0HRk\nU84seybze89XaJBjUnAQkQJn61a/t0R0NNxxB8ydC7VrB11VeHp3ybu0G9uO5jWbM/PmmZxW5rSg\nS5Iwp64KESlQ5szxO1ru3w+ffQbt2gVdUXhyzjFw9kAGzh5I38i+vNL2FU4ooq8EOT61OIhIgZCS\nAk8+CVddBbVqwdKlCg1Hk5SSRO9PezNw9kAGtRjE8GuHKzRIhun/KSKS7y1fDrfd5ld/fOwxeOQR\n7Wh5NLsO7OKGj25gxpoZjO44mm4Xdgu6JMln9FdLRPKtffvg2Wfh+ef92gwzZ8KVVwZdVfjauGsj\n1469ltXbVzO562RanNUi6JIkH1JwEJF86csv/doM69bBAw/AQw9ByZLHf15htXzLclqPaU1Kagpf\n9/yaCypp9SvJGo1xEJF8JSEBunb1e0yccYbfAvvJJxUajuWrdV9x2buXUbZEWRbcskChQbJFwUFE\n8oXUVBgxAurWhalT/fLRM2dCvXpBVxbexi0bx9WjryaiSgRf9/yaqmWrBl2S5HMKDiIS9r7/Hpo0\n8QMgO3aEFSugRw9tTnUszjli5sfQeXxnbjj3BiZ3nUy5kuWCLksKAAUHEQlbe/bAgAEQEQF//AGz\nZsG770LFikFXFt5SUlO4e8rd3PvFvTx4xYOM7jia4kWLB12WFBAaHCkiYWnSJOjXDzZuhMcfh/vv\nhxIlgq4q/O08sJOeE3vyyYpPeP3a1+l7Sd+gS5ICRsFBRMLKhg1w110wfjxcfTVMm6blojNi3Y51\nvLzwZd6Kf4sUl8InnT/hujrXBV2WFEAKDiISFlJS4PXX/bTKUqVgzBi/KZXGMRzbot8XETM/hvE/\njuekEifRr0E/7mh4B2eUPSPo0qSAUnAQkcAtWeIHPn7zjf9z0CA4+eSgqwpfKakpTFw5kZj5Mcz9\ndS61Tq7FsNbDuLn+zZxY/MSgy5MCTsFBRAKze7dfInrYMDj3XL+L5WWXBV1V+Np1YBcjl45k6IKh\nrNmxhqbVm/JJ509od047ihYpGnR5UkgoOIhIICZO9Fteb9vml42+5x4oVizoqsLTbzt/4+WFLzMi\nbgS7D+6m8/mdGXfDOC45/ZKgS5NCSMFBRPLUr79C//4+OLRtC6++CjVrBl1VeIrbEEfMghjGLRtH\nmWJl6BPZh/4N+3NmuTODLk0KMQUHEckTycnw8su+a6JsWfjoI+jUSYMf00t1qXy28jNiFsQwZ90c\napavyUstX6Jn/Z6cVOKkoMsTyfwCUGbWxMw+NbPfzSzVzP6R7vzI0PG0j0nprilhZq+Z2VYz22Vm\n483stOy+GREJT4sWQYMGcN990LOn3wb7+usVGtLac3APw78ZTp1X69DhPx1ITk1m/A3j+bn/z9x5\n6Z0KDRI2stLiUAZYCrwD/Pco10wGegCH/rNwIN35oUAboBOwE3gN+BhokoV6RCRMJSbCww/D8OFQ\nvz4sXOgDhPxlw64NvLroVd5Y/AaJBxK5/tzrGd1xNI2qNgq6NJEjynRwcM5NAaYAmB313wsHnHNb\njnTCzMoCvYAuzrnZoWM9geVm1tA5tyizNYlIeHHOL+B0112wcye89JIf13CCOkf/tHTTUmLmx/Dh\nDx9S8oSS3BpxK/0v7U+N8jWCLk3kmHLrr3EzM0sAtgMzgEecc3+EzkWGXnf6oYudcyvNbD3QGFBw\nEMnH1qzxS0VPngwdOvhxDWdqLB/gxy9M/nkyMQtimLFmBtXKVeO5q5+j98W9tQGV5Bu5ERwm47sd\n1gC1gEHAJDNr7JxzQGXgoHNuZ7rnJYTOiUg+lJQEMTEwcKDfhOqTT6B9+6CrCg97k/Yy+tvRDFkw\nhJXbVtLwjIb85/r/8M96/+SEImqGkfwlx/8f65wbl+bXZWb2PbAaaAbMzM69o6OjKVfu8FQeFRVF\nVFRUdm4rItk0b55f8fHHH+Huu314OFELGLJp9yaGfzOc4d8M5499f9CxXkfebf8ujas25ug9vSJZ\nFxsbS2xs7GHHEhMTc/Q1cj3qOufWmNlWoDY+OGwCiptZ2XStDpVC545qyJAhRERE5F6xIpIp27fD\nAw/AiBF+0OPixXDxxUFXFbzvE75nyIIhjPl+DMWKFKP3xb2589I7qXVKraBLkwLuSP+Yjo+PJzIy\nMsdeI9eDg5lVBSoAG0OH4oBkoAUwIXRNHaAaMD+36xGR7HMOYmMhOhr27/eLOPXtC0UL8arHzjmm\nrp5KzPwYpv0yjaplq/L0VU9za+StlC9ZPujyRHJMpoODmZXBtx4camc7y8wuAv4IPR7Hj3HYFLru\neeAnYCqAc26nmb0DxJjZdmAX8DIwVzMqRMLfqlVw++3w5Zdwww0wdCicfnrQVQVnf/J+PvjuA4Ys\nGMKPW34kokoEY/45hhvOvYFiRbWGthQ8WWlxuATf5eBCj5dCx0cB/wYuBLoD5YEN+MDwmHMuKc09\nooEUYDxQAj+9s18WahGRPHLgAAweDM88A1WqwKRJ0KZN0FUFZ/Oezbz+zeu89s1rbN27lX/U+Qev\nX/s6Tao10fgFKdCyso7DbI694mTrDNzjANA/9BCRMDd7tu+KWLXKr/746KNQunTQVQXjxy0/MmT+\nEEZ/N5oiVoSe9Xtyd6O7ObvC2UGXJpInNA9IRI5q61a4/3547z2/3fWSJXD++UFXlfecc0xfM52X\n5r/ElFVTqHJiFR6/8nFuu+Q2Til1StDlieQpBQcR+RvnYNQo37qQkuJnTfTuDUUyvbtN/nYg+QCx\nP8QSMz+G7zd/T/3K9Xm/w/t0Pr8zxYsWD7o8kUAoOIjIYVas8N0Ss2dD165+uehKlYKuKm9t3buV\nNxa/wauLXiVhTwLtzmnHsNbDaFajmcYvSKGn4CAigJ9W+eyz8NxzUL06TJsGV18ddFV5a+XWlQxd\nMJRR347C4bj5opu5u9Hd1K1YN+jSRMKGgoOI8OWXforlunV+QaeHHoKSJYOuKm8455i1dhYxC2L4\n/KfPqVSmEg81eYi+l/SlYumKQZcnEnYUHEQKsYQEuOceGDsWrrwSPvsM6haSf1wfTDnIuGXjiJkf\nw5JNS7jgtAsY2X4kUedHUeKEEkGXJxK2FBxECqHUVHj7bRgwwK/2OHIk3HwzFIbu+z/2/cGIuBG8\nsugVNuzaQOvarfmi2xdcfdbVGr8gkgEKDiKFzPff+8GP8+ZBz55+UaeKhaBF/udtPzNs4TBGLh1J\nSmoK/7rwX9zd6G7OO+28oEsTyVcUHEQKib174ckn/SyJ2rVh1izfPVGQOef4av1XxMyP4dOVn1Kx\ndEXuv+x+/t3g35xW5rSgyxPJlxQcRAqBSZOgXz/YuBEef9wv6lSiAHfjJ6UkMf7H8cQsiGHxhsXU\nq1iPEdeNoOsFXSlVrFTQ5YnkawoOIgXYhg1w993w0Ud+auW0ab61oaDasX8Hb8W9xcuLXua3nb9x\nzVnXMLnrZFrWakkRK2SrV4nkEgUHkQIoJQXeeOOvaZVjxkBUVMEd/Lhlzxae+/o5RsSP4EDyAbpe\n2JXoRtFcWOnCoEsTKXAUHEQKmCVL4Lbb4Jtv/J+DBsHJJwddVe7YeWAnMfNjeGn+SxSxItx16V30\na9CPKidVCbo0kQJLwUGkgNi9Gx57DIYNg3PPhblz/cZUBdH+5P28/s3rPPPVM+w+uJv+DfvzwBUP\nUKF0haBLEynwFBxECoCJE6F/f7+b5aBBEB0NxYoFXVXOS05NZtTSUQycPZANuzbQ6+JePHblY1Qt\nWzXo0kQKDQUHkXzsxx/9OIaJE6FtW3j1VahZM+iqcp5zjo+Xf8wjMx5h5baV3HjejTx11VOcU+Gc\noEsTKXQUHETyocWL/YZUEybAmWf6WROdOhW8wY/OOb785UsemvEQizcsplWtVoztNJaIKhFBlyZS\naGl+kkg+4Zzf6rpVK2jQAH74Ad55B1atguuvL3ihYeFvC2nxfgtaftCSE4qcwMybZzKl2xSFBpGA\nqcVBJMw55xdwevZZv0z0RRfBf/7jWxiKFg26upy3bPMyHpn5CJ+s+ITzTzufiV0mct0512kfCZEw\noeAgEqZSUmD8eD/Y8dtv/QyJ//0P2rQpeK0LAGt3rOWJWU8w+rvRVCtXjfc7vM9NF9xE0SIFMB2J\n5GMKDiJh5uBBGD0ann8efv4ZWrb0+0o0bVowA8PmPZt5Zs4zvL74dU4pdQrDWg+jT2QfihctHnRp\nInIECg4iYWLvXr/V9QsvwG+/wT//CWPHwiWXBF1Z7kjcn8hL818iZn4MRYsU5fErH+euRndxYvET\ngy5NRI5BwUEkYDt2wPDhMGQIbN8OXbvCgAF+EaeCaF/SPl775jUGfT2IvUl7ubPhnQy4YgCnlDol\n6NJEJAMyPavCzJqY2adm9ruZpZrZP45wzZNmtsHM9prZNDOrne58CTN7zcy2mtkuMxtvZtrjVgqV\nzZv9GgzVq/vtrm+4wc+QGDWqYIaG5NRk3op7i7NfOZsHvnyAG869gVX9V/H8Nc8rNIjkI1mZjlkG\nWAr8G3DpT5rZAOAOoA/QENgDTDWztB2WQ4FrgU5AU+B04OMs1CKS7/z6K9x1F9SoAa+84veTWLPG\ntzrUqBF0dTkv1aUybtk4zht+Hn0+70OT6k1Y3m85b7R7gzPKnhF0eSKSSZnuqnDOTQGmANiR50fd\nBTzlnPs8dE13IAHoAIwzs7JAL6CLc2526JqewHIza+icW5SldyIS5n76yQ94fP99KFsWHngA7rgD\nTimg/9h2zvHF6i94cPqDLNm0hDa12/Bhpw+5uMrFQZcmItmQowtAmVlNoDIw/dAx59xOYCHQOHTo\nEnxgSXvNSmB9mmtECoylS6FzZ6hbFyZP9uFh3Tq/IVVBDQ3zf53PVaOuovWY1pQuVpo5PeYwqesk\nhQaRAiCnB0dWxndfJKQ7nhA6B1AJOBgKFEe7RiTfmzvXL9o0aRKcdRa88QZ07w4lSwZdWe75YfMP\nPDzjYT5d+SkXVrqQz6M+p+3ZbbV4k0gBolkVIjnIOfjiCx8Y5syB886DDz7wLQ4nFOC/bWu2r+Hx\nWY/zwXcfUPPkmoz55xi6nN+FIqZV7UUKmpz+T9kmwPCtCmlbHSoBS9JcU9zMyqZrdagUOndU0dHR\nlCtX7rBjUVFRREVFZbdukWxJTfUbTj37LMTH+70kPvkErrsOihTg786E3Qk8Pedp3ox7kwqlK/Ba\n29foHdFbizeJBCQ2NpbY2NjDjiUmJuboa5hzf5sYkfEnm6UCHZxzn6Y5tgF4wTk3JPR7WXyI6O6c\n+yj0+xb84MgJoWvqAMuBRkcaHGlmEUBcXFwcERHa4EbCR1ISxMb6ZaFXrIDmzf0Uy+bNC+Yqj4fs\n2L+DF+e9yJAFQyhetDgDLh9A/4b9KVO8TNCliUg68fHxREZGAkQ65+Kze79MtziYWRmgNr5lAeAs\nM7sI+MM59yt+quUjZrYKWAs8BfwGTAQ/WNLM3gFizGw7sAt4GZirGRWSX+zbByNHwuDBfqDjddf5\n3xs1Crqy3LU3aS+vLnqV575+jv3J+7nr0rv4v8v/j5NLnRx0aSKSR7LSVXEJMBM/CNIBL4WOjwJ6\nOecGm1lp4E2gPPAV0MY5dzDNPaKBFGA8UAI/vbNflt6BSB7audMPcoyJgS1b/NiFTz+FCy8MurLc\nlZSSxLtL3uXJOU+yec9m+kT04ZGmj1DlpCpBlyYieSwr6zjM5jjTOJ1zTwBPHOP8AaB/6CES9rZt\ng2HD/IJNe/ZAjx7wf/8HtWsf96n52qHFmx6d+Sir/1jNTRfcxMBmA6l1Sq2gSxORgBTgcd4i2ff7\n77514c03/YyJ226De+6BqlWDrix3OeeYsmoKD814iKWbltLunHZ8fOPHXFipgDetiMhxKTiIHMHq\n1X78wnvvQalSEB0Nd94Jp54adGW5b+76uTw4/UG+Wv8VTao14eueX3N5tcuDLktEwoSCg0gaP/zg\nZ0h8+CFUrOg3n7r9dr9EdEH3XcJ3PDzjYT7/6XMuqnQRk26aROvarbV4k4gcRsFBBFi40AeGiROh\nWjV4+WXo1cu3NhR0v2z/hcdmPsbY78dy1slnEdsplhvPu1GLN4nIESk4SKHlHMyc6Rdtmj4d6tTx\nXRM33QTFigVdXe7buGsjT895mhHxIzitzGm8fu3r9Lq4F8WKFoI3LyJZpuAghU5qKnz+uQ8MCxfC\nxRfDRx9Bx45QtGjQ1eW+7fu2M3juYIYtHEbJE0ryTPNnuKPhHZQuVjro0kQkH1BwkEIjORnGjfNd\nEj/8AE2a+N0qW7Uq2Ks8HrI3aS8vL3yZ5+c+z8GUg9zT+B7uu+w+ypcsH3RpIpKPKDhIgXfgAIwa\n5bez/uUXaNMGhg/3waEwSEpJ4u34t3lqzlNs3buV2yJv4+GmD1P5RG1GKyKZp+AgBdbu3TBiBLz0\nEmzcCJ06wfjxvmuiMEh1qXz4w4c8OvNR1mxfQ7cLuzGw2UBqnlwz6NJEJB9TcJACZ/t2ePVVv9Jj\nYiJ06wYDBkDdukFXljecc0z6eRIPzXiI7xK+o32d9kzsMpHzTzs/6NJEpABQcJACY9MmGDLEd0Mk\nJ8Mtt8B990H16kFXlne+WvcVD05/kLm/zuXK6lcyr9c8Gp/ZOOiyRKQAUXCQfG/tWnjhBXjnHShe\nHO64A+6+GypVCrqyvLN001IenvEwk36eRESVCKZ0nULLWi21eJOI5DgFB8m3li+H556DMWOgfHl4\n9FHo18//XJA559iydwsrtq5gxdYVTF8znXHLxnFOhXMYd/04Op3bSYs3iUiuUXCQfOXAAfjsM79Q\n06RJcPrp8OKLcOutUKZM0NXlrKSUJH7Z/sufAWHFthV//rxj/w4AilpR6lSsw1vXvUWP+j04oYj+\nSotI7tJ/ZSTsOQeLF/uwEBvrBz82bOhnTPzrX1CiRNAVZs/2fdtZsXUFK7et/CskbF3B6u2rSU5N\nBqBsibLUq1iPOhXr8I9z/kHdinWpW7EutU6pRfGixQN+ByJSmCg4SNjasAE++MAHhuXLfetCnz5w\n881Qr17Q1WVOSmoK6xLX/RkKVm5d+WcLwuY9mwEwjOrlq1O3Yl1a1279ZzioW7EulcpU0ngFEQkL\nCg4SVvbv9xtNvfcefPGFH+zYoYOfLXH11eG/JPTug7t9KEjTvbBy60p+2vYTB1IOAFC6WGnqVKhD\n3Yp1aV6j+Z/h4OwKZ2vZZxEJewoOEjjnYMECv7rjhx/6tRcaN4bXX4cbbwy/wY7OOX7f9fth3QqH\nHr/v+v3P604/6XTqVqxLk2pNuDXiVupU9GGhatmqGrwoIvmWgoME5tdfYfRoHxh++gmqVvWzIm6+\nGc45J+jqYF/SPn7+4+e/WhC2/dXNsCdpDwDFixbn7FPOpm7FuvSo34O6FetSp0Id6lSsQ9kSZQN+\nByIiOU/BQfLU3r0wYYIPC19+CSVL+qWgX3sNrroq77sinHNs3rP5r7EHaQYort2xFocDoGLpitSt\nWJeIyhHcdP5NPiBUrEON8jU0k0FEChX9F09ynXMwd64ftzBuHOza5TeYeustuOEGKJsH/zBPSkli\n9fbVh3UrHAoJaac2nnXyWdStWJfrz73+z7EHdSrUoULpCrlfpIhIPqDgILlm3Tp4/33furB6tV/6\nOToauneHWrVy5zX/2PfHEQcnpp/aeCgUaGqjiEjmKDhIjtqzBz7+2LcuzJzpF2W6/np4+21o2hSK\n5NCYwITdCcRtjPvb4MQte7cAf01trFOhjqY2iojkIAUHybbUVPjqKx8Wxo/321k3a+Z/79QJTjwx\n+6+RlJLEvF/nMXX1VKasmsKSTUuAw6c2tqjZQlMbRURyWY4HBzN7HHg83eEVzrlz01zzJHALUB6Y\nC9zunFuV07VI7vrll7+6ItauhbPOgvvv910RNWpk//5rd6xl6qqpTFk9hem/TGfXwV2cWvpUWtVu\nxT2N76FJtSacWe5MTW0UEclDudXi8APQAjjUHpx86ISZDQDuALoDa4GngalmVs85dzCX6pEcsmsX\nfPSRDwtz5vjWhBtvhB494IorIDs9AHuT9jJ77ew/WxVWbltJUSvKZWdexgNXPECrWq24uMrFCgoi\nIgHKreCQ7JzbcpRzdwFPOec+BzCz7kAC0AEYl0v1SDakpvrxCqNG+fEL+/ZBixZ+DYaOHbO+uZRz\njuVblzNl1RSmrJrCnHVzOJBygGrlqtG6VmsGtRhE85rNKVeyXM6+IRERybLcCg5nm9nvwH5gPvCg\nc+5XM6sJVAamH7rQObfTzBYCjVFwCCs//+zDwujRsH49nH02PPSQ31iqWrWs3XPH/h1M/2U6U1ZN\nYerqqfy681dKnlCSK6tfyXNXP0erWq2oW7GuBi+KiISp3AgOC4AewEqgCvAEMMfMzseHBodvYUgr\nIXROApaY6NdaeO89mDfPr7HQpYtfzbFx48x3RaS6VOI3xv/ZqrDgtwWkuBTqVaxHp3qdaF27NU2r\nN6VUsVK58n5ERCRn5XhwcM5NTfPrD2a2CFgH3AisyM69o6OjKVfu8GbrqKgooqKisnPbQi8lBaZP\n92FhwgQ4eBCuuQbGjvUbTJXK5Hd6wu4Evlj9BVNWT+GL1V+wde9WypYoS4uaLRh+7XBa1WpF9fLV\nc+W9iIgUZrGxscTGxh52LDExMUdfw5xzOXrDI76IDw/TgLeB1UB959x3ac7PApY456KP8vwIIC4u\nLo6IiIhcr7ewWLHir66I33+HunX9IMdu3eCMMzJ+n6NNlYysEkmrWq1oXbs1jao2oljRYrnzRkRE\n5Kji4+OJjIwEiHTOxWf3frm+joOZnQjUBkY559aY2Sb8jIvvQufLApcCr+V2LQLbt8N//uNbFxYu\n9DtPRkX5wNCgQca7ItbuWPvnOIW0UyVb1mrJPY3voWWtlpxW5rTcfCsiIhKA3FjH4QXgM3z3xBnA\nQCAJ+DB0yVDgETNbhZ+O+RTwGzAxp2sRLzkZvvjCty5MnAhJSdC6tR/LcN11fqOp4zk0VfJQWEg7\nVXLA5QNoXbu1pkqKiBQCudHiUBUYC1QAtgBfA42cc9sAnHODzaw08CZ+AaivgDZawyHnLVvmWxY+\n+AA2bYLzzoOnn4auXaFKlWM/V1MlRUTkSHJjcORxRyo6557Az7aQHLZtG8TG+taFxYvhlFPgppt8\nV0RExLG7IjRVUkREjkd7VRQASUkwZYpvXfjsM79g07XX+sWarr0WSpQ48vOONlWybsW6miopIiJH\npOCQj337rW9ZGDMGNm+Giy6CwYN9C8NpRxmXeKSpkicVP4mrz7paUyVFROS4FBzymS1b/PoK770H\nS5dCxYp++uTNN0P9+n+//mhTJSOqRNAnoo+mSoqISKYoOISx7dthyRKIi4P4eP/nzz9DsWLQrh0M\nHAht2vjf0zreVMlrzrqGSidWCuZNiYhIvqbgECa2bvXh4FBAiI/321aD30Sqfn0fEh5+2I9bqFjx\nr+dqqqSIiOQVBYcAJCQc3ooQH+83kQI46SQ/+6F9e4iM9D+fcw4ULfrX851z/Ljl6FMln23xLC1q\nttBUSRERyXEKDrnIOdiw4fCAEBfnjwGcfLIPBp07/xUSatWCIkdoGPh95+/MXjebmWtm/jlVskTR\nEjSr0YxBLQbRunZrTZUUEZFcp+CQQ5zzrQZpQ0J8vG9dAN+1EBnpBzEeCgk1ahx9XYXfdv7GrLWz\nmL12NrPWzWLVH6sAOPfUc+lUrxOtareiafWmlC5WOm/eoIiICAoOWeIcrFnz9+6Gbdv8+UqVfDjo\n08cHhMhIqFr12Isv/Zr4K7PWzvJhYd1sVm9fDcB5p55Hq1qtGNRiEE2rN9X+DyIiEigFh+NITYVV\nqw4PCUuWwI4d/nzVqj4c3Hmn/zMiAk4//fj3XZ+4/rCg8Mt2PxLygtMuoE3tNjSr0Yym1ZtyaplT\nc/HdiYiIZI6CQxopKbBy5eGtCEuWwK5d/nz16r714L77/J8XX+xbFzJi3Y51Piis82Fh7Y61AFxY\n6ULand2OK2tcSdPqTalYuuKxbyQiIhKgQhsckpPhxx8PDwlLl/L/7d17bJX1Hcfx9/eUm4oMByro\nRBTFoYUWKywItZ3ituwPiXfxFiSEkG2J8pcuLtE53WLcyLKLxmSZxjhRk21xJjqFuTLmjWw4BSyo\nQLmX6yjaQmT63R+/h/b0eNqeYnt+Pef5vJInTZ/n6en3+fbb029/z+VHW1vYPmFCaA7uvbdjJGHU\nqMJfv+lgU/uIQkNTA1tatgBQdXoVcy6YQ93ZoVEYdWIvXlRERCSyVDQOn34aZorMPt3w3ntw5Ei4\n7mDixNAYXHNN+Dh1KowcWfjru3tHo5CMKGxt2YphVI2p4uqvX03d+Dpqx9WqURARkZJWdo3DkSOw\nZk3nkYQ1a0LzkMnApEmhObj55vCxujo8O6E33J3NBzd3GlHYdmgbhlE9ppprJ11L/fh6asfVcsoJ\np/TPgYqIiERQ0o1DW1uY6Cm7SVi3LpyGqKiAiy4KpxvmzQsfp0wJT2HsLXdn4383tt8a2dDUwPZD\n28lYhuox1Vx/4fXUj69n1rhZahRERKSslVTj8M47sHJlR6PQ2Bjuehg8GCZPhunTYdGiMJIwZQoM\nG3Z838fd+ejAR6zYsqJ9RGHHxzvIWIaLx17MjRfd2N4ojBzWi3MaIiIiJa6kGocFC2Do0DB9dG0t\n3DdnfGsAAAeBSURBVHVXaBIqK2HIkON/XXfnwwMftt8a2dDUwM6Pd5KxDDVja5hbObe9UdBjnEVE\nJM1KqnF49tlwAWPubJC95e58sP+DTo3Crk92UWEV1JxRw62Tb6VufB2zxs1ixNARfRO8iIhIGSip\nxuH884+vaXB3Nuzf0OmBS82fNFNhFVxyxiXcXnU7dWfXMXPcTDUKIiIi3SipxqFQ7s76fevbb49c\n0bSC3a27qbAKpp05jXlV86gfX8+lZ13KyUN7eUuFiIhIipVF4+DuNO5r7DSisKd1D4Myg5h2xjTm\nT53f3igMHzI8drgiIiIlqyQbB3fn/b3vdxpR2Nu2l0GZQUw/czoLpi5obxROGnIc91+KiIhIXpnY\nAfTGc2uf47rnr+O0n59G5WOVLH5lMc2fNLOwZiHLblvGwbsP8vr813noioe4csKVZds0LF26NHYI\nA4LyECgPHZSLQHnooFz0vaiNg5l938w2m9lhM3vLzKZ1t/+SN5ewp3UPi2oWsfy25Ry85yAr71jJ\ng5c/yOxzZ5dto5BLvwiB8hAoDx2Ui0B56KBc9L1opyrM7EbgF8BCYBWwGHjFzCa6+758X9NwRwMz\np88sYpQiIiKSLeaIw2LgcXd/yt3XA4uANmB+V19wwqATihWbiIiI5BGlcTCzwUAN8Ldj69zdgeXA\njBgxiYiISM9inaoYDVQAu3PW7wYuyLP/MIDGxsZ+Dqs0tLS0sHr16thhRKc8BMpDB+UiUB46KBed\n/nYe5wxOnVn4R7+4zGwssAOY4e5vZ61/GLjM3Wfk7H8z8IfiRikiIlJWbnH3Z77si8QacdgHfAac\nnrP+dKA5z/6vALcATcCRfo1MRESkvAwDxhP+ln5pUUYcAMzsLeBtd78z+dyArcCv3P2RKEGJiIhI\nt2I+OXIJ8KSZ/ZuO2zFPBJ6MGJOIiIh0I1rj4O7Pm9lo4AHCKYr/AN92972xYhIREZHuRTtVISIi\nIqWnpOaqEBERkbjUOIiIiEjBBnTjYGb3mdnnOcv7sePqb2ZWa2Z/MbMdyTFflWefB8xsp5m1mdky\nMzsvRqz9radcmNkTeWrkpVjx9hcz+6GZrTKzQ2a228z+bGYT8+xX1nVRSB5SVBOLzOxdM2tJljfM\n7Ds5+5R1PUDPeUhLPeQys3uSY12Ss/5L18SAbhwSawkXT45JlllxwymKkwgXi34P+MJFKGZ2N/AD\nwgRh04FWwgRhQ4oZZJF0m4vEy3SukbnFCa2oaoFfA98AZgODgVfNrH0Cl5TURY95SKShJrYBdwMX\nEx7h/xrwgplNgtTUA/SQh0Qa6qFdMtP0QuDdnPV9UxPuPmAX4D5gdew4Iufgc+CqnHU7gcVZn48A\nDgM3xI43Qi6eAP4UO7YIuRid5GNWmuuiizyksiaSY98P3JHWeugiD6mqB2A4sAG4HPg7sCRrW5/U\nRCmMOJyfDFNvNLOnzeys2AHFZGbnEDrm7AnCDgFvk94JwuqTYev1ZvaomX01dkBFMJIwAnMAUl0X\nnfKQJVU1YWYZM7uJ8CycN9JaD7l5yNqUpnr4LfCiu7+WvbIvayLmA6AK8RYwj9A9jQXuB/5hZpXu\n3hoxrpjGEN4o800QNqb44UT3MvBHYDMwAfgZ8JKZzfCkpS43yVNWfwn8092PXfOTurroIg+Qopow\ns0rgTcIjhT8Grnb3DWY2gxTVQ1d5SDanqR5uAqqBS/Js7rP3iAHdOLh79nO115rZKmALcANh+ElS\nzt2fz/p0nZmtATYC9YRhunL0KHAhMDN2IJHlzUPKamI9UAV8BbgOeMrMLosbUhR58+Du69NSD2b2\nNUIjPdvdj/bn9yqFUxXt3L0F+AAouyuDe6EZMAqfICxV3H0zYRK1sqwRM/sN8F2g3t13ZW1KVV10\nk4cvKOeacPf/ufsmd3/H3e8lXAx3Jymrh27ykG/fcq2HGuBUYLWZHTWzo0AdcKeZfUoYWeiTmiip\nxsHMhhN+2N2+UZSzpOibgSuOrTOzEYSrzN/o6uvSIum6R1GGNZL8sZwDfNPdt2ZvS1NddJeHLvYv\n25rIIwMMTVM9dCEDDM23oYzrYTkwmXCqoipZ/gU8DVS5+yb6qCYG9KkKM3sEeJFweuJM4MfAUWBp\nzLj6m5mdRGiQLFl1rplVAQfcfRthOOpHZvYRYarxnwDbgRcihNuvustFstxHOH/ZnOz3MGFUqk+m\njx0ozOxRwi1kVwGtZnbsv4YWdz821XzZ10VPeUjqJS018VPC+futwMnALYT/ML+V7FL29QDd5yFN\n9ZBc99fpOUdm1grsd/fGZFXf1ETsW0d6uK1kaXJQhwlF8QxwTuy4inDcdYRbzD7LWX6ftc/9hFtr\n2gi/AOfFjrvYuSBcCPVXwhvCEWAT8Bhwauy4+yEP+XLwGXB7zn5lXRc95SFlNfG75PgOJ8f7KnB5\nmuqhpzykqR66yM1rZN2O2Vc1oUmuREREpGAldY2DiIiIxKXGQURERAqmxkFEREQKpsZBRERECqbG\nQURERAqmxkFEREQKpsZBRERECqbGQURERAqmxkFEREQKpsZBRERECqbGQURERAr2f8qZkHMYQztp\nAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ "<matplotlib.figure.Figure at 0x7ff59814def0>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "plt.plot(N,k_jacobi,label='Jacobi-Verfahren')\n",
+ "plt.plot(N,k_gauss_seidel,label='Gauß-Seidel-Verfahren')\n",
+ "plt.show()"
+ ]
+ }
+ ],
+ "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.5.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/10_Poisson1d_function.ipynb b/10_Poisson1d_function.ipynb
new file mode 100644
index 0000000..40593fa
--- /dev/null
+++ b/10_Poisson1d_function.ipynb
@@ -0,0 +1,168 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Wirkung der 1d-Poissongleichung auf einen Vektor\n",
+ "Es soll eine Funktion erstellt werden, die die Wirkung der 1d-Poissongleichung auf einen Vektor x simuliert. Die Idee dafür ist es am Anfang und Ende von x eine zusätzliche 0 einzufügen (sogenanntes zero padding) und dann entsprechende Ausschnitte auf x zu addieren."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[0 1 2 3 0]\n"
+ ]
+ }
+ ],
+ "source": [
+ "from numpy import pad, array\n",
+ "\n",
+ "a=array([1,2,3])\n",
+ "print(pad(a,1,mode='constant'))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Nun die Funktion, die die Wirkung der 1d-Poissonmatrix auf einen Vektor berechnet:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "def poisson1d(x):\n",
+ " h=pad(x,1,mode='constant')\n",
+ " y=2*x-h[:-2]-h[2:]\n",
+ " n=max(x.shape)+1\n",
+ " return y*n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Wir testen an zwei einfachen Beispielen, um uns von der Richtigkeit zu überzeugen:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[ 0 0 16]\n",
+ "[ 10. 0. 0. 0. 0. 0. 0. 0. 10.]\n"
+ ]
+ }
+ ],
+ "source": [
+ "from numpy import ones\n",
+ "\n",
+ "print(poisson1d(a))\n",
+ "\n",
+ "b=ones(9)\n",
+ "print(poisson1d(b))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Als Hinweise für den 2d-Fall: "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[0, 0, 0, 0],\n",
+ " [0, 1, 2, 0],\n",
+ " [0, 3, 4, 0],\n",
+ " [0, 0, 0, 0]])"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "c=array([[1,2],[3,4]])\n",
+ "pad(c,1, mode='constant')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[1 2]\n",
+ " [3 4]]\n",
+ "[1 2 3 4]\n"
+ ]
+ }
+ ],
+ "source": [
+ "d=array([1,2,3,4])\n",
+ "n=max(d.shape)\n",
+ "e=d.reshape((int(n**(1/2)),int(n**(1/2))))\n",
+ "print(e)\n",
+ "print(e.ravel())"
+ ]
+ }
+ ],
+ "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.5.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
--
GitLab