assignment_2_352210.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before you submit this assignment, make sure everything runs as you expect it to.  \n",
    "* **restart the kernel** (in the menubar, select Kernel -> Restart)\n",
    "* **run all cells** (in the menubar, select Cell -> Run All)\n",
    "\n",
    "Make sure to fill in any place that says \n",
    "```python\n",
    "# YOUR CODE HERE\n",
    "```\n",
    "and fill in your name (and of your fellow students, if you submit together), and according matriculation number below. If you are more than one student, separate by a comma **,**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "NAME = \"Kittiphon Boonma\"\n",
    "MATNUM = \"352210\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assignment 2: Finite Difference models\n",
    "##### (30 Points)\n",
    "\n",
    "After the introduction to scientific computing with Python in the exercise - and the fundamentals of fluid mechanics and porous media flow in the lecture - we are now ready to bring those two things closer together to get closer to solving actual reservoir engineering problems with numerical methods! \n",
    "\n",
    "During this exercise, we will first have a look at several methods to solve the types of algebraic systems that we obtain for 1-D Laplace and Poisson equations. Starting from the matrix we developed at the end of the last exercise, we will move from direct solvers to iterative solvers. Finally, we will have a look at the implementation of different boundary conditions and sources and sinks.\n",
    "\n",
    "If you follow these examples carefully and implement the code correctly for the 1-D example that we consider here, then you are well prepared for the assignment!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='https://fonts.googleapis.com/css?family=Marvel' rel='stylesheet' type='text/css'>\n",
       "<link href='https://fonts.googleapis.com/css?family=Lato' rel='stylesheet' type='text/css'>\n",
       "<link href='https://fonts.googleapis.com/css?family=Source+Code+Pro' rel='stylesheet' type='text/css'>\n",
       "<link href='https://fonts.googleapis.com/css?family=Open+Sans' rel='stylesheet' type='text/css'>\n",
       "<link href='https://fonts.googleapis.com/css?family=Raleway' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    div.cell{\n",
       "        width:900px;\n",
       "        margin-left:1.5%;\n",
       "        margin-right:1.5%;\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Marvel', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Marvel', sans-serif;\n",
       "    }\n",
       "    h3 {\n",
       "        font-family: 'Marvel', sans-serif;\n",
       "        margin-top: 10px;\n",
       "        margin-bottom: 2px;\n",
       "    }\n",
       "    h4 {\n",
       "        font-family: 'Mavel', sans-serif;\n",
       "    }\n",
       "    h5 {\n",
       "        font-family: 'Marvel', sans-serif;\n",
       "    }\n",
       "    div.text_cell_render{\n",
       "        font-family: 'Open Sans', 'Raleway', 'Helvetica Neue', 'Lato', sans-serif;\n",
       "        line-height: 120%;\n",
       "        font-size: 110%;\n",
       "        width: 775px;\n",
       "        margin-left: auto;\n",
       "        margin-right: auto;\n",
       "    }\n",
       "    div.cell.code_cell{\n",
       "        background-color: rgb(240,240,240);\n",
       "        border-radius: 5px;\n",
       "        width: 775px;\n",
       "        margin-top: .7em;\n",
       "        margin-left: 9.5em;\n",
       "        margin-right: auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "        font-family: 'Source Code Pro';\n",
       "        font-size: 100%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 100;\n",
       "        font-size: 48pt;\n",
       "        line-height: 100%;\n",
       "        color: #0168B3;\n",
       "        margin-bottom: 0.7em;\n",
       "        margin-top: 0.7em;\n",
       "        display: block;\n",
       "    }\n",
       "    .text_cell_render h2 {\n",
       "        font-weight: 100;\n",
       "        font-size: 28pt;\n",
       "        line-height: 100%;\n",
       "        color: #799096;\n",
       "        margin-bottom: 0.7em;\n",
       "        margin-top: 0.7em;\n",
       "        display: block;\n",
       "    }\n",
       "    .text_cell_render h3 {\n",
       "        font-weight: 100;\n",
       "        font-size: 24pt;\n",
       "        line-height: 100%;\n",
       "        color: #bd131c;\n",
       "        margin-bottom: 0.7em;\n",
       "        margin-top: 0.7em;\n",
       "        display: block;\n",
       "    }\n",
       "    .text_cell_render h4 {\n",
       "        font-weight: 100;\n",
       "        font-size: 20pt;\n",
       "        line-height: 100%;\n",
       "        color: #d64911;\n",
       "        margin-bottom: 0.7em;\n",
       "        margin-top: 0.7em;\n",
       "        display: block;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 330;\n",
       "        font-size: 15pt;\n",
       "        color: #799096;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb(189, 19,28)\n",
       "    }\n",
       "    /*div #notebook{\n",
       "        background: #fff;\n",
       "        width: 1000px;\n",
       "        padding-left: 1em;\n",
       "        margin: auto;\n",
       "    }*/\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                            extensions: [\"AMSmath.js\"]\n",
       "                        },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', \n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "            });\n",
       "</script>\n",
       "        \n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# use this cell to import further needed libraries. \n",
    "#But generally, these should suffice.\n",
    "\n",
    "import numpy as np\n",
    "import math as mt\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "import colormaps as cmaps\n",
    "plt.register_cmap(name='viridis', cmap=cmaps.viridis)\n",
    "\n",
    "\n",
    "from IPython.core.display import HTML\n",
    "css_file = 'nre2_style.css'\n",
    "HTML(open(css_file,'r').read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the tridiagonal matrix\n",
    "##### (10 Points)\n",
    "\n",
    "During the last lecture, you have seen the meaning of the tridiagonal matrix that we developed at the end of the last exercise: the matrix can be used to solve a problem that can be described with the Laplace equation in 1-D (for Dirichlet boundary conditions).\n",
    "\n",
    "Here the matrix as a reminder:\n",
    "\n",
    "$A = \\begin{bmatrix}\n",
    "-2 & 1 & 0 & 0 & \\cdots & 0\\\\\n",
    "1 & -2 & 1 &  0 & \\cdots & 0\\\\\n",
    "0 & 1 & -2 & 1 & \\cdots & 0\\\\\n",
    "\\vdots & & & & & \\vdots \\\\\n",
    "0 & \\cdots & 0 & 1 & -2 & 1 \\\\\n",
    "0 & \\cdots & 0 & 0 & 1 & -2\n",
    "\\end{bmatrix}$\n",
    "\n",
    "And the system of algebraic equations that we want to solve:\n",
    "\n",
    "$$A u = b$$ \n",
    "\n",
    "where $u$ is the vector with the (unknown) state variables that we want to know, and $b$ encapsulates the boundary conditions.\n",
    "\n",
    "The straight-forward way to solve this equation is to use the inverted matrix $A^{-1}$. \n",
    "\n",
    "**Exercise**:\n",
    "\n",
    "* write function to create the tridiagonal matrix **(2 Points)**\n",
    "```python\n",
    "def tridiag_121(n):\n",
    "    \"\"\"\n",
    "    Create a tridiagonal matrix of 1 -2 1 form of size n x n \n",
    "    \"\"\"\n",
    "```\n",
    "* extend it to solve the following problem:  **(3 Points)**  \n",
    "solve for temperatures (in a steady-state temperature field without sources and sinks) on 10 nodes, considering $b$ with values of of $T_0 = 10^o$C and $T_1 = 20^o$C (These are Dirichlet boundary conditions, which we will get to know later)\n",
    "```python\n",
    "def inv_tridiag_121_b(n,b0,b1):\n",
    "    \"\"\"\n",
    "    Make a matrix inversion of a tridiagonal matrix of size \n",
    "    n x n with the boundary conditions b (array of size n)\n",
    "    returns the array u\n",
    "    \"\"\"\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 10\n",
    "dx = n\n",
    "b0 = 10.\n",
    "b1 = 20.\n",
    "max_iter = 10000\n",
    "epsilon_tol = 1e-6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def tridiag_121(n):\n",
    "    \"\"\"\n",
    "    Create a tridiagonal matrix of 1 -2 1 form of size n x n\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    A = np.zeros(shape=(n,n))\n",
    "    i,j = np.indices(A.shape)\n",
    "    A[i==j] = -2.\n",
    "    A[i==j-1] = 1.\n",
    "    A[i==j+1] = 1.\n",
    "    return  A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numpy.testing import assert_array_equal\n",
    "assert_array_equal(tridiag_121(10)[7],[ 0.,  0.,  0.,  0.,  0.,  0.,  1., -2.,  1.,  0.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def inv_tridiag_121_b(n,b0,b1):\n",
    "    \"\"\"\n",
    "    Make a matrix inversion of a tridiagonal matrix of size n x n with the \n",
    "    boundary conditions b (also of size n x n)\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    A = tridiag_121(n)\n",
    "    b = np.zeros(n)\n",
    "    b[0] = -b0\n",
    "    b[-1] =-b1\n",
    "    return np.dot(np.linalg.inv(A), b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numpy.testing import assert_array_almost_equal\n",
    "assert_array_almost_equal(inv_tridiag_121_b(10,10.,20.),[ 10.90909091,  11.81818182,  12.72727273,  13.63636364,\n",
    "        14.54545455,  15.45454545,  16.36363636,  17.27272727,\n",
    "        18.18181818,  19.09090909])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, create a plot of the results: What do you observe? **(Plot: 3 Points, Text answer: 2 Points)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEZCAYAAAB1mUk3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcHHWd//HXezgCkSNcyhGSCAICggkDyGUYL0DxXmSR\nVQkoh+wKEa/V1Y0HeOAqAVyVVQggggewuPkph0eGexWGGYTAkrhAQI5IBDSQlQDz+f1R1aHSdM/0\ndLq6errez8djHtNVXcfn01Vd367vt6q+igjMzKzceooOwMzMiufCwMzMXBiYmZkLAzMzw4WBmZnh\nwsDMzHBh0DKSpkoaltTUZyrpSElXtTquUda5o6RBSX+R9E/tXHedeOZI+kH6eltJf5WkouNqlKR5\nkr7YxvV9WdJJ7VpfjfUfIOnuJue9U9LMVse0JiQtlzStBcu5VNLBax5Re7kwqJLu4DdKelLSMknX\nS+ptcPaGbtqoVXBExMURcUhTQTfvk8BvImLjiPhWm9ddTwBExIMRsVHkcCOMpPskvb7Vy20nSZsD\n7wfOyWn5w5K2G2maiLghInZuZvkR8aqIuK656PIRERtGxP0tWNTXgNNasJy2cmGQIWlDYD5wJrAJ\nsA3wBeCZVq+K5KBX9K/eqcDCem82e5bTDuPpjGFN1cl1FvCLiGj1vlkxYiEsaa2c1pubdsUcEbcA\nG0raox3ra5WO/bIXZEcgIuInkXgmIn4VEXdC8qWU9FlJ90t6VNL5kjaqtaDqX59pFciF6eC16f8n\n06qQ10g6StL1men3k/Q7SU9I+q2kfTPvLZD0RUk3pPNfJWnT9L0Jkn6QntVU5t2iRny/Bl4H/Hu6\njFek1RzflvRzScuBPkkbSbpQ0p/SnP4ls4yj0hi+ma7rD5L2Tcc/kH5GH6j3YUuaJqk/raa6Gtg8\n895qZ09pzqem63saeHka27mSHpb0oKQvZQ+cko6VdFea352SpqfbYAowPx3/8RpxbSZpfprTnyVd\nm3nvlWksT0i6Q9Lb6uR2l6S3ZIbXSj/D6enwPukZ6BNKquoOrNq+q+VaYxVv5oX9CEkHpp/BJyQt\nlfSQpHdIerOke9L94dOZ6feSdFO6/ocknS1p7fS9a0l+qPw+/Yzek1n+JyU9ApxXGZfOs136WVXy\n2zrNt2ZVkDLfDyXfjR9LuiBd3x1KD6Tp+n5aNe+ZkuamrzeS9P1a+0DV/rkMmCNp+3SfezKN75LM\ncledDWn0/f56SV+X9Lik/5VUfVZ/LXBordw7VkT4L/0DNgQeA84HDgEmVb1/DLCI5Bf1ROAy4ML0\nvanA80BPOnwf8PrMvHNqTKvM+0cB16WvNwEeB44kKbCPSIc3Sd9fACwGtgcmpMNfTt87DvhZOl7A\nDGCDOvkuAI7JDM8DngD2SYcnABcC/5nmOxW4Bzg6E/NK4APpur4ELAHOBtYB3gT8FZhYZ/03AV9P\np31tOm29z3MBcD/wyvQzWTuN69vAeiQFyX8Dx6bTvwd4ENgjHd4O2DazbV43wn7w5XS5PcBawP7p\n+LXTz/1T6evXpTHvkPn8vpi+/hxwUWaZhwIL09fbAMuAg9PhN6TDm9XJda0aMf4J6M0MHwg8C/xL\nGvOH0mkuSrfdLsAKYGo6/R7A3ul2m0JyhnhSZnnDwMtrLP/L6faakI57IDPNB4E7gfWBq4GvjfAZ\nr/p+kHw3VgAHp/F8Gbg5fW8K8BTwknS4B3gY2CsdHmkfOCqN+cR0vvWAi4FPp++vC+yXiel5YLv0\n9Wj7/TMkxwMBJwAPVeX3UeDSoo9pYzr+FR1Ap/0BOwHnAQ+QHOh+BmyRvvcr4ITMtDum0/TQXGHQ\nk3k/Wxi8D/jvqrhuAj6Qvl4AfCbz3odJqgwAjgZuAHZrINdahcH5meGedKffKTPuOJJ2hkrM92Te\ne1Wa1+aZccuA3Wuse9v0s1s/M+6H9T6jNNbPZ6Z9KfA3YEJm3BHAr9PXVwEfqZP3atumxvtfSA8E\n21eNPwB4uGrcxcC/Zj6/SmGwPUlBsV46fBHw2fT1J4ELqpZzFfD+WrnWiXElsGNm+EDgadIfGMAG\nJAf0PTPT3Aq8vc7yTgYuywwPkx4YM8v/G7BO1bgHqpZzBfB7YCg77UjbgOS7cU3mvZ2BpzPD1wHv\nS1+/CVicvn5ZnX0gu3/eX7XeC4DvAtvUiGmY5EdDI/v9osx766f76ksz4z4E/Gq072An/bmaqEpE\n3BMRx0TEFJKD29bA3PTtrUl++VYsIfmF+LIWh1G9nsq6tskMP5p5vYLkyw/wA5JfZT+S9EdJX9XY\n6kofzLzenCS/B0aIY2nm9f8BRMSyqnEb8GJbA09ExP9VLbvR2KaS/EJ9JD1Vf4LkS16pEtsW+N9R\nllfP6em81yip+vpUJuYHq6at/jwAiIj/Be4C3iZpfeDtJIVdJfbD07grse8PbJlZRPV6qj1Bciab\n9edIj0Sk24Lk7IDMuA0AJO2QVoU9IulJkgbPzRnZYxHx7CjTfB/YFTi7gWmzqvfn9fRCm9UlwHvT\n1+8lKYAhOWuotQ9k86j+HD9BcrD/XVoddXSNWBrZ71fFm+7DYvX9fEPgyRrL7lguDEYQEYtIqoxe\nlY56mOSLXDGV5DR0KS/2NMkpZkX2ix6M7GFgWtW4KcBDo8xHRDwXEV+KiF2B/YC3kVTjNCob2zKS\n/KpzHjWOBjwCbJIeKCumjCG2B0l+FW4WEZtGxCYRMSkids+8v30Dy3nxmxFPR8THI2J7koP4KZJe\nR7JdqmMcabv8iKSq7x0kVUT3ZWK7MI27EvuGEfH1RmMk+fW94yjTjOQ7wN0kZz+TSKqXRmuUHzEm\nSS8h+eF0LvB5SZPWIL6sn5K0X20DvIsXCoPR9oEXxRwRf4qI4yJiG5LqnW/rxVdNtWK/3xm4fQzT\nF86FQYaknSSdku50SNqW5JfIzekklwAfVdLwuQHJr6kfRcRwZRGZxQ0BR0haW9KewGGZ9x4jOSWt\nd7D6BbCDpCPShse/J9m55jeQQ5+kV6W/qp4i2amHR5mtpjSvnwCnSdpA0lSSutAfjBRCg8t+gKTa\n4guS1pF0AEnB1dCyIuJR4BrgDEkbKrFdpsHy+8DHMw2R26fbE5LCu+5lk5IOlVTZNsuB50g+w98C\nT6eNmmtL6gPeSrJf1PIj4CCSaryLM+MvIjljOEhSj6T10sbYrevFVMMvgL4xTF9tQ+CvEbFC0ivT\nGLMeZYTPqI6zgN9FxHFpfGty2euqbZ+eaV5LUg13b0Tck44fbR948UKlwyrfb5Jf7sNUfT+a3O+r\nHQhcOYbpC+fCYHXLgdcAv1VyNc1NJL/AKlecnEeyQ1xHUo2wAsje9JP9FfI54BUkDb9zeKGKoHJa\neRpwY3p6u3c2iIh4nOQg83GSXykfBw6NiCdqrKfalsClwF9IGgUXUH8nrl5OreWeRJLnvSR5XxQR\n80ZYfyPLrDgS2Af4M8nndcEI89ZazgdIGgHvIvmcf0p6BhYRl5J8xhdL+itJG8Cm6XxfAT6Xfvan\n1FjuDsCv0n3gRuDfI+LatNrjbcBbSLbLt0jq+RfXijE9WN2c5vjjzPg/kpwtfIbkh8ESkm3cU2s5\ndVwIvFnShBGmGWlbfBz4h/SzOYek4Mr6PHBh+hkdxigkvZ2k4DsxHXUKMEPSe+vMMlqO1e9fTNLQ\n/sOq8XX3gTr2Ivl+/5WkfeOkeOHeguw6m97vJe0FLI+IW0eYvuNUGpvyWbg0mWSnfRlJ6fu9iDhL\n0iYkX46pJFdNHB4Rf8ktELMuJOlU4E8RcVbRsdgLJF0KfD8i2vpEgTWVd2GwJbBlRAyl1SoDJL+I\njiZp7Do9bZzbJCL+ObdAzMxsRLlWE0XEoxExlL5+iqTBajJJgVCpErgAeGeecZiZ2chyPTNYbUXJ\nA6D6Sa7MeTAiNsm893hEbFp7TjMzy1tbGpDTKqJLgZPTM4SxNDKamVnO1s57BUqed3Ip8IOI+Fk6\neqmkl0XE0rRd4U915nUhYWbWhIgY08Mc23FmcB5wV0ScmRn3XyRPXYTk1u6fVc9UUcRt2XPmzCn8\n1nDn7Jyds3Nu9q8ZuZ4ZSNof+AfgDkmDJNVBnyF53vdPJB1Dco314XnGMVb3339/0SG0nXMuB+ds\n9eRaGETEjSRPUKzljXmu28zMGuc7kGuYNWtW0SG0nXMuB+ds9bTt0tJmSIpOjs/MrBNJIjqwAXnc\n6e/vLzqEtnPO5eCcrR4XBmZm5moiM7Nu42oiMzNriguDGspYx+icy8E5Wz0uDMzMzG0GZmbdxm0G\nZmbWFBcGNZSxjtE5l4NztnpcGJiZmdsMzMxaaXh4mMHBQQBmzJhBT0/7f3O7zcDMrEALBweZ3dvL\nkpkzWTJzJrN7e1mYFgydzoVBDWWsY3TO5eCc8zM8PMw5xxzD3KEh3r1iBe9esYK5Q0Occ8wxDA8P\ntyWGNeHCwMysBQYHB+lbtGi1g2oPcOCiRauqjTqZ2wzMzFpgYGCAJTNn8u4VK1Ybf9nEiUy77jp6\ne3vbFovbDMzMCjID6B8eJlshNAxcu+OOzJgxo6CoGufCoAbXq5aDcy6HtuT80EP0HHoox3/608x+\n9au5bOJELps4kZNf/WqOP++8Qq4oGqtc+0A2MyuFbbaBP/yBXTfYgLmf/eyqNoIzC7q0tBluMzAz\n6zJuMzAzy9vDDxcdQS5cGNTgetVycM7l0LKcV66EOXNgjz3g8cdbs8wO4sLAzGw0g4Ow114wMAC3\n3Qabblp0RC3nNgMzs3pWroTTToPvfAf+7d/g/e8HjakqvhDNtBn4aiIzs3qWL4f774ehIdh666Kj\nyZWriWpwvWo5OOdyWKOcN9sMLrig6wsCcGFgZma4zcDMLGkb+N734LjjYJ11io5mjfk+AzOzsapc\nKXTllfDUU0VHUxgXBjW4XrUcnHM51M25ct/AwQfDxz4G8+fDJpu0NbZO4quJzKx8HnsM3vhG2Hbb\nUlwp1Ai3GZhZVxqxL+IIuOoqOOSQcXHfwFi5zcDMjAb6IpbgzW/uyoKgWS4ManC9ajk45+5U3Rfx\npuOsL+KiuDAws64y3vsiLkquhYGkcyUtlfT7zLhXS7pZ0qCk30naM88YmtHX11d0CG3nnMuhlDkX\nHcA4kfeZwTzg4KpxpwNzImIGMAf4es4xmFmJjPe+iIuSa2EQETcAT1SNHgY2Tl9PAh7KM4ZmlKFe\ntZpzLoeuz3nZMnre8Y7V+iL+woQJ46ov4qIUcZ/BR4GrJX0DELBfATGYWTfafHNYvJhd119/VV/E\ny269lc8de6wLglEUURh8GDg5Iq6QdBhwHvCmehPPmjWLadOmATBp0iSmT5++qt6z8ivHw2s+3NfX\n11HxtGO4Mq5T4mnXcDb3Tognr+HrrrsOgOOPP74j4slzuL+/n/PPPx9g1fFyrHK/6UzSVGB+ROye\nDj8ZEZMy7/8lIjauM69vOjOz2pYsgalTi46iI3XqTWdK/yoeknQggKQ3AIvaEMOYVP+CKgPnXA5d\nkXPlmUJ77QWPPDLq5F2RcxvkWk0k6WKSK7s2k/QAydVDxwJnSVoL+BtwXJ4xmFkXGRyEWbNeeKbQ\nVlsVHVHX8LOJzKzzjdO+iIviPpDNrDutXAmPPuonjObI11rVUMY6RudcDuM25w02gHPOaaogGLc5\nt5kLAzMzc5uBmXWQlSvhW9+CD38Y1l+/6GjGrU69tNTMbHSVvoh/8xt4+umioykdFwY1lLGO0TmX\nQ0fmXKsv4s03b9niOzLnDuSricysOH/5C8yc6b6IO4DbDMwsNyP2Q1yxYAH09fm+gRZym4GZdYxR\n+yGueN3rXBB0ABcGNZSxjtE5l0O7cq7uh/jdBfZDXMbt3AwXBmbWcnX7If6f/3E/xB3KbQZm1nID\nAwMsmTmTd69Ysdr4yyZOZNp119Hb21tQZOXgNgMz6wgzZsygf/Jk90M8jrgwqKGMdYzOuRzalXPP\n009z/FNPMXvKFC6bOJHLJk4srB/iMm7nZvg+AzNrvQ03ZNf77mPu2muvaiM4s96lpdYR3GZgZtZl\n3GZgZu23eDH4R9u458KghjLWMTrncmhpzpVnCu2/P9x3X+uW22Jl3M7NcGFgZmNXecLowEDyTKHt\ntis6IltDbjMws8a5L+JxwX0gm1n+nnzSTxjtQq4mqqGMdYzOuRzWOOd114UzzxxXBUEZt3MzXBiY\nmZnbDMyshpUr4Ywz4IQTYOONi47Gxsj3GZjZmqtcKXT99fDMM0VHY23iwqCGMtYxOudyGDHnWn0R\nv/SlbYstL2Xczs3w1URmBitWwH77weTJvlKopNxmYGaJm2+GffbxfQNdoJk2AxcGZiXQUMf01jVy\naUCW9B5JG6avPyvpckl7NBvkeFDGOkbn3L2yHdP/fP/9mb3HHrU7pu9SZdnOa6qRnwefi4jlkg4A\n3gicC3wn37DMrBWqO6af+cwzzL39ds458si2d0xvnW3UaiJJgxExQ9JXgDsi4uLKuNyDczWR2Rpx\nX8TllNd9Bg9JOgf4e+AXkiY0OJ+ZdQKfAVgDGjmoHw5cDRwcEU8CmwKfyDWqgpWxjtE5d6cZO+9M\nf8Sqjun7KV/H9GXYzq0wamEQESuAPwEHpKOeAxbnGZSZtUbPxIkcf/31zJ4+ncsmTuTaCRMK65je\nOlsjbQZzgD2BnSJiR0lbAz+NiP1HXbh0LvBWYGlE7J4Z/xHgRJKC5ecR8c915nebgVkL+NLScsnl\nPgNJQ8AM4LZKo7Gk32cP7iPMewDwFHBhZXpJfcBngLdExHOSNo+IZXXmd2Fg1qi77oKddoK11io6\nEitYXg3IK9MjcqQreUmjC4+IG4AnqkZ/GPhqRDyXTlOzIChSGesYnfM4VnmmUF8fLFo04qRdk/MY\nlDHnZjRSGPwkvZpokqRjgV8B31uDde4IzJT035IWSNpzDZZlVm7VfRHvvHPREdk41dDjKCS9CTgI\nEHB1RPyy4RVIU4H5mWqiO4DfRMTJkvYCfhwRNXvTdjWRWR3PPgunnuq+iK2m3PpATg/+DRcAo3gQ\nuDxd7i2ShiVtFhF/rjXxrFmzmDZtGgCTJk1i+vTp9PX1AS+c/nnYw6Ub7umhf/Fi+Pa36TvssOLj\n8XChw/39/Zx//vkAq46XY9VIA/Jy0vYCYF1gHeDpiNiooRVI00jODHZLh48DtomIOZJ2BH4ZEVPr\nzFvImUF/f/+qD7wsnHM5OOdyyOXMICI2zKxAwDuAfRoM6GKgD9hM0gPAHOA8YF5aXfQM8IGxBGxm\nZq3X1COs/WwiszZZuRJOPx0+9CHYcsuio7FxIpczA0nvzgz2kNyA9rcxxmZmYzU4CLNmwbbbJoWB\nWY4aubT0bZm/g4HlJFVFXavSMFMmzrmD1OqLuEVnBR2bc47KmHMzGmkzOLodgZgZSUGw776w1Vbu\ni9jaqm6bgaSzeeEqoheJiJPyCioTg9sMrHwGB2H6dN83YE1rdZvBrWsYj5lVaeiBcSV5tLR1lrpt\nBhFxwUh/7Qyy3cpYx+ic85fti3jJzJmF9EXs7Wz1NHI10RbAp4BdgPUq4yPi9TnGZdZVsn0RV36B\nvfP225l9xBHMvftuP1LaCtfIHcjXAD8GPg6cABwFPBYRn8o9OLcZWJdwX8TWTnk9wnqziDgXeDYi\nro2IYwCfFZiNlfsitg7WSGHwbPr/EUmHSppB0g9y1ypjHaNzzteMV72KfolscVBEX8TezlZPI08t\nPVXSxsDHgLOBjYCP5hqVWZfpmTAh6Yv4Qx/iwLQDmv4dduAE90VsHaKRNoMtIuKxNsVTvW63GVhX\ncV/E1g559YG8CLifpBH58oio7sYyNy4MbFy6446kL+J11y06EiupXBqQI2JH4LPArsCApP8n6X1N\nxjgulLGO0Tm3QOWZQm94Ayxc2Nplt4i3s9XT0DlqRPwuIk4B9gYeB7r6pjOzMavui9h3Eds400g1\n0UbAu4AjgO2B/wR+EhEDuQfnaiLrdM8/D1/8ovsito6SVx/ItwNXAF+MiJubisysW/X0wDrr+Amj\nNu41Uk20XUR8tEwFQRnrGJ1zkyT47GfHTUHg7Wz1NNKA7HoaM7Mu11QfyO3iNgPrGCtXwle+krQJ\nbLdd0dGYjSiXS0sl7d/IOLOuVblS6JZbYP31i47GLBeNtBmc3eC4rlHGOkbnXEOtvoi32qotseXF\n29nqqXs1kaR9gf2ALSSdknlrI2CtvAMzK9Tzz8NrXwtbbOErhawURuoD+UCgj6QPg+9m3loOzI+I\nxbkH5zYDK9Jdd8HOO/u+ARt38no20dSIWLJGkTXJhYHlwQ+Ls26XV+c250v6TfVfkzGOC2WsYyxL\nztl+iH++//7M7u1l4S23FB1W25RlO2eVMedmNHIH8sczr9cD/g54Lp9wzPJT3Q/xpsDMoSFmz5zJ\n3CuvpKevr9gAzQrU1H0Gkn4XEXvnEE/1elxNZC1Ttx/idddl2o030rvnngVFZtZauTybSFK2i8se\noBfYeIyxmXWutdd2I7GVXiNtBgPAren/m0m6v/xgnkEVrYx1jGXIecaMGfTvsMOqfoj7KaYf4iKV\nYTtXK2POzRj1zCAiXt6OQMzy1tPTw/Hz5jH76KM5cPFi7nz+eS575SvdD7EZjV1auh5wInAAEMD1\nwHcj4m+5B+c2A8uBLy21bpfXfQY/IbnR7KJ01JHApIh4T1NRjoELA1sjQ0Owww7wkpcUHYlZW+V1\nn8GrIuKDEbEg/TuWpD/krlXGOsauyrnyTKGDDko6p6+jq3JukHO2ehopDG6TtE9lQNJrSBqUzTpP\ndV/E++wz+jxm1lA10d3ATsAD6agpwD0kN55FROyeW3CuJrJGRcDnP+++iM3Irw/kQ5qMB0nnAm8F\nllYXGpI+Bnwd2DwiHm92HWZAcuDfZBM/YdSsSY1UE50aEUuyf9lxo8w7Dzi4eqSkycCbgEIegDea\nMtYxdkXOs2ePqSDoipzHyDlbPY0UBqs1Fktam+Qu5FFFxA3AEzXeOgP4RCPLMDOz/I3Un8Gngc8A\n6wMrgEr900rgPyLi0w2tQJpK0v/B7unw24G+iDhF0n1Ab71qIrcZ2IusXAmnnQaHHQa77VZ0NGYd\nqaWXlkbEVyJiQ+DrEbFRRGyY/m3WaEFQI8D1SQqYOdnRzSzLSqhypdBtt8FmmxUdjVlXaaQB+UpJ\nM6tHRsR1Taxve2AacLskAZOBAUl7R8Sfas0wa9Yspk2bBsCkSZOYPn06femjhit1ga0erozLa/md\nOFyde9HxrDa8335w6qn0n3UWnHgifaedBtIaL3/u3Llt2Z86aXhoaIjZs2d3TDztGK6M65R48vr+\nnn/++QCrjpdj1cilpfMzg+sBewMDEfH6hlYgTSOpJnrROX1aTbRHRNRqVyismqi/v3/VB14WHZtz\nRNIX8SabwDnntPRKoY7NOUfOuRxyeRxFjZVsC8yNiL9rYNqLSfpR3gxYCsyJiHmZ9+8F9nSbgY3o\nD3+A7bf3fQNmDWpXYSBgYUTsMqYZm+DCwMxs7HJ5NpGksyWdlf59i+Sppbc1G+R4kK1rLIsich4e\nHmZgYICBgQGGh4eTK4XaWPh7O5dDGXNuRiMNyNnnED0HXBIRN+YUj5XEwsFBzjnmGPoWLQLggsmT\nOf7559n17LPhzW8uODqz8mm0P4NXpIN/aEc/Bpl1u5qoCw0PDzO7t3dVx/SQ9Dg2e8oU5t57Lz1r\nrVVkeGbjXkuriSStLel04I/ABcCFwIOSTpe0zpqFamU2ODhI36JFq+18PcCBy5YxODRUVFhmpTZS\nm8HXgU2Bl0dEb0TsQXKfwCTg39oRXFHKWMfonMvBOVs9IxUGbwWOjYjllRER8Vfgw8Bb8g7MuteM\nGTPo33HHVR3TQ/k6pjfrNCM9m2hRROw41vdayW0G3avSgHxg2oDcv8MOnDBvHru6MDBbYy29z0DS\nFcDlEXFh1fj3AYdHxNubjrTR4FwYdIfBQZg6FTbddLXR7pjeLB+tvs/gH4F/lNQv6Rvp37XASSRV\nRV2rjHWMueS8ciX867/CwQfDwoUverunp4fe3l56e3sLKQi8ncuhjDk3o+59BhHxEPAaSa/nhT4N\nfhERv25LZDa+DQ7CrFkwZYp7HzMbB8b8OIp2cjXROBQBX/gCfPvb7ovYrCB59YFs1jgJttnGZwNm\n44xb7GooYx1jS3M+9thxURB4O5dDGXNuhgsDMzNzm4E1qdIX8SGHwL77Fh2NmWXk8ghrsxcZGoK9\n94aBgeT+ATMb91wY1FDGOsaGcl65EubMgYMOglNOgfnzx0XbQD3ezuVQxpyb4auJrHGHHgoTJvhK\nIbMu5DYDa9ySJclNZL5vwKyjtaUP5HZyYWBmNnZuQG6RMtUxVvohPuecc5J+iCFpG3j++WIDa4My\nbecK52z1uDAosYWDg8zu7WXJzJksPflkZvf2svCSS2CvveDyy4sOz8zayNVEJVW3H+K11mLu979P\nz1FHuW3AbJxyNZE1rG4/xOuuy+Buu7kgMCsZFwY1lLGOsb/yokSFQCm3s3O2OlwYlJT7ITazLLcZ\nlJj7ITbrTr7PwEY2NAQvfelqdw+7H2Kz7uMG5BbpujrG7DOF7r57tbcq/RAvX768dAVB123nBjhn\nq8fPJup2Q0NJX8STJ/uZQmZWl6uJutmpp8JZZ7kvYrOScR/ItrrttvPZgJk1pFyVxA3qmjrGI49s\nuCDompzHwDmXQxlzboYLAzMzc5vBuFfpi3jffZP+iM2s9Hxpadlk+yLeffeiozGzcSzXwkDSuZKW\nSvp9Ztzpku6WNCTpMkkb5RlDMzq+jjGHvog7PuccOOdyKGPOzcj7zGAecHDVuGuAXSNiOrAY+HTO\nMXSfww5LzgaGhuADH/Alo2a2xnJvM5A0FZgfES+qx5D0TuDvIuL9deZ1m0EtjzwCW27pQsDMahqP\nbQbHAFcWHMP4s9VWLgjMrKUKu+lM0r8Az0bExSNNN2vWLKZNmwbApEmTmD59On19fcALdYGtHq6M\ny2v5tYaHh4f53ve+B8Cxxx5LT08P/b/8JUTQd9BBua+/Ove819cJw3Pnzm3L/tRJw0NDQ8yePbtj\n4mnHcGW/wwerAAAJ40lEQVRcp8ST1/f3/PPPB1h1vByziMj1D5gK/L5q3CzgRmDCKPNGERYsWNDW\n9d15223xkenT47KJE+OyiRPjI9Onx50XXxyx++4R8+a1JYZ259wJnHM5lDHn9Ng5pmN1O9oMppG0\nGeyWDh8CfAOYGRF/HmXeyDu+orkvYjNrtY5rM5B0MXATsKOkByQdDZwNbAD8UtJtkr6dZwydzn0R\nm1knyLUwiIgjI2LriJgQEVMiYl5E7BARUyNij/TvxDxjaEa2rrEwbS4EOiLnNnPO5VDGnJtR9NVE\npee+iM2sE/jZRB3AfRGbWSu5D+TxYmgINtwQtt9+1Sj3RWxmrdJxDcjjVW51jNlnCi1evNpblb6I\ne3t7CykIyliv6pzLoYw5N8M9nbWL+yI2sw7maqJ2+MpX4Iwz3BexmbWF+0DuVLvs4rMBM+tobjOo\noeV1jO94R8cXBGWsV3XO5VDGnJvhwsDMzNxm0DKVvoh33RUOP7zoaMysxHxpaVGyfREfcEDR0ZiZ\njZkLgxoarmPMoS/iopSxXtU5l0MZc26GryZaE+9/Pzz9tK8UMrNxz20Ga2LZMthsM983YGYdxc8m\nMjMzNyC3yovqGFeuhBUrComlXcpYr+qcy6GMOTfDhUHG8PAwAwMD3HPPPQwPpz0MVK4Umjev2ODM\nzHLkaqJUpU+BvkqfAq94Bce/5jXsesUVfqaQmY0rbjNoUt1O6TfaiLkLF9IzeXLuMZiZtYrbDJpU\n3Sl9P2mn9M89x+DSpcUF1kZlrFd1zuVQxpyb4cLAzMxcTQQjVBNNn87cgQF3QWlm44rbDNaAO6U3\ns27hNoM1sOuMGcwdGGDaddex7Jvf5MzbbitVQVDGelXnXA5lzLkZfjZRRqVT+uXLl7tqyMxKxdVE\nZmZdxtVEZmbWFBcGNZSxjtE5l4NztnpcGJiZmdsMzMy6jdsMzMysKS4MaihjHaNzLgfnbPW4MDAz\nM7cZmJl1G7cZmJlZU3ItDCSdK2mppN9nxm0i6RpJ90i6WtLGecbQjDLWMTrncnDOVk/eZwbzgIOr\nxv0z8KuI2An4DfDpnGMYs6GhoaJDaDvnXA7O2erJtTCIiBuAJ6pGvwO4IH19AfDOPGNoxpNPPll0\nCG3nnMvBOVs9RbQZvDQilgJExKPASwuIwczMMjqhAbnjLhe6//77iw6h7ZxzOThnqyf3S0slTQXm\nR8Tu6fDdQF9ELJW0JbAgInauM2/HFRRmZuPBWC8tbUfnNkr/Kv4LmAV8DTgK+Fm9GceajJmZNSfX\nMwNJFwN9wGbAUmAOcAXwU2BbYAlweES4hcfMrEAdfQeymZm1Ryc0IHcUSYdI+h9JiyR9quh48iZp\nsqTfSFoo6Q5JJxUdUztI6pF0m6T/KjqWdpG0saSfSro73d6vKTqmPEn6qKQ7Jf1e0g8lrVt0TK3W\nyht7XRhkSOoBvkVyo9yuwHslvbLYqHL3HHBKROwK7Av8YwlyBjgZuKvoINrsTOAX6QUbrwbuLjie\n3EjaGvgIsEd68crawBHFRpWLlt3Y68JgdXsDiyNiSUQ8C/yI5Ca5rhURj0bEUPr6KZIDxDbFRpUv\nSZOBtwDfLzqWdpG0EfDaiJgHEBHPRcRfCw4rb2sBL5G0NjAReLjgeFqulTf2ujBY3TbAg5nhP9Ll\nB8YsSdOA6cBvi40kd2cAn6AD73HJ0cuBZZLmpdVj/yFp/aKDyktEPAx8A3gAeAh4MiJ+VWxUbdPU\njb0uDAwASRsAlwInp2cIXUnSocDS9Gyo+rLnbrY2sAfw7xGxB7CCpDqhK0maRPILeSqwNbCBpCOL\njaowDf3ocWGwuoeAKZnhyem4rpaeRl8K/CAi6t730SX2B94u6V7gEuB1ki4sOKZ2+CPwYETcmg5f\nSlI4dKs3AvdGxOMR8TxwObBfwTG1y1JJLwNIb+z9UyMzuTBY3S3AKyRNTa88OILkJrludx5wV0Sc\nWXQgeYuIz0TElIjYjmT7/iYiPlB0XHlLqw0elLRjOuoNdHcD+gPAPpLWkySSfLu1wbzejb0wyo29\nWe24A3nciIjnJf0TcA1JQXluRHTrDgSApP2BfwDukDRIckr5mYi4qtjILAcnAT+UtA5wL3B0wfHk\nJiJ+J+lSYBB4Nv3/H8VG1XrZG3slPUByY+9XgZ9KOob0xt6GluWbzszMzNVEZmbmwsDMzFwYmJkZ\nLgzMzAwXBmZmhgsDMzPDhYF1IEnPp8/PuVPSoKRTMu/1Spqb47oPlLRvZvh4Se9rYjmvlvTmzPDb\nJH2yRTHe18A0CyRNGW06swrfdGad6On0+TlI2hy4RNJGEfH5iBgABtZk4ZLWSh9RUEsf8BRwM0BE\nnNPkaqYDewJXpsuZD8xvclnVfHOQtZzPDKyjRcQy4Djgn2DVL/f5mdeD6VnEgKSXpOM/lXZoMijp\ny+m4BZLOkHQLcJKkzSVdKum36d++kqYCJwCz02XuL2lO5cxE0klppzBD6Z2fSNpL0k3p+m+QtEN6\nh+8XgcPT5bxH0lGSzk7nmSrp1+lyfpk+Upv0iaJnSrpR0h8kvbvOx/JYJv8FmQ5rfpCZ5s9AvQLP\n7EV8ZmAdLyLuU9Iz2RaVUen/jwEnRsTNkiYCf5N0CPA2YK+IeCZ9emXFOhGxF4CkHwLfjIibJG0L\nXB0Ru0j6LrA8Ir6ZTvfGzPyfAqZFxLNp/wCQPO/mgIgYlvQG4CsRcZikfwV6I+KkdDlHZeI+G5gX\nERdJOjodflf63pYRsb+knUmeMXN5jc8j20PZdGAX4FHgRkn7RcRNEXFYY5+uWcKFgY0XtR41fSNw\nRnpgvzwiHkoP3vMi4hmAiHgyM/2PM6/fCOycPsQMkkccTxwlhtuBiyVdAVyRjpsEXChpB5KDfSPf\nqX154eD/A+BrmfeuSOO+W1Ijz6H/XUQ8AiBpCJgG3NTAfGarcTWRdTxJ2wHPRcRj2fER8TXgg8D6\nwA2SdhplUU9nFwu8JiJmpH9TImLFKPMfStIt6h7ALUq6Sf0SyZNPdyM5I1mvgZRGqvN/pirG0WSn\nfx7/wLMmuTCwTrTqIJhWDX2HpCpl9Ymk7SJiYUScDtwK7AT8EjhaaS9ekjaps45rSPpBrizr1enL\n5cBGNeeAKRFxLUmnMBsBGwAb80KfF9mngI60nJuA96av3wdcX2e6snS8Yx3AhYF1ovUql5aSHLSv\niogv1phutqQ70uqRlcCVEXE1SV37rZJuI2lXgBf/Gj8Z2FPS7el6jk/HzwfeVWlArsynpAOgiyTd\nTnI105lpH8KnA1+VNMDq36cFwC6VBuSqdZ9EUmANkTw+vFIoVcc41quGfJWRNc2PsDYzM58ZmJmZ\nCwMzM8OFgZmZ4cLAzMxwYWBmZrgwMDMzXBiYmRkuDMzMDPj/JGY1GoKUKkUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa62dba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# YOUR CODE HERE\n",
    "plt.plot(inv_tridiag_121_b(10,10.,20.),'ro--')\n",
    "plt.xlim([-1,10])\n",
    "plt.title('Solutions from direct solver (matrix inversion)')\n",
    "plt.xlabel(\"Discretisation 'n'\")\n",
    "plt.ylabel(\"Output values\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot show a linear trend. This linear interpolation is similar to the characteristic arose from physical problems such as Diffusion Problem where we have two heat sources of different temperatures on each of the two sides of a block --> Laplace del^2 T = 0 --> meaning that the edge values take on the values of boundary conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We have seen the previous example before and discussed it in the solution notebook to the last exercise. We will now have a look at other ways to solve the equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterative solvers\n",
    "\n",
    "We didn't discuss iterative methods in the lecture yet (you'll have to wait until next Monday for that), but the basic idea is pretty simple: instead of solving the system of equations at once (as with direct solvers), we start with an initial guess of the solution and then use the algebraic equations to iteratively refine our guess until we reach a solution that is close to the true result.\n",
    "\n",
    "So, let's first go back one step from the matrix equation above and look again at the underlying equations. We will start with the discretisation of a 1-D Poisson equation (more on this equation as well in the next lecture!) into a central finite difference scheme, considering constant $\\Delta x$.\n",
    "\n",
    "The partial differential equation is:\n",
    "\n",
    "$$\\frac{d^2 u}{d x^2} = q$$\n",
    "\n",
    "And the central difference discretisation:\n",
    "\n",
    "$$\\frac{u_{i-1} - 2 u_i + u_{i+1}}{\\Delta x^2} = q_i$$\n",
    "\n",
    "We can solve the discrete equation for $u_i$:\n",
    "\n",
    "$$u_i = \\frac{1}{2} \\left( u_{i-1} + u_{i+1} - \\Delta x^2 q_i\\right)$$\n",
    "\n",
    "If we now consider iterative solution methods, then we basically want to use this expression to obtain a new estimate for $u$ at a time $n+1$ from the values at the previous iteration step $n$. There are two standard ways of doing this that we are now implementing in this exercise.\n",
    "\n",
    "In addition to defining the iteration methods to solve the problem, we need a way to determine if we converge towards a solution. As we discussed during the lecture, in a general numerical modelling context, we can not actually determine the true convergence because we do not know the true answer. So, what we usually do instead is to check if the solutions do not change anymore from iteration to iteration. This check can be performed in several ways (see lecture). We will here use the L2-norm for convergence, calculated as:\n",
    "\n",
    "$$\\epsilon = \\sqrt{\\sum_{i,j} \\left(u_{i,j}^{n+1} - u_{i,j}^n \\right)^2}$$\n",
    "\n",
    "We stop the iterations when we reach a certain tolerance level so that\n",
    "\n",
    "$$\\epsilon < \\epsilon_{tol}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Jacobi method\n",
    "##### (10 Points)\n",
    "\n",
    "The straight-forward way to solve for $u_i$ at a new iteration step is to use the values from the previous iteration, i.e.:\n",
    "\n",
    "$$u_i^{n+1} = \\frac{1}{2} \\left( u_{i-1}^n + u_{i+1}^n - \\Delta x^2 q_i\\right)$$\n",
    "\n",
    "This method is called the Jacobi method.\n",
    "\n",
    "**Exercise**: solve the 1-D Laplace problem of before (note: no sources/sinks, so $q$ = 0!) with the Jacobi method with the following steps:\n",
    "\n",
    "- Write a function for a single iteration **(2 Points)**  \n",
    "```python\n",
    "    def jacobi_iter(u,q,dx):\n",
    "```\n",
    "- Write a function to update a value $u_i$ **(3 Points)**\n",
    "\n",
    "These steps will take some work, so here the base construct that you could use:\n",
    "```python\n",
    "def solve_laplace_jac(b0,b1,dx,max_iter,epsilon_tol):\n",
    "    # initialize u with all ones\n",
    "    u = \n",
    "    # Set boundary conditions b0 and b1 at beginning and end of u\n",
    "    # set parameters\n",
    "    # no sources / sinks, so set q to 0\n",
    "    # do a for loop for jacobi iterations\n",
    "        # copy results from previous iteration\n",
    "        # perform iteration\n",
    "        # check for convergence by calculating epsilon\n",
    "        # if epsilon is smaller than epsilon_tol, break for loop\n",
    "    return u,i # i being the current iteration\n",
    "```\n",
    "- Apply this function __solve_laplace_jac__ iteratively until a convergence criteria is met\n",
    "- The function should return u and the current iteration i_jac:\n",
    "    - *u,i_jac = solve_laplace_jac(b0,b1,dx,max_iter,epsilon_tol)*\n",
    "- Plot the results **(3 Points)**\n",
    "- Compare them to the direct solution **(2 Points)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# First step: define the function for a single iteration:\n",
    "\n",
    "def jacobi_iter(u, q, dx):\n",
    "    \"\"\"Function to perform the Jacobi iteration for all values in the array u\n",
    "    \n",
    "    u = 1-D np.array with values of previous iteration\n",
    "    q = 1-D np.array with source/sink terms (not required for Laplace solution!)\n",
    "    dx = float : node spacing (not required for Laplace solution!)\n",
    "    \"\"\"   \n",
    "    # Hint: do not forget to create a copy of u before you overwrite the values!!\n",
    "    \n",
    "    # YOUR CODE HERE\n",
    "    u_iter=u.copy()\n",
    "    for i in range(1,dx-1):\n",
    "        u[i] = 0.5 * (u_iter[i-1] + u_iter[i+1] - (dx**2 * q[i]))\n",
    "    return u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 30.        ,  32.22221993,  34.44444006,  36.66666086,\n",
       "         38.88888218,  41.11110451,  43.33332743,  45.55555125,\n",
       "         47.77777545,  50.        ]), 253)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def solve_laplace_jac(b0,b1,dx,max_iter,epsilon_tol):\n",
    "    \"\"\"\n",
    "    Function to solve the 1D problem in iterations until \n",
    "    convergence is met\n",
    "    Initialize u within this function!\n",
    "    b0 = boundary condition 1\n",
    "    b1 = boundaray condition 2\n",
    "    dx = discretization (= n)\n",
    "    max_iter = maximum iteration number\n",
    "    epsilon_tol = convergence limit\n",
    "    \"\"\"\n",
    "    \n",
    "    # YOUR CODE HERE\n",
    "    # initialize u with all ones\n",
    "    u = np.ones(dx)\n",
    "    \n",
    "    # Set boundary conditions b0 and b1 at beginning and end of u\n",
    "    u[0] = b0\n",
    "    u[-1] = b1\n",
    "    \n",
    "    # set parameters\n",
    "    \n",
    "   \n",
    "    # no sources / sinks, so set q to 0\n",
    "    q = np.zeros(dx)\n",
    "    \n",
    "        \n",
    "    # do a for loop for jacobi iterations\n",
    "    for iter in range(max_iter):\n",
    "\n",
    "        # copy results from previous iteration\n",
    "        u_previous = u.copy()\n",
    "                  \n",
    "        # perform iteration\n",
    "        u_current = jacobi_iter(u, q, dx)\n",
    "        u = u_current\n",
    "               \n",
    "        # check for convergence by calculating epsilon\n",
    "        epsilon = np.sqrt(np.sum((u - u_previous)**2))\n",
    "                      \n",
    "        # if epsilon is smaller than epsilon_tol, break for loop\n",
    "        if epsilon < epsilon_tol:\n",
    "            break\n",
    "        \n",
    "    return u, iter\n",
    "\n",
    "####################################################################\n",
    "\n",
    "solve_laplace_jac(30.,50.,10,10000,1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "assert_array_almost_equal(solve_laplace_jac(30.,50.,10,10000,1e-6)[0],[ 30.,32.22221993,34.44444006,36.66666086,\n",
    "                                                                        38.88888218, 41.11110451,  43.33332743,\n",
    "                                                                        45.55555125, 47.77777545,  50.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now: create a plot and compare the result for the direct solution you obtained before:\n",
    "*Note* Your function should return two parameters, $u\\_jac$ and the current iteration $iter\\_jac$, so be called:\n",
    "```python\n",
    "u_jac,iter_jac = solve_laplace_jac(b0,b1,dx,max_iter,epsilon_tol)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEZCAYAAACervI0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl4VOX1xz9nQJYgOygiq0tcqJI4YFU0oG3VqrhgrdYi\nIhaxSjUqihtFbVXUaqFaKmABpa4tilLbSq0EVLRImGDJT42iBBDEJQbRiEDm/P64d8IQZpI7ySx3\nkvN5nnkyd3vvd+69ed/3nve854iqYhiGYTQ/ApkWYBiGYWQGawAMwzCaKdYAGIZhNFOsATAMw2im\nWANgGIbRTLEGwDAMo5liDUAzREQuFJF/ZVpHBBEZKiLrG3H8TSIyM5maPJzzOBEpE5GvROTMdJ47\njp45InKH+/14EXknaluuiIREZIuIjBeRNiKyUEQqReTpzKlOLiKyWETGZFpHNmENQCNwK9K3RGSr\niHwsIi+KyJBM66oPVX1CVU/NtI5aeJqQEquxUNW7VfWy1MiKyx3AH1S1g6q+kOZz14mqvqaqh0Wt\nugF4RVU7qupDwE+A7kBnVT0/3fpEJCwiB6T7vMaeWAPQQETkWuAB4LfAPkAf4I/A8Ezqqg8RaZFp\nDY1E8NhYpJi+wP/F2ygikkYt9dEXKK21XKYNmAWapOfHD/fPE03g/6VuVNU+CX6ADsBWYEQd+7QC\npgIfAxuA3wN7uduGAuuB64HN7j5nAT8G3gM+B26KKmsy8FfgKeArYAVwZNT2icAH7rbVwNlR2y4G\nXsNprD7H6bleDLwatU8YGAeUARXAQ1HbAsD9wGfAGuBKd/9AnN890f29XwHvACd6vB7rauk5IGp5\njqs7B6gCdrrX/yugh3t95kXtf6Z7HSqAV4BDo7Z9BFwHrAK+BJ4EWrnbugIL3fVfAEvi/MYPXA1V\nroa9gMU4nYHXgG+AA4D9gOfdssqAX9S6p88A89wyVgEHAze6z0Q58MM6nq98oBjY4j4XTwJ3RD9f\n7vf/uFq/dc/zBPAdsN1dvsTdbwxOg/YF8E+gT637cYX7G9a46w4FFrn7vwOcV+t+PQT83T3HG0B/\nd9sSt7yv3W3nxfhtBwJFQCXwKfBk1LbjgOXuPfovcGzUtsXu72jlbj88als39351c5fPAELufq8B\nR9R6Rm5w78m3OP8DMZ/rbP9kXEA2foBT3H+gmJWgu88dwDKcSqUr8Dpwu7ttKLADuAVoAfzCfdD/\nglPJHe4+rH3d/Se7/7TnuPtfB3wItHC3nwvs634/z/3niixf7J7rCvdBbu2uWxqlNQy8ALQHerta\nTna3XY5Tme4HdAT+DVTH+u1ALrAu6tx9ov7x67se0Q1ANTEagFj7Rl2fx6I0fA2c5F6r64H3gZbu\n9o+AN4F9gU44ld5l7ra7gOnudWoBDKnj/n4UXQngVD5rcSrGANASp7J7EKeBGOhe12FRmquAH7r7\nP+re05uinokP45x7L/dcV7n7novzPMa8Rq62MbGul7t8Fk7lnutquRl4vdbz8ZJ7vVrjPKPrgFE4\nb2QDcToIh0bdr8+AoFveX4AnapXXv45r+wRuBwinMj/O/d4Zp1G/0C33Ane5c+3fCTwC/CaqzCuA\nf7jf83Ea2UGu/ovc+xnpkHwErAR6ur837nOd7R8zATWMrsDnqhquY58LcSq4L1T1C+B2nActwnbg\nLlWtxunBdQOmqmqVqv4fTsU0MGr/YlV9zt3/AaANcAyAqs5X1c3u97/iVHhHRx37sapOV9Wwqn4X\nR+/dqrpVVdfj/CPluevPA6ap6iZV3QJMqeM3V+P8w35PRFqq6jpV/cjj9YimMeaTnwJ/V9VX3Gv1\nO6AtTs8xwjRV3ayqlTg9/shv3YHT0PVX1WpVfb2ec9XWOVdV33Wfix7uOSeq6g5VXYVTKY2K2v9V\nVX3Z3f+vOM/AlKhnoq+IdIhx3mNwGrQ/uDrnA2/Vo7UuxuHc/zJXyxQgT0R6R+1zl6pWus/PGcBH\nqvqYOqwC5uM8KxGeU9Vit7zH2XWNI9R1j3fg/Pb9VXW7qi5z15+OY7p6wn2WnwLeJbbZ9UngZ1HL\nF7o6AMYCD6vqClf/PJwO1jFR+09T1Y3u763ruc5qrAFoGF8A3USkruvXE6fXEKHcXVdThrrdCZzX\nTHB6iESt2ztquWbg0z1uQ6Q8ERnlenl8KSJfAgNwKpM9jq2DzVHfq6LO3bPW8XHLUtU1QCFwG7BZ\nRJ4QkR5R5dR1PZJFT7fsiCbF0bx/1D7xfut9OGauRSLygYhMTPDc0demJ1ChqlVR68rr0PEtTqci\n+pkQdn8Gosv+uNa68hj7eaUvME1EKkSkAuf51lpaN9Ta/5jI/u4zdyHOW1WET6K+R19jL1yPUzct\nF5H/icgl7vrd7q1L7WsaYTHQVkQGi0hfnM7Ugij919XS34vdn8ea3xvnud4vgd/jW6wBaBhv4PQY\nzq5jn49xHrQIfYGNjThnTW/MHWDsBWwUkT7ATOAKVe2sqp1xBvyie1iNGXTb5J4rQp+6dlbVp1T1\nBHb99nvcvxvxfj2qcMwMEXpEfa/vt9Q+DzjXbkOMfXdDVb9W1QmqeiDOOMK1InJifcfF0bYR6CIi\n7aLW9WHPirshbGLPSq/O+1IP64BxqtrF/XRW1b1V9c2ofaJ/23qgqNb+HVR1fCM07DqR6qeqepmq\n7o9jgpzueg1tBPrV2j3mNXXfPJ7BaZh+hvNW+E2U/jtj/N5ol1itVV7t57quN+GswRqABqCqX+HY\nUf8oImeJSFsRaSkiPxaRyIPxFHCriHQTkW7AJJwBv4YSFJGzXa+Ea4BtOLbsdjg21c9FJOD2lr7X\niPPU5hngahHpKSKdcAbHYuL6m58oIq1wTFzfutrAeSX3ej1CwIXu7zkVx6YdYTPQNY5pJKL3dFdH\nSxGZgHOt3qjvh4rI6SJyoLu4FWfwtC4zX1xUdQPOmMfdItJaRI4ELqVxz0CEN4CdIvIr9zeOYHeT\nX6LMAG4WkcMBRKSjiPykjv3/DuSKyEj3/HuJyCAROcTj+T7BGSSPiYj8REQiDVwlzj0IA/8ADhaR\nC0SkhYicDxyGY8aLxZPA+TiNwBNR62cBl4vI0e752onIabUa62g9dT3XWY01AA1EVR8ArgVuxTHd\nrMMZaIq8Zv4Wx1vnbRxvghXAnXUVWc/y8zgP85fAz4FzXPvvOzheOm/i/GMNwPFqSOjn1LE8C8fb\n420cr5MXgZ1xxj9a4/SMPsPprXXHGdSExK5HIU4P/Euc3ttzNcJU38P5x/7QfX2PfjtAVcuAkThe\nKJ/h2I2Hq+rOOL81moOBl0VkK84g9R9VdUmcfeu7X7ja++Nci/nAJFVdXMf56zuHs1J1BzACuATH\nXHOeW35C5USVtwDnvj0lIpU49yh6nkjt3vDXwMk4g7Ab3c8UnPvvhduAx9z7F6uhGQz8V0S+wvl/\nukpV16pqBc74wwQcj7YJwOmq+mUcnctxPLL2w/FsiqwvxhkHeMg1eZXhOEbE/L3U/VxnNbLL5Jii\nE4isxXFVCwM7VPVoEekMPI3zOrUW+Kk7wGjEQEQmAweq6qh6d069llOBP6lq/0xrMQyjcaTjDSCM\n4/qWr6qR19QbgZdV9RAcP+0m0Zo2RdywAT92X7n3xzF9PZtpXYZhNJ50NAAS4zxn4fg94/6tazDV\nyCyC47JZgWMCKsVpBAzDyHLSYQL6EGcgpxqYoaqPiMiXrrdKZJ8KVe2SUiGGYRjGbrRMwzmGqOom\nEemO41/9Ht4G0AzDMIwUkvIGQFU3uX8/E5EFOO5qm0VkX1Xd7HpxfBrrWBGxhsEwDKMBqGq9M+pT\nOgYgIjkisrf7vR2O69j/cOLOjHZ3uxjHxTEm6oN4GdGfyZMnZ1xDNmjyqy7TZJqanK5QiOojj+RX\nHTpQTWLmlFQPAu8LvCYiIRw/9YWqughnduiPXHPQD8iiWXVr167NtIQ98KMm8Kcu0+QN0+SdjOna\nvh0mT4aTTyZw3XWMW7yYwrw85ufk1H+sS0pNQOoETKodBAp1JnT8MJXnNgzDaLKUlMDo0dCrl/O9\nZ08GAFOLiwmFQjBokKdibCZwgowePTrTEvbAj5rAn7pMkzdMk3fSqiuq18+118LChdBzVwy7QCBA\nMBj0XFzK3UAbg4ion/UZhmGkjehe/8yZu1X8tRERNNODwKmiX79+iIh9svjTr1+/jDw7RUVFGTlv\nXZgmb/hRE6RBVz29/saQjnkASae8vBx7M8huxFcpcw3Dp8Sw9SeTrDQBua83GVBkJAu7h4ZRB9u3\nw513wp/+BL/7HVx0ESTQafJqAsrKNwDDMIwmS4p7/dFk5RiAYTQUP9qRTZM3/KgJkqgrhbb+eFgD\nkAF++ctfcueddeWG8SeBQIAPP/ww0zIMo+lRUgJHHw3Fxc73UaMSMvk0mExPo67r48jbk3jrq6ur\ndcWKFbpixQqtrq6OuU99NLaMvn37atu2bbVDhw7auXNnHTJkiD788MMaDocbpCcR5s6dq8cff3yd\n+5SWlurJJ5+sXbp00c6dO+ugQYP0n//8p6fyA4GArlmzJhlS495Dw2hWfPed6q9/rdq9u+qjj6om\nqZ5w/7/qrWObzBtAKFRKMFhIQUE5BQXlBIOFhEKlaS9DRHjxxRfZsmUL5eXl3Hjjjdxzzz1ceuml\nno6vrq5O6HzRqGq93jXDhw/nlFNOYfPmzXz66af84Q9/oEOHeOl19yzfMIwkkalefzReWolMffD4\nBlBdXa15eb9SqFZQ9+Os89qLT0YZqqr9+vXT//znP7utW758uQYCAS0tLVVV1dGjR+ukSZNUVbWo\nqEh79eql99xzj/bo0UNHjRqlqqoLFy7UvLw87dSpkw4ZMkTffvvtmvLWr1+vI0aM0O7du2u3bt30\nV7/6lb7zzjvapk0bbdmype69997auXPnPbR9/vnnGggEdMuWLXH1z5w5Uw866CDt2rWrnnXWWbpx\n48aabSJS8wYwevRovfLKK/X000/X9u3b6zHHHKMffvih5+sU796mmsWLF2fkvHVhmrzhR02qDdCV\nol5/NDSnN4BQKERZ2TB2H9IIUFY21ImLkaYy4jF48GB69erFq6++GnP7J598QmVlJevWrWPmzJmE\nQiEuvfRSZs2aRUVFBePGjePMM89kx44dhMNhzjjjDPr378+6dev4+OOPueCCCzj00EN5+OGHOfbY\nY9m6dSsVFRV7nKdr164cdNBB/PznP+f555/n0093j8L9yiuvcPPNN/O3v/2NTZs20adPHy644IK4\nv+vpp5/m9ttvp7KykgMPPJBbbrmlUdfJMJo8fuj1R9EkGoB4VFU5MZFEnM9tt8Xe77bbnP2qqlKn\npWfPnjErZYAWLVpw++23s9dee9G6dWtmzZrF5ZdfzqBBgxARLrroIlq3bs2bb77J8uXL2bRpE/fe\ney9t2rShVatWHHfccZ51LF68mP79+zNhwgR69uzJ0KFDWbNmDQBPPPEEl156KQMHDmSvvfbi7rvv\n5o033mDdunUxyzrnnHMIBoMEAgF+/vOfU1JSkviFSTPDhg3LtIQ9ME3e8KMm8KgrTR4+4XCY4uJi\nz/s3iQYgPz+f3NwinPzzEcLk5S2hujq/xqBTVwNQXZ1PXt6eZeTmLiE/P7/RGj/++GO6dImd9bJ7\n9+7stddeNcvl5eXcf//9dOnShS5dutC5c2c2bNjAxo0bWb9+PX379iUQaNit69mzJ3/4wx94//33\nKS8vp127dowaNQqAjRs30rdv35p927VrR9euXfn4449jltWjR4+a7zk5OXz99dcN0mQYTZo09fqj\nxzC90iQagEAgwOzZ48jLKyQnZz45OfMZOPBqZs8e57miTEYZ8XjrrbfYuHEjJ5xwQszttQdue/fu\nzS233EJFRQUVFRV8+eWXfP3115x//vn07t2bdevWEQ6H6y2nPvbff3+uvPJKVq9eDTiNQ3n5rofn\nm2++4YsvvqBXr14Jletn/OhLbpq84UdNUIeuNPr1h8NhxoyZQUnJVKqqRng+rkk0AAD5+QMoLp7K\n0qX9WLq0HytXTiM/f0Day4hm69at/P3vf+dnP/sZF110EYcffrin48aOHcvDDz/M8uXLAaci/sc/\n/sE333zD0UcfzX777ceNN95IVVUV3333HcuWLQNg3333ZcOGDezYsSNmuZWVldx2222sWbMGVeXz\nzz9n9uzZHHvssQD87Gc/Y86cObz99tt899133HzzzRxzzDH07t27wdfAMJolabb1xx7DrJ8m0wDA\nrljYEbt0psoYPnw4HTt2pE+fPtx9991MmDCB2bNnez4+GAwya9Ysxo8fT5cuXcjNzeXRRx+t0bdw\n4ULef/99+vTpQ+/evXnmmWcAOOmkkxgwYAA9evRgn3322aPcVq1asXbtWn70ox/RsWNHjjzySNq0\nacOcOXMA+MEPfsBvfvMbRowYwf77789HH33EU089VXN8Im8Y3/ve93jyySc9758u/GhHNk3e8Jum\niL29ffv2u97IMzCbtzFYMDgjI9g9NLKZ0lCIGWPGMKysDICi3FzGTZzIgClTPMXrTxbhMDz+OCxY\nEGbNmkJWrZqK069vwvkADKOh+NGObJq84RdN4XCYGWPGMLWkhBFVVXSpqmJqSQkzRo4kXFiYtl7/\nokVw1FEwfTpcc02AOXN2jWF6xaKBGoZhJEAoFGJYWVmtGUMwtFUrQkccQTDFfv2rVsH118PatXD3\n3TBiRGR4wRnDDIVCXlMCmwnIyAx2D41spbi4mPKCAkbUmjg0PyeHfkuXJpSTtyH85S/w1VcwdixE\neY/vhuUDMAzDSAH5IjwaDnM2u2zoYWBJbi7nJGHOUH2MHJm8smwMwGhW+MWOHI1p8kbGNbkePoFT\nT2XcTTdROHAg83NyuL11a64eOJBxs2c3es5QNNu2wc6dSSsuJtYAGIZh1Ectv/4Bv/41U1eupN/S\npfSYNo1pK1cyIEm9/3AY5s2DQw6Bl19OSpFxsTEAIyPYPTSygkbm5k2URYvghhugbVu47z44/viG\nlWNjAIZhGI0hjbl5KyrgggtiefakFjMBZSF1pWZ84oknOPXUU5N+zkcffTRuLKNsIuN25BiYJm+k\nTVOCs3mToatTJ6etKS2Fc89NX4ToJtUARKZmFxcXxwyWlq4yTjzxxIRCPyRKXSEZLrzwQv71r3/F\n3f7nP/+Zww47jI4dO7Lffvtxxhln8M033zT6vIbRJMhQvP5AAC68ML5bZ8rOm97TpY7SUIjCYJDy\nggLKCwooDAYpTTCRSzLKSAcNtZ0vWbKEW265haeffpotW7bwzjvvcP755ydZnb/xWzwZME1eSamm\nRsTwSUTXtm3w1lsN1JgCmkQDUHtq9ojI1OwxYzz34pNRRm0qKysZPnw4++yzD127dmX48OFs3Lix\nZvuXX37JmDFj2H///enatSsjRuwK4zpr1iwOPvhgunXrxtlnn82mTZt2K/vFF1/kwAMPZJ999uGG\nG26oWV+XqWbFihUcd9xxHHnkkQB06tSJiy66iHbt2gHw1VdfMWrUKPbZZx/69+/PnXfeGfe3BQIB\nZsyYQW5uLl26dGH8+PGJXyDD8ANp6PVHe/ZMn57UohtFk2gA4k7NLitLKCVkY8uojROjewzr169n\n3bp15OTkcOWVV9ZsHzlyJN9++y3vvPMOn376Kddccw3gLTXjggULWLlyJStXruT555/fzeQUz1Tz\n/e9/n5deeonbbruNZcuWsX379t22jx8/nq1bt7J27VqKiop47LHHaiKFxuLFF1+kuLiYVatW8cwz\nz7Bo0aKEr1G6ada27QRoFpqSFLmzPl3RMXsefxzq+JdKO02iAYhL7ZyQdX1SkBOyS5cunHPOObRu\n3Zp27dpx0003sXTpUgA2bdrESy+9xIwZM+jQoQMtWrSo6bl7Sc1444030rFjR3r16kVhYaGn0MvH\nH388zz77LKFQiDPOOINu3bpx3XXXoaqEw2GefvpppkyZQk5ODn379uW6665j3rx5ccu76aabaN++\nPb179+bEE0/MipSQhgGkzdZ/000wfjxMmgTLljXcrTNVNIkGID8/n6Lc3FrJHGFJXh751dXU5ISs\n45NfXU1RXt6eZeTmNjgl5Lfffsu4cePo168fnTp1YujQoVRWVqKqbNiwgS5dutChQ4c9jvOSmjE6\nS1ffvn13My3VxSmnnMLzzz9PRUUFzz//PHPnzuWRRx7h888/Z+fOnfTp02e3cuOlgwQnAU2EbEkJ\n2exs2w2kyWpKYrz+mPkAanHVVen37EmEJtEABAIBxs2eTWFeHvNzcpifk5Pw1OxklFGb+++/n/ff\nf5+33nqLysrKmt6/qtK7d28qKir46quv9jjOS2rG9evX13xft24dPRvwEJ944omcdNJJrF69mm7d\nutGyZcvdzlteXs7++++fcLmG4UuS2OuPzr9bUFBOMFhIKFS6x3777Zd+z55EaBINAMCA/HymFhfT\nb+lS+i1d2qCp2ckoI5qtW7fStm1bOnToQEVFBbdFZaXv0aMHP/7xj7niiiuorKxk586dvPrqq4C3\n1Iz33XcflZWVrF+/nmnTpu0xRhCLF154gaeffprKykoAli9fzpIlSzj22GMJBAL89Kc/5ZZbbuHr\nr7+mvLyc3//+91x00UUN/v1+pFnYtpNAk9KU5CxdtfPvVlV1oaRkKiNHzmiww0imaDINAPgnJSQ4\nA7HXXHMNVVVVdOvWjeOOO47TTjttt33mzZtHy5YtOfTQQ9l3332ZNm0a4C0141lnnUUwGOSoo45i\n+PDhjBkzpl5NnTt3ZtasWeTm5tKxY0dGjRrFxIkTaxqPBx98kJycHA444AAKCgoYOXIkl1xySdzf\nF4/XXnstpmnLMNJOCmz9sfPvBvjww6ENdhjJFBYLKAUEg0EmT57MmWeemWkpvsXv99DIclIYw6e4\nuJghQ8r57rsRu63PyZnP0qX9Up4PwAteYwE1qTcAP1BaWsq7777b4IFjwzAaSYo9fLp0yQeKoJbL\nSG7ukqz7v7cGIInceOONnHrqqdx777272esN/9CkbNspJCs1JdnWH4/+/QMsW7Yr/27r1rczcODV\nzJ49Lqn5ANKBRQNNIlOmTGHKlCmZlmEYzY80Ru4EOOqoXfl3V6z4nLFjJ2Vd5Q82BmBkCLuHRlJI\nka1/2zZ46CFYs8YpOtuwMQDDMJoEcSP0psDWH4nZc+ih8PrrcPXVjRTvc6wBMJoVWWnbzgB+0RQd\noffFIUOcCL3Ll6fE1v/yyxAMOjF7/vIXeO45pyGoD79cq4aQlWMAffv2tdj0WU50qAvDiEV0hN4A\n0AUoKCmhcOhQpp50EoEk2/pXr4Zbb01fNi4/kJVjAIZhNH2Ki4spLyhgRK0gjfNbtaLf668THDQo\nQ8r8j40BGIbRNGnZslFd9K1bk6gly7EGIEH8aO/zoybwpy7T5A0/aMofMICi9u1rplsV0bgIvdu2\nOY5CBxwA776bPJ1+uFYNJS0NgIgERCQkIi+4y5NFZIOIrHQ/yc9ibhhG9lJSQuCYYxh30EEUHn44\n83NyWNK6dYMi9EZn43r9dXj1VW+Du82BtIwBiMg1QBDooKpnishkYKuqPlDPcTYGYBjNiRh+/WHV\nmiBr+fn5CVX+H3wAP/kJtG0L993nv4QsqcLrGEDKvYBEpBdwGnAncG30plSf2zCMLCLObN6ASIMD\nrPXs6XiMnn128/HsSYR0mIB+D1wP1O7KjxeREhF5REQ6pkFHUvCjvc+PmsCfukyTN9KqyWMMn4Zo\nysmBc85JbeXvx/vnlZS+AYjI6cBmVS0RkWFRm6YDd6iqishvgQeAS2OVMXr0aPr16wdAp06dyMvL\nq0kNF7nw6VwuKSnJ6PljLUfwi57IciRHsF/02P3z4fIHHzDsoYegVy+Kpk+Hbt0Y5tbWiTxPX34J\n8+cXcdBBzfP+FRUVMXfuXICa+tILKR0DEJG7gJHATqAt0B54VlVHRe3TF1ioqkfGON7GAAyjKZJA\nDJ9wOBx3DCASs+eee5ywDbfemhb1vscX8wBU9WZV7aOqBwAXAK+o6igR6RG12whgdSp1GIbhIxKI\n4RMv924szx6r/BMnU/MA7hWRt0WkBBgKXJMhHQlT+7XPD/hRE/hTl2nyRko0JRivP17u3TFjZnDe\neWGmT4fHH/cesydV+PH+eSVtsYBUdQmwxP0+qp7dDcNoSjQgXn+83LtlZUO5++4Qp5wSNM+eRmKx\ngAzDSB2NiNdfXFxMQUE5VVX+zb3rV3wxBmAYRjOmkfH6+/bNp337IppC7l2/Yg1AgvjR3udHTeBP\nXabJG43S1MjcvJGYPYcdFuD448fxve/5O/euH++fV7IyH4BhGD6lEbl5w2FnUPfWW+GooyIxewYQ\nDmd/7l2/YmMAhmE0niTk5l26FCZObF4xe1KF1zEAawAMw2gc0b3+mTMblaVL1WL2JAMbBE4RfrT3\n+VET+FOXafKGJ02NtPXHoq7K34/XCfyryws2BmAYRuI00NZfUQF33QXt2ztth5FZzARkGIZ3Gmjr\n37YNHnwQ7r0Xzj3Xqfz32y8NepspvskHYBhGE6GBvf5582p79qRWpuEdGwNIED/a+/yoCfypyzTV\nTTgcpri4mBkzZhAOuxOwGmnrX7MmOTF7/HSdovGrLi/YG4BhGACUhkLMGDOGYWVlbK6upvDhhxk3\ncSIDpkxpkF9/hNtuS75WIznYGIBhGITDYQqDQaaWlNSYBcJAYYsWTH3kEQIXX1yvrX/LFuiYNbn9\nmjbmBmoYhmdCoRDDyspqxd2Eoa1aETriiDor/4oKmDDBic2/ZUvKpRpJxBqABPGjvc+PmsCfukyT\nN4oiX+qo+CMxew45BL7+GkKh1L4B+PE6gX91ecEaAMMwyAeKwuFacTdhSW5uzMiby5btno3r4YfN\nrTMbsTEAw2jORPn1l44fz4xnn2Xo++8DUHTwwVw+Zw4DYjQAGzfChx9azB6/YrGADMOom1DI8evv\n3bsmhs/OnTt56qmnALjgggto2dIcBbORpA0Ci8h5ItLe/X6riDwrIkclQ2Q24kd7nx81gT91mSZ2\n+fWfcgp0QZ6OAAAgAElEQVRcd12NX38oVMrgwdcyblwOv/jFhwwefC0vvlhKWVl65cXDj/cO/KvL\nC17GACap6lYROR74IfBn4E+plWUYRkoIhWDw4D2ydNVOwP7ddwWUlEzlrLNm8MYb4frLNbKSek1A\nIhJS1XwRuRv4n6o+EVmXcnFmAjKM5FBPDJ94+Xfbtp3Pq69a/t1sI5mxgD4WkRnAj4B7RKQ15j1k\nGNlDtK0/wdm8Fpu/aeOlIv8p8BJwiqpWAl2A61Oqysf40d7nR03gT13NSlMcW38s8vPzyc0tYlcC\ndue7nxKw+/HegX91eaHeBkBVq4BPgYjD107g/VSKMgyjkcSx9ccjEAgwe/Y48vIiCdiX+DIBu5Fc\nvIwBTAYGAYeoaq6I9AT+qqpDUi7OxgAMIzE8xOsvL4c//9kJ0la7bg+Hw4RCIcB5K7DKPztJZiyg\nc4AzgW8AVHUj0L5x8gzDSDr19PojMXuOcp24d+7cs4hAIEAwGCQYDFrl3wzwcoe3u91wBRCRdqmV\n5G/8aO/zoybwp64mqSna1j9hwh62/m3b4L77dsXsWb0a7rgDWrVKoaYU4EdN4F9dXvDiBfSM6wXU\nSUTGAmOAWamVZRiGJyIePn36xPXwefJJJ3aPZeMyauMpFISI/Ag4GRDgJVX9d6qFuee1MQDDiEW0\nrf/++2HkyLiDvKrmztncsFhAhtFUie71z5jRoCxdRtMmmbGAtorIV+5nm4hUi8hXyZGZffjR3udH\nTeBPXVmtqbat/4UXair/8nLH4Wf27DRrSiN+1AT+1eUFL/MA2qtqB1XtALQFzgWmp1yZYRi7iHj4\nrFzp2Ppd985oz54DDoDzzsu0UCObaJAJyGIBGUaaiGPr37kTpk6Fe+6Bc891XgwsIYsRIWmxgEQk\nOjpUAGdS2LZGaDMMwwt1ePgEArB5MyxdCocdljmJRnbjZR7A8KjPKcBW4KxUivIzfrT3+VET+FNX\nVmiqw9YfIRBwfPtTVflnxXXyCX7V5YV63wBU9ZJ0CDGM5kYk7MJ7771HQUGBM/M2Rq//yy+hc+dM\nqzWaInHHAETkQdzZv7FQ1atSJSpKg40BGE2S0lCIGWPGMMxNt1V00EGM+/73GbBgQY2tv3ydcOut\nsGKFM3u3RYsMizayhkbPAxCRi+s6UFUfbaA2z1gDYDRFwuEwhcEgU0tKamywYaCwQwemlpZSmdOL\nu+6COXNg/HjHCtTeom8ZCdDoeQCq+mhdn+TKzR78aO/zoybwpy4/aAqFQpzw7rs1/3xFOP+Ix2/f\nzn2Pb94tZs/tt2em8vfDdaqNHzWBf3V5wYsXUHdgInA40CayXlVPSqEuw2iyhMNhdny3Y4/127/b\nSf/+YfPsMdKGl3wAi4CngQnA5cDFwGeqOjHl4swEZDQ1tm9nw5VXctcjs3mI8G4moJOkH/f99xkG\nDx6cQYFGUyCZ+QC6quqfgR2qukRVxwDW+zeMRHFn8+793nv8s+XNDCWPx8jhMXIYykBWtrrCYvAb\nacXL0xZ5V90kIqeLSD5OXuBmiR/tfX7UBP7UlRFNUX79n4+ewPjei9nAFl7jLS5mKRfzAK+xggMP\nW2/5d+vAj5rAv7q84KUB+K2IdASuwzEDPQJck1JVhtFUcHv92/+7kjtGlHDIby/iwINasHjxOPLy\nriUnZy2tW3/CwIHXWP5dI+14GQPorqqfpUlP7XPbGICRnUTF8Pn2t/dz4OSRnHmW7Bazx/LvGqki\nafkARKQMWIszEPysqn6ZFIUesAbAyEpixOv/7DPo3j3TwozmQtIGgVU1F7gVGAAUi8jfRWRkEjRm\nJX609/lRE/hTV0o11RHDp67Kv9ldpwbiR03gX11e8PTOqarLVfVa4GigAkhoIpiIBERkpYi84C53\nFpFFIvKeiLzkjjEYRvYSCvHt9wbz/jO7x+s3DD/jxQTUATgHuAA4EHgOeEZViz2fROQaIAh0UNUz\nReQe4AtVvVdEJgKdVfXGGMeZCcjwN9u3U3n9nQRm/Ilb2tzP4XeN5PJfitX9RkZJ5hjAR8ACnEr/\njQYI6QXMAe4ErnUbgHeBoaq6WUR6AEWqemiMY60BMHzLlqIQ3/x0NG9X9mH1+BmMu72nxewxfEEy\nJ4IdoKrXNKTyd/k9cD27RxbdV1U3A6jqJ8A+DSw77fjR3udHTeBPXUnR5Nr6W55+Ci99bwID177A\nhAcaXvk32euUZPyoCfyrywte8gE0uAsuIqcDm1W1RESG1XWaeBtGjx5Nv379AOjUqRN5eXkMG+YU\nFbnw6VwuKSnJ6PljLUfwi57IcklJia/0JOX+vf8+wx56CPr0Yfnc6fTv3o39ekqj9EXww/Xx87If\nn6doMqmnqKiIuXPnAtTUl15oUE5gz4WL3AWMBHbiJJRvjzOGMAgYFmUCWqyqe4S/MhOQ4Rvi5OY1\nDD+SNBOQiAzxsi4WqnqzqvZR1QNwBpFfUdWLgIXAaHe3i4HnvZRnGOkmFIJfHhPi68MGw0rz8DGa\nFl7GAB70uC4RpgA/EpH3gB+4y1lB7dc+P+BHTeBPXV41lZfD6Au388oJk7m/9BTa3Bo7N286NaUT\n0+Qdv+ryQtwxABE5FjgO6C4i10Zt6gAknJxOVZcAS9zvFcAPEy3DMFLN1q1OEpYVs0I81XY03U7o\nQ8s/l6Sk4jeMTFNXSsihwDCcHAAPR23aCixU1fdTLs7GAIw0U1W5nddPu5OTyv5Ei9+brd/ITpI5\nD6CvqpYnTVkCWANgJAtPgddixPAxjGwkmfMA5orIK7U/SdCYlfjR3udHTeAfXaFQKcFgIQUF5QwZ\n8iLBYCGhUClffOHuUEcMn3Tgl+sUjWnyjl91eaHeeQA4OQAitAHOxXHrNAzfEw6HGTNmBiUlU3H6\nO10oKSmgoKCQE06Yyj/uXLWr119itn6jedGgeQAislxVj06BntrnMROQ0SiKi4spKCinqupsIOSu\nzafdXs/w7qhX6PXCAvPrN5ocXk1A9b4BiEh0+scATlA3i95pZA0tqj/keIJcRhkA8+jFFTu+Yu+y\ng63XbzRrvIwBFAMr3L9v4KSGvDSVovyMH+19ftQE/tA1cOBAjgncwxJKuIgqelPFvyjj8b2+Zu//\n/McXlb8frlNtTJN3/KrLC15iAfVPhxDDSAWrVq1irH69W08nAJzXYier3n6bYDCYKWmGkXG8uIG2\nAa4AjscJ2vYq8LCqbku5OBsDMBLk5ZfhuOMgJ8dZLn7zTcpPOIERO3f3W5ifk0O/pUutATCaJMl0\nA30MJx3kg8BD7vd5jZNnGMklFIKTT4YrroC1a3etzL/sMopycghH7RsGluTmkp+fn36hhuEjvDQA\n31PVS1V1sfsZi9MINEv8aO/zoyZIj67ycic222mnwdlnQ2kpHH7QLr/+wIQJjFu8mMK8PObn5HB7\n69ZcPXAg42bPjj0ZLAP48f6ZJu/4VZcXvMwDWCkix6jqmwAi8n2cQWHDyChr1sDRR8P48VBWhpOQ\nJTKbt3fvGg+fAcDU4mJCoRCfr1jBpLFjfVP5G0Ym8TIG8A5wCLDOXdUHeA9nMpiq6pEpE2djAEY9\nfPEFdO3K7vH6f/c7C9lsNGuSNg8AODUJegwjJXTtSsxev2EY9ePlPfi3qloe/Ylel2qBfsOP9j4/\naoLk6Vq0CKZPj7EhOobPddfBwoX1Vv5+vFamyRt+1AT+1eUFLw3AbgO+ItISZzawYaSUiGfP+PEx\n6vVQCAYPhuJip9c/apSZfAwjQerKB3ATcDNOLt8qIPLftR2Yqao3pVycjQE0S8rL4dZbHZ/+SZNg\n7FjYay93o9n6DaNekpkP4O50VPZxzm0NQDPksstgv/2cyMzt20dtiLb1z5xptn7DiEMyJ4L9U0QK\nan+SoDEr8aO9z4+aoOG6Zs500jLWVP4NsPUnW1MqMU3e8KMm8K8uL3jxAro+6nsb4GicwHAnpUSR\n0WxQ9WC9MQ8fw0gZCecDEJHewFRVPTc1knY7l5mAmiiLFsHNN8PTT8OBB8bYwWz9htFgkjkPoDYb\ngMMacJxhEArBxIlOvJ6774YDDoizk/X6DSPl1DsGICIPisgf3M9DONFAV6Zemj/xo73Pb5rC4TDF\nxcXMmDGDcNgJw7Zpk+OpedppcM45Tsyec8+t1alPoq0/Hn67VmCavOJHTeBfXV7w8gYQHfdnJ/Ck\nqr6eIj1GlhMKlTJmzAzKyoZRXb2Zhx8uZPbscey77wAOPBD++Mdanj27DrRev2GkGa/5AA5yFz9I\nRx6AqHPbGEAWEQ6HCQYLoxKwA4TJyyukuHhq7ABsZus3jKTTaDdQEWkpIvfi2PwfxckLsF5E7hWR\nveIdZzRfQqEQZWXDoFb+rbKyoYRCoVgH2Gxew8ggdY0B3Ad0AfqralBVjwIOBDoBv0uHOD/iR3uf\nXzS9+SZ8+230mqLYO6bB1h8Pv1yraEyTN/yoCfyrywt1NQBnAGNVdWtkhap+BfwSOC3VwozsIRKz\nZ+rUfPr2LYJa+bdyc5fsyr5lvX7D8A11xQIqU9XcRLclExsD8D9XXQV//Sv8+tfwi1/A6tWlXHLJ\nw7z3Xm8AcnPXMXfuL8kfcLDZ+g0jTSRjHsD/icgoVX2sVsEjgXcbK9BoGpx/vlOvRzx7WrGdE3iV\nW3gPgKVyCK3eHQKjLzQPH8PwGXWZgK4ErhSRIhG53/0sAa7CMQM1S/xo78ukpiFDdlX+4XCYGWPG\nMG3VKs7bto3u27YxbdUqZlx0EeFrrkmrrT8edv+8YZq841ddXojbAKjqx6r6feAOYK37uUNVj1bV\nj9Mjz/AD4TAsWOD8rYtQKMSwsrJaPkAwtFUrQkccYSYfw/AZCccCSic2BpB5Fi2CG26Atm3hhReg\ne/f4+xYXF1NeUMCIqqrd1s/PyaHf0qUEg5ZHyDDSQTLDQRvNkOhsXJMmwbJldVf+APkiFIXDtXyA\nYElu7i4vIMMwfIM1AAniR3tfsjW98ooTs+fss+PE7KmN69cfOPVUxt10E4UDBzI/J4fbW7fm6oED\nGTd7duxZwBmgOdy/ZGCavONXXV5oSDRQo4lTUABlZXFi9tSmpMSJ4dOrF5SUMKBnT6beeiuhUIjP\nV6xg0tixvqn8DcPYHRsDMBqGxfAxDN+SynwARhMgHIbHH3eSrV9wQYIH1+r1Z9q10zCMhmHv5gni\nR3tfopoWLYJgEKZPh759EzgwEsPn5JPh2mvr9etvCtcqHZgmb/hRE/hXlxfsDaAZUTsb14gRCVht\nrNdvGE0OGwNoJqg6Xj2nnAJjxzqmH0+Yrd8wsg6vYwDWABjxie71z5xpvX7DyBJsIliK8JO9L1bu\nXXB6+40iQVt/PPx0rSKYJm+YJu/4VZcXrAHIUkKhUoLBQgoKyrn66s0Eg4UUF5cybx7k58PWrfWX\nEZOSEjj6aIvXbxjNADMBZSHxcu+2aVPIwIFT+d3vAhx/fIKFmq3fMJoMNg+gCRMv9244PJSHHgox\naFCCQdfMw8cwmiVmAkoQf9r7igBo2TLBTnuSbP1xVfnwWpkmb5gm7/hVlxdS2gCISGsR+a+IhETk\nfyIy2V0/WUQ2iMhK93NqKnU0NfLz88nNLaLO3Lv1YbZ+w2j2pHwMQERyVLVKRFoAr+NkFPsxsFVV\nH6jnWBsDcNm2Db79Fjp3dpZDoVLGjJlBWdlQAA4+uIg5cy4nP39A3QWZrd8wmjy+cQNV1Uh2kNY4\nYw6RGt1qHQ+EwzBvHhxyiBO7J0J+/gCKi6eydGk/li7tx8qV0+qv/K3XbxhGFClvAEQkICIh4BPg\n36r6lrtpvIiUiMgjItIx1TqSRTrtfYsWwVFHOTF7Hn/cSc4STSAQIBgMsnXr1rpDLqfY1h8PP9pG\nTZM3TJN3/KrLCyn3AlLVMJAvIh2A50TkcGA6Tn5hFZHfAg8Al8Y6fvTo0fTr1w+ATp06kZeXx7Bh\nw4BdFz6dyyUlJSk/3wknDGP4cHj77SIuuwwmTRqGSPz9I8Tc/sEHDHvoIejVi6Lp06FbN4a5vf5U\nX6+SkpKUlu/X+5focgS/6PHrsh+fp2gyqaeoqIi5c+cC1NSXXkjrPAARmQR8E237F5G+wEJVPTLG\n/s12DGDJEjjuuARi9tTGbP2G0WzxxTwAEekG7FDVLSLSFvgRMEVEeqjqJ+5uI4DVqdSRjQwdWv8+\n4XCYUCgEOJ5BNWYg8+s3DMMDqR4D2A9YLCIlwH+Bl1T1H8C9IvK2u34ocE2KdSSN2q99jWHbNliw\noGHHloZCFAaDlBcU8OKQIRQGg5QuX54RW388knmtkoVp8oZp8o5fdXkhpW8Aqvo/4KgY60el8rx+\nJ5KN69ZbnUHe006DVq0SOT7MjDFjmFpSQgDoAhSUlFA4dChTTzqJgPX6DcPwgMUCSjOLFsENN0Db\ntnDffSQeswcoLi6mvKCAEVVVu62f36oV/V5/neCgQUlSaxhGNuKLMQBjd2bNcir9hLNxeSXhWBCG\nYTRnLBZQgjTG3vfzn0NpKZx7buPq6fwBAyhq374mEEQRTlCIJbm53kNBpAE/2kZNkzdMk3f8qssL\n9gaQRnJyklBISQmB0aMZd9BBFHbtytC1a1ldXc38Qw/l8tmz654QZhiGEYWNASSZbdvgwQedwd0f\n/CCJBcfw6w+rxnYDNQyjWWNjAGkm2rMnPx/OOiuJhcfx6w+IEAwmGPvfMAzDxbqMHomXfxf2jNmz\nYAHk5ibhpB5j+PjVBulHXabJG6bJO37V5QVrADwQK/9uKFQKOCafu++GSZNg2bKGuXXGxCJ3GoaR\nYmwMoB7i5d/NyyukuHhq8u3uFsPHMIxG4pt8ANlOvPy7ZWVDawZgk4b1+g3DSCPWACRMUfKLbGS8\nfr/aIP2oyzR5wzR5x6+6vGANQB2Ew7B6dT47dhTRqPy7dWG9fsMwMoSNAcTh1VfhqqugTRsYO7aU\nBx9sQP7dujBbv2EYKcLrGIA1AHFYtAi2bt0Vsydu7P2GEO3XP3OmRe40DCOp2CBwIzn55N1j9njO\nv1sXKcrN61cbpB91mSZvmCbv+FWXF5r9TOCKCsfMk5Q4PXVhWboMw/AZzdYEFInZc++98OijTlKW\nlGC2fsMw0ozFAopD7Wxcr74Khx6aopNZr98wDB/TrMYAtmyBYHBXzJ7nnku88vdk70uRrb9RmjKA\nH3WZJm+YJu/4VZcXmtUbQMeOMHUqFBSk0ApjvX7DMLKEZjsGkHTM1m8Yhk9o1m6gFRXw97+n8YQ2\nm9cwjCykSTUA27Y5SdcPOQT+/e/klh0zH0Cabf3x8KsN0o+6TJM3TJN3/KrLC02iAQiHYd48p+Jf\ntszx7Jk2LXnll4ZCFAaDlBcUsPnqqykMBil96inr9RuGkdU0iTGAiRNh6VKn95+0hCwu4XCYwmCQ\nqSUlUdkAoLBFC6Y+8giBiy+2it8wDF/RrOYBTJoE7dqlph4OhUIMKyurlQ0AhrZqReiIIwha5W8Y\nRpbSJExAe++dvk54UeSLjyp+v9og/ajLNHnDNHnHr7q84PsGIDLgWlEBEybAO++k9/z5QFE4XCsb\nACzJzU1OPgDDMIwM4fsxgCOP/BUnnTSOv/xlAOeeC7ffDvvum4aTR/n1l44fz4xnn2Xo++8DUHTw\nwVw+Zw4DrAEwDMOHNJl8AFBNx46FLFs2lcMPT9MLSyjkzObt3bsmXn9S8wEYhmGkkCY0ESzAjh1D\n+fbbJCdgj0XEr/+UU+C663bz609KPoAU4VcbpB91mSZvmCbv+FWXF5qEF1BSiO71WwwfwzCaAVlh\nAsrLK6S4eGpqet4Ww8cwjCZGk5kHMHDg1cyefXlqKn/r9RuG0YzxlzE7BitXTiM/f0ByC63D1l8f\nfrT3+VET+FOXafKGafKOX3V5wfdvAEnv+Vuv3zAMA8iCMYCk6TNbv2EYzYQmMwaQFKzXbxiGsQe+\nHwNoFI2w9cfDj/Y+P2oCf+oyTd4wTd7xqy4vNN03AOv1G4Zh1EnTGwMwW79hGM2c5jkGYL1+wzAM\nzzSNMYAU2Prj4Ud7nx81gT91mSZvmCbv+FWXF7L/DcB6/YZhGA0ie8cAzNZvGIYRk6Y9BmC9fsMw\njEaT0jEAEWktIv8VkZCI/E9EJrvrO4vIIhF5T0ReEpGOngpMo60/Hn609/lRE/hTl2nyhmnyjl91\neSGlDYCqfgecqKr5QB7wYxE5GrgReFlVDwFeAW6KV0YkJzChEAweDMXFTq9/1KiMmHxKSkrSfs76\n8KMm8Kcu0+QN0+Qdv+ryQsq9gFS1yv3aGsfkpMBZwKPu+keBs+MdX5ifT+lll2W01x9NZWVlxs4d\nDz9qAn/qMk3eME3e8asuL6R8DEBEAkAxcCDwR1V9S0T2VdXNAKr6iYjsE+/4qW+/TeHatUwtLSXQ\nq1eq5RqGYTQb0vEGEHZNQL2Ao0VkAM5bwG67xTs+AAzduZPQ5s0pVOmdtWvXZlrCHvhRE/hTl2ny\nhmnyjl91eSGtbqAiMgmoAn4BDFPVzSLSA1isqofF2N+/PqqGYRg+xosbaEobABHpBuxQ1S0i0hZ4\nCZgCDAUqVPUeEZkIdFbVG1MmxDAMw9iDVDcAR+AM8gbcz9OqeqeIdAGeAXoD5cBPVTV7R1IMwzCy\nEF/PBDYMwzBShy+DwYnIqSLyroiUuSaijCMifxaRzSLydqa1RBCRXiLyioiUuhPtrvKBppiT//yA\niAREZKWIvJBpLRFEZK2IrHKv1/JM6wEQkY4i8lcRecd9tr6fYT257vVZ6f7d4pNn/RoRWS0ib4vI\n4yLSygearnb/7zzVB757A3DdRsuAHwAbgbeAC1T13QzrOh74GnhMVY/MpJYI7gB6D1UtEZG9cdxt\nz/LBtcpR1SoRaQG8Dlylqhmv3ETkGiAIdFDVMzOtB0BEPgSCqvplprVEEJG5wBJVnSMiLYEcVf0q\nw7KAmvphA/B9VV2fQR09gdeAQ1V1u4g8Dbyoqo9lUNMA4ElgMLAT+Cdwuap+GO8YP74BHA28r6rl\nqroDeApn4lhGUdXXAN/8k4Izh0JVS9zvXwPvAPtnVlXcyX8ZRUR6AacBj2RaSy0EH/0fikgH4ARV\nnQOgqjv9Uvm7/BBYk8nKP4oWQLtII4nTYc0khwH/VdXvVLUaWAqMqOsA3zx4UewPRN/cDfigUvM7\nItIPJ9zGfzOrpMbUEgI+Af6tqm9lWhPwe+B6fNAY1UKBf4vIWyIyNtNigP7A5yIyxzW5zHQ9+PzC\n+Ti93IyiqhuB+4F1wMdApaq+nFlVrAZOcGOt5eB0eHrXdYAfGwAjQVzzz9+Aq903gYxSa/Lf90Xk\n8EzqEZHTgc3u25K4H78wRFWPwvlnvdI1NWaSlsBROLP2j8KZt+MLF20R2Qs4E/irD7R0wrFM9AV6\nAnuLyIWZ1OSafu8B/g38AwgB1XUd48cG4GOgT9RyL3edEQP39fNvwDxVfT7TeqJxTQeLgVMzLGUI\ncKZrb38SOFFEMmarjUZVN7l/PwOewzGBZpINwHpVXeEu/w2nQfADPwaK3WuVaX4IfKiqFa655Vng\nuAxrQlXnqOogVR0GVOKMp8bFjw3AW8BBItLXHVW/APCL14bfeo8As4H/U9VpmRYCzuS/SHhv13Tw\nIyCjg9KqerOq9lHVA3Cep1dUdVQmNYEzWO6+vSEi7YCTcV7jM4Ybo2u9iOS6q34A/F8GJUXzM3xg\n/nFZBxwjIm1ERHCu0zsZ1oSIdHf/9gHOAZ6oa3/fJYRR1WoRGQ8swmmg/qyqfriwTwDDgK4isg6Y\nHBkoy6CmIcDPgf+5NncFblbVf2VQ1n7Ao663RmTy3z8yqMfP7As854Y8aQk8rqqLMqwJ4Crgcdfk\n8iFwSYb14Nq0fwhclmktAKq6XET+hmNm2eH+nZlZVQDMdyfa7gCuqG8A33duoIZhGEZ68KMJyDAM\nw0gD1gAYhmE0U6wBMAzDaKZYA2AYhtFMsQbAMAyjmWINgGEYRjPFGgDDd4hItRuHZrUb/vfaqG1B\nEZmawnMPFZFjo5bHicjIBpQzUER+HLU8XERuSJLGjzzss9idDGQYcfHdRDDDAL5x49BE0oo+KSId\nVPU2VS3GCXvdYESkhTt9PxbDcMJ+vwGgqjMaeJo8YBBOSF5UdSGwsIFl1cYm7xhJwd4ADF+jqp/j\nzP4cDzU99IVR3yOJQordcAqIyEQ3SUdIRO5y1y0Wkd+LyFvAVW7Iir+Jk7zmvyJyrIj0BS4HCt0y\nh4jI5MgbiIhc5SZIKXFnhiMig0VkmXv+10TkYHcG7R3AT91yzhORi0XkQfeYviLyH7ecf7uhqnEj\ncE4TkddF5AMRiRfK97Oo379YdiVvmRe1zxfUEwjMMOwNwPA9qvqRG2K6e2SV+/c6nOnub7ihAraJ\nyKnAcGCwqn7nRm2MsJeqDgYQkceBB1R1mYj0Bl5S1cNF5GFgq6o+4O73w6jjJwL9VHWHOHHzwYn/\ncryqhkXkB8DdqvoTEfk1TqKXq9xyLo7S/SAwR1X/IiKXuMvnuNt6qOoQETkMJwbWszGuR3SGrjzg\ncJzQ26+LyHGqukxVf+Lt6hrNGWsAjGwhVhC+14Hfu5X5s6r6sVthz1HV7wBUtTJq/6ejvv8QOMwN\n5AVOON+cejSsAp4QkQXAAnddJ+AxETkYp4L38j91LLsq/Hk4IXwjLHB1vyMi+3goa3kkoqiIlAD9\ngGUejjMMMwEZ/kdEDgB21g4DrKr3AJcCbYHXROSQeor6JrpYnLSC+e6nT1Qms3icDjyEEx75LTfg\n3W9woosegfPm0cbDT6rLhv9dLY31Eb1/NdapMxLAGgDDj9RUfK7Z5084ZpLddxI5QFVLVfVeYAVw\nCE4yjEvcUNSISOc451gEXB1V1kD361agQ8wjoI+qLsFJkNIB2BvoyK58FdFRM+sqZxlOaGOAkcCr\ncfbzW+hxo4lhDYDhR9pE3EBxKup/qeodMfYrFJH/uaaP7cA/VfUlHNv5ChFZiTNOAHv2uq8GBonI\nKmDMsW4AAACJSURBVPc849z1C4FzIoPAkePESbzzFxFZheOFNM0NtXsvMEVEitn9/2kxcHhkELjW\nua/CaaRKcMJ5Rxqi2hoT9fYx7yAjISwctGEYRjPF3gAMwzCaKdYAGIZhNFOsATAMw2imWANgGIbR\nTLEGwDAMo5liDYBhGEYzxRoAwzCMZoo1AIZhGM2U/wcinqIXnAJAGAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa728e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# YOUR CODE HERE\n",
    "n = 10\n",
    "dx = n\n",
    "b0 = 30.\n",
    "b1 = 50.\n",
    "max_iter = 10000\n",
    "epsilon_tol = 1e-6\n",
    "\n",
    "#Direct Solution\n",
    "direct = plt.plot(inv_tridiag_121_b(n,b0,b1),'bo--', label='Direct Soln.')\n",
    "#plt.xlim([-1,10])\n",
    "\n",
    "#Jacobi Method\n",
    "u_jac,iter_jac = solve_laplace_jac(b0,b1,dx,max_iter,epsilon_tol)\n",
    "jacobi = plt.plot(u_jac,'ro-', label='Jacobi Soln.')\n",
    "\n",
    "\n",
    "plt.title('Comparing solutions from different solvers')\n",
    "plt.legend(loc = 'upper left')\n",
    "plt.xlabel(\"Discretisation 'n'\")\n",
    "plt.ylabel(\"Output values\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first obvious observation would be the difference in gradients of the two methods.  In Jacobi iteration, the start and end values fixed and only the nodes inbetween were updated by the loops whereas in the direct solver, all the nodes, including the boundary nodes, were updated as the iteration progressed.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Gauss-Seidel method\n",
    "##### (10 Points)\n",
    "\n",
    "You probably noticed one aspect when you wrote the function for the Jacobi method above: when you iterate along the unknown values $u_i$ you actually have already updated information about $u_{i+-1}$ because you calculated it in the previous step! So: why not use this newer information instead of the one from the last time step?\n",
    "\n",
    "This is exactly the idea of the Gauss-Seidel method, a simple extension to the Jacobi method. Instead of using the previous value for $u_{i+1}$, we simply use the already updated value:\n",
    "\n",
    "$$u_i^{n+1} = \\frac{1}{2} \\left( u_{i-1}^{\\color{red}{n+1}} + u_{i+1}^n - \\Delta x^2 q_i \\right)$$\n",
    "\n",
    "This is a simple adaptation, but leads to a faster convergence of the problem!\n",
    "\n",
    "**Exercise**: \n",
    "\n",
    "- Implement a function to perform one Gauss-Seidel iteration **(2 Points)**  \n",
    "\n",
    "```python\n",
    "def gauss_seidel_iter(u,q,dx):\n",
    "```\n",
    "\n",
    "- Adapt your *solve_laplace_jac* function from above to use the Gauss-Seidel method as function *solve_laplace_gs* **(3 Points)**  \n",
    "    These steps will take some work, so here the base construct that you could use:\n",
    "```python\n",
    "def solve_laplace_gs(b0,b1,dx,max_iter,epsilon_tol):\n",
    "    # initialize u with all ones\n",
    "    u = \n",
    "    # Set boundary conditions b0 and b1 at beginning and end of u\n",
    "    # set parameters\n",
    "    # no sources / sinks, so set q to 0\n",
    "    # do a for loop for jacobi iterations\n",
    "        # copy results from previous iteration\n",
    "        # perform iteration\n",
    "        # check for convergence by calculating epsilon\n",
    "        # if epsilon is smaller than epsilon_tol, break for loop\n",
    "    return u,i # i being the current iteration\n",
    "```\n",
    "- Apply the function **solve_laplace_gs** iteratively until a convergence criteria is met\n",
    "- The function should return u and the current iteration i_gs:\n",
    "    - *u,i_gs = solve_laplace_gs(b0,b1,dx,max_iter,epsilon_tol)*\n",
    "- Plot the results **(3 Points)**\n",
    "- Check if my claim that it converges faster is actually true! **(2 Points)**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gauss_seidel_iter(u, q, dx):\n",
    "    \"\"\"Function to perform the Gauss-Seidel iteration for all values in the array u\n",
    "    u = 1-D np.array with values of previous iteration\n",
    "    q = 1-D np.array with source/sink terms (not required for Laplace solution!)\n",
    "    dx = float : node spacing (not required for Laplace solution!)\n",
    "    \"\"\"  \n",
    "    # Note: copy of u not required (as we operate always on updated array!)\n",
    "    \n",
    "    # YOUR CODE HERE\n",
    "    \n",
    "    for i in range(1,dx-1):\n",
    "        u[i] = 0.5 * (u[i-1] + u[i+1] - (dx**2 * q[i]))\n",
    "    return u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 40.        ,  42.22222075,  44.44444185,  46.66666338,\n",
       "         48.88888537,  51.11110781,  53.3333306 ,  55.55555365,\n",
       "         57.77777683,  60.        ]), 135)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def solve_laplace_gs(b0,b1,dx,max_iter,epsilon_tol):\n",
    "    \"\"\"\n",
    "    Function to solve the 1D problem in iterations until \n",
    "    convergence is met\n",
    "    Initialize u within this function!\n",
    "    b0 = boundary condition 1\n",
    "    b1 = boundaray condition 2\n",
    "    dx = discretization (= n)\n",
    "    max_iter = maximum iteration number\n",
    "    epsilon_tol = convergence limit\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    # initialize u with all ones\n",
    "    u = np.ones(dx)\n",
    "    \n",
    "    # Set boundary conditions b0 and b1 at beginning and end of u\n",
    "    u[0] = b0\n",
    "    u[-1] = b1\n",
    "    \n",
    "    # set parameters\n",
    "    \n",
    "   \n",
    "    # no sources / sinks, so set q to 0\n",
    "    q = np.zeros(dx)\n",
    "    \n",
    "        \n",
    "    # do a for loop for jacobi iterations\n",
    "    for iter in range(1,max_iter):\n",
    "\n",
    "        # copy results from previous iteration\n",
    "        u_previous = u.copy()\n",
    "                    \n",
    "        # perform iteration\n",
    "        u_current = gauss_seidel_iter(u, q, dx)\n",
    "             \n",
    "        # check for convergence by calculating epsilon\n",
    "        epsilon = np.sqrt(np.sum((u - u_previous)**2))\n",
    "                      \n",
    "        # if epsilon is smaller than epsilon_tol, break for loop\n",
    "        if epsilon < epsilon_tol:\n",
    "            break\n",
    "        \n",
    "    return u, iter\n",
    "\n",
    "#########################################################\n",
    "\n",
    "solve_laplace_gs(40.,60.,10,10000,1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "assert_array_almost_equal(solve_laplace_gs(40.,60.,10,10000,1e-6)[0],[ 40.,42.22222075,44.44444185,46.66666338,\n",
    "                                                                      48.88888537,51.11110781,53.3333306 ,\n",
    "                                                                      55.55555365,57.77777683,  60.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot $u$ from the Gauss-Seidel Iteration as you did for the Jacobian.\n",
    "*Note* Your function should return two parameters, $u\\_gs$ and the current iteration $iter\\_gs$, so be called:\n",
    "```python\n",
    "u_gs,iter_gs = solve_laplace_gs(b0,b1,dx,max_iter,epsilon_tol)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEZCAYAAACervI0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd4FNX6xz8nQOihIy2EIqBwVUIARarlig28wFUQATGK\n2EVQwYJgQaxXbDQFwUJRsWDnd68EEFBM2EWJAjYCAqIIgcQAgez7+2Nm4yZkk0mym0yS9/M882Rn\n5sw5352dnDPznnfe14gIiqIoSsUjorQFKIqiKKWDDgCKoigVFB0AFEVRKig6ACiKolRQdABQFEWp\noOgAoCiKUkHRAaACYowZboz5tLR1+DHG9DXG7CzG8fcYY+aGUpODNs82xmwzxhwyxgwsybaD6HnF\nGPOQ/bmXMeb7gH3tjTEeY8xBY8wtxphqxpgPjDGpxpilpac6tBhjVhpj4ktbR1lCB4BiYHekXxtj\n0owxu4wxHxljepa2roIQkUUicmFp68iFoxdS8hosRGS6iFwfHllBeQh4TkSiRGR5CbedLyLyhYic\nGrDpbuBzEakjIi8A/wYaAfVEZGhJ6zPG+IwxbUq6XeVEdAAoIsaY8cB/gEeAxkBL4EVgQGnqKghj\nTKXS1lBMDA4HizATA3wXbKcxxpSgloKIAZJzrW+TIrwFGqLrxw2/nyPKwf9L/oiILoVcgCggDRic\nT5lIYAawC/gVeAaoYu/rC+wE7gL22mUuAy4CtgL7gHsC6poCvAUsAQ4BicDpAfsnAj/a+zYD/wrY\ndzXwBdZgtQ/rzvVqYE1AGR8wFtgG7AdeCNgXATwN/AH8BNxsl48I8r0n2t/3EPA9cI7D87Ejl542\nAeuv2LprABnAcfv8HwKa2OfntYDyA+3zsB/4HDglYN8vwARgE3AAWAxE2vsaAB/Y2/8EVgX5jj/a\nGjJsDVWAlVg3A18AfwFtgKbA+3Zd24Drcv2mbwKv2XVsAtoBk+xrIgU4P5/rKxZIAg7a18Vi4KHA\n68v+/D9b62G7nUXAUSDTXr/GLhePNaD9CXwCtMz1e9xkf4ef7G2nACvs8t8Dl+f6vV4APrTbWA+0\ntvetsutLt/ddnsd3awskAKnA78DigH1nAxvs3+groEfAvpX294i093cM2NfQ/r0a2uuXAh673BfA\nabmukbvt3+Qw1v9Antd1WV9KXUBZXID+9j9Qnp2gXeYhYB1Wp9IAWAs8aO/rCxwD7gMqAdfZF/rr\nWJ1cR/tijbHLT7H/aQfZ5ScAPwOV7P1DgJPsz5fb/1z+9avttm6yL+Sq9rbVAVp9wHKgNhBta7nA\n3ncDVmfaFKgD/B+Qldd3B9oDOwLabhnwj1/Q+QgcALLIYwDIq2zA+Xk1QEM6cK59ru4CfgAq2/t/\nAb4ETgLqYnV619v7HgVm2uepEtAzn9/3l8BOAKvz2Y7VMUYAlbE6u+exBogz7PPaL0BzBnC+XX6h\n/ZveE3BN/Byk7Sp2W7fZZYdgXY95niNbW3xe58tevwyrc29va7kXWJvr+vjMPl9Vsa7RHcAorCey\nM7BuEE4J+L3+AOLs+l4HFuWqr3U+53YR9g0QVmd+tv25HtagPtyud5i9Xi/39wReBh4OqPMm4GP7\ncyzWINvV1j/S/j39NyS/ABuBZvb3DXpdl/VFTUBFowGwT0R8+ZQZjtXB/SkifwIPYl1ofjKBR0Uk\nC+sOriEwQ0QyROQ7rI7pjIDySSLyrl3+P0A14CwAEVkmInvtz29hdXjdA47dJSIzRcQnIkeD6J0u\nImkishPrH6mzvf1y4FkR2SMiB4HH8vnOWVj/sP8wxlQWkR0i8ovD8xFIccwnVwAfisjn9rl6CqiO\ndefo51kR2SsiqVh3/P7vegxroGstIlkisraAtnLrXCAiW+zroond5kQROSYim7A6pVEB5deIyH/t\n8m9hXQOPBVwTMcaYqDzaPQtrQHvO1rkM+LoArfkxFuv332ZreQzobIyJDijzqIik2tfPpcAvIvKq\nWGwClmFdK37eFZEku743+Psc+8nvNz6G9d2bi0imiKyzt1+CZbpaZF/LS4At5G12XQxcGbA+3NYB\nMAaYLSKJtv7XsG6wzgoo/6yI7La/b37XdZlGB4Ci8SfQ0BiT3/lrhnXX4CfF3pZdh9i3E1iPmWDd\nIRKwrVbAevbEp33cr/76jDGjbC+PA8aYA0AnrM7khGPzYW/A54yAtpvlOj5oXSLyEzAOmArsNcYs\nMsY0Cagnv/MRKprZdfs1CZbm5gFlgn3XJ7HMXCuMMT8aYyYWsu3Ac9MM2C8iGQHbUvLRcRjrpiLw\nmjDkvAYC696Va1tKHuWcEgM8a4zZb4zZj3V9Sy6tv+Yqf5a/vH3NDcd6qvLzW8DnwHPshLuw+qYN\nxphvjTHX2Ntz/LY2uc+pn5VAdWNMN2NMDNbN1HsB+ifk0t+CnNdj9vcNcl03LcT3cS06ABSN9Vh3\nDP/Kp8wurAvNTwywuxhtZt+N2ROMLYDdxpiWwFzgJhGpJyL1sCb8Au+wijPptsduy0/L/AqLyBIR\n6c3f3/1x++9unJ+PDCwzg58mAZ8L+i652wHr3P2aR9kciEi6iNwpIm2x5hHGG2POKei4INp2A/WN\nMTUDtrXkxI67KOzhxE4v39+lAHYAY0Wkvr3UE5FaIvJlQJnA77YTSMhVPkpEbimGhr8bEvldRK4X\nkeZYJsiZttfQbqBVruJ5nlP7yeNNrIHpSqynwr8C9E/L4/sGusRKrvpyX9f5PQmXGXQAKAIicgjL\njvqiMeYyY0x1Y0xlY8xFxhj/hbEEuN8Y09AY0xCYjDXhV1TijDH/sr0S7gCOYNmya2LZVPcZYyLs\nu6V/FKOd3LwJ3G6MaWaMqYs1OZYntr/5OcaYSCwT12FbG1iP5E7PhwcYbn+fC7Fs2n72Ag2CmEb8\nei+xdVQ2xtyJda7WF/RFjTGXGGPa2qtpWJOn+Zn5giIiv2LNeUw3xlQ1xpwOXEvxrgE/64Hjxphb\n7e84mJwmv8IyB7jXGNMRwBhTxxjz73zKfwi0N8aMsNuvYozpaozp4LC937AmyfPEGPNvY4x/gEvF\n+g18wMdAO2PMMGNMJWPMUOBULDNeXiwGhmINAosCtr8E3GCM6W63V9MYc3GuwTpQT37XdZlGB4Ai\nIiL/AcYD92OZbnZgTTT5HzMfwfLW+QbLmyARmJZflQWsv491MR8ArgIG2fbf77G8dL7E+sfqhOXV\nUKivk8/6S1jeHt9geZ18BBwPMv9RFevO6A+su7VGWJOaULjzMQ7rDvwA1t3bu9nCRLZi/WP/bD++\nBz4dICLbgBFYXih/YNmNB4jI8SDfNZB2wH+NMWlYk9QvisiqIGUL+r2wtbfGOhfLgMkisjKf9gtq\nw9oocgwYDFyDZa653K6/UPUE1Pce1u+2xBiTivUbBb4nkvtuOB24AGsSdre9PIb1+zthKvCq/fvl\nNdB0A74yxhzC+n+6TUS2i8h+rPmHO7E82u4ELhGRA0F0bsDyyGqK5dnk356ENQ/wgm3y2oblGJHn\n9yX/67pMY/42OYapAWO2Y7mq+YBjItLdGFMPWIr1OLUduMKeYFTywBgzBWgrIqMKLBx+LRcCs0Sk\ndWlrURSleJTEE4APy/UtVkT8j6mTgP+KSAcsP+1yMZqWR+ywARfZj9zNsUxf75S2LkVRik9JDAAm\nj3Yuw/J7xv6b32SqUroYLJfN/VgmoGSsQUBRlDJOSZiAfsaayMkC5ojIy8aYA7a3ir/MfhGpH1Yh\niqIoSg4ql0AbPUVkjzGmEZZ/9VacTaApiqIoYSTsA4CI7LH//mGMeQ/LXW2vMeYkEdlre3H8ntex\nxhgdGBRFUYqAiBT4Rn1Y5wCMMTWMMbXszzWxXMe+xYo7M9oudjWWi2OeiAviZQQuU6ZMKXUNZUGT\nW3WpJtVU3nRtXLyY06tVo1XlyszGCujllHBPAp8EfGGM8WD5qX8gIiuw3g79p20OOo8y9Fbd9u3b\nS1vCCbhRE7hTl2pyhmpyTqnpyszkmf796T98OBNGjuSDL79kYZMmjI2MdFxFWE1AYgVMyh0ECrFe\n6Dg/nG0riqKUW7xeGD2abjVq4E1MpFmXLgB8sWsXHo+Hp7p2dVSNvglcSEaPHl3aEk7AjZrAnbpU\nkzNUk3NKVFdmJkyZAhdcAOPH02vt2uzOHyAiIoK4uDjH1YXdDbQ4GGPEzfoURVFKDPuunxYtYO5c\naBY8mK4xBintSeBw0apVK4wxuugS9qVVq1Zhv54TEhLC3kZhUU3OCbeuzPR0pvTtyzO9esH48fDB\nB/l2/oWhJN4DCDkpKSnok4FSEhjjptS+SkXDs2QJo6+5huioKOauXg0B5p5QUCZNQMYYHQCUEkGv\nNaU0yExPZ9ollzBrzRqeuu46Rs6ejYlwbrCxr9sC717K5BOAoihKucXrZdz557OjUqUcHj7hoEzO\nAShKecKNtm3V5JyQ6Qrw8Hl82jQ+2LMnrJ0/6ABQKtx4441Mm5Zfbhh3EhERwc8//1zaMhSl/OH1\nQvfukJQEXi+1x44tlMmnyJT2a9T5LZa8Ewm2PSsrSxITEyUxMVGysrLyLFMQxa0jJiZGqlevLlFR\nUVKvXj3p2bOnzJ49W3w+X5H0FIYFCxZIr1698i2TnJwsF1xwgdSvX1/q1asnXbt2lU8++cRR/RER\nEfLTTz+FQmqZIdi1piih4Ghamuy7806RRo1EFi4UCVE/YV+3Bfax5eYJwONJJi5uHH36pNCnTwpx\ncePweJJLvA5jDB999BEHDx4kJSWFSZMm8fjjj3Pttdc6Oj4rK6tQ7QUiIgV6rQwYMID+/fuzd+9e\nfv/9d5577jmiooKl1z2xfkVRQoNnyRK6NWrEzOXLrSeAUaOgpL3OnIwSpbXg8AkgKytLOne+VSBL\nQOzF2ub0Lj4UdYiItGrVSv73v//l2LZhwwaJiIiQ5ORkEREZPXq0TJ48WUREEhISpEWLFvL4449L\nkyZNZNSoUSIi8sEHH0jnzp2lbt260rNnT/nmm2+y69u5c6cMHjxYGjVqJA0bNpRbb71Vvv/+e6lW\nrZpUrlxZatWqJfXq1TtB2759+yQiIkIOHjwYVP/cuXPl5JNPlgYNGshll10mu3fvzt5njMl+Ahg9\nerTcfPPNcskll0jt2rXlrLPOkp9//tnxeSorBLsGQ8nKlSvD3kZhUU3OKayuo2lp8kCfPtLIGFk4\nZoz4imityA8q0hOAx+Nh27Z+5JzSiGDbtr54PJ4SqyMY3bp1o0WLFqxZsybP/b/99hupqans2LGD\nuXPn4vF4uPbaa3nppZfYv38/Y8eOZeDAgRw7dgyfz8ell15K69at2bFjB7t27WLYsGGccsopzJ49\nmx49epCWlsb+/ftPaKdBgwacfPLJXHXVVbz//vv8/nvOKNyff/459957L2+//TZ79uyhZcuWDBs2\nLOj3Wrp0KQ8++CCpqam0bduW++67r1jnSVHKO/67/qQtW/AmJjJq7tySsfUHoVwMAMHIyICuXa2n\nKmNg6tS8y02dapXLyAiflmbNmuXZKQNUqlSJBx98kCpVqlC1alVeeuklbrjhBrp27YoxhpEjR1K1\nalW+/PJLNmzYwJ49e3jiiSeoVq0akZGRnH322Y51rFy5ktatW3PnnXfSrFkz+vbty08//QTAokWL\nuPbaaznjjDOoUqUK06dPZ/369ezYsSPPugYNGkRcXBwRERFcddVVeL3ewp8YhX79+pW2hBNQTc5x\npMv28Fl53XVW5M4wefj4fD6SkpIcly8XA0BsbCzt2ydg5Z/346Nz51VkZcVmG3TyGwCysmLp3PnE\nOtq3X0VsbGyxNe7atYv69fPOetmoUSOqVKmSvZ6SksLTTz9N/fr1qV+/PvXq1ePXX39l9+7d7Ny5\nk5iYGCKKeNfQrFkznnvuOX744QdSUlKoWbMmo0aNAmD37t3ExMRkl61ZsyYNGjRg165dedbVpEmT\n7M81atQgPT29SJoUpVwT4OEzftu2sN31B85hOqVcDAARERHMnz+Wzp3HUaPGMmrUWMYZZ9zO/Plj\nHXeUoagjGF9//TW7d++md+/eee7PPXEbHR3Nfffdx/79+9m/fz8HDhwgPT2doUOHEh0dzY4dO/D5\nfAXWUxDNmzfn5ptvZvPmzYA1OKSk/H3x/PXXX/z555+0aNGiUPUqhcON/u2qyTlBdeWK3BnKGD65\n8fl8xMfPweudQUbGYMfHlYsBACA2thNJSTNYvboVq1e3YuPGZ4mN7VTidQSSlpbGhx9+yJVXXsnI\nkSPp2LGjo+PGjBnD7Nmz2bBhA2B1xB9//DF//fUX3bt3p2nTpkyaNImMjAyOHj3KunXrADjppJP4\n9ddfOXbsWJ71pqamMnXqVH766SdEhH379jF//nx69OgBwJVXXskrr7zCN998w9GjR7n33ns566yz\niI6OLvI5UJSKiGfJElZ27Jjt1x9uD5+85zALptwMAPB3LGy/Xbq06hgwYAB16tShZcuWTJ8+nTvv\nvJP58+c7Pj4uLo6XXnqJW265hfr169O+fXsWLlyYre+DDz7ghx9+oGXLlkRHR/Pmm28CcO6559Kp\nUyeaNGlC48aNT6g3MjKS7du3889//pM6depw+umnU61aNV555RUAzjvvPB5++GEGDx5M8+bN+eWX\nX1iyZEn28YV5wvjHP/7B4sWLHZevyLjRtq2aCsZvb69du3b2E7k/cmf/4cP586KLwnrXHwo0GJyi\n5INea0peJHs8zImPp+/Wrbxx7BjNO3Wi79ChPPzQQ1bkzk8+CXsYBwCfD954A957z8dPP41j06YZ\nWPf15TgfgKKUJ9xo21ZNwfH5fMyJj2eG14scPsynx4+TsWkTo+69lztGjCiRGD4AK1ZY0aFnzoQ7\n7ojglVf+nsN0ikYDVRRFKQQej4fM5GROA44Bh4H/AnWN4YfGjcPu179pE9x1F2zfDtOnw+DB/ukF\naw7T4/HgMCWwmoAUJT/0WlNyk5SUxPbevfEdPswEYCcQDVwRGcmwtWvp6rT3LSKvvw6HDsGYMRDg\nPZ4DzQegKIoSBmKNYaEIvYFUoCPWIPBd06Z0KQHTz4gRoatL5wAUpZRxi207ENV0In4Pn4d792bs\nPffwTJMmjI2M5IrISHrGxNBh4MBivzMUyJEjcPx4yKrLEx0AFEVRCiAwhs+YNWvo9MADfLFrF8PW\nraPJc8/x0c8/88xzz4WkLZ8PXnsNOnSA//43JFUGRecAFCUf9Fqr2BQ3N29hWbEC7r4bqleHJ5+E\nXr2KVo/OASiKohQHr5f7LriA740Je27e/fth2LC8PHvCi5qAyiD5pWZctGgRF154YcjbXLhwYdBY\nRkrxKG3bdl5UaE0BMXweevTRAv36Q6Grbl0YPRqSk2HIkJLLC1OuBgD/q9lJSUl5BksrqTrOOeec\nQoV+KCz5hWQYPnw4n376adD98+bN49RTT6VOnTo0bdqUSy+9lL/++qvY7SpKuSBXbt7q111XIvH6\nIyJg+PDgbp1ha7dkmwsfyR4P4+LiSOnTh5Q+fRgXF0dyIRO5hKKOkqCoNulVq1Zx3333sXTpUg4e\nPMj333/P0KFDQ6xOKSxui3EDFU9TZno6v02YUKTInYXRdeQIfP11EUWGgXIxAAS+mj04I4PBGRnM\n8HqZEx/v+C4+FHXkJjU1lQEDBtC4cWMaNGjAgAED2L17d/b+AwcOEB8fT/PmzWnQoAGDB/8dxvWl\nl16iXbt2NGzYkH/961/s2bMnR90fffQRbdu2pXHjxtx9993Z2/Mz1SQmJnL22Wdz+umnA1C3bl1G\njhxJzZo1ATh06BCjRo2icePGtG7dmmnTpgX9bhEREcyZM4f27dtTv359brnllsKfIEVxAX4Pn+fe\nfz9skTsDPXtmzgxp1cWiXAwAHo+Hftu25UrmCH23bStUSsji1pEbK0Z3PDt37mTHjh3UqFGDm2++\nOXv/iBEjOHz4MN9//z2///47d9xxB+AsNeN7773Hxo0b2bhxI++//34Ok1MwU82ZZ57JZ599xtSp\nU1m3bh2ZmZk59t9yyy2kpaWxfft2EhISePXVV7MjhebFRx99RFJSEps2beLNN99kxYoVhT5HSgW3\ntxeCUGsKjNw5YeRIpm3bVqTInQXpCozZ88YbkM+/VMnjJHFwaS04TAqfmJgoy2rUkIBs7iIgb4Mk\n5toWbEkEWZbH9rdr1JDExMQCUjDnpF+/fjJv3rwTtns8Hqlfv76IiOzevVsqVaqUZ4L2a6+9ViZO\nnJi9np6eLlWqVJGUlBQRsZKzr1ixInv/zJkz5fzzzxcRkQULFkjv3r2Davv0009l4MCBUq9ePald\nu7aMHz9efD6fZGVlSWRkpGzZsiW77Jw5c+Scc87Js15jjKxbty57/YorrpDHH388/xNTBgl2DYYS\nNyY7L++aNi5eLKdXqyaXNG4su5KSilVXfromTRJp107k7bdFfL5iNVMoqEhJ4WNjY0lo3z5XMkdY\n1bkzsVlZjoaA2KwsEjp3PrGO9u2LnBLy8OHDjB07llatWlG3bl369u1LamoqIsKvv/5K/fr1iYqK\nOuE4J6kZA7N0xcTE5DAt5Uf//v15//332b9/P++//z4LFizg5ZdfZt++fRw/fpyWLVvmqDdYOkiw\nEtD40ZSQRaei2duLSkg02R4+SWPGFDs3b175AHJz220l79lTGMrFABAREcHY+fMZ17kzy2rUYFmN\nGtx+xhmMnT+/UCkhi1tHbp5++ml++OEHvv76a1JTU1m9ejVgPXVFR0ezf/9+Dh06dMJxTlIz7ty5\nM/vzjh07aFaER9dzzjmHc889l82bN9OwYUMqV66co92UlBSaN29e6HoVxZUEePhct3VrsXLzBubf\n7dMnhbi4cXg8ySeUa9q05D17CkO5GAAAOsXGMiMpiVarV9Nq9Wqe3biRToW8cw9FHYGkpaVRvXp1\noqKi2L9/P1MDstI3adKEiy66iJtuuonU1FSOHz/OmjVrAGepGZ988klSU1PZuXMnzz777AlzBHmx\nfPlyli5dSmpqKgAbNmxg1apV9OjRg4iICK644gruu+8+0tPTSUlJ4ZlnnmHkyJFF/v6KMyqCvT0U\nFFlTiHPz5s6/m5FRH693BiNGzCmyw0hpUW4GAHBPSkiwJmLvuOMOMjIyaNiwIWeffTYXX3xxjjKv\nvfYalStX5pRTTuGkk07i2WefBZylZrzsssuIi4ujS5cuDBgwgPj4+AI11atXj5deeon27dtTp04d\nRo0axcSJE7MHj+eff54aNWrQpk0b+vTpw4gRI7jmmmuCfr9gfPHFF3mathSlpPEsWcKHp5wS0ty8\neeffjeDnn/sW2WGktNBYQGEgLi6OKVOmMHDgwNKWohQTt19rSt4ExvB5/sYbGfrCCyEzwiclJdGz\nZwpHjw7Osb1GjWWsXt2KuLi4kLRTHJzGAipXTwBuIDk5mS1bthR54lhRlOIRGLnTm5jI0BdfDOkM\nbP36sUAC5HIZad9+VZn7v9cBIIRMmjSJCy+8kCeeeCKHvV5R8qNc2dvDSIGaMjOZe+ml2X794crN\n27p1BOvW/Z1/t2rVBznjjNuZP39sSPMBlAQaDTSEPPbYYzz22GOlLUNRKh5eL4weTa86dcIeuROg\nS5e/8+8mJu5jzJjJZa7zB50DUJR80WvN5WRmwrRpMGsWPPUUjBwZEnPPkSPwwgvw009W1WUNnQNQ\nFKVcECxCr3g8OSJ3hsLDxx+z55RTYO1auP324qp3NzoAKEopUybt7SWEP0Lv9t69Gde9O7d36YJ3\n7Vqm9O3LhJ49Q+LX7+e//4W4OCtmz+uvw7vvWgNBQbjlXBUFnQNQFMWVBEbofQdIAgZs2sT5vXtz\nZsOGvPTFF1aUtRCxeTPcf3/JZeNyAzoHoCj5oNda6ZGUlMRLPXqw5tgxMoEfsUwWdYAbJk3i0enT\nS1egi9E5ACUsTJ8+neuvvz7o/tatW/P5558XWE9KSgoREREhf3U+XPUqpcM/K1dmKvC7vd4EiI+M\nZPCQIUWuMy0tBMLKCToAhIElS5Zw1llnUatWLZo0aUKPHj2Y5SJXgmPHjjFhwgSio6OJioqiTZs2\njB8/3tGx99xzD3Pnzg2JjoLCSfTs2ZO6devSsGFDevfuTVJSUrHrdSNutCG7QVNsp06sskOKCBAD\npAHfNW1KlyKYfo4csRyF2rSBLVtCp9MN56qolMgAYIyJMMZ4jDHL7fUpxphfjTEb7SVkWcxFhEmT\nJhXrsb04dTz99NPccccdTJw4kb179/Lbb78xe/Zs1q1bx7Fjx4qsKZQ8+uijbNy4kcTERA4dOkRC\nQkKR/qHCRVpaGgMGDOD222/nwIED7Nq1iylTplC1atXSlqaUFF4vEWedxdiTT+aZRo0YGxnJNZGR\n9IyJocPAgYXyuQ/MxrV2LaxZ42xyt0LgJGlAcRfgDuB1YLm9PgUY7+C4/JId5Mlbb70ltWvXlrff\nfjtomYIoah0HDx6UmjVryrvvvptvuY8++khiY2MlKipKWrZsKVOnTs3el5CQIC1atMhRvlWrVvK/\n//1PREQ2bNggXbt2laioKGnSpIlMmDBBRESOHDkiI0aMkAYNGkjdunWle/fu8vvvv+fZ/qWXXirP\nPvtsUH27d++WIUOGSKNGjaRNmzby3HPPZe+bOnWqjBgxInv91VdflZiYGGnYsKFMmzYth1afzyfT\np0+Xtm3bSsOGDWXo0KFy4MABERHZvn27RERESFZW1gntJyYmSr169YLq8/l88vDDD0tMTIycdNJJ\ncvXVV2cn1sldb79+/WTy5MnSs2dPqV27tvTv31/+/PPPoHXnJr9rTQk9R9PSJOX220UaNRJZuFDE\nTlaUmJgoiYmJeV4v+fHDDyJnnCFy1lkia9aESbQLwWFCmJLo/FsA/wf0yzUATHBwbH5fLgezZ8+W\njh07Srt27QSQdu3aSceOHWX27NmOT1px6/j000+lSpUqBV6kq1atks2bN4uIyLfffitNmjSR999/\nX0SsASA6OjpH+cBOtUePHvL666+LiMhff/0lX331lYhY2bsGDhwoR44cEZ/PJxs3bpS0tLQ823/k\nkUekZcsR8yaQAAAgAElEQVSWMnPmTPn2229z7PP5fBIXFyePPPKIHD9+XH755Rdp27ZtdgayqVOn\nysiRI0VEJDk5WWrVqiVffPGFZGZmyvjx46VKlSrZWmfMmCE9evSQ3bt3S2Zmptxwww1y5ZVXikj+\nA8ChQ4ekYcOGcvXVV8snn3ySPWj4mTdvnrRr1062b98uf/31lwwePDhbU14DwMknnyw//vijHDly\nRPr16yf33HNPvr9PIDoAlBz+LF13tG4tsmtXSOr86y+Rd94p2WxcbsBNA8BbQGegb64B4BfAC7wM\n1AlybH5fLgc+n0/efPNNiY6OFiyT4QnLlClT8qxvypQpeZaPjo6Wt956S3wOr57XX39dmjZtmmPb\n2WefLXXr1pXq1avLmiC3IOPGjZPx48eLSMEDQN++fWXq1Kmyb9++HGXmz58vPXv2lG+++aZAnT6f\nT2bOnCm9evWSatWqSfPmzWXhwoUiIvLVV19JTExMjvLTp0+X+Ph4Eck5ADz00EPZHbqINSBFRkZm\naz311FPl888/z96/e/fu7AEyvwFARGTLli1yzTXXSHR0tFSpUkUGDhyY/URz3nnnyaxZs7LLbt26\nNWi9/fr1k2nTpmWXnTlzplx00UUFniM/JTEAlPf0iwVxNC1NHujTRxoZIwvHjBFfkGvCjedJxJ26\nnA4AYX0PwBhzCbBXRLzGmH4Bu2YCD4mIGGMeAf4DXJtXHaNHj6ZVq1YA1K1bl86dOwdrC2MMqamp\ndOzYkZ07d/LKK68wxIG3wNSpU7OTtbz99tvEx8cTHR3Nzp07s+t1QoMGDdi3bx8+ny/bRrl27VoA\nWrZsme2Z8tVXX3HPPfewefNmMjMzyczM5PLLL3fUxrx585g8eTKnnHIKbdq04YEHHuCSSy5h5MiR\n/PrrrwwbNoyDBw8yYsQIpk2bxvr167noooswxhATE8O3336LMYYbb7yRG2+8kaNHjzJv3jzi4+M5\n88wzSUlJYdeuXdSvXx+wbhB8Ph99+vQ5Qcvu3btzBL2rUaMGDRo0yF5PSUlh0KBB2edCRKhSpQp7\n9+4t8Ht26NAhO9H9tm3buOqqqxg3bhxvvPHGCSkzY2JiOH78eNB6mzRpkkNjUVNX+if7/KkJQ7Ue\n7vrdvP7D55/zwpNPEh0VxczZs2nYvn12lq7c5b1eb9D6DhyAZcsSOPnkkv8+fkrzfCYkJLBgwQKA\n7P7SEU5GiaIuwKPADuBnYA+QDryaq0wM8E2Q4/Mb3U7g0Ucflbffflt8Pp+8/fbbMn369MINm8Ws\nIzU1VWrVqiXvvPPOCftatGghq1atEhGRtm3byrPPPiuZmZkiYj0B+O+qv/76a2nQoEH2ccePH5ea\nNWtm31UH8vbbb0u1atUkIyMjx/aUlBTp2LGjzJ8/37H2hg0byjvvvCPr16+X9u3bBy0X+ATw4IMP\n5vsEcMopp+RIHB9IQU8AuXnhhRfk9NNPF5G8nwAiIyODPgHMmzcvu2zu5PYFEexaU4rJ0aMiDzwg\nS2vXzveuX0TynQM4fFjkySdFGjYUefjhcIsuO+CGpPAicq+ItBSRNsAw4HMRGWWMaRJQbDCwORTt\n3XPPPQwZMgRjDEOGDGHSpEklWkedOnV44IEHuOmmm1i2bBnp6emICF6vl4yMjOxy6enp1KtXjypV\nqrBhwwYWLVqUva99+/YcOXKETz75hOPHj/PII4+QmZmZvf+NN95g37592e0ZY4iIiCAhIYHNmzfj\n8/moVasWVapUCeop8eyzz7Jq1SqOHDlCVlYWCxcuJD09nS5dutC9e3dq167NE088kb0/OTmZxMTE\nE+r597//zYcffpjt4fTAAw/4B24Axo4dy7333suOHTsA+OOPP1i+fHn2/sCygWzdupX//Oc/2Qnp\nd+7cyeLFi+nRowdgpcx85pln2L59O+np6dx3330MGzYsx5OG4mICcvNesWVLvrl5g+Xezcuz5/77\nS/h7lAecjBKhWMg5B/Aq8A3WHMB7wElBjslvdHMtixYtku7du0vNmjWlcePGctZZZ8nLL78sx44d\nExGRZcuWSUxMjERFRcmAAQPk1ltvzb6rFhFZuHChNG3aVE466SR5+umnpXXr1tl31SNGjJDGjRtL\n7dq15R//+IcsX75cREQWL14sHTp0kFq1akmTJk1k3LhxQe+u586dK3FxcVK3bl2pV6+enHnmmfLx\nxx9n79+zZ49ceeWV0qRJE6lfv7706NEju/3AJwARywuoZcuW0rBhQ3n00UdzaPX5fPLMM89Ihw4d\nJCoqSk4++WS57777RCT/J4Bdu3bJFVdcIc2bN5datWpJixYt5MYbb8ye1PZ7AUVHR0vjxo1l1KhR\nkpqamme955xzTr5PAP5J7GCUxLXmRhtyWDTZd/2BHj75kZWVJZ073yqQJSACKwWsbYMHZ7nGs8eN\nvx8OnwA0FISi5ENJXGsJCQnZdl23EGpNniVL2DpxIsNOOw3mznUUvC0pKYk+fVLIyPCnXkwA+lGj\nxjKWLWtF//5xrojZ48bfT0NBKEoZwW2dB4ROU2Z6OlP69qX/8OH4Lr64mJE7/9bUqJF7Ara58fdz\nig4AiqKEhdy5eYfPmlWoXjsmJpbatRMoD7l33YoOAIpSyrgxlkyxNGVmsuCyy4qcm9cfs+fUUyPo\n1Wss//iHu3PvuvH3c4rmA1AUJXTYuXn71K9f6Ny8Ph+88YblzdOliz9mTyd8vrKfe9et6CSwouSD\nXmsOCUFu3tWrYeJEePJJ6NUrTDorCE4ngXUAUJR80GutYMTjwVxzDbRo4djDJ2hd4p7J3bJMufYC\niomJyQ7RoIsu4VwCQ06ECzfakJ1o8nv4jDn77JDl5jX5dFluPE/gXl1OKJMDwPbt20vsBbbcy8qV\nK0ut7bKkya26Cqtp+/btpX25u5JAD5+H1q6FUaMc3brv3w933gkPPlgCIpUCKZMmIEVRSofM9HSm\nXXIJs9as4anrrmPk7NlBwzgEcuQIPP88PPEEDBkCU6ZA06YlILiCYowzE5B6ASmK4gyvl+cvvpik\nrKxCefi89lpuz54w61QcUyZNQKWJG+19btQE7tSlmvLH5/ORlJTEnDlzssOXk5lp3bJfcAG3T5tW\naL/+n36y3Dvffbd4nb+bzlMgbtXlBH0CUBQFgGSPhznx8fTdupXXjx5l86xZ3DBpEp0ee8zy8PF6\nqVyESV471YbiQnQOQFEUfD4f4+LimOH18g4QDzwC/FipEjNefpmIq68ucJL34EGoU6ck1CoF4XQO\nQE1AiqLg8XjITE7mNGACkAY8ALzr8zH5hx/y7fz9nj0dOliDgFJ20AGgkLjR3udGTeBOXaopOP0q\nVaIj8Ku9HgUMrVKFQYMG5VneH7OnQwdITwePJ7xPAG45T7lxqy4n6ByAoiiYrVu55ehRYoAaQANg\nP/Bd06Z0yWPCd906uPJK9ewp6+gcgKJUZOwYPp/MmIH3/PNZvm4dvfbv58xKlZjXuDGnDBzIM889\nd8Jhu3fDzz9rzB634nQOQAcARamoeDwwejRER2fH8Dl+/DhLliwBYNiwYVSurEaCskjIJoGNMZcb\nY2rbn+83xrxjjHHuBFzOcKO9z42awJ26VBN/+/X37w8TJmTH8PF4kunWbTxjx9bguut+plu38Xz0\nUTLbtpWsvGC48bcD9+pygpNJ4MkikmaM6QWcD8wDZoVXlqIo4WDjokW83KYNJCVZsfvtGD4+n4/4\n+Dl4vTPIyBjM0aN98HpncNllc1i/3ldwxUqZpEATkDHGIyKxxpjpwLcissi/Lezi1ASkKCEhMz2d\nRy6+mNlffMEz11/PVbnSM56YgN2ievVlrFnTiri4uJKWrBSDUMYC2mWMmQP8E3jcGFMVdR9VlDLD\nxkWLGH3ttcRERRU6S5eDAJ9KGcZJR34F8BnQX0RSgfrAXWFV5WLcaO9zoyZwp64KpSkzkzcGD+bC\nESO4a9QolucTwyc2Npb27RP4OwG79dlNCdjd+NuBe3U5ocABQEQygN8Bv8PXceCHcIpSFKWYeDzQ\nrRvnHDqENzGRkXPm5Bu2OSIigvnzx9K5sz8B+ypXJmBXQouTOYApQFegg4i0N8Y0A94SkZ5hF6dz\nAIpSOBzk5k1JgXnzrCBtuft2n8+Hx+MBrKcC7fzLJqGMBTQIGAj8BSAiu4HaxZOnKEqo8SUlQbdu\nJ3j4+PHH7PFbgY4fP7GOiIgI4uLiiIuL086/AuDkF860b8MFwBhTM7yS3I0b7X1u1ATu1FUeNflz\n8w7v2dPq4XPl5j1yBJ588u+YPZs3w0MPQWRk+DSFAzdqAvfqcoITL6A3bS+gusaYMViRYl8KryxF\nUZzgXbKE0ddcQ4uoKOauW/f37X0AixdbsXs0Zo+SG0ehIIwx/wQuAAzwmYj8X7iF2e3qHICi5EGO\n3LxjxjBy1qygk7wi6s5Z0QhpTmC7wy+RTl9RlALweJg/YABJx4458uvXzl8JhpNYQGnGmEP2csQY\nk2WMOVQS4tyIG+19btQE7tRVpjUFxPC5Pldu3pQUy+Fn/vwS1lSCuFETuFeXE5y8B1BbRKJEJAqo\nDgwBZoZdmaIof2P79bNxI3i9RFx9NSYiIodnT5s2cPnlpS1UKUsUKRy0xgJSlJIhMz2dH+66i07L\nlsHTT8OIEWAMx4/DjBnw+OMwZIj1YNC0aWmrVdxCyOYAjDGB0aEisF4KO1IMbYqiOMDv4RPXqBHz\nvN4crp0REbB3L6xeDaeeWooilTKNk/cABgQs/bHyRV8WTlFuxo32PjdqAnfqKgua/H79FwwfzvhR\no3h5+/YcnT9YA8CTT4av8y8L58ktuFWXEwp8AhCRa0pCiKJUNPxhF7Zu3UqfPn2IiIjgm6VLGTV6\nNC0CInceOAD16pW2WqU8EnQOwBjzPPbbv3khIreFS1SABp0DUMolyR4Pc+Lj6bt1K28cO0bzU0/l\nhrPOYv+bb/LL0KGMnDWLHTsjuP9+SEy03t6tVKm0VStlhWLnBDbGXJ3fgSKysIjaHKMDgFIe8fl8\njIuLY4bXyztYr9bPB1ZHRTEjOZnUGi149FF45RW45RbLy6e2Rt9SCkGxg8GJyML8ltDKLTu40d7n\nRk3gTl1u0OTxeDiyeTOnAfdiTardC7yTlsa/b5uVI2bPgw+WTufvhvOUGzdqAvfqcoITL6BGwESg\nI1DNv11Ezg2jLkUpt/h8PmKOZ5EJfG5vOwIMFkPPoZcxbZp69iglg5N8ACuApcCdwA3A1cAfIjIx\n7OLUBKSUMzLT05nYpw+zPR7igdeAaGAnEE1DFmz4mG7dupWuSKXME8p8AA1EZB5wTERWiUg8oHf/\nilJIvEuW0L1RI75LSaFhpdv5iCZcTSR3U53mxPBjpW4ag18pUZxcbcfsv3uMMZcYY2Kx8gJXSNxo\n73OjJnCnrlLRlJnJsiuu4ILhwxkzaBQN++/lN+MjhZ28wDpG8wxb+JGOp52s+XfzwY2awL26nOBk\nAHjEGFMHmIBlBnoZuCOsqhSlvGDH8Dl730FGXZXIA5/N4eR2lVm5ciydO4+nRo3tVK36G2eccYfm\n31VKHCdzAI1E5I8S0pO7bZ0DUMomAbl5Dz/yNG2njGDgZSZHzB7Nv6uEi2K/BxBQ0TZgO9ZE8Dsi\nciAkCh2gA4BSFslKTKTStddCy5YwZw40a8Yff0CjRqWtTKkohGwSWETaA/cDnYAkY8yHxpgRIdBY\nJnGjvc+NmsCdusKpyR/DZ6A/N+/y5dkxfPLr/CvaeSoqbtQE7tXlBEfPnCKyQUTGA92B/UChXgQz\nxkQYYzYaY5bb6/WMMSuMMVuNMZ/ZcwyKUmbxLllC1waNWJ20hZfWr7eys2gqLsXlODEBRQGDgGFA\nW+Bd4E0RSXLciDF3AHFAlIgMNMY8DvwpIk8YYyYC9URkUh7HqQlIcTWZ6elMPu8S5m9YQ+fqYxj0\n5CxuvClC+36lVAnlHMAvwHtYnf76IghpAbwCTAPG2wPAFqCviOw1xjQBEkTklDyO1QFAcS0HEzy8\nNvAyXks/ygXXfcLdT3fRmD2KKwjli2BtROSOonT+Ns8Ad5EzsuhJIrIXQER+AxoXse4Sx432Pjdq\nAnfqCokmOzdv5Uv6UyPuEd7duYeH5xa98y+35ynEuFETuFeXE5zkAyjyLbgx5hJgr4h4jTH98msm\n2I7Ro0fTqlUrAOrWrUvnzp3p18+qyn/iS3Ld6/WWavt5rftxix7/utfrdZWekPx+P/xAvxdegJYt\n2bBgJm0aNaRZ84hi6fPjhvPj5nU3Xk+BlKaehIQEFixYAJDdXzqhSDmBHVduzKPACOA4VkL52lhz\nCF2BfgEmoJUickL4KzUBKW4hMz2dzePH0+W993Lk5lUUNxIyE5AxpqeTbXkhIveKSEsRaYM1ify5\niIwEPgBG28WuBt53Up+ilDQeD1x+muXh85+PPwavVz18lHKDkzmA5x1uKwyPAf80xmwFzrPXywS5\nH/vcgBs1gTt1OdWUkgIjh6Yz8cy+JCQP546Ro3htx44TcvOWpKaSRDU5x626nBB0DsAY0wM4G2hk\njBkfsCsKKHRyOhFZBayyP+8Hzi9sHYoSbtLSrCQsCbPf5MjRq4mpF8WmT63cvIpS3sgvJWRfoB9W\nDoDZAbvSgA9E5Iewi9M5AKWEyUjNZO3F06j/3fMkD72ckbNmYSI0Ro9StgjlewAxIpISMmWFQAcA\nJVQ4Crzm8cDo0Tli+ChKWSSU7wEsMMZ8nnsJgcYyiRvtfW7UBO7R5fEkExc3jj59UujZ8yPi4sbh\n8STz5592Aduvn/79T4jhUxK45TwFopqc41ZdTijwPQCsHAB+qgFDsNw6FcX1+Hw+4uPn4PXOwLrf\nqY/X24c+fcbRu/cMHh31Jh+MG8fkbt0sDx+961cqEEV6D8AYs0FEuodBT+521ASkFIukpCT69Ekh\nI+NfgMfeGkuNyou4/rRneMPr4akxYxg1e7a6dirlBqcmoAKfAIwxgekfI7CCumn0TqXMUCnrZ3oR\nxxi28ibH2EdT9h7fw3cpUXgT1cNHqbg4mQNIAhLtv+uxUkNeG05RbsaN9j43agJ36DrjjDM4K+Jx\nVuGlBof5P47zHTtpVKkSH+7Z44rO3w3nKTeqyTlu1eUEJwlhWotIG/tvOxG5QES+KAlxilJcNm3a\nRKtjqZwG3AtkAo2A3b5Mpj74YOmKU5RSxokbaDXgJqAXVtC2NcBsETkSdnE6B6AUkv/+F84+G2rU\nsNaTvvyS7b164cvKYgKwE4gGroiMZNjatXTt2rUU1SpKeAilG+irWOkgnwdesD+/Vjx5ihJaPB64\n4AK46SbYvt3admzDBmKvv55VNWsCkAp0tP9+17QpXVxg/lGU0sTJAPAPEblWRFbayxisQaBC4kZ7\nnxs1QcnoSkmxYrNdfDH861+QnAwnt7Ry857fqxdmwgTGrlzJM02aMDYykisiI+kZE0OHgQPzfhms\nFHDj76eanONWXU5w8h7ARmPMWSLyJYAx5kysSWFFKVV++gm6d4dbboFt26B2bSs37+hrrqFFVBSL\nv/wS06ULnYAvdu3C4/GQmJjIR2PGuKbzV5TSxMkcwPdAB2CHvaklsBXrZTARkdPDJk7nAJQC+PNP\naNDAitc/7ZJLmLVmDU9ddx0jZ8/WGD5KhSVk7wEAF4ZAj6KEhQYNAI+HzwYPJikjQ/36FaUQOLlF\nekREUgKXwG3hFug23Gjvc6MmCJ2uFStg5sw8dgTE8Ll06lQ+cODX78ZzpZqc4UZN4F5dTnAyAOSY\n8DXGVMZ6G1hRworfs+eWW/II0ePxQLdukJQEXi/m6qvV5KMohSS/fAD3YL07Ux3IAPz2pExgrojc\nE3ZxOgdQIUlJgfvvt3z6J0+GMWOgShVrX2Z6Oom33cbZH34ITz2l6RkVJQ9CmQ9gekl09kHa1gGg\nAnL99dC0qRWZuXbtv7f7PXza1a3Lm4mJmObNS0+koriYUL4I9okxpk/uJQQayyRutPe5URMUXdfc\nuVZaRn/nn5lu+fVfMHw440eO5M1du4rc+bvxXKkmZ7hRE7hXlxOceAHdFfC5GtAdKzDcuWFRpFQY\nRAq23ny3bBnDR4ygRZRG7lSUUFPofADGmGhghogMCY+kHG2pCaicsmIF3HsvLF0KbdvmUSAzE6ZN\n48cXXmDdkCHq168ohSCU7wHk5lfg1CIcpyh4PDBxohWvZ/p0aNMmSKHRoyE6mpO//ZaTNUuXooSF\nAm+pjDHPG2Oes5cXsKKBbgy/NHfiRnuf2zT5fD6SkpKYM2cOPp8PgD17YNQoK2bPoEFWzJ4hQ3KZ\ngAJz806YAB98EPIUjW47V6CanOJGTeBeXU5w8kztTwbjTwgzUURGhFWVUmYJTMB+++17sxOwi1im\nnm3b4MYb/3br9ONdsoS7mjdHEhOt3LyjRql7p6KEGaf5AE62V38siTwAAW3rHEAZwufzERc3LiAB\nO4CPzp3HkZQ0I88AbBrDR1FCT7HdQI0xlY0xT2DZ/Bdi5QXYaYx5whhTJdhxSsXF4/GwbVs/cl5W\nEWzb1hePx3NCee+SJXRv1IikLVvwJiYyau5c7fwVpQTJ77/tSaA+0FpE4kSkC9AWqAs8VRLi3Igb\n7X1u0fTll3D4cOCWhLwLZmbyv5Ejs/36ncTwCRVuOVeBqCZnuFETuFeXE/IbAC4FxohImn+DiBwC\nbgQuDrcwpezgj9kzY0YsMTEJgC9gr4/27VcRGxv7d+Fu3ei1b5/e9StKKZNfLKBtItK+sPtCic4B\nuJ/bboO33oIHHoDrroPNm5O55prZbN0aDUD79jtYsOBGYju1g2nTYNYsjeGjKGEmFO8BfGeMGSUi\nr+aqeASwpbgClfLB0KFWv+4P2xBJJr1Zw72yhUXHj9OCTsg33WD0cIiOtjx81K9fUVxBfs/eNwM3\nG2MSjDFP28sq4DYsM1CFxI32vtLU1LPn352/z+djTnw8z27aRMTRo6zIymLvN99w0ejRZI0bFxa/\n/sKiv58zVJNz3KrLCUEHABHZJSJnAg8B2+3lIRHpLiK7Skae4gZ8PnjvPetvfng8HjKTkzkNmAAc\nBj4BKhnDAz/+qCYfRXEZhY4FVJLoHEDps2IF3H03VK8Oy5dDo0bByyYlJfFDr14sO3KEd7CmgqOB\nKyIjGbZ2LV27di0h1YpSsQlnLCClApA7Zs/gwQXfwMcaw6M+H9uBGkBLYCfwXdOmdNEonoriOtT/\nrpC40d4Xak2ff27F7PnXv4LE7MmNHcMn4sILeei++6japAk3REZyRWQkPWNi6DBwYJ5vAZcGFeH3\nCwWqyTlu1eUEfQJQTqBPHytmT2A2rqB4vVbkzhYtwOulU7NmfHH//Xg8HhITE/lozBjXdP6KouRE\n5wCUIpGZns6aG2/kvM8+U79+RXEZOgeg5IvPB2+8YUXlHDascMd6ly5l9OjRtK5Th3M2biSiRYvw\niFQUJazos3khcaO9r7CaVqyAuDiYORNiYpwfl52b98orGT9yJO/s3p1v518ezlVJoJqc4UZN4F5d\nTtAngApEUTx7/Gx5912GDR+uuXkVpRyhcwAVBBHLq6d/fxgz5sSELEGxc/PufuEF/vfvfzNi1iwN\n3qYoLsfpHIAOAEpwAj185s4t9TAOiqI4o9gJYZS8cZO9L6/cu2Dd7RcLf27eCy6A8eOLHMPHTefK\nj2pyhmpyjlt1OUEHgDJKXrl3k5KSee01iI2FtLSC68gL79Kl3NCsGT7Nzaso5R41AZVBguXerVZt\nHGecMYOnnoqgV6/C1am5eRWl/KDvAZRjguXe9fn68sILHrp2jStUfX6/fvXwUZSKhd7iFRJ32vsS\nAKhcuZDWmsxM1sXHZ/v1hzo3rxvPlWpyhmpyjlt1OSGsA4Axpqox5itjjMcY860xZoq9fYox5ldj\nzEZ7uTCcOsobsbGxtG+fQL65dwvC64Xu3Tlz716+2bhRc/MqSgUk7HMAxpgaIpJhjKkErMXKKHYR\nkCYi/yngWJ0DsDlyBA4fhnr1rHWPJ5n4+Dls29YXgHbtEnjllRuIje2Uf0W2X7/m5lWU8otr3EBF\nJMP+WBVrzsHfo2uv4wCfD157DTp0sGL3+ImN7URS0gxWr27F6tWt2Ljx2QI7/7/Wr4fu3SEpST18\nFEUJ/wBgjIkwxniA34D/E5Gv7V23GGO8xpiXjTF1wq0jVJSkvW/FCujSxYrZ88YbcMstOfdHREQQ\nFxdHWlpaviGX/TF8uvfuXaK5ed1oG1VNzlBNznGrLieE3QtIRHxArDEmCnjXGNMRmImVX1iMMY8A\n/wGuzev40aNH06pVKwDq1q1L586d6devH/D3iS/Jda/XG/b2evfux4AB8M03CVx/PUye3A9jgpf3\nk9f+Hz//nBeefJIWUVE8OHMma1q1op991x/u8+X1esNav1t/v8Ku+3GLHreuu/F6CqQ09SQkJLBg\nwQKA7P7SCSX6HoAxZjLwV6Dt3xgTA3wgIqfnUb7CzgGsWgVnn12ImD25UL9+Ram4uOI9AGNMQ+CY\niBw0xlQH/gk8ZoxpIiK/2cUGA5vDqaMs0rdvwWV8Ph8ejwewPIOyzUBeL99ecQXeQ4fUr19RlKCE\n+5awKbDSGOMFvgI+E5GPgSeMMd/Y2/sCd4RZR8jI/dhXHI4cgffeK9qxyR4P4+Li2N67N+O6d+f2\nLl1I3rAhO4ZP3P33836I/foLSyjPVahQTc5QTc5xqy4nhPUJQES+BU7ogURkVDjbdTv+bFz3329N\n8l58MURGFuZ4H3Pi45nh9fIOkATcvmkTc/r2Zca55xLh9WrkTkVRCkRjAZUwK1bA3XdD9erw5JMU\nOmYPQFJSEi/16MGaY8fIBH4E2gGHgVH33MO0Rx8NrWhFUcoUrpgDUHLy0ktWp1/YbFx58c/KlWl7\n7Bj32etHgKGRkQwaPDgUUhVFqQCoW0ghKY6976qrIDkZhgwpXuffqV07Ho+I4GGgEhADpALfNW1K\nFya9lAoAAA79SURBVBdN+LrRNqqanKGanONWXU7QJ4ASpEaN4tfx4/Ll/HvoUOpUr067atU4Ny2N\nWsCXTZvSYeDAfF8IUxRFCUTnAELMkSPw/PPW5O5554WwYjuGz4EXX+STyy/nyhdfRCBvN1BFUSo0\nOgdQwgR69sTGwmWXhbDygNy89b75huG2h48B4uIKF/tfURTFj94yOiRY/l04MWbPe+9B+/YhaNRh\nbl632iDdqEs1OUM1OcetupygA4AD8sq/6/EkA5bJZ/p0mDwZ1q0rmltnXniXLmXESSdx7OuvNXKn\noihhQecACiBY/t3OnceRlDQj5HZ3jeGjKEpxcU0+gLJOsPy727b1zZ6ADRXepUvp3qgRSVu24E1M\n1CxdiqKEFe1dCk1C6KvMzMRz/fVFzs3rVhukG3WpJmeoJue4VZcTdADIB58PNm+O5dixBIqVfzc/\n7Ny8nXftInnTJr3rVxSlxNA5gCCsWQO33QbVqsGYMck8/3wR8u/mh+bmVRQlTDidA9ABIAgrVkBa\n2t8xe4LG3i8CB9esoc6tt0KLFjB3rkbuVBQlpOgkcDG54IKcMXuc5t/ND39u3i79+pF5220hy83r\nVhukG3WpJmeoJue4VZcTKvwAsH8/ZGSEvx3PkiV0sz181nz9NZHx8WryURSlVKmwJiB/zJ4nnoCF\nC62kLOEg0K//6TFjGDFrlk7yKooSVjQWUBByZ+NaswZOOSVMjXm9/DRsGJsPHtTcvIqiuI4KdSt6\n8CDExf0ds+fddwvf+Tuy9wXE8Dn13ntZtnt3WDt/t9og3ahLNTlDNTnHrbqcUKGeAOrUgRkzoE+f\nMJrfAyJ3orl5FUVxMRV2DiDUZKan8/4113D5qlXq168oSqlSod1A9++HDz8sufb8Hj6vrl7N0a++\n0sidiqKUCcrVAHDkiJV0vUMH+L//C23deeUD8Pv19x8+nDtHjWL5nj1Ubd06tA07wK02SDfqUk3O\nUE3OcasuJ5SLOYBwe/YkezzMiY+n79atvH70KJtnzeLS+HjunjiRllFR6uGjKEqZpFzMAUycCKtX\nW3f/oUrI4sfn8zEuLo4ZXi/vAPHAfOB/ERH0uf56hr34ovr1K4riKirUewCTJ0PNmuExu3s8HjKT\nkzkNOAakAfcCh0WoX7++dv6KopRZykXvVatWeOdc/1m5MlOBI/b6EWBolSoMGjQofI0WArfaIN2o\nSzU5QzU5x626nOD6AcA/4bp/P9x5J3z/fQkL2LKFu48eJRNIBWLsv981bUoXtfsrilKGcf0cwOmn\n38q5547l9dc7MWQIPPggnHRS+NsOjOFzx6BBLF+3jl7793NmpUrMa9yYUwYO5Jnnngu/EEVRlEJS\nbvIBQBZ16oxj3boZdOxYMg8snsWLGR0fT8uoKOZ88gnNunQJaT4ARVGUcFKOXgSL4Nixvhw+HNoE\n7HmSmcnWm2+m/1VXZfv1+907Q5EPIFy41QbpRl2qyRmqyTlu1eWEcuEFFBI8Hhg9mg7R0WzdvJl6\nHTuWtiJFUZSwUiZMQJ07jyMpaUZ47rw1N6+iKOWMcvMewBln3M78+TeEpfP/c+VKGowbB9HRGrlT\nUZQKh7uM2XmwceOzxMZ2Cmmd/hg+seefT8YttxQqN68b7X1u1ATu1KWanKGanONWXU5w/QAQ6jt/\nz+LFdGvUiI1btvDl119TY8wYNfkoilIhcf0cQKj0aW5eRVEqCuVmDiAkeDzsGT6cLZqbV1EUJZvy\nfQvsz83bvz8x99zD0l27it35u9He50ZN4E5dqskZqsk5btXlhPL7BGD79auHj6IoSt6UuzmAzPR0\nFo8Ywai1azFPP61+/YqiVDjKUSgI5/g9fN5ev56M9es1N6+iKEo+lIsBIDs3b0AMn5onnxyWttxo\n73OjJnCnLtXkDNXkHLfqckKZnwP49dNPuWTQIM3NqyiKUkjK7hyAHcMnc+ZMPrzySgbNmKF+/Yqi\nKJT39wACPHwiN21isHr4KIqiFJqw3jIbY6oaY74yxniMMd8aY6bY2+sZY1YYY7YaYz4zxtRxVGGA\nXz8TJhQqhk+ocKO9z42awJ26VJMzVJNz3KrLCWEdAETkKHCOiMQCnYGLjDHdgUnAf0WkA/A5cE+w\nOvw5gT2LF3NRw4Yc+vJLy6+/lDx8vF5vibdZEG7UBO7UpZqcoZqc41ZdTgi7CUhEMuyPVe32BLgM\n6GtvXwgkYA0KJ3DL6acjVauyzOPh6TFjqD1rFpSirT81NbXU2g6GGzWBO3WpJmeoJue4VZcTwj4A\nGGMigCSgLfCiiHxtjDlJRPYCiMhvxpjGwY7/ODkZU7kyGzdsoEXXruGWqyiKUmEI+620iPhsE1AL\noLsxphPWU0COYsGOPwIcz8pi1rvvhlGlc7Zv317aEk7AjZrAnbpUkzNU0/+3d/axWpZ1HP98BRsi\nQ60sUt5iqUEzjZcI0UWDGuak0cysXI71oit30FqxuZXNVolbErOtctopEIgkpFgqUJ5ZgIPDq6I4\nK0xeLKYZ5stCoG9/XNcDt6dzDgcmXDc9v8929jzPfe77ur/nPs99/e7rd13X9+o5ddXVE47rMFBJ\n3wBeBT4PTLC9W9IAoM328E72r+8Y1SAIghrTk2GgxzQASHorsM/2i5JOAZYBt5Ly/y/YnilpBnCG\n7U77AIIgCIJjw7EOAOeTOnlPyj8LbX9H0puBXwKDgGeAK22fuD0pQRAEJyC1ngkcBEEQHDtq6Z0g\nabKkJyU9lVNExZF0t6Tdkh4traWBpIGSHpL0eJ5o11IDTZ1O/qsDkk6StEHSb0praSDpr5I25+u1\ntrQeAEmnSbpX0tb83RpbWM+5+fpsyK8v1uS7fqOkLZIelTRP0ptqoGl6vu96VB/UrgWQh40+BUwE\nngXagatsP1lY18XAy8Ac2+8tqaVB7kAfYHuTpH6k4bYfq8G16mv7VUm9gFVAi+3ilZukG4FRQH/b\nU0rrAZC0DRhl+5+ltTSQ9DPgYdutknoDfW3/q7As4GD9sBMYa3tHQR1nASuBd9t+TdJC4Le25xTU\n9B5gATAG2A88AFxne1tXx9SxBfB+4E+2n7G9D/gFaeJYUWyvBGpzk0KaQ2F7U37/MrAVOLusqi4n\n/xVF0kDgo8BdpbV0QNToPpTUH7jEdiuA7f11qfwzk4C/lKz8K/QCTm0ESdIDa0mGA2ts77V9APgD\n8PHuDqjNF6/C2UD1n7uTGlRqdUfSUJLdxpqySg6mWjYCfwdW2G4vrQmYBXyNGgSjDhhYIald0hdK\niwHeCTwvqTWnXO7MI/jqwidJT7lFsf0s8H1gO7AL2GP7d2VVsQW4JHut9SU98Azq7oA6BoDgCMnp\nn0XA9NwSKEqHyX9jJY0oqUfSZcDu3FpS/qkL422PJN2sX86pxpL0BkaSZu2PJM3bqcUQbUknA1OA\ne2ug5XRSZmIIcBbQT9KnS2rKqd+ZwArgfmAjcKC7Y+oYAHYBgyufB+ZtQSfk5uciYK7tX5fWUyWn\nDtqAyYWljAem5Hz7AuBDkorlaqvY/lt+fQ64j5QCLclOYIftdfnzIlJAqAOXAuvztSrNJGCb7Rdy\numUxcFFhTdhutT3a9gRgD6k/tUvqGADagXdJGpJ71a8C6jJqo25PjwA/BZ6wPbu0EEiT/xr23jl1\n8GGgaKe07ZtsD7Y9jPR9esj2Z0tqgtRZnltvSDoV+AipGV+M7NG1Q9K5edNE4ImCkqp8ihqkfzLb\ngQ9I6iNJpOu0tbAmJJ2ZXwcDU4H53e1fuwVhbB+QdD2wnBSg7rZdhws7H5gAvEXSduDmRkdZQU3j\ngc8Aj+Wcu4GbbD9YUNY7gJ/n0RqNyX/3F9RTZ94O3JctT3oD82wvL6wJoAWYl1Mu24BphfWQc9qT\ngC+W1gJge62kRaQ0y778emdZVQD8Kk+03Qd86XAd+LUbBhoEQRAcH+qYAgqCIAiOAxEAgiAImpQI\nAEEQBE1KBIAgCIImJQJAEARBkxIBIAiCoEmJABDUDkkHsg/Nlmz/+5XK70ZJ+sExPPcHJY2rfL5W\n0tVHUc4Fki6tfL5c0tffII1P92CftjwZKAi6pHYTwYIAeCX70DSWFV0gqb/tb9leT7K9Pmok9crT\n9ztjAsn2+xEA2z85ytNcCIwmWfJieymw9CjL6khM3gneEKIFENQa28+TZn9eDwef0JdW3jcWClmf\n7RSQNCMv0rFR0nfztjZJsyS1Ay3ZsmKR0uI1aySNkzQEuA64IZc5XtLNjRaIpJa8QMqmPDMcSWMk\nrc7nXynpnDyD9hbgylzOJyRdI+mOfMwQSb/P5azIVtVkB87ZklZJ+rOkrqx8n6v8/W06tHjL3Mo+\n/+AwRmBBEC2AoPbYfjpbTJ/Z2JRfv0qa7v5Itgr4t6TJwOXAGNt7s2tjg5NtjwGQNA+43fZqSYOA\nZbZHSPox8JLt2/N+kyrHzwCG2t6n5JsPyf/lYtv/kTQR+J7tKyR9k7TQS0su55qK7juAVtv3SJqW\nP0/Nvxtge7yk4SQPrMWdXI/qCl0XAiNI1turJF1ke7XtK3p2dYNmJgJAcKLQmQnfKmBWrswX296V\nK+xW23sBbO+p7L+w8n4SMDwbeUGy8+17GA2bgfmSlgBL8rbTgTmSziFV8D25p8ZxqMKfS7LwbbAk\n694q6W09KGttw1FU0iZgKLC6B8cFQaSAgvojaRiwv6MNsO2ZwOeAU4CVks47TFGvVIslLSv4vvwz\nuLKSWVdcBvyQZI/cng3vvk1yFz2f1PLo04M/qbsc/t4OGg9Hdf8DxENdcAREAAjqyMGKL6d9fkRK\nk7x+J2mY7cdt3wasA84jLYYxLVtRI+mMLs6xHJheKeuC/PYloH+nR8Bg2w+TFkjpD/QDTuPQehVV\n18zuyllNsjYGuBr4Yxf71c16PPg/IwJAUEf6NIaBkirqB23f0sl+N0h6LKc+XgMesL2MlDtfJ2kD\nqZ8A/vepezowWtLmfJ5r8/alwNRGJ3DjOKWFd+6RtJk0Cml2ttq9DbhV0npefz+1ASMancAdzt1C\nClKbSHbejUDUUeORjvaJ0UHBERF20EEQBE1KtACCIAialAgAQRAETUoEgCAIgiYlAkAQBEGTEgEg\nCIKgSYkAEARB0KREAAiCIGhSIgAEQRA0Kf8FFGHIEDbI6eQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xb58df98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# YOUR CODE HERE\n",
    "n = 10\n",
    "dx = n\n",
    "b0 = 30.\n",
    "b1 = 50.\n",
    "max_iter = 10000\n",
    "epsilon_tol = 1e-6\n",
    "\n",
    "#Direct Solution\n",
    "direct = plt.plot(inv_tridiag_121_b(n,b0,b1),'bo--', label='Direct Soln.')\n",
    "#plt.xlim([-1,10])\n",
    "\n",
    "#Jacobi Method\n",
    "u_jac,iter_jac = solve_laplace_jac(b0,b1,dx,max_iter,epsilon_tol)\n",
    "jacobi = plt.plot(u_jac,'ro-', label='Jacobi Soln.')\n",
    "\n",
    "\n",
    "#Gauss-Seidel\n",
    "u_gs,iter_gs = solve_laplace_gs(b0,b1,dx,max_iter,epsilon_tol)\n",
    "gauss_seidel = plt.plot(u_jac,'k*--', label='Gauss-Seidel Soln.')\n",
    "\n",
    "plt.title('Comparing solutions from different solvers')\n",
    "plt.legend(loc = 'upper left')\n",
    "plt.xlabel(\"Discretisation 'n'\")\n",
    "plt.ylabel(\"Output values\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print in the following code cell the number of iterations required for conergence by the two solvers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of iteration required for convergence (Jacobi): 253\n",
      "The number of iteration required for convergence (Gauss-Seidel): 134\n"
     ]
    }
   ],
   "source": [
    "# YOUR CODE HERE\n",
    "print('The number of iteration required for convergence (Jacobi): %d'%iter_jac)\n",
    "print('The number of iteration required for convergence (Gauss-Seidel): %d' %iter_gs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Describe and interpret the results, and explain why Gauss-Seidel converges with less iterations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Despite having produced the same results (as seen in the graph above), the solutions reach convergence almost half the time quicker with Gauss-Seidel method (134 iterations) than with Jacobi method (253 iterations). \n",
    "\n",
    "The difference between the Jacobi and Gauss-Seidel rule appears in the moment of updating. The calculation of a new nodal value depends on the values of the surrounding nodes.\n",
    "\n",
    "Jacobi method uses the value of the previous iteration cycle, whereas Gauss-Seidel mehtos uses the newly found value. \n",
    "\n",
    "Thus, in Jacobi iteration, the updating occurs when all nodes have been recalculated (have to store a whole set of values/data), while in Gauss-Seidel iteration each node is updated directly when its new value is found (only storing a few relevant value in the process). Therefore, Gausess-Seidel method converges with a lot less iterations compared to Jacobi method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}