Ex10-Diffusion-equation-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 10: Diffusion problems\n",
    "\n",
    "Up to now, we were concerned with steady-state problems and solutions to the Laplace and Poisson equation. However, in subsurface reservoir problems, we are a lot more often concerned with the evaluation of properties that change with time! \n",
    "\n",
    "From now on, we will investigate several of these transient problems and investigate how (a) they relate to the problems we considered previously, and (b) present some new numerical challenges.\n",
    "\n",
    "Today, we will look at problems that can be described with a diffusion-type equation: \n",
    "\n",
    "$$\\frac{\\partial u}{\\partial t} = \\kappa \\frac{\\partial^2 u}{\\partial x^2}$$\n",
    "\n",
    "We will talk in the lecture about how diffusive systems generally behave over time. Typical examples of diffusive problems in the context of reservoir engineering are, for example, conductive heat transport:\n",
    "\n",
    "$$\\frac{\\partial T}{\\partial t} = \\kappa \\frac{\\partial^2 T}{\\partial x^2}$$\n",
    "\n",
    "where $\\kappa$ is thermal diffusivity, or groundwater flow in a confined aquifer:\n",
    "\n",
    "$$\\frac{\\partial h}{\\partial t} = \\frac{T}{S} \\frac{\\partial^2 h}{\\partial x^2}$$\n",
    "\n",
    "where $T$ is transmissivity, and $S$ the (depth-averaged) storativity.\n",
    "\n",
    "If we want to solve these equations numerically, we now have to find discretisations in time and space. We will look at several options to do that in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#from IPython.core.display import HTML\n",
    "#css_file = 'nre2_style.css'\n",
    "#HTML(open(css_file, \"r\").read())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "plt.rcParams['figure.figsize'] = [10., 7.]\n",
    "plt.rcParams['font.size'] = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explicit scheme to solve the diffusion equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We discussed in the lecture different schemes to solve the diffusion equation \n",
    "\n",
    "$$\\frac{\\partial u}{\\partial t} = \\kappa \\frac{\\partial^2 u}{\\partial x^2}$$\n",
    "\n",
    "with the finite difference method. The simplest - and most intuitive - way is to use a first-order Taylor expansion for the time derivative:\n",
    "\n",
    "$$\\frac{u_i^{k+1} - u_i^k}{\\Delta t} = \\kappa \\frac{u_{i-1}^k - 2 u_i^k + u_{i+1}^k}{\\Delta x^2} $$\n",
    "\n",
    "If we now separate the variables for a new time step $k+1$ and the current time step $k$, we obtain the simple explicit formulation for the 1-D diffusion equation:\n",
    "\n",
    "$$u_i^{k+1} = \\frac{\\kappa \\Delta t}{\\Delta x^2} \\left( u_{i-1}^k - 2 u_i^k + u_{i+1}^k \\right) + u_i^k$$\n",
    "\n",
    "For simplification, we substitute $\\frac{\\kappa \\Delta t}{\\Delta x^2} = \\beta$ and obtain:\n",
    "\n",
    "$$u_i^{k+1} = \\beta u_{i-1}^k + (1 - 2 \\beta) u_i^k + \\beta u_{i+1}^k $$\n",
    "\n",
    "As this equation provides a function for an update of one node value from all the known node values at the previous time step, we can define a simple python function to solve the diffusion problem.\n",
    "\n",
    "We now need an initial distribution of our property $u$ in the domain. Let's use a simple function for an initial distribution of our property:\n",
    "\n",
    "$$u(x, t=0) = u_{max} \\exp \\left[ - \\left( \\frac{x-0.5}{\\sigma} \\right)^2 \\right]$$\n",
    "\n",
    "Let's have a look at a diffusion problem in a domain with the following properties:\n",
    "\n",
    "Variable | Value  | Explanation\n",
    "------------- | -------------\n",
    "$L$  | 1. | Domain size\n",
    "$\\kappa$ | $10^{-6}$ | Diffusivity\n",
    "$dt$ | 25. | Time step\n",
    "$nt$ | 200 | Number of time steps\n",
    "$nx$ | 100 | Number of cells\n",
    "$u_{max}$ | 10. | Maximum value of initial perturbance\n",
    "$\\sigma$ | 0.05 | Half-width of initial perturbance\n",
    "\n",
    "Write a code to solve this problem and create a plot with the initial and final distribution of the property:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0xa78dc50>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEZCAYAAACXRVJOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGm5JREFUeJzt3XuUnHWd5/H3J+kkJCGBdMiFJBAJImAwaABhd0R60V0u\nXgBnmUEYUMeds+vsMOquLMw6u4Q5jOOcs+INnLOMrq4zw+IBV0FHRpSxnYVRGQkJAQIMAUMupDuE\nhBDCJZfv/vFUk0rbl6ruqvo9l8/rnDqp6n76eb55Tld9+nd7HkUEZmZmE1IXYGZm+eBAMDMzwIFg\nZmY1DgQzMwMcCGZmVuNAMDMzwIFgBSPpYUnvbMW2kn4g6fIG9/W0pLMbrbNRkvZLWtLq/ZqNRVfq\nAsyaEREnjWVbSdcCx0bEFXXfP7/F5Y2FFwJZbriFYJaWUhdgNsCBYIVS33Uj6VpJ35L0vyXtlLRG\n0vLB20o6B/ivwG9LelHSg7Xv/0TS79aeL5F0j6TnJPVL+mtJMxuo5+2SnpWkuq9dJGl17flpkv5R\n0nZJmyR9WdKQLfP6emqvPyTp/9W9PkHS3ZK2SVor6eJmz5/ZSBwIVnTvA24BDgO+B9w0eIOI+CHw\nGeBbETEjIt42xH5U22Y+cCKwCFgx2sEj4n5gF1A/vvBB4K9rz/cBnwC6gX9R2+73G/h/vX4IAEnT\ngLtr+z0CuAS4SdIJTezLbEQOBCu6eyPih5FdlOuvgGVj2UlErIuIeyJib0RsAz4PnNXgj98KXAog\naQZwfu1rRMTKiLg/Ms8ANzex33rvBZ6OiG/W9rUa+L+AWwnWMh5UtqLbUvd8N3CIpAkRsb+ZnUia\nC3wROBM4FJgIPN/gj98C3CfpPwAfAB6IiA21/R4H3ACcCkwle8890ExtNYuBMyQN1KRajX81hn2Z\nDcktBKuK0WbzfAbYDyyNiMOB36HBAd+IWAusJ2sZfJAsIAb8BbCWbIbT4cCnR9jvS8C0utfz655v\nAHojorv2mBURMyPiPzZSo1kjHAhWNsN92PYBb6gf/B1kBtlYwIuSFgJXNXncW4CPk7Uwbhu0350R\nsbvW3/+xEfaxCviApKmS3gh8tO573wfeJOl3JHVJmiTpVI8hWCs5EKxoRvtLP4Z5fhtZWGyT9Msh\nvn8dcAqwg2xw+ttNHvdW4J3APRFR39X0KeAySTuB/1nbbrj9fh7YQ9YN9nUODEwTEbuAf0M2mLy5\n9vgsMHmUuswapnbeIEfS18gGw/oiYlnta7OAb5H1if4K+K2IeKFtRZiZWUPa3UL4OnDOoK9dA/w4\nIo4H/h74ozbXYGZmDWhrCwFA0mLge3UthMeAsyKiT9J8soEy94OamSWWYgxhbkT0AUTEFmBughrM\nzGyQPAwq++JeZmY5kGJhWp+keXVdRv3DbSjJYWFmNgYR0fSFEzvRQhAHzw2/E/hw7fmHgDtG+uGI\n8COCa6+9NnkNeXmkPBef+Uxw+eUHHnfdVd1zkbeHz8WBx1i1NRAk3QL8I9mCmmckfYRs7vS/lvQ4\n8K7aa7Pc27sXrr8eenrg3e+GWbPgq19NXZVZ67S1yygiLh3mW+9u53HN2uGxx2DRIvjd2gWqTz8d\nzjsvbU1mrZSHQWVrQE9PT+oSciPVuVi5EpYvP/D6uONg61bYvj1JOYB/L+r5XIxf29chjIekyHN9\nVi2f+AQsXAhX1V3l6Mwz4brr4OyW323ZbOwkETkdVDYrhcEtBMher1yZph6zVnMgmDVg/35YtQre\nNuheaw4EKxMHglkDnnwSZs+G7u6Dv+5AsDJxIJg1YKjuIoATT4QNG+DFFztfk1mrORDMGjBcIHR1\nwUknwerVna/JrNUcCGYNGC4QwN1GVh4OBLNRRDgQrBocCGajWL8epk6FefOG/r4DwcrCgWA2ipFa\nB5CNITz5JLz8cudqMmsHB4LZKEYLhClT4PjjYc2aztVk1g4OBLNRjBYIAKec4m4jKz4HgtkIIuCB\nB0YPBI8jWBk4EMxG8Oyz2WUrFi0aeTsHgpWBA8FsBI8+CkuXgka5buTSpdm2vjivFZkDwWwEGzfC\nUUeNvt2MGTB5ctp7I5iNlwPBbASNBgJk223c2N56zNrJgWA2go0bRx8/GLBokQPBis2BYDYCB4JV\niQPBbAQbNjQXCBs2tLces3ZyIJiNwGMIViUOBLNh7N4NL70ERxzR2PbuMrKicyCYDWPTJli4cPQ1\nCAMcCFZ0DgSzYTQzfgAHxhC8OM2KyoFgNoxmZhgBzJwJEybACy+0ryazdnIgmA2jmQHlAR5YtiJz\nIJgNo9kWAngcwYrNgWA2DAeCVY0DwWwYzQ4qgxenWbE5EMyG4TEEqxoHgtkQXnkFdu6EOXOa+zl3\nGVmRORDMhrBpEyxYkE0jbYYDwYrMgWA2hLGMH4DHEKzYHAhmQxjL+AHAYYdlK5V37mx9TWbt5kAw\nG8JYppxCdt0jdxtZUSULBEmflPSwpIck/Y2kyalqMRtsrIEADgQrriSBIGkBcCWwPCKWAV3AJSlq\nMRuKA8GqqCvhsScC0yXtB6YBmxPWYnaQsQ4qgweWrbiStBAiYjPwOeAZYBOwIyJ+nKIWs6GMdVAZ\nvDjNiitJC0HS4cAFwGLgBeB2SZdGxC2Dt12xYsXrz3t6eujp6elQlVZVr74K27fD3Llj+/lFi+A7\n32ltTWYj6e3tpbe3d9z7USS4m4ekfwucExG/V3t9OXB6RPzBoO0iRX1WbU8/DT09sH792H7+oYfg\nsstgzZqWlmXWMElERIP3+jsg1SyjZ4AzJB0iScC7gLWJajE7yHgGlMGDylZcqcYQ7gduBx4EVgMC\nbk5Ri9lgGzaMffwAYNYseO012LWrdTWZdUKydQgRcV1EnBgRyyLiQxGxJ1UtZvXG20Lw4jQrKq9U\nNhtk40ZYuHB8+/DUUysiB4LZIM8+m13pdDyOPBK2bGlNPWad4kAwG6SvD+bNG98+5s3L9mNWJA4E\ns0EcCFZVDgSzQRwIVlUOBLM6r7wCu3dnU0fHw4FgReRAMKvT359dskJNr/E8mAPBisiBYFanFd1F\n4ECwYnIgmNVpVSDMmQPbtsH+/ePfl1mnOBDM6rQqECZNgpkzs1AwKwoHglmdLVtaEwiQ7ceL06xI\nHAhmdVrVQgCPI1jxOBDM6jgQrMocCGZ1HAhWZQ4Eszp9fTB/fmv25UCwonEgmNVxC8GqzIFgVvPa\na/Dii9Dd3Zr9ORCsaBwIZjX9/dmCsgktelc4EKxoHAhmNa3sLgIHghWPA8GsptWBMHcubN3qy1dY\ncTgQzGpaHQhTpsD06bB9e+v2adZODgSzmlYHArjbyIrFgWBW40CwqnMgmNU4EKzqHAhmNQ4EqzoH\ngllNOwJh/nwHghWHA8Gsxi0EqzoHghmwdy/s2AFHHNHa/ToQrEgcCGZkC8i6u2HixNbu14FgReJA\nMKM93UXgQLBicSCY0dp7KdcbCISI1u/brNUcCGa09sY49Q45JHvs2NH6fZu1mgPBjPZ1GYG7jaw4\nHAhmOBDMwIFgBjgQzCBhIEg6TNJtktZKekTS6alqMXMgmEFXwmN/EfhBRFwsqQuYlrAWqzgHglmi\nFoKkmcCZEfF1gIjYGxE7U9RiBg4EM0jXZXQM8Jykr0taKelmSVMT1WIVt28fPP88zJnTnv07EKwo\nUgVCF7AcuCkilgO7gWsS1WIV99xzcPjh0NWmDlQHghVFqjGEjcCGiPhl7fXtwNVDbbhixYrXn/f0\n9NDT09Pu2qxi2tldBA4Ea7/e3l56e3vHvR9FojX1kn4K/F5EPCHpWmBaRFw9aJtIVZ9Vx49+BJ/9\nLNxzT3v2/9JL2VVUd+8GqT3HMKsniYho+rct5SyjPwT+RtIk4CngIwlrsQprdwth+vTsKqovvggz\nZ7bvOGbjlSwQImI1cFqq45sNaHcgwIFuIweC5ZlXKlvldTIQzPLMgWCV50AwyzgQrPIcCGYZB4JV\nngPBLONAsMpzIJhlHAhWafv3w9at7btsxYC5cx0Iln8OBKu055+HQw+FKVPaexy3EKwIHAhWaZ3o\nLgIHghWDA8Eqrb+/c4HQ39/+45iNhwPBKq1TLYSZM2HPnux6RmZ55UCwSutUIEjuNrL8cyBYpXUq\nEMCBYPnnQLBKcyCYHeBAsEpzIJgd4ECwSnMgmB3gQLBKcyCYHeBAsMqK6Nw6BHAgWP41dMc0Sf99\nqK9HxJ+0thyzznnhBZg8GaZO7czxHAiWd43eQvOluueHAO8F1ra+HLPO6WR3ETgQLP8aCoSI+Fz9\na0n/A/hhWyoy6xAHgtnBxjqGMA1Y1MpCzDqt04Ewa1Z26YpXXuncMc2a0egYwhogai8nAnMAjx9Y\noXU6EKTsvgj9/XD00Z07rlmjGh1DeG/d871AX0TsbUM9Zh3T6UCAA91GDgTLo0bHENa3uxCzTuvr\ng+XLO3tMjyNYnnkdglVWyhaCWR45EKyyHAhmB3MgWGU5EMwO5kCwSopwIJgN5kCwStq1K5sGeuih\nnT2uA8HyzIFglZSidQAOBMs3B4JVkgPB7Nc5EKySUgXC7Nmwcyfs2dP5Y5uNxoFglbRlS5pAmDAB\n5sxxK8HyyYFglbR5MyxcmObYCxZkxzfLGweCVdKmTekCYeHC7PhmeeNAsEravDn7Sz0FtxAsr5IG\ngqQJklZKujNlHVY9biGY/brULYSPA48mrsEqyC0Es1+XLBAkLQLOB76aqgarppdfzu5cNnt2muMv\nXOhAsHxK2UL4PHAVB+7EZtYRmzfDkUdml65IYcECdxlZPiUJBEnvIbvr2ipAtYdZR6TsLgJ3GVl+\nNXoLzVb7DeD9ks4HpgIzJH0zIq4YvOGKFStef97T00NPT0+narSSSjmgDNDdfaDbatq0dHVYefT2\n9tLb2zvu/SgibY+NpLOA/xwR7x/ie5G6PiufG26AZ56BL3whXQ1LlsDdd8Mb35iuBisvSURE0z0v\nqWcZmXVc6hYCeOqp5VPyQIiInw7VOjBrl9RjCOBxBMun5IFg1mluIZgNzYFgleMWgtnQHAhWKRH5\nCAQvTrM8ciBYpezYAZMmdf5eyoN5cZrlkQPBKiUPrQNwl5HlkwPBKiUPA8pwIBC8zMbyxIFglZKX\nFsL06TBlCmzfnroSswMcCFYpeWkhgKeeWv44EKxS8tJCAI8jWP44EKxSNm/OVwvBgWB54kCwStm0\nKV8tBHcZWZ44EKxS3EIwG54DwSpj3z7o74f581NXknELwfLGgWCV0d+f3Zxm0qTUlWTcQrC8cSBY\nZeRp/ADcQrD8cSBYZeRpyinAvHmwdSvs3Zu6ErOMA8EqI0+L0iDrujriCOjrS12JWcaBYJWRtxYC\neHGa5YsDwSojT1NOB3hg2fLEgWCVkbdBZfDAsuWLA8Eq41e/gsWLU1dxsMWLs7rM8sCBYJWwf3/2\nwXvMMakrOdixx8JTT6WuwizjQLBK2LQpW5Q2bVrqSg62ZAmsW5e6CrOMA8Eq4amnsr/G82agheA7\np1keOBCsEtaty/4az5tZs2DCBNi2LXUlZg4Eq4innspnIIDHESw/HAhWCXntMoIsqBwIlgcOBKuE\nvHYZgQeWLT8cCFYJeW4huMvI8sKBYKW3cyfs3g1z56auZGhuIVheOBCs9AYGlKXUlQzNLQTLCweC\nld66dfntLgJYtCi7m9urr6auxKrOgWCll+cppwBdXXDUUb6mkaXnQLDSy/OA8gBPPbU8cCBY6eV5\nyumAY4/1wLKl50Cw0st7lxG4hWD5kCQQJC2S9PeSHpG0RtIfpqjDym/vXtiwAd7whtSVjMwtBMuD\nrkTH3Qv8p4hYJelQ4AFJd0fEY4nqsZLasAHmzYMpU1JXMjK3ECwPkrQQImJLRKyqPd8FrAVydrdb\nK4O8TzkdMBAIvgy2pZR8DEHSG4C3Ar9IW4mVURHGDwBmzIDp06GvL3UlVmWpuowAqHUX3Q58vNZS\n+DUrVqx4/XlPTw89PT0dqc3KoSiBAAdaCfPnp67Eiqa3t5fe3t5x70eRqI0qqQv4PnBXRHxxmG0i\nVX1WDhdfDL/5m3DJJakrGd1ll8G558Lll6euxIpOEhHR9MVaUnYZ/S/g0eHCwKwVithCMEsl1bTT\n3wAuA86W9KCklZLOTVGLlVdEcQaVwVNPLb0kYwgRcR8wMcWxrTq2b89Cobs7dSWNWbIE/vIvU1dh\nVZZ8lpFZuzzyCJxwQn4vez3YiSfCo4966qml40Cw0lq5Ek45JXUVjZszBw491Fc9tXQcCFZaK1fC\n8uWpq2jO8uVZ3WYpOBCstBwIZs1xIFgp7d6dzdhZujR1Jc1xIFhKDgQrpTVrsgHlvF/UbrDly+GB\nBzywbGk4EKyUithdBLBgQTYravPm1JVYFTkQrJSKGgiSu40sHQeClVJRAwEcCJaOA8FK57XXYO1a\nWLYsdSVj40CwVBwIVjqPPJJdBmLatNSVjI0DwVJxIFjpFLm7CLL7P+/aBf39qSuxqnEgWOkUPRAG\nBpYffDB1JVY1DgQrnaIHArjbyNJwIFip7N0LDz0Eb31r6krGx4FgKTgQrFQefxwWLoSZM1NXMj5v\ne5sDwTrPgWClUobuIoDjjssGlbdvT12JVYkDwUrl3nvh7W9PXcX4TZwIp54K992XuhKrEgeClcb+\n/XDnnfC+96WupDXe/364447UVViVOBCsNH7xC5g9O+tuKYMLLsgCbt++1JVYVTgQrDS++1248MLU\nVbTOkiUwfz787GepK7GqcCBYKUTAd74DF12UupLWuvDCLOjMOsGBYKWwdi28/HI5ZhjVu+iiLBB8\nwxzrBAeClcJAd5GUupLWOvnkbLHdww+nrsSqwIFgpfDd75avuwiygHO3kXWKA8EKb+NGWLcOzjwz\ndSXtceGF2fiIWbs5EKzw7rgD3vMemDQpdSXt8Y53wDPPwPr1qSuxsnMgWOGVcXZRva6ubLGdu42s\n3RwIVmgPPJDdIe2cc1JX0l4f/jB86Uvw6qupK7EycyBYYUXApz4F111X3NtlNuqss+DNb4abbkpd\niZWZIscTnCVFnuuztL73PbjmGli9OutWKbu1a+Gd78wu8d3dnboayzNJRETTk7AdCFZIe/bAW94C\nN9wA55+fuprO+djHYOrU7P9tNhwHglXKV76SDSbffXf5FqONpK8Pli7NLuR37LGpq7G8ciBYZezY\nASecAH/3d8W/VeZY/OmfwoMPwm23VSsMrXFjDQQPKluhbNuWzSi69NJqhgHAJz8JTzwBV13laxxZ\nayULBEnnSnpM0hOSrk5VhxXHpk3ZoGpPD3zuc6mrSWfaNOjtze6m9tGPZtc6MmuFJIEgaQJwI3AO\nsBT4oKQTUtRSFL29valLSOqJJ7IVu1dcAeed11v5rpLubvjxj2HNml4uvhh2705dUXpVf4+0QqoW\nwtuBf46I9RGxB7gVuCBRLYVQ1V/2++/PuofOOAM+/Wm4+urqnovBpk+Hc87p5bDD4Jhj4I//GDZv\nTl1VOv69GL9Us7cXAhvqXm8kCwmrqIjsfgZPPw2rVmWPn/4Utm6FK6/MZhUdfnjqKvOnqwu+8Y1s\nbcKXvwwnnZRd5O+007IxlmXL4Mgjy3udJ2ut3C/nKcsN08fr8cezyzTkWf0A58DziAOP/fuz+wPv\n25ddgmHgsWNHNlgswdFHZx9kJ5+czaY5+2yYODHN/6dIjj8ebrwRrr8e7rorW6x3442wZg3092dr\nF7q7YcYMmDIle0yenJ3bCROyfwe64Qb/O/h5XhXhPZJ3SaadSjoDWBER59ZeXwNERPz5oO08h8LM\nbAwKsw5B0kTgceBdwLPA/cAHI2Jtx4sxMzMgUZdRROyT9AfA3WQD219zGJiZpZXrlcpmZtY5uVip\n3MgiNUlfkvTPklZJKu0a1dHOhaRLJa2uPe6V9JYUdbZbowsXJZ0maY+kD3Syvk5q8P3RI+lBSQ9L\n+kmna+yUBt4fMyXdWfucWCPpwwnK7AhJX5PUJ+mhEbZp7nMzIpI+yELpSWAxMAlYBZwwaJvzgL+t\nPT8d+HnquhOeizOAw2rPzy3juWjkPNRtdw/wfeADqetO+DtxGPAIsLD2+ojUdSc8F38E/NnAeQC2\nAV2pa2/T+XgH8FbgoWG+3/TnZh5aCI0sUrsA+CZARPwCOEzSvM6W2RGjnouI+HlEvFB7+XOyNR1l\n0+jCxSuB24H+ThbXYY2ci0uBb0fEJoCIeK7DNXZKI+cigBm15zOAbRFRyot7RMS9wPYRNmn6czMP\ngTDUIrXBH3KDt9k0xDZl0Mi5qPfvgLvaWlEao54HSQuACyPiL4ACzJIfs0Z+J94EdEv6iaR/knR5\nx6rrrEbOxY3AmyVtBlYDH+9QbXnU9Odm7hem2dAk/SvgI2TNxir6AlDfh1zmUBhNF7AcOBuYDvxM\n0s8i4sm0ZSVxDvBgRJwt6VjgR5KWRcSu1IUVQR4CYRNwdN3rRbWvDd7mqFG2KYNGzgWSlgE3A+dG\nxEhNxqJq5DycCtwqSWR9xedJ2hMRd3aoxk5p5FxsBJ6LiFeAVyT9A3AyWX97mTRyLj4C/BlARKyT\n9DRwAvDLjlSYL01/buahy+ifgDdKWixpMnAJMPhNfSdwBby+ynlHRPR1tsyOGPVcSDoa+DZweUSs\nS1BjJ4x6HiJiSe1xDNk4wu+XMAygsffHHcA7JE2UNI1sALGM63oaORfrgXcD1PrL3wQ81dEqO0sM\n3zpu+nMzeQshhlmkJunfZ9+OmyPiB5LOl/Qk8BLZXwGl08i5AP4b0A18pfbX8Z6IKNWFARs8Dwf9\nSMeL7JAG3x+PSfoh8BCwD7g5Ih5NWHZbNPh7cT3wjbqpmP8lIp5PVHJbSboF6AFmS3oGuBaYzDg+\nN70wzczMgHx0GZmZWQ44EMzMDHAgmJlZjQPBzMwAB4KZmdU4EMzMDHAgmJlZjQPBzMwAB4JZUySd\nWrs50WRJ02s3pHlz6rrMWsErlc2aJOlPgKm1x4aI+PPEJZm1hAPBrEmSJpFdaO1l4F+G30RWEu4y\nMmveEcChZHfkOiRxLWYt4xaCWZMk3QH8H+AYYEFEXJm4JLOWSH75a7Miqd2e8rWIuFXSBOA+ST0R\n0Zu4NLNxcwvBzMwAjyGYmVmNA8HMzAAHgpmZ1TgQzMwMcCCYmVmNA8HMzAAHgpmZ1TgQzMwMgP8P\nbjFu0SVuZgwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x3d510f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# your code here:\n",
    "xvals = np.linspace(0,1,100)\n",
    "plt.plot(xvals, 10. * np.exp(-((xvals - 0.5)/0.05)**2))\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('u')\n",
    "plt.title(\"initial value\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explicit equation in matrix form\n",
    "\n",
    "The explicit solution of the equation is easily possible with a simple loop scheme that we have seen above. However, in order to show you the similarities to the types of problems we solved before, we will now have a look at another representation of the equation. This will also lay the groundwork for the subsequent, more complex, schemes.\n",
    "\n",
    "We start with the equation that we solved before:\n",
    "\n",
    "$$u_i^{k+1} = \\beta u_{i-1}^k + (1 - 2 \\beta) u_i^k + \\beta u_{i+1}^k $$\n",
    "\n",
    "If you look at the equation on the right hand side and think about the solutions to the Laplace and Poisson problems that we analysed before, you will see a similarity: we can write it in matrix form and obtain a tridiagonal matrix, again!\n",
    "\n",
    "We define: \n",
    "\n",
    "$$A = \\begin{bmatrix}\n",
    "1 - 2\\beta & \\beta & 0 & 0 & \\cdots & 0\\\\\n",
    "\\beta & 1 - 2\\beta  &  \\beta &  0 & \\cdots & 0\\\\\n",
    "0 &  \\beta & 1 - 2\\beta  &  \\beta & \\cdots & 0\\\\\n",
    "\\vdots & & &  \\ddots & & \\vdots \\\\\n",
    "0 & \\cdots & 0 &  \\beta & 1 - 2\\beta  &  \\beta \\\\\n",
    "0 & \\cdots & 0 & 0 &  \\beta & 1 - 2\\beta \n",
    "\\end{bmatrix}$$\n",
    "\n",
    "And writing the variables $u$ in a column vector,\n",
    "\n",
    "$$u^k = \\begin{bmatrix}\n",
    "u_1^k\\\\\n",
    "u_2^k\\\\\n",
    "u_3^k\\\\\n",
    "\\vdots\\\\\n",
    "u_{i-1}^k\\\\\n",
    "u_i^k\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "We can write the equation in matrix form:\n",
    "\n",
    "$$u^{k+1} = A u^k$$\n",
    "\n",
    "## Boundary conditions\n",
    "\n",
    "One thing that we neglected so far are boundary conditions. We will first have a look at Dirichlet boundary conditions. If you write down the equation for the first node, you can see that the implementation of boundary conditions will be very similar to the Laplace equations before:\n",
    "\n",
    "$$u_1^{k+1} = \\beta\\,u_{0}^k + (1 - 2 \\beta)\\,u_1^k + \\beta\\,u_{2}^k $$\n",
    "\n",
    "We now simply assign the boundary value $b_0$ for $u_0^k$ and, accordingly, for the last value $u_n$ and obtain a column vector for the boundary conditions as:\n",
    "\n",
    "\n",
    "$b = \\begin{bmatrix}\n",
    "b_0\\\\\n",
    "0\\\\\n",
    "0\\\\\n",
    "\\vdots\\\\\n",
    "0\\\\\n",
    "b_1\n",
    "\\end{bmatrix}$\n",
    "\n",
    "If we now add this vector to the previous matrix equation, we obtain the explicit equation for the diffusion equation with dirichlet boundary conditions as:\n",
    "\n",
    "$u^{k+1} = A u^k + \\beta b$\n",
    "\n",
    "**Your task**: write the code to solve the problem described above with a matrix method (with Dirichlet BC's $b_0 = 0$ and $b_1 = 0$).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.525\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8HOWd+PHPs321K+1q1btkSVZx7w0bYWxMSShxKCEJ\nLfVIclzKXXKXXICUS3Jpl3KX+3EJJARCQguB0Iwxpti49241S7J639Wuts7vj7EtV7At2fKK7/v1\n0ovd2ZlnnpmX+e6z3/nOM0rTNIQQQowNhtHugBBCiJEjQV0IIcYQCepCCDGGSFAXQogxRIK6EEKM\nIRLUhRBiDHnfoK6U+p1Sqk0pteO4ZclKqRVKqf1KqVeVUq4L200hhBBn42xG6o8Ay05a9g1gpaZp\nZcAq4F9HumNCCCHOnTqbm4+UUgXAC5qmTT7yfh9wuaZpbUqpTGC1pmnlF7arQggh3s/55tTTNU1r\nA9A0rRVIH7kuCSGEOF8jdaFU5hoQQohLgOk8t2tTSmUcl35pP9OKSikJ+EIIcR40TVPnus3ZjtTV\nkb+jngfuOvL6TuBv79Mx+dM07r///lHvw6XyJ+dCzoWci/f+O19nU9L4J2AtMF4p1aCUuhv4IbBU\nKbUfuPLIeyGEEKPsfdMvmqbdfoaPloxwX4QQQgyT3FF6EVVVVY12Fy4Zci6GyLkYIudi+M6qTn1Y\nO1BKu9D7EEKIsUYphXYBL5QKIYSIAxLUhRBiDJGgLoQQY4gEdSGEGEPO945SIeJX22YY7Bp67ymH\nxPzR648QI0hG6uKDJRaFQ6+AKQHMDoj4oeXd0e6VECNGRurig8XfBlY35Fymvx/sht2PgKaBOufq\nMSEuOTJSFx8svkZIzBt6b00GYhDqG7UuCTGSJKiLDxZfEzhzh94rpb/3NY1en4QYQRLUxQeLt+nE\nkTqAM09fLsQYIEFdfHCEfRAJgC3lxOWJuXpaRogxQIK6+ODwNoEzB9RJ/+wd2eBvh1h4dPolxAiS\noC4+OE6+SHqU0QL2VBhoufh9EmKESVAXHxzeky6SHs+ZK3l1MSZIUBcfDFpUH4mfKagn5kkFjBgT\nJKiLDwZ/G1hdYLKd/nOnXCwVY4MEdfHB8F6pF9BvQopFISg3IYn4JkFdfDCc6SLpUUodScHIaF3E\nNwnq4oPh/UbqIBdLxZggQV2MfUdvOrKnvvd6MlIXY4AEdTH2+dvAkXnqTUcnc2Tp68qD0kUck6Au\nxr5grz7d7vsxWsFghsjAhe+TEBeIBHUx9p1tUAd9vWDvhe2PEBeQBHUx9klQFx8gEtTF2BfsO8eg\nLrXqIn5JUBdjn4zUxQeIBHUxtsUiEB4AS+LZrS9BXcQ5CepibAv16wFdGc9ufQnqIs5JUBdjW7AX\nLK6zX9/i0reRWnURpySoi7HtXPLpoM/iqIwQ8V+4PglxAUlQF2PbuQZ1kBSMiGsS1MXYdi7ljEdJ\nWaOIYxLUxdgmI3XxASNBXYxtEtTFB8ywgrpS6stKqV1KqR1KqceVUpaR6pgQwxaL6tPuWpLObTsJ\n6iKOnXdQV0plA18CpmuaNhkwAbeNVMeEGLZQP5idYDjLGvWjrG4ISVAX8Wm46Rcj4FBKmYAEoHn4\nXRJihAR79YdNnyur1KqL+HXeQV3TtGbgp0ADcBjo1TRt5Uh1TIhhO598OoDJDiiIBka8S0JcaKbz\n3VAp5QZuAAqAPuBppdTtmqb96eR1H3jggWOvq6qqqKqqOt/dCnH2zjeow1Be3ZQwsn0S4gxWr17N\n6tWrh92O0s7zJ6ZS6qPAMk3TPnPk/SeBOZqmffGk9bTz3YcQw1LzHCTmQ/r0c992/xOQNhU8FSPf\nLyHOglIKTdPUuW43nJx6AzBXKWVTSingSmDvMNoTYmSNxEhdiDgznJz6BuBpYCuwHVDAQyPULyGG\nT4K6+AA675w6gKZpDwIPjlBfhBg5WhRC3nOvUT/K6ob++hHtkhAXg9xRKsamkBfMDjCc57hFRuoi\nTklQF2PT+daoH3W0Vl2IOCNBXYxNw8mnAxjtgAYRqVUX8UWCuhibzmfK3eMpJVPwirgkQV2MTef6\nGLvTsUheXcQfCepibAoPo/LlKEui3o4QcUSCuhibQl6wOIfXhtmptyNEHJGgLsamsA/MicNrw5Ko\ntyNEHJGgLsaeWFSvWjE7hteOJVFG6iLuSFAXY0/YByYHqGH+8zZLUBfxR4K6GHvCI5BPB70NuVAq\n4owEdTH2hHx66mS4zA49jaNFh9+WEBeJBHUx9oS9w79ICqCM+kMywgPDb0uIi0SCuhh7Ql69HHEk\nSFmjiDPDmnpXiEtS2AeO7DN/HhnU/46yJILBePp1paxRxBkJ6mLsCXkh+Qzpl8ggrP/5UGVMLALJ\n42DCx06/vpQ1ijgjQV2MPaH3yKk3vgPJJVB5s/4+GoJ1P4P+JkjKPXV9KWsUcUZy6mLsOVNJY8gH\nTe9C0ZVDy4wWKLwCalecvi0paxRxRoK6GFu097ib9NBqyJgKds+Jy7NmwmA3dFefuo05Uf8yECJO\nSFAXY0t4QC9DVCdd+Az0QOtWKKw6dRuDEYqW6qN1TTvxM5mpUcQZCepibDlTOWP965Az98w3JaVP\n0kf5nbtPXC4ljSLOSFAXY0vIe2rg9ndC137IX3jm7ZQBxi2DupUnjtbNToj4QYtdmP4KMcIkqIux\nJXyaKQI690LaRDDZ3ntbT6le8hjoGlpmMILJLneVirghQV2MLacrZ+w+ACnj339bpcAzXl//eGbJ\nq4v4IUFdjC0nlzNGgtDfCO5xZ7e9Zzx0nRTULZJXF/FDgroYW0InPfGotw4Sc8FkPbvtk4uh7xBE\nw0PLpKxRxBEJ6mJsCZ90ofRsUy9Hme3gzIS++qFlUtYo4ogEdTG2nFzS2H1AT6mci5NTMFLWKOKI\nBHUxdmgxvfzwaFD3d0E0Ao6Mc2sn5aSLpTJSF3FEgroYO8IDevnh0Wl0u/dDSqle1XIunFkQ9sNg\nj/5ecuoijkhQF2PHyeWMXQfPPfUC+o1IntKhFIxMvyviiAR1MXYcX84YDesXO5NLzq+tlPHQfVB/\nbXZAZEDuKhVxQYK6GDuOL2fsq9dz6Wb7+bWVXAK9tfpDNAwmMNr0lIwQlzgJ6mLsOL6csfugnkI5\nXxYn2FOhr+HIe7lYKuKDBHUxdhxfztjXAO7C4bXnKtDvRgUpaxRxY1hBXSnlUko9pZTaq5TarZSa\nM1IdE+KcHR2px6Lga4HEnOG1l5QL3ib9tYzURZwY7kj9F8BLmqZVAFOAvcPvkhDn6WhOfaANbO73\nn5Xx/STm6s8uBSlrFHHjvIO6UioJWKhp2iMAmqZFNE3rH7GeCXGujs6l7j18+odInyu7R38w9dF2\nJf0i4sBwRupFQKdS6hGl1Bal1ENKqfMsNRBimLSYXnZoduij68QRCOpKHRmtH5YHUIu4MZygbgKm\nA/+tadp0wA98Y0R6JcS5Cvv1skODSc+Dj0RQB0jK0duT9IuIE6ZhbNsENGqatunI+6eBr59uxQce\neODY66qqKqqqqoaxWyFO4+iNR9GQ/vg6Z+bItJuYC80bIGuGjNTFBbV69WpWr1497HaUdvLT089l\nY6XeBD6jadoBpdT9QIKmaV8/aR1tOPsQ4qz0HoSWdyH7cjjwAsz64unXG+yHDQ9BoGdoWeYkmHyL\nPj3AyYL9sPGXMO9fYNOPYPa3zn0uGSHOg1IKTdPO+R/bcEbqAP8IPK6UMgO1wN3DbE+I83O08qW/\n6cwXScMB2Ph/kDUFCi/Tl8WisOUPsOd5qLzh1IBtTQJl0ts3mCESAHPChT0WIYZhWCWNmqZt1zRt\nlqZpUzVN+4imaX0j1TEhzsnR9MuZKl+iYdj0CCQXQckSffoAsx2sTph5D3RVQ83rp2/7aF5dLpaK\nOCB3lIqxIeTT7/o8XeWLFoOtj+sBfMKNp47GzXaY/VloWA8N605t+2gFjNkpF0vFJU+Cuhgbwj4w\nWCHUDwlpJ37WtFHPjU+5/fR5cwBbEsz5LOz7Owye9IPz6J2l5kR9P0JcwiSoi7Eh7NXLGp3ZQw/J\nAH2UXvMGlF0Dxve5hORIg9xZUPvmicsTc8HbrNfAS/pFXOIkqIuxIeTTK1tOzqe37tTTKylnOa96\n0eXQtAFCA0PLzHY9n64pSb+IS54EdRH/NE0fQfs7T5zES9OgehUUX3n2ZYh2t17iWP/OicsTc/TK\nFxmpi0ucBHUR/6JBQOkzMx4/Uu/cD7EwZFSeW3vjFsOhNRAJDi1LyoVBr4zUxSVPgrqIf2EfmBx6\ncLd5hpYfG6Wf4z9zZ5qerjm+EiYxBwZ75UKpuORJUBfxL+TVbxByZA6lWXrq9btGs6eeX5vFi6F2\nNUQj+ntHJgS69OoaIS5hEtRF/Av7QOPE+V7q3oaiRSdWwpwLVy4kZkLrdv292Q4mux7ko6Fhd1mI\nC0WCuoh/YR9Eo0NBPRKEjn2QM2147ebMgOZtQ++dmWCwSQpGXNIkqIv4F/Lq87o4jgT19j2QXKCX\nIQ5HxkTortHbhiPtG+RhGeKSJkFdxL+QV68rd2To71u265N2DZfZpl8wbdulv3dmQiwmI3VxSZOg\nLuLfYLc+c6LJCpFB6Dyo15qPhKypQykYZ6ae2pGgLi5hw516V4jRF+gFR7r+um03eIrOPD1uNAT7\nX4fAcVUsacWQP/P062dMgF3P6L8E7Kn69oO9I9t/IUaQjNRF/Ds65wvoo+qsM5Qx9rfC6z+FvhZw\nZ+t/rizY8yps+OOJNxsdZbJC6nh9ugGDUZ9f3d9x4Y5FiGGSkbqIb7GIPld6Uq5+QbO7Bqbefup6\n9eth+3Mw6cNQNO/EaQMK58DWp2Hlj2Hu3eDOOXHb7Kn6jUj5c8GeAoHuC3tMQgyDjNRFfAv7IIY+\nUm/bpV/YNNtPXKfuXdj7KlR9CcbNP3UeGJMVZn0cKq6CN38Nvs4TP0+vgN4GCPr0vLpMFSAuYRLU\nRXzzd+kTd9mTT5968XXCjudh/qfBlf3ebRXMhoqleipGiw0tN1ogrRxad0BSHoQHR/44hBghEtRF\nfPM16aWHkRB01504eZcW0wN0xdL3D+hHlVbpufP9K09cnjVFD+quQv2u0lhkpI5AiBElQV3EN18r\nWBOh66B+w5HJNvTZ/pV6gC6tOvv2lAFmfwIOvAE9jUPLU8dDTwMYrXr6ZqB1xA5BiJEkQV3EN38X\n2JKhY7+eIjmqpxH2v6EH6HOdpTHBA1M+Ahse1S/Cgv5rwJUL3bV6Dr6vYeSOQYgRJEFdxLdgr/5M\n0o59kFamL9M0vZpl8vV6gD4f+TP1KXgPHvdou7QyfT+WBPA1D7/vQlwAEtRF/NJiEPLrc6lrGjiP\nTBPQcRCCXr1U8XwpBRM/BAdW6fl60H8JdOwDqwsGpFZdXJokqIv4NdirB19/nz6KPlqquOcVvTzx\nXNMuJ3NlQ+o4qD3yaLukLP0GJZMTBnuG17YQF4gEdRG/fK1gMEBf01A+vbMGBrogf9YZN9OCA2je\nzqG/SPjM+6hYBvtX6bl1ZYDUMohE9LtYo++xnRCjRO4oFfHL1wJaFHobYUapvuzoKP00D8fQBr1o\n219Eq35Xz4vrSyEWRU39MKp0Aerk7ZLz9L+6d6Fkkf6L4NDbYDRDoBOcWRf2GIU4RxLURfzyNYNm\ngKRM/S7S7nrobzsll65FI2i7X0PbvQJVNBvDR76LsicNfd5ZT2zTM2h7VmKY+VFU3uQT91O5DNY+\nrN+NmjYedj4FiQ7wtUlQF5ccCeoifg20QUwNpV72vArlS8Aw9M9ai8XQ3v4dWtCP4bpvoJIyTmlG\npRZiWPYVOLyL2LuPoQa9GEoXDK3gKdS/OOrXw7gF4EiFSL++fyEuMZJTF/EpFoHBPgiH9ZSIt10f\nqRfNO7aKpsXQ3v2jHtCXfPG0Af0opRQqdxKGq76MtuU5tPrNJ65QsVS/IUnT9C+RcEhuQBKXJAnq\nIj75O/W7O2MxcOVA3VoomKPnugFN09A2PoXW24Jh8b2oI8vfj3JlYljyj8TW/QmtadfQB6kl+n+7\naiG9EiJR/UKtEJcYCeoiPg20gabAmarXq9dvgHHHjdJ3vYrWsh/Dkn9Ema3n1LRKycOw+F5ibz+M\n1ll/ZKHSfwXUrgV3PsQ0vRb+dHOwCzGKJKiL+DTQpqdA3LnQvBOSMiBRT69o3Y1ou1ZgWPIllPUM\nT0B6Hyq9GDXnNmLvPIJ2tHSxcLa+r0gQ7E4w2GCgfaSOSIgRIUFdxCdfK4QC4Bmnj56L5gOgxaLE\n1vwBNfMjKEfysHahimZBYjra9pf0BdZEyKyAhk36L4RYVC6WikuOVL+I+NR/GBT6FAE9jbDgswBo\nu1aA1YEqWXDazWL+PoIH16EFB44tM7qzsBTPPCXvrpTCMO/jxJ7/LlrBdFRKnv7lsf1ZKCiH7ia5\nWCouOTJSF/EnGoJBL5jN0LIPCmaB0YzW24K2+zUM8+9AnfR0o3BrNd7X/pfeJ/6VaFcTKOORPwPB\nfe/Q8+hX8a9/hqjvxEfVqQQ3asZyYmt+jxaLQHqpPheMwaaP1PubLuKBC/H+ZKQu4s9AO2AEYwwa\ntsKiL6BpMT3tMu16lDPl2KpaNMLA248RbtyJbdJSHAs/icHmOKXJaE8zg7tW0ffkt3Es/ATW0rnH\nPlMl89DqNqLteg01+Rr9gml7Ddgc0N9yMY5YiLM27JG6UsqglNqilHp+JDokxPvyteoXSZUFHB5w\nZcGhLfrt/mWLjq0WC3jpf+HHxAZ6cN36PexTrz5tQAcwJmfjWPgJkq7/Ov51T+Ff9zTakUfaHU3D\naLtXoA0emf2xvRaMRggH5Zml4pIyEumX+4A9I9COEGenuwaMJhiMQMFs/eLolr9hmH4j6sjMjJHu\nw/Q98x1MGcUkXnMfBov9fRrVmVLzcC2/n3DLAbyv/AotrJcsqsRUVNEstJ2vgN0FKQUQDuiljVIB\nIy4hwwrqSqlc4FrgtyPTHSHOQm+D/qBpnxdyp6HVrNMDbbb+fNKorxvvCz/BPvMGHPNuQRnO7Z+5\nISGJpOv/BWWy4n3tN2ixIyP2ydeiVa9FG+iB/NngHwSU/jQkIS4Rwx2p/xz4Z0Abgb4IcXYGugEj\nJKagmSxo2144MkpXaOEg3pd+gW3SEmzll533LpTRhHPxp9FCAfzrntKXJbhRpQvRtv8dciZDKAJ2\nN3RXj9CBCTF8532hVCl1HdCmado2pVQVeoGZEBdWyKdXnwz0QeZ4tANvQ3IOKqMETYvhe/0hTCm5\n2KZde8qmsaCfgepN+Gu2EAsFji03JaXiKJuLPa/yhKl3ldFE4tVfou+Z72JMzsZWsRA1aRmxZ7+F\nNvEqlDMJsIBXyhrFpWM41S8LgOuVUtcCdiBRKfWopml3nLziAw88cOx1VVUVVVVVw9it+EBr363P\nwjjQi5ZagrbmSQxL7gMgsP5ZYgEvSUs/f0JJY7Ctnp53nyHQuAd7XiWO0lkYHe4jn2qEOhrofvvP\nRL1dOMrmkTxvOUa7EwCDzUnStffR99wPMSalYc4pR1UuQdv6AsqTA52toPn0OWjOMc0jxPFWr17N\n6tWrh92O0rThZ06UUpcDX9U07frTfKaNxD6EAGDrH/RyQjRiKRPB14uh6nOEGncz8MbvcN38AIYj\nc6Vr0Qg96/5K//bX8cxfjqN8PsYzVL8AhHvb6Nv8EgPVm0i98h4cJTOOfRZq2KW3f+v3UEYjsWe+\niWHibNTBDWBWMP9e/e5WIUaIUgpN0845AyJDCxFfepsgEkOzm9FqNqEmXY0WCjCw+mEcVXcfC+jB\n9kMcfuybhNrryb3jByRNXfqeAR3A7M4g9cq7Sb/2i3S9+RjtL/03saAfAEv+RMxF0/CvfQJltqHK\nq9BaGsGdBsoErTsu+KELcTZGJKhrmvbm6UbpQoyoWBQGeiA0iBYYAFcmKqWAgXefwpxTiSV/EgCD\nTftoefoHuGZeR8aNX8PkPLc5YOx5FeTe8QOUyULzk98jGvAC4Jh3C+HD+wgd2qEH9fZDaLYEiMSk\nAkZcMmSkLuJH32GIRtEyytF6BjBMvIpw0x7C9dtIWHAbAIGG3bQ+/1+kX/sFEicsOmW6AE3TCHQ0\nMnD44LG/kLf7lF0ZzDZSl34ae8EkWp78HpGBPpTZhvOKexhY/Xs0ZUTlTkTr69QfnNHfqv9XiFEm\n0wSI+NG+G6IxsLkAA1p6Cb4nv43j8jswWB3467bT/vJvyPjwfdjzKo5tFgsH6d2/kf6arfTXbMdg\nsWKyJ+ofahrB3jbMTg9JJVNxlc7AmV+pPwlJKTwLb6PXZKHlL98h6+ZvYs6txFwwGf/aP+OYtITY\nyz9BK5mJat+lP7jDkTY650aIIySoi/jRugs0iLUcQOUUENj8AqbMEiyFUxlsPkD7y/9D5o1fxZY9\nHtAfZ9e98x1aVv8Ja0oO7rLZZC26BWvyiY+102JR/C019Fdvo+HF/4fZmUzu0jtIyCpGKUXy/OVg\nNNHy9A/Iuf1BEubfSt8T3yQyfh4GuxnNkICKRKHzgAR1MepGpPrlPXcg1S9iJGgxeP6f0SxJxHp7\n0CbOpn/DGty3fZ9YLMbhP/07qUs+haN4OgC+xn00vvIwymgkd8kdOPMr3mcHR3cTpWvbKlrefBJn\n4QRyl96J+UhOvuO13xH19ZBx41cIHVxHYPsKkjIMaG1+DMqPyi6D+V+4YKdAfLBI9YsY2/qb9Xx6\n1IDKGUeg5gC2KcvAmkDb8z8jaepVOIqno2ka7RtepPapn5A5/0bK7v6Psw7oAMpgJHX6Uiq/8Ess\nSWns+903GDh8AIDUxXcSC/rpWfs0ltK5KGUgFDKA1Q7mBOisl7y6GHWSfhHxoWUnWjSK1tdBNCWF\nSG8XjuuW0fnabzG50nHPvp5YJETDiw8RaK2j7O7vn5JmGeztoLd+N/2NB4iGBo8ttyQm4y6oxFVQ\ngcmmP/7OaLGTc+XHceSOp+bPPyT7yo+TOvVKMq6/j8OP/TuW1HwSLrsd34s/xjx9CbGajRhVCPxd\n4Ei9qKdGiONJ+kXEh9d/SKyzES1rBt76LdgnLyVIAv07VpH9sfuJhcNUP/F9rO508j98L0aLDYDQ\nQB9Na/9Oy5bXiYVDuAsrceWXY7IP1awP9nTQe2gP3qaDJKTnkTvvQ2RMWYTBqI95BjsPU/OXH+Eq\nnU7O0jsJtR+i5ZkfkH3b/QTf/F8MSSnY2uswGCKo6TdC6ZJROUVibDnf9IsEdXHp02Joz36ZmD9E\nuGIZwT0rsC/5As3P/ZLs2+7HkOCi+rEHcRZOJOfKT6KUYrC3g4a3nqV122rSJy4gd/6HcWTkn1Li\neLxYJExPzXYOvfVXAt0t5F92I9mzr8ZothAZHKD68e/izCsnZ+md9G9fiW/3W2TMmkv/m38jqXI+\nhoaNGDIL4Yp/vnjnRoxZEtTF2NVxAO31XxA12PAOajhybHS0GnFWXoazchHVj3+XhJxScq+6C4Dm\njSuoffUPZM1cSt6CG7Amec55l/2NB6hb9WcC3a1UfvSfSMobTyTg4+BjD5JUPJWsqo/R9uyPsCW7\nsQY7iZpSSGjZjsFpR33k5/AeXx5CnA0J6mLsWv8w0erNBJOyiCgb0eBhBmNppF3/T9T++QfY0vLI\nu+YzhPq72fvsLwn5eqm8+Ss4MwuONRHye+ltPEhPw0H6DtcQDQWPfWZxuknOL8WdV4oruwij2QLo\nNyq173ibA3//P3JmL6PwiluJhQIc+OMDJJfPJW36Eg7/8eukTywjUF2P02nBFB3EcNVXIaX4op8m\nMbZIUBdjlvbc14n29dIfNWGb+SHaV/+ZnLt+QsMrj2C0JlBw/b301u9h159+RM6cayi84hYMRhOx\nSJjmHWupW/sSPYcO4M4txp1fiju3BJN16ElIg31d9DQepLfhAAPdbeRMXUjR/GtJLihDKUWwv4t9\nz/6asN/L5Du+hVKKA7//dzLm34DNGKVnzZ/xTLuO2KH1JIS8GMbPQc25a/ROmBgTJKiLsWmwn9jT\nX2dQg4inhP62Jlz5GXhNhXjrd1H6yfvp2PUuB154iAm3fQ1PyVQioUEOvv40te+8QGJGPkULriN7\n8vxjI/D33l0Ph9avoH7ty5isNsqv/jjZU/SHbdStfJzWbW8y9a4HMCg48Id/p2j5l/G/+WuM+Qsw\n1KzFaVGYnU7UR358oc+MGOMkqIsxSdv5AtHNL9AfjEDZYoIte7BmpNO8az/j7/kBLZtfp+ndvzPl\nzm/jyCykeds77HjuITyF5VReeweJGXlDbWka/p4OelsOEQ2Hji23OBLx5IzDkuAcWjcWo23fZnb9\n7bdYE91MWf4PJGUV6vn6FX9k4sf/FRUZpOGF31A8JYv2fW14xk+Fph04TRqG5T/SH7EnxHmSoC7G\npNhfv0Ggt5ewxUZvrx/P9HnUv/06xZ/8Di3b3qL7wGYm33U/sZjG5sd/QtDXx5Tl/0Ba6RQAvJ0t\n1G5YxeFdG+huqsFoMuPOLsJktR3bR6C/h97mOqyOJFLySymccTkF0y7DbEsgFo1S987f2fvK4xTM\nvYoJ191JT+1O9jz5Mybc9jXCHYfo3vgcWbOWMVC3G1t/E4kWI+YZ16Gm3DRap02MARLUxdgT8hN5\n4sv0D0YI2l1YimbTsukVcuZcSW/ARPvONUz/zH/QdWg/mx//KcWX38D4K28BTaN63QoOvP0i/W1N\nepCevhBPfin2RPdpd6XFYng7W+io3UPtxjdoPbCdnMqZVC5ZTmbpZILeXjb/6WeEBvqZfde/Eerr\nYOfjP2Dix/+Vvrd+CwmpmEOD2GI9mIMhnMluDMv/8yKfMDGWSFAXY46262X87z5NwOLAPxghklKM\nKdQKKaU079rCtE9/n+p3/s6hdSuYdcfXSSmeQM26lWx74Q84UzKZeNXN5FTOwmDSbyKKhEP0tjXR\n29pINBw+th+rw4knq4DElAzUkUfSBX391G58g52vPoErI4/pN9xDamEZB994hoOrnmHG7V/BbDax\n58mfMWnxIg5v2kTqjKsJbH4ep8VEkkVhvv2nKFvSqJw7Ef8kqIsxJ/r0v9Db3o4vqmEpmUvPoYMk\npRhpqO0R8ngYAAAgAElEQVRkyj3fY8dfHyISDDDrzn+lr62JNY/+BKsziek33ENW+TQioSANuzdR\ns/ltmvZswdvdTlJqJu7MPMyW49Ivvl56mhsY9HvxZBdSNGUu42YsJKOonFg0wsF3XmL7i4+RNq6C\nubffh7+zhQ2//wHjr/worvQs9j/zX5RfdjnNWzaSnpeL5u3GEfWTOONa1IyPjuIZFPFMgroYU7RB\nL/4//hP9MSNRFaW7P0bqwpupffm3TLr7+2z76/+R4MlgykfvZevzf6B24xvMu/0+CqYvpL3+AFte\nfoLaLWtIKyilePplFEyeQ3JWPkaT+Yz7DPoH6Dh0gNqta6nd8jbh4CATLr+OKUuXY01wsP3vf2T/\nWy8y59Z7ySydxJrffJPsKQtIcZmoe/slCi9bTv+et7CF+0kwgNuTjOX2/7qIZ02MJRLUxZgSXfcY\nPZtX4o1oBKIxEqdeQ/2aFxk/OZede3pJL5tO5pTLeOt3/0FqUQVzbvsibbX72PTi4/S0NDLt6luY\nsOha7Ilu/N4+Dh/YRVdrIz2th+lpO0w0clz6JcGJJzMXT2YuqTkFZJdMwGQ209lUy/YVz7D/3ZWU\nzrmCmdfdTjQY4O1HfkhiWjazb/48mx/7MUlp6bgNXfT0G0lyJmDyN2OOxkgyaiTd/kNUUsZ7HKkQ\npydBXYwp/t/fS3e/l6DdQzQaobN/kJzpC9m3+nkKL78VrA42//W3zP/El0nMLGD1oz/H193B7Bvu\nYPzcK2k7VM2+9W9Ss30dHY115JROICU7H09mLskZOZgs1mP7Cvj66W5toqe1ibZD1XS3NlJQOZ2S\naXOpmLsYk8nE9teeYftrz1Jx2dXMuv6T7HrlCeq3vE3VZ77FgVcexTx4mITMSaBFiDXtICnJjT0a\nIKViFuYlXxrFMynilQR1MWbE2qrpfOo7eENRfGGNsMWJs3gmDVveJL98HN0RDy37trDoU99kz5oV\n7H3nFWbfeCcVC65h97uvs+nVZwj4+pmwYAklU+eRXzEFg8FIf3cHXW3N9Ha0njhStzvwZGTjycgm\nIdGF39tH3Y4NHNyyln0b3qRo4kxmXr2crIIS1jz1EPXb17Hg1s9js1pY/8SvmHXLvbS9+QjWrEng\n68Jh9BMLhHAwSEpiAomffgil5NEF4txIUBdjhv+Z+2k/VI0vZiRmNhK1p9HZ2UtmfgaNdQ1YknMZ\nf/kNrHz4P8kpm8rsG+9i+1uvsP7vT5BfOY2ZV32E3LLJ1O3ZSvWOTdTs3ERT9T4SEl14MrJITss6\nZaTe095CV+thlEFRVDmNkskzKZk0k7ScfHaveY1NK54lEgqx+PbP40nPYtUjP8aRnMbM6z7G2j/8\nmKICN36/EUd6PoN1m0lOTcMU6CHRZCDjms9hHH/ZKJ5REY8kqIsxIRYK0P7/PkfvYJigzYOvu5VB\nayrJ4yZxaOsbZBZXoLkL2Lnqb1Td8WW6O9p459k/MG7ybBbd/Ck6WhrZuPIF9mx4m9ySimPBOb98\nIkoZ6elso6ejjWgkcmyfNnsCyemZuD2p+Pp7qN21heodmzi4fSOBgX5mXHEdMxdfR6Cvm5WP/TdK\nKa742Odor97Fnrde4vJPfIn9Lz5Ecs44wgMBEp02Ai01JNosOE2Qnp6G685fjuJZFfFIgroYE/yr\nH+HwxlfpH4zgixrQDApb/lQa9++kIN9FXcsgJnsSExbfyGt//DWezDwuW343B3ds4s3nHic5LZOZ\nV36IylkLaWmsp3rXNqp3b6Ph4F78vn7cqekkp2ZgtgzNAxMY8NHd3oqvvw9PeibjKiZRMmEqJZOm\nYbFY2bL6JTavepEkTxpX3nI3JgOsevx/SS8oZlrVtbz9+C8pKSvC39ODxZWF1teGyxokFjVijwVI\ntpvJuesnGJOzRvHMingjQV3EPU3TaPn1XbT3+/GFomjOVKLhQVq7Byiau4ydK/9C6YIP0ef1c3DL\nGq742D/Q3FDH2peeonzGfOZfdwutTY1sXP0K+7ZuoHD8BIonTqVk4lSyC0uIxTR6uzro6eo4Kaee\ngCctgyRXMqFggEP7d1O9aysHd26lr7uTGYuWMvPypURDAV5/8hEG/T6uWH4n/u42Nr/2VxYtv5um\nzSswamES0woJDXhR3XW4swox+VpxmA3kFFfivvmB0Tu5Iu5IUBdxz7/lRepe/QP94Rhho51AMEx3\nv4+cGVeyZ80KJpans2F3J8VT5+FIy+KNpx9l8oLFlM1cyMbVK9jy9krKp81hxuVLSUxOpe7AXqr3\n7KB6z04625pJTknDk55JckoaJvNQvXpgYICezna6O9sY8PaTN248JZWTKKmcTG5BEfX7drLxjVfp\nbm/hihtuI6+4lDXP/wn/gJcrbvoEO1a9gMWkkZEELe0BsguKCbcdQItESbKCzaBIc1op+OL/YbQ5\n3+MMCDFEgrqIe4d+eRetXX0MRDQGYma83n6Sx03gcH0dWWWT2LPuDWZ86E7WrXwBV0o6FXOuYO2K\nF+hub6XqwzeTmJrBjo3vsnnNatyeFCbOmMu4iokkuj0Eg0F6u7vo7mynp7OTaHQop2612fGkpeFJ\nTcfl9mBQ0N7cSPXuHWxb/w5J7mRmLbqS4rJKDm7fyIZVrzB1wRWUTpjM2889Rk5xGZkZaexb8ypT\nlyyn+t0V5GalEouEUQEfDosBt91E4eQ5eG74l1E8wyKeSFAXcc23+00OPPtLugdCeAfDBLBit8bo\njTiIWhIgEsBsilF/uI9ZV93EprdfJxgIsOCam2g+3MRbL/+NvOJSZlx2BXZHEtX797Bn+1bqqveR\n7EmloHg8qekZeFLTcKekYj5ubvWAf4Cezg66OztoazlM7f69GE1GSsonMHHaLDKys2lvamDtypew\nWK0suvp6tEiIN5//C5PnXk5qWhqbXn+e2ZVpVDcMUDRhOu37NpKV7gIsmEP9OC0G0hJtlH/1UQxm\n25lPhBBHSFAXce3ATz7B4a4++oJR/KEoBncmwb42OgNGssdPZt/WNRQU5dPhN9HR3MTMK65lz/bN\nNNQe4LKlH8aW6GLL+jXs3bGVsolTmDxzDo5EN8FwiMb6ehrqa+np6qS7s4Oe7i6i0eixfVutNjyp\nqSSnpJKWnkFR8XjSMzMxGhTthxvZvO4dujvamb3wCkrKK2is3sfWtW8yb/HVJDkdrH/9JaYvXEzP\ngXVYXJnYLDaMRAh3NZKSU4TJ24rNZCTJZqRk2nwyP/L1UTzTIl5IUBdxq3fTi+x67iG6fUEGohA2\nWokZLbS2tOApnszhuoNkZ6eyfcd+piy6mpr9e+nv62Xq/MtpbDjEtvVrmTb3MvLGjaevr49tm9ZT\nX1tNSVkl40rH40lNx2yxogHhSJjgYPCE/RuMBqwWKyaTCS0WIxgYoLW5iYP79tLWcpgpM2ZTMXES\nRmDXlvW0NDUw/4qrsJlNrH/jFSbPXoBZaRzavYHJ02dxcPdOSion4z20C1tCEk5TFJvJiF1FyXDZ\nmfr1RzFJbl28DwnqIi5pWozN37uF5m4fvlCUQCRK1Oahp7uDqNGCxZVJc1MT6ckWeoMm+rx+SifP\nZNum9RjNJiqnzqazs5O1b71B8fhyissrMJostDQ3c2D/XurraklLz6CgsIhkTwoOhxOz1YrRMHSH\nZyQSZXDQz4DXS0d7G7U1B4lGYxSXllJWMQF3UhIDPi/bNrxLLBZjzmWLMCnY8NZKissrycrKZuua\nN5hSUUhXZy9OTwaRgV5SnQpff4AkuwWHxYBdxXAnmCmpmEjxPTLXunhv5xvUTReiM0KcraYX/of2\n3gH6AmECMQOhSIxALEh3v4/k7HEcqqvDk1PAnr1bKJ40m321W7E21JNTUs76NW8TM24nu2AclTPm\nsG7NOxxqOkxZ5USciS7KJ00nM7eIuro61q7fRH9fH0kuF263G/Px1S+BQXq6uxgcHCQ1NY2i4hJy\nc3NJTHQSjkTZu28vO7ZuxuVyM23mLHwDfnZv3YTd4cLh8rBj6xYsiR58IY2m1nYqM4voaKjDbkkh\npEHvQJBYBKJmE5qmcejgXrI7GrGn5b3HmRHi/MhIXYyaiN/L6gdvo6U3gDcYZlAzE4rF6A3EGPT3\n4cgspr62ltyScg7s3EhW8UQCoQh1NdVMmDqD3n4v2zZvonLSVDxp6XR1d7Nj23aUUsycPZuJkyZT\nXFJCSUkpBYWFuJOTMRqNZ+xPMBiks6OD2toaag4e5ODBA2zfto1tWzaTlZ1NWXkFdpuFhtpq2tta\nmDV7HpHQIAd272DazNl0Ha7FboridGcR9HuxG2M4LRAMRkgyxXAk2LEbYngSrBTlZDDr3x67iGdb\nxBtJv4i4s/Enn+FAdQ09/jC+UJRgTGMQC129XoyGmJ6G6evDmeTC19OKN2IlIzef7du2kZWThys5\nhQMHD+L1erly6VUsWXoVc+bNIzc3D6XO+f+FM4pEIuzbu5d33n6Tla++yrp31zK+rJyMjDTqqvdj\nt9koyM+nZs8Oxud76AsoPB4PA93tZCaZCYQVCcYoSQl27IYITouJ9KQE5t32D2QvkOeYitOToC7i\nSvu21ax66Lu09Q/qUwJENMLhEEFjIp1dHdgT7PSHDESUkXBwEKs5Rm1zH0Wl5TQ2HSYUCrP85pu5\n8ablTJk2DYPh4s2CGAgEePONVTz7zNO88tKLlBSXYLOZaaqvoSIviebuMAX5BXS3NpGbbMEfNmA3\nargSzNiNBhIsBtw2E7meRK750V8xWaTEUZxKgrqIG9FImOe+/CGaerx0DwTpD4TwhRWRWJSIOZH2\n9lbsiUn0DoTBbGXA20ei005NYyc33LSc2z/5SWbPmXtRA/mZ+P1+Xn35JR79/cNs37qVkmwHhzsG\nKBo3jq7Ww+SlOvEHwiRYLbjMUazWBBymKIkWEylOK5VlpVz+b4+M9mGIS9BFD+pKqVzgUSADiAH/\np2naKVPRSVAXJ1v5vbvYtfcA7f0B+gJhfMEI3mCEiKbAlkh7Wwsmh4uBwQhRDNhNMT79iRu55fP/\njsvlGu3un1FtTTWP/Og+Hn9xE5lZ2fR0tpKX7sbb14/b7SbJGEYZTThUjKQEC64EC5lJdpbd/RWK\nLpc0jDjRaAT1TCBT07RtSiknsBm4QdO0fSetJ0FdHLN/xeO89Ptf09oXoMcfwheM4I+CdyCANxTD\nZHfS19NJSNkZX1bOfV/5KksLe1FpkyFlwmh3//1VP4vPmMFjL2/lV//1c8wG8Pd3k56egSk6iIkY\nDosJp9WEy24iNdFOnsfBJ37+HHaXZ7R7Ly4ho55+UUo9B/xK07TXT1ouQV0A0N/eyCNfuZWGrgG6\nfINHRukhBjUTwcEA3YMa0ZjGjIp8vv6D/2HW7Dn6hjt+A8U3gCN7dA/gbDS+AWiQt5hwOMyzT/6Z\nH3z7a4QNCYQG+kmyGUl0OEgwajgtRlx2C+kuG+OzU7njNytG9AKviG+jWqeulCoEpgLrR6I9MfbE\nIhF+/7XbqWnvp6P/6Cg9Sv9gCO9gBA3F5Omz+dY9S5k9ezYUHAnomgaDPWBNHt0DOFs2D/RWA2A2\nm7n1459keamXv2wZ4D9++r94B7oZiAyQYNRItFkIhCKEojGisQ6eu/9ubvrO70e3/yLuDTuoH0m9\nPA3cp2ma73TrPPDAA8deV1VVUVVVNdzdijjz0Bc+xL7DXbT0+ekdCNEbCOMdDBNRJrJSXPzoPx5k\n0Y13wcZfgbtwaMPIACgDmOyj1fVzY0uGYPcJi0wpxXx8mZFbPrWP33/3C3z/oWcYjEUJhGMEggYC\n4SjBSITwxi2k/vGnLPzkV0ep82I0rV69mtWrVw+7nWGlX5RSJuDvwMuapv3iDOtI+uUD7tFv3MGb\na9fT3Ounyxekxx/EH9ZIsBj5/n/+lI+mVmO46kGIReDdH8Fl3wLDkfFGXx00vQET7hndgzhbkQBs\n/TnM/Ib+ZQTQtR8a3oJpn4H2vfRtf5GfvXqIh377MEaDIsluJtlhJdVppTAtkc998z+pXHTt6B6H\nGHXnm34Zbk3Yw8CeMwV0IZ7+z39hxVvvUtPez6EuH009A3gDEe699Wp2vPA7bqmaiCG9DIxm6G+A\nxNyhgA7gb4GEzNE7gHNlsoPRDsGeoWWuAvAe1r+0UkpxxXp48HvfZ/1//yPzZk2n2xeksXuAQ50+\n9jX38Yv7v0rj3i2jdwwirp13UFdKLQA+DixWSm1VSm1RSl09cl0T8e6vv3yQJ558in0tvdR2eOno\nH2RORREbnnuYby4rwl5RBW27IGOivkFv/YmpF4CBVnDEUVAHvb8DrUPvTTawp4C3GYwmSCuD9r3k\nz7mGpx64iz99+9M4ExJo7vVT09bH7uZuvvXp22ip2Tt6xyDi1nkHdU3T1miaZtQ0baqmadM0TZuu\nadorI9k5Eb+e+Nm3+Z+HHmLn4W4OdfkwGk38/tuf4+l/XEZeQR54CsDqgK5qyKjQN+qrB1fhiQ35\nW+NrpA56f/2tJy5zFepfWqB/ibXtgoJZ0LaPKz56N5t/8mnuvWYmfYEwNW1etjV28qXbrqWpWgK7\nODejf0ueGHN+/W/38l+//g27D/fQ2R/k5rllrHnwFpbNmYwavxB1aD0UzYfOA+DOB3MCRMPga4Gk\n42YujIUh0A0JaaN3MOfj5JE66L9A+ur11+kV0F2j59yzJ0GgE4vVyje+/CVe/49PUZyRTFO3n831\nHdx1/RL2bV13sY9AxDEJ6mJEfe3j1/KLhx+nut1Los3Kn792E9/53C24Ji5CNe1A5U4EX6cezFqP\nS714myAhHUzWocYCHXqJoMF8+p1dqk43UncXQt8h0GJgtoO7ADr2w7gFqLp1qInLMHXXk5/m5skH\nP8u3blmEbzDKjqZuPnbjdbz+7J9G5VBE/JGgLkaEpmlcO7OUx19cTZd3kFvmjeeZb95OZck4VG8z\nVpMBVTwX1bIDCo/UoLfvOS6fXjc28ukAVjdEQxAeGFpmSdR/kQy06e8zJ0HrTkgpAoMB5XSjAn3Y\ny+ZhjQS4YVo+T3/lQ5RmuKjv9HHPZz/LD//lC6NzPCKuSFAXw9bd2U5lrod1exuwmoz84o5FfGrx\nJOxaEEeSG/vExai6DaiKxXBoIxTN09MPCSlgd+uN9B46fT49HoO6UmcerR/Lq0+Ajr0Qi0LRfFT9\nOtSkazC2H8Bkd5JUMo3cvHx+es9VfOmqKQwEI/z0Nw9z7bxJF/toRJyRoC6G5U8P/y+VxYW09fqZ\nX5rBL+9eTFlRLp6sApJKZ6L1HsZq1FCFM1C9jeDKBmcaNG2GnOl6I9GwXs54upF6vF0kPeq0efVx\n0FOjv7a5IClH/7VSMBtadqPyJqH6W7FXXoahowaLyYg7NYOrZlXw008sJC3RxrvbDzAu3UVHR9vF\nPyYRFySoi/N29WWz+MKX7sNkUNy7ZAL3XDmVtIwMkpJTMYZ8mL2tJMz4MBx8BzVxGdS8A+PmQySo\nV39kT9Mb6q0DZ6aenjhKi4G/LY6DetapI3XPeOit1b/EAHJmQNNGvQooawKqYRNq4jKMzbsxpxeR\nmF+ONdSPw6SRl+bm326cydJJufT7g1SOK+BXP/nBxT8uccmToC7O2Z49u8jxJLJh6w5KM5L48tVT\nmFyYicflxK6FsGqDuMvnAgpzNIjKmYiKBMDfBTlToXUHeMaBNVFvsHMvpFScuJNgLxitJwb6eJJw\nmpG6xQGODP1LDCBrMnTXQtALpVfAwTdRxfOgsw575eXE6rfgGDeFlIJy3K4k0j0uPjSjhM9WlZNo\nM/Ptbz/AnMllhMPhi3984pIlQV2ck/s+/xkWzJqBFouybGI2ty0oJzc9mZQkBy5nAinFk7CnFxGr\n20jC3OWwdxVq8jVwYBWUVoHBCE2bIHem3qCmQdc+SC0/cUfxepH0KHuqfldp7KSAm1IOXUdqz002\nPbfevBU8+eBMRbXuRk1YiqH2XawVl2MzgNbThCunmLTsPLLS06gsyuGeqkom53moqTtEToqLFS+/\ndPGPUVySJKiLs9J46BAluZn88Y+PkpOcwPJZ45hTlk9uajJpqamk5hTgyhtPrO0gjtRMTNllmLrq\nUNnlKLMV2vbptemBHuhvhvRKvWFfix7oT65F97dCQtbFP9CRYjCBLRX87ScuT62Azn36lxnoX25N\nm/TXZYth/+uo8YvQ2qqxF04m0nIA94xrMPnaMUWCuDwe0tNSKchKZ9nUQpZNzCXBYuDW5Tdx49VL\nZNQuJKiL96ZpGg9+61+ZUlmG39fPzKI0rp5SSGluOllpKaSmZ5CcnoUxMojJ20by3I8Q3r8G+9Rr\n0PauQk27EQ6u1gO62aZfIM2aos/1AvooPaVCrxg5XryP1OH0F0sT0vQvsaPLU0og5IP+FsishGgE\n1duAmvIhtB0vkjDnI9CwHVtmMSmFZVhCXhLMihS3i9z0FKYUZnLNlAKKMxJ5Z8075Gek8NYIzPQn\n4pcEdXFGdbW1TCgp5Jf/9XMyXDYWlWUzsySLouw0MtNSSE1LI8FmweBtw1NUiSUlF615D7bJV2Go\nWaPXpVsToH4DlF6uj04Pb4TcWUM76dx7auoF4nN6gJOdrqxRKf1LrPNICkYZIGemfsFUGY6M1leh\nxi+AgR4siR40LUpiXhnh5v2kVMzB7UkhwRAh2WknI8VFcVYKC8tzmFWUhkVp3HDdMj720ZsIBoMX\n/5jFqJOgLk4Ri8V48N//jRmTK+nt7qQ8x828kizG56WSl5FKekoK7qQkLNEADouBtKlXEG7YSVLp\ndLT+dmwlM9HqNqImXwu1ayB7gl6P3tug78Cdr/832K9PA3ByfXp4AKJB/SaeeHa6kTroX2Kdxz31\nMXcGNG/Ra9bzZ0JPI8rbjmHGTWibn8Vx+Z0MbnmBlIW3EmvZh9XpJqN0MokJVlwJVtLcieSlupmU\nn8a80gwKUp2sfPUVinIyeP21Vy/e8YpLggR1cYLdu3ZSWVLEL3/+M9ITbUzLT2VKfjrjctLISUsh\n1ZNMos2MNRYgraAUZ2YRkUPb8Fx2K4Mb/4qj6m60bS+gJlyFMtvg4JswfrHeeNMGyJk1lGrp2g+e\nUj0dcbyjpYzx/mi3hAz9WLTYictdhRDo0r/UAJwZYE/Wpw0wmqFkkX5hOX8amCwY+1uwli9Ca9qN\nPWc8iZ5Uot1NpJdMIi0jE1eClVSXk5xUN+OyPEwrSGdibjJmItx84w3cfMOH8Hq9F/3wxeiQoC4A\nCAaDfPW+L7Jw7mx8vZ2Mz3QxuSCV8bmp5GelkpWWgsftJsGg4TArsifOJebrxm41Y8spQ3XWYima\njslkhI5aVOVivS49OQ/cuXreuHk75M0e2umZUi8DLXrpX7wz2fWSzMETn4SEwah/mXXtH1pWMB/q\n39JflyyC5t0oXzuGmcvRtvwN+7RriXQeIql4OqHD+8mceRXhw/twpqSTmT+OJIcNt8NGlsdNfkYy\n5TkpTClIIy/VyVtvrKK0IJen/vLExTt2MWokqAtWrXyNsqJ8/vjI78hIslGZ46E8N5Vx2WnkpHvI\nTE0h2ZWETQvh9rjImDCbQN02MqZWEe5owFUxl3DDLhJmf4TYuj+hpt+E0mKw7zWY9GF9J4fWQtYk\nsCXp76MhvV7bU3Zqh/rrILHg4p2ACykxXz+ek6WWD+XVQb8Ry9cOfU1gSdBz6zv/jsoohdRC2LsS\nZ9VdBN59krRln8O39WUy5l6HKdCD2ajIKi4n2ZWEO8FChsdNbkYK4zI8VGR7KMtJxm7S+Pyn7mHx\nwvk0NTVetMMXF58E9Q+wzo4Obr7xem6+6Xq0kJ+STDfluSmMy/KQm+oi0+MmPTUFd1IipvAAabn5\nuLPHMVi7leyqW/FueYn0qz+L/+3HcSz6BNSsA0sCqngu7F8FmRX6tAD/v73zjpKruvP85773Kudc\n1blbrVaOKCGBJCSSsEww2CRjsI2N7TUOc86EPWd3zM6s1zthx+Mw4xl7vGtsPOAxyQIMJoogCRBC\nOedWd6tzdVd1V3zv3f3jtZAQzSAMaommPudUV737blf96veqvu/W7/3u7xplOLIOGpeffPH+feCr\nsioWnopRskoGBJrG1BdnjeBEGNj/zvbwJOukphesbUWDhovg0Fpru3kZ9B2G/iMoC25E7nkRzRPE\nVjcd48AGQguvRT+8GV/tFLxeLzLbR6p5KqFgAL9TIxYKUJ2M0lgVozkVZmp1hNqohz07tjJn2mT+\n9nvfRdf1MXNDhbGjIuofQwzD4Mc/+D7TWibwyovPURX00FIVpiEWoDrsIxHwEgv6icciBHxe1GKW\n6onTcLhcMNBBfNFq8jueJ3zRTZR2voCtejK2WD1y2xMoi25FFLNw8CWY9gnrBdvesC6O+k4JqXRs\nhNS8dxqXOWJNsf+oLDT9XgQnWO/p9ElINpcVgunacrKt7kLo2WPl8mt2mLYKtq1BuIOI2Z/E3HAf\n7sU3Um7fg8MfRvOFcblcKIog3jSVYudBUi0zCQUCeDRBLBwiGQtTHQtRnwjRkorQnAoRdNv5++99\nl5mTJ7Jh3boxdUeFs09F1D9mbHz9NS6YMZW/+sv/hs+h0JQIMiEVoiERpiYeJhH0EfG7iSfiBPwB\nlPwgNZNnQTmPx2HDVz8Veo/iiNVjt9swug/jXnIL5mv3I6auRAQSsOspq0iVJ2xdJDy8FpqWnzQi\n32+t2Rmb/k4DBw5AsHmMvDEGaG5wJaxfH6dTtQDaXz85EcnmgpoFcHgktt6wCPKD0LUHMWmZ9Yvn\n6GZ8l32F3Cv3EV5yA/nDbxKbtphC5yFq568k17qDVMsMggEfTmEQDYdIhIPEAx6SQQ91ER/NiSB1\nUS/Dg31cfdXl3PKZG+jqHCVLp8JHkoqof0w43tHBzZ++nk9ctpJ0Tye1ES8TkmEa4iHq4mGq4xES\n0RBBr5NYNEo4Gscc7KZ26gXoQ/1EUnUIoeD1hzCG+ggvvJrc+t/gvfyriOO7IdNtFe0a6oFjm2DK\n5R5Vl9cAACAASURBVNYLd+2yRt3hU8IpHRshMefkBKQTSGmFKoITx84xY0HoXUIwoSYwS5A5Jcbd\neLGVs17OWxdUp6+G7WsQApQLb0O++QiqP4pr9iryL/+KxNXfZmD9b6m55CYyO16madn1DB/ZRmrS\nbAJuFw4hicfixEJ+YgEv1bEw9YkITakYE5IhEgE3a595illTW/ib732XUqk0dn6pcFaoiPo4J5/P\n87/++n8wZ/oUXnzuaWIBFw3xALVRP1UhLzG/i3jYTzwaxmdXiQQDJOqaKPUeo2HOYkrpDuKNUyj0\ntpGYuZShXS8TX/0Nhp7/Oa55V6N6gpiv3Y9y4WcRigabH7RSGE8U6zr0gjVKP5GeaBrQ+SZUzX+n\nsYU+kLqVCjieCDaPLupCgdQC6Hj9ZJsrBLEp0LrB2q6ZbcXbD61HRGoRzRdibrgPx6wrEDYnxuFN\nxK/8Cum1v6R+9V30vbaGCStvZPjQm9TMXIRL6DjtGolUNSGvk1goQCoepjoWpDrsoz7qoz7mx+/Q\n+P7ffo/pk5p5fM3vkCd+PVT4yFER9XGKaZr8+le/ZGpzIz/6/t8TcGrUx/w0xoM0JCLUp+Iko0Ei\nPjfxcBi/x4VbE9RMnEaxu5X6mQvJdR4iNWUeuWN7qF1xC/0v/zvxa/+E4hu/Q/WGcUxfgfny/0U0\nX4hItkDrG5AfgEkrLSO6d0Np2Frl5wS9u8AVAU/8nUafGKV/1PPTT8edtCZTnZ7aCJCaa/mknD/Z\nNmEFHHoRyjnLF/NugR1PQC6NmHMNDPfD3hfxrryT4p51qNIgMH81mfW/peG6b9C7/hEmrLyJgd3r\naViwAiWXxufzE08k8doF8UiIWNBHzO+iKhKwRu5VEepiQfT8EF+64zaWLV7I5k2bxs5HFT40KqI+\nzpBS8uzTf2D+rOn86bfuRkOnPh6kMRWiJhIg5nMRC3hIRkOEfW6ioQCRWAyRH6Rh+lwKPW1UTZxK\nsf84tbOXktn3BvWrvkjvUz8hdvmXke270HuP4r30S7DzGdCLiNlXQyELWx+BebdaI0tTh12/gylX\nv31yUcdGqF44uvHjMfQC1oj83Ubrdq9VZ71r88k2fwqS02Hf09Z2IAUTl8Kbv0EoGsqyLyO3PI7I\nZ/Ct+gbDL96Lp2YKjuoWhjY+Rv0nv0bv+kdouerzpLe/RNOFl1PqPkqysQW/y4HHoZGIRQj53CSj\nYWqTUZIBD0m/k5qwh7jfxdGD+1h12SXcfMOnOHhgFLsrnLdURH0csWH9OpYvWcTtt9xEuqeTRMBF\nVcBN3OekOhKiNhUjFvQSDXqJhkN4NKhuaIJ8hlRdE7KYI15Ti1nKk5o8l/7ta2n45FfoefKfiCz/\nLDYFCtuewXfVtyDdjtz5DMrSOxGKClsehIYFVglZgCOvWMvVJaaeNDDXZ1VljE17p/FG0Vp82t84\nNs4aa94ttRHeecEUoGWVVTogO3IBc/Jl1gj92CaEP45YeBPmiz9FDSXxLLudoad+RGjBtWi+MLnt\nz1J31ZfoeuHXlrDveJmJF61iuHU39dPnoRSyRGJxogEfQa+LeCRI2OckFQnSWBWnMRWhNuon4nXw\n8tpnWbJgHl+58wu0t7WdfT9V+MBURH0c8OqG9VyxcjnXX72aowf3EfZo1IS9NCQjNFTHiQe8hH0u\n4pEQQY+TRCKOXZhEozFcDgdeh0YwWYtLNdHsDiJVdQzu20jT9d+m76l/IXThDTij1Qyt/QW+Vd9A\n0WyYL/4byuLbEN4ItG+D9DGYepVlUDELB5+Hqde83dCO1yA51xrJn07mCHirrRrj45FAE2Rb35na\nCBBsBEwYPHKyzeGF5kutXztSWj6bfytseRgKWZSmBYhEC3L9fdibLsA5YwVDT/2Q6Mo7EDYHhX0b\nqL/mbjqfv4+Jl99Kes+r1M9eRLGvg/pJM9AHe6hpaMatGkTCISI+N4lIgFQ0RNhtJxn00JiK0pSK\nEPE6ePzRB7lg5jS+9fWvVcT9PKci6h9RpJS88vJLXLZ8KddfvZo927cQ8diojweZUJMkHvQSctup\nikWIBrzEQj5Cfi8hr4tINI5WzFI3ZQ7FzkPUzr2YfOt2ws2z8Ljd5I8fpumau+l94scE5l2Fu7qZ\nzBPfx7PsdtRwFeZz/4Romo+omw25NLz5G5h3s5VbDbD391blQe8pcfPCABzfBDUXjv6G0vvHVyrj\n6WguK7Z+qnCfQAioWwaHnn77aL1+CRQGoXuntR2uh4aFsPE+kCZi4Y3IbDdy8xqccz6BFqsn++QP\niV/+ZaRRprD7RZo+/ad0vvgbmpZ9imJ/F8FwCLvTRcDrwe3xEAmFcamCZDKB3+UgFgkS9jqpTcZJ\nhQOE3Taqo34aU1GifhcPPvBr5s2azt1fvYvW1qNj4roK74+KqH/EMAyDNY8+ykUL53HzDdexb9c2\nwm6NCVUx6pMRwm4biYCX2kSUsMdBLBwg4LKRiCfQ9Dw1jRMppTtpmDmfwUPbmHzpDfRtfpamlTdj\ndB7AyGepu+w2uh79OwLzPoGnfhqZx/4ez5JbsDfOwVz7M+vn/9zrQC/Bup9ZS7HFRmLhA8esC6QT\nL3u74YeehupF4Byl8qJpwMDe8RlPP5VgC6R3j74vOcfKQ+/ZfrJNUa1fO7vWnFzXdPpqa9btjscR\nmgPl0ruRhzci976EZ/kdqP442ad+SHzVfwEJmVcfovnGP6dnw6PEm6biiVZjpo9RO+di8u37qZs2\nFzPTQ6qmAbdqEgkGCflcxMJ+In43Mb+XhlSCqNdB3O+kIRUl5nfz6IMPsHDOLG6/9Wa2b9t29n1X\n4YypiPpHhGw2y7/+5J+ZOXkid3/1SxxvPULQqTKhOk5dIkzAoVATj5IMW/HRWDhA2O8hFPAT8DgJ\nhSO4FJNgLEHA40ZFUjNlFt1vPsvkq++ib8MjuFMTSC28iq413ye68vN4aiefFPTmBcgNvwbTQCz5\nHAJg0/3gi5/MdinnYfN9lhCdOv1/8BikD0Hd0tHfXN92cEbAPUpGzHgiPhv6dkFplIqJQoHmq+Dg\nUycFHCA2yZqNu/MRa1tR4cIvQOsmaH0D4fShXP5N5NbHoXUrnks+j+qPM/SHHxG74i7skWr6n/03\nJnzmLxg+ths1n6b2wtX0vPEkk1dcS+bobhqnzUUWh0lWVeFQJMlEAr/LSTToIxbyEg8HCLns1Fcl\nSQbcRL126pNR4kE3zz/9e65csYyrLlvJ0089iWma73xvFcaUiqif5+zft5c//fa3mNRYx//+q++g\nF4YIOASNVVHqkxECdkFtMkbU5yYa8BAL+YgEfPhcdhLxOHZ0apuaKaWP0zR7MZkDm6mduxQ50IFA\n0rziMxx/5hdUXXILvliK3j/8K8lr/gS7x09mzd+dFPRNDyH7j6FccpeVj77vOch0Wel2Qlhhg60P\nQKzFKk51AinhwBPQdBlojne+QWlC+8tQ/S6CP56weSE2C46vH31/qAm8VdB22tT9GZ+G/sNwbCSf\n3eGDJV+CzQ9Ztdd9MZSVX8dc/ys4vuctYc8+/n8ILboO7+Ql9Kz5B2ovvQ1HIMbg5meY+qm76du5\nnmRdE05fALI9pJqnY2S6qaprxKbnicfj+Jw2oiE/Eb+LqmiYgEOlLhUnGfQQ9dhprk4S9TnYuXUT\nX/7C7UxvaeYnP/4RA+n02fVlhXelIurnIYVCgf944H5WLl3C8iWLeOC+XxALegl7bFSFfUyoSeLX\nJLXJGBGf24qZh31Egj78LjuJRAKbUSBVU48c6iOSqsPvtCOLwzRfeBk9bzxJ3UXX4FSh57XHmXDj\nn1E+spnM1mdJ3fgdGOwk+4cf4135ZewT5iHX3Yvs3Idy2TesGumtb8C+F2DJnSfj6IfWWvXBp5x2\ncbR7m5XemJzzjvcJQP/IjNPxmvVyOqnF0L3ZWghkNCasgmOvWBebT6A54ILbYc/jMNhutQVr4ILP\nwLqfQrYbEa1HWfFVzJd+jjy0Ec8ln8dWO5XMw3+Nd8Jswss/S9fv/gF/opbEkutof/Kn1M27BIc/\nRP7oNiYtv4bskZ00TJ6FTRH4PS4CgQB+h0Y4GCTocRIN+wl5HFTFI/htCo3VKcIejUTQy8TaJD4b\nmKUcf/Pdv2JKcyNfvP02NqxfV5nINMaIs+1wIYSsHNT3RkrJqxvW8+tf3ssjDz2I3W4Do0zA48Tn\ncaMYJZLxGOX8EIl4HFnKE43FMAtDxBJJStkBUnWNDHW1UjdlFtn2g0TrmnFqCkNdx2hcdDn9O17G\nl2ogOX0RnWsfIDR1MZHpF9H71E9w1k4hvOxWCq8/TKl1O/5V30TxhjFf/ClIibL8LmsB6UPrYdeT\ncPHXrPxpgL6DsPlXsOSb1ozIE+gF2PhDmPLpkQyP09+0Cdv/BWovhVDLmPj5vODQY1ad9dqVo+8/\n8KRVf37KDW+fiNX+Jux7Ci769snw1uEN1sSkpV+DQBUy3YH57A8QUy9FmXYZxf2vMfzyfXiW34ES\nrqHrsR9gC8bxzb+Wtj/8HIDg7Es59Nz9OMNVKL4orRufJzJlPscP7kZodkybk+6jB3FGqmg/tBdP\ntIr2o4dxBqN0tB3D5gnQ3d0NmoOB7BDZfBEDld6BQYRqx+f387k7vsCNt9xCQ8PH5OT9ISCEQEr5\nvmfiVUT9HCKlZOvmzTz04G/5j/t/TaGQxywXCQcDCKNMNORHxcTnceG023BpKm6XE6ddxeWw4/d6\n0QvDJKpqGepupXridAq97fijSYLROP37t1C/4FLkcC9D7YdoXPEZ8ke2kes8TO2qL2Mc30Nm+wtE\nLrkNd3ICQ8//DGF3473sqwi9gLn2XxGBFGLJbVbIZf9aa4S+7OvgjVlvYrANNv4bzLoJYqcseGEa\nsO2XVq56y9WjO6B/D7SthRl3jb9ZpP8ZhX7Y8TOY/Y3Rq1HqRXjzX6xSCjWL375v56OQaYd5X7QW\n8gYrvr7lIbjoKxCuQw73Yz79A0RqMmLeDRj9bWSf+jG2+lm4FlxHev1DDB/YSHjprQwP9NK17mEi\nF1xOSZe0v/oEsZnLGOzroffQLkITZ9O2+01s3iAlA/qOt+KKVNF+eD+eSJLjba3YPUG6ujrRXF66\nu3vQXG66evpQ7E56+vrR0RguFCnrJhNbJnHLZ2/j6muvo6q6+uz7+iNMRdQ/Iui6zqsb1vPEY2t4\n+Le/JZcbRhplAl4PwiwTi4YxCjkSsSjl/BDJRJLi0ADJVBXFbD+xZBWUcng8Xnx+P/n+Lqqbp5Dr\nOoY3EiMUTZI+tI3U9EU4VBg4uIWaRavRZIn+rWuJL1yNL9VA/9pf4Ug0EV52K/qhjeTffBz3/Gtx\nTF+BPPwG8vXfIKZfbi1LJ03Y+QS0bYGlX7eqLwL0H4JN98KMG95eCgBg3xprybYZn3vncnVgxdp3\n/AyqLoLI1HfuH+8ceAScYahZNvr+fNoS9smfgsgpC4lIE3Y8Yl18XnCnNSMVrLkCm+63ctlT05HF\nYcz1v4RMD8rSO5HuALl191Pu2Iv3ki9iIOh55udo/ij++VfT/frvGWrbR3ThJ0m37qVv/2Yi05bQ\n13aEbPcxfHVT6Ni3DZsngKna6Ty6H2+sls72o2hOL0PDwwznC5hCoz+dRrE76e3rR2gO0pkMulTJ\nDA+TKxnopqRY0pnY0sKnb7yZVZ/4BBNbRlks5WNORdTPY453dPDC88/x+8cf44XnnkXTVPRSgYDP\ni2KUiYSDmKUCsUgYvZAjGo1iFnOEQiFsCqiYhKNRCuleErX1GLlBnC434WiCzPHDxBsm4XTYyBzb\nR9X0RdiEzsDBbaTmrsDldNK/9QX8zXOITFvM0JanKQ90El1xBzanm9yG34Ci4V1xJ4rTg3ztAWTv\nEZRldyIi9TDcB6/da8V1F9x2cuWinj2w5X6YfYuVoXEqbeutGZIXfOXdJxP1boP2l2Dm16zMj48b\n+R7Y+f9g+pfAGRq9z2ArbP8VzP4ieJMn26WEvU9YFTAX3gXOgNXee9A6VlUzYeY1SEVDHliHfONh\nxOzViEnLKbduZejFe7HXzsA5dzXZva8y+MbjeKdchK16CsdfeRizXCIw7SL6Du8k07afYMt8Bns6\nSbcdxF87if6udvJDgzjDSTqPHMDmCVAo6aT7enAEInS2teLwhejvT2MgyBXL5PMFShKyQznKKAwN\n58jrJsWSgSnBHwhw1Sev5oorV7Hk4qV4PJ6zfwzOcyqifh7R1dnJhvXreP7ZZ3j+2Wfo6enBblNB\nmvg9TjQg4POAXiYSDqAXckTCEcxSnmDAj8NuwyzliSeSFAb7CUVjOB02Spk0sdoGZGEY9BKxumb0\nwR6QBrEJ0zCzvRTTXSRmLcWGwcCudQQnLSA4YTb5vesodh8ltOg6XKkm8hsfxehvx73gOmwT5sHe\nl5Dbn0I0XICYdz1CtVsldLc8BJMuhZZLLPGVplXv++ALcMEdED4tRtq1FQ78HubeBa7w6A7KHIF9\n/wFTPgee5Oh9Pg4cfxW6NsL0L1p110ejaysc+gPM+gK4o2/fd+A5q5rjzJsgOjJxq5SzRuzZblh4\nuxVnz3RhvnIvlPMoF1yPGW2kuOVJCjtfwDHlYuyTl5LZ9hzZnS/jnboUNVpHz5vPUMr245+ymKG+\nLnp3v46/cQYlAzr3vIknXosuVLoP7cEZSVIql+k73oYrnCTd10vZMMHmoLvzOE5fiP7+fgwpKOgG\nQ0PDSM1Of3oAYXOSHhjEVDQGh3JoNjuFYonJU6Zw6eVXsmz5cuYtWPixFPmKqJ8jCoUCO7dv5/XX\nXmXtC8+xedMmBgYGsGsqSAO3y4mGid/nQTF0/F4PmpA4bBo+jxujmCcSi0KpgMNhJxAIUMj2E4km\nsGsqpeE00WQtqjQoDw0QqWlAQ5LvO060cTIOh53htn34qhrxJ+ow+o9R7O8iPG0xrkCY/L7XMMsF\nAnOuwBmIUtz9EnrPEVxzV+OYfBEc24Lc8hiEqlEu+BQiWAXd+2HHY1Zsd/6tEBqp55Ltgm3WyJ6Z\nnwHPKSJj6tYFvr49MP1Wa6m60ch1w+57ofn68bNk3Qfh6NMwdMw6wSm20ft0bLSEfeLVkJj59n1d\nu2DHgxCfClNWW7+MpIQjr8L2NZCaDtNWIV0haN2Cuelh8IRQZn8S0xMh/8YaSgc3Ym9eiDZhHkP7\nNjK0+xWc1ZPRUi1k2vaTPbgFb+NMdM1J+tAODEPHHq1jKN1LtrsNZ6yWfC7HQOcxnOEkuXyOwb4e\nHIEo2cFBcvk8mstLX3cXisNFoaQzOJhBcbhIDwwgVDtD+QKFsk7JhKHhHLoUDOeLqJqNQrFEXX09\niy5czNLllzBv/nyaJjSjKOP7F945EXUhxJXAP2KlRv5cSvk3o/QZF6IupaSt7Rh7du1i65bNvPbq\nBnZs3053Vxd2mw1p6jjsNmyqwKGpeFx2VGnidbvRhInTbsftcoBewu/3owmJIg2CoTCylMduU/EH\nQxi5IZxOB75gGD2XwaYq+KMJlHIRI5chUFWPw26j2NuGKxjFG6tCKRcodB7CW9OCJ1GHGO4j37oT\nV900PE1zUPKDlPatQzg8OGesxJ5ohoPrkQfWQ7gGZeZV1gLH3fth77Mw1AvTroK6C6zReTFrFeg6\nuh5aroT6C98eMsmnYef9Vv705Bveue7oCUoZ2PFzqFsJ0Zmj9/m4IU048DBIAyZ++t1DUdl2y8fh\nidD8ibfXzynnrRoxfQeg5QpIzQFVs0bt+56Hgy9D3XxoXor0hJH71yF3Pguqhpi0DJmcTHHfBgq7\n1qIGk9ia5lMuF8nuWoeRy+Csn4kuVDJHd1PK9OGsaqFULjPYcRCEhuqLkssOku3pwB6MUzZM0l3t\nqC4fUrWR7ulE2BxIxUZfbzeK3YluQjqdRrW7yBVKDOeGQbMzNJRDl1DQDXKFIoZUGMrn0U1BoVRG\nNyUSBSlN6hsamTNnLgsWXci0GTOYMmUqwdC7hLI+goy5qAshFGAfsBLoADYCN0kp95zW7yMj6oZh\n0NHRTuvRo+zbs4ed27exa9cuDh86SE9PDwJQBAgkNpuGTQGXw4ZdU7ErArfTgSokDpuVnaJKE6fT\ngctuR5hlenNlJlZHoFzE4/HhdNgwi3ncHg8utxtZzKGpCr5gBAUTc3gQTziK0+2F0jCyOIw3Vo3D\n5UYODyKLWdzJRpweP6I4RLm3FWeqGWesFk1I9LadUC5ha5qLIzkBZbgP2rZBtgfRvBjRcjFCEXD0\nDTj6Gtg90LQEGhdZopHpgMMvQ+d2qJpt1fl2nxJSyfdZOdVdW6FhBdQsGT2LRZrQsxXanofkIqha\nwtq1a1m+fPlYHdrzmrXPP8fyVLtV7Kv+CvDWjN5RL8CehyHTBrWLITX/7RO6eg/Aweesyo71S6Bu\nkVUYrJCxxP3I6+CNQsNCZNUs6D+GuXctdOyGZAtUz0A3oXRsB+VjO1Cj9YhYA+VCjnzHAfRsH7Z4\nI7rQyKd7yPccQwslkZqTfHaA4b4uNH8EE41cJk1+eAjNG6RULpPp70HYXaDYyGTS6LoEm4Ph4WGK\npRJStTOcy3OkK00y4mc4V0Q3THQJuXyBsoSiblAoltClIF8sUywblE1JqWwgFAXdMLDbHVRVVTFh\n4kSmTZvOtBkzaWqaQH1DA+FIBPERyrI6F6K+CPiOlHLVyPZfAPL00fq5FnUpJblcjr6+Xnp7emg7\n1srRI0c4euQwrUeO0tHRTk9PN4ODGYrFIooi3hJvVVHQFIFNU3BoKnZNwW5TcaiqNSK3azhUFVUB\np90Sd1WAc0ToNQXsmorDbsemwHPbW/nk/InY7RpCL2PTNNxuD6oCslzA5fLgdHtQMJDFYVxePw63\nF0VI5PAgNocTpz+MZrcjynkY6sceSmD3hVFtNpR8FnOwEy1SixaMY3N5UPIDiN7DoNoQNTMRkTrQ\nVETvAejaZwlu7RyrUJQ3Cukj0L3Hqt9ilqFuMdQvOpllURqGgYPQvQMGDlllY6sXgcM/ivNNK35+\n9GnrJFF/BfhqAbjnnnu45557xuhTcH5zzz33cM93/tI68R17HgINUHOJlR0zGpk2OPYypA9aVS8j\nk8Bfd3J5wMxxOPyidTL2Ja2VlOKTwJOA7n1WaKZrr3W845OQwRpkLgNdB5DtO8EdgEg9ZcWOPpxB\n72/HHE4jQjUYmgO9mKeU7kIvDiM8EQxFo5TPURzsxVRs4PCgGwb5TJpSuYzi8KAbklx2kLJhgOag\nrBvkhrNIoWIIhUKhRKlU5Jmd7Vw8uYZ8oYAB6AYUSiVKummdcMplCmUDw5QUdYOibqAbJoWyQals\nUDZNS+wNiWFKdMNEAuaIBrndbkKhEIlEkuqaGurqG6hvbKS+roFUdTXRaJRwJILdbh+LQ/+f8seK\n+ig1UM+YauCUxRVpAxb8sU8mpaRcLlMqlSgWCmSzWbJDGbKZLIPpfgYGBhgYGGAw3U9/Xx8DA2ky\ng4MMDg4ylM0ynMuRz+UoFAuUSmX0so5uGBgjtSgUYTlJCFCEQBECVbFuNlUQdCjYPW5smoJdVbFp\nYuRewaGqaKqCXRPYVBWbqmDXFDRVwaaqbz3WFCv0ommq9byKwKapaKol8H6XjZTXjs2moikOVGGi\nqiY2mx2bLYCGRJjDaDYbdm8ITVVQZM46Afhc2JwuVKWMohdR9DyKy46mmqhGDhUVVdVR3XZEoRsy\nJYQeAIcLkZoA+QFEx0Yr5u2LWTVbmhcDJgx1weZfWmEWf5WVbz7zBkvIiwPQuQXyvdZamvl+CDZY\nQjL5+pMjRdOActaqa5LrhMHDlqDb3FCzHMLTPl656O8XoUB8jpXe2bHOSvfU3BBoBH+DVRvH7rfa\n/DUw7WbrWBx/wyqWNtxlCbs3aa0sVTvPWkZwqBv6D8Lmf7fCZN6YdYynXmYVBhtOIw6+hBjqAdNA\nphpAcyL1PPZCL7b8IJQHMD0uDFnAKJcwTB2HQ0U3BaaeRap2DLtADwYoF3IY5UGkasMX8FLWy5SK\nBXTTwB90oZtQLuuUyyaGzY2BQlk3KdsFuq7ic2jUBJ2UdY2ybqJLiWG6KOm6tW2alA1JaUTMy4ZJ\nSbfuy7ol5iVDp6zLt9qKb/UzKOkGQ/3dDPR0smvbZkwJhpRIU2JKkEiktD6qqqKgaRo2TcPhsON0\nuXC53Hg9Hnx+P/5AgEAgQDAUJhwOEwiFCAaCBMMRAoEAPp8fr9+Px+3G4XTicDhQ1VHSez9kPoio\nnzFRr5MTY3XrhCkZuWO0MbwY+XNCAoQQp7QJFGE9Vk4T6RNC7dEEmkPDptqxqcrIzRJom6JgUwV2\nTUFVVWxCoNlU7KpAU6y+mirQNBVNEdgVBVU9OWLXVAVNiLdEXFMVFEXBNvJYUwWqGLkfGelb/68Q\n8jhojHlQhBhpF9gUFaGAKk70B0VRUBUTRUirn3CgKKAoBYQQKIpAsWsjJ6qcteyZIaxQikNBqALk\nIOQzkMc6owHYBYgcDLdC/pglJIpg5MnBqUGpC9qOQ4cKQrXKAKg2UO1gc4CrFoQB/duhd5M1mjdK\nYBTA5rGExxm1ZojWXwGOwFn7XI1LVAfUrrBOhLkuyByG3u3WybWUsXytuazjodism9MBjlooFyC9\nF3rLVgVNs2ydbKVpHUuXHfQB6E1bbaYcuTdBk2CaiKFWkCAk1pdVSKRbRTFLaOWiNeVfStAkpsc2\nIoYFTMPE1EB6VHRTwTBNDNPAtAlMpxPTNCmbEsOQmFLFMG3WSNq0RtTWzSTktjMh5rW2DUvEdRMM\nc+SxYaIbEkPKt0T9xIhcNyRl08QYEf6ybqCboBsGZcNEN0fajBHBN6yc+ZJujvSRlE6065KSeeL5\nTcxSjqFCjkx/L8dNiYnENK3B6AlXyRFBs7blaZr3dsRbf07cibfGPKfu+2P4oOGXe6SUV45sv2v4\n5Y83r0KFChU+vox1TF0F9mJdKD0OvA7cLKV8l4LRFSpUqFDhbPNHh1+klIYQ4uvA05xMaawIDqGa\nCwAAA6JJREFUeoUKFSqcQ8765KMKFSpUqDB2fGhTsoQQVwoh9ggh9gkh/vxd+vxQCLFfCLFFCDH7\nw3rt84338oUQ4hYhxNaR2ytCiBmjPc9HnTP5TIz0my+EKAshPjWW9o0lZ/j9WC6E2CyE2CGEeGGs\nbRwrzuD74RdCrBnRie1CiDvOgZljghDi50KILiHEu64J+L51U0r5gW9YJ4cDQD1gA7YAk0/rswp4\nYuTxQuDVD+O1z7fbGfpiERAYeXzlePTFmfjhlH7PAY8DnzrXdp/Dz0QA2AlUj2xHz7Xd59AX/xX4\n3gk/AH2Adq5tP0v+uAiYDWx7l/3vWzc/rJH6AmC/lPKolLIMPACctgQO1wC/BJBSvgYEhBCJD+n1\nzyfe0xdSylellIMjm69i5fyPN87kMwFwN/Ag0D2Wxo0xZ+KLW4CHpJTtAFLK3jG2caw4E19IwDfy\n2Af0SSn1MbRxzJBSvgL8Z2v/vW/d/LBEfbSJSKcL1el92kfpMx44E1+cyp3Ak2fVonPDe/pBCFEF\nXCul/AkfKDP3vOdMPhMtQFgI8YIQYqMQ4rYxs25sORNf/BiYKoToALYC3xwj285H3rdujsnkowqj\nI4S4BPg81k+wjyP/CJwaUx3Pwv5eaMBcYAXgATYIITZIKQ+cW7POCVcAm6WUK4QQE4BnhBAzpZRD\n59qwjwIflqi3A3WnbNeMtJ3ep/Y9+owHzsQXCCFmAj8FrpRSjsel18/ED/OAB4RVZSkKrBJClKWU\na8bIxrHiTHzRBvRKKQtAQQjxEjALK/48njgTX3we+B6AlPKgEOIwMBl4Y0wsPL9437r5YYVfNgLN\nQoh6IYQduAk4/Yu5BvgcvDUbdUBK2fUhvf75xHv6QghRBzwE3CalPHgObBwL3tMPUsqmkVsjVlz9\na+NQ0OHMvh+/Ay4SQqhCCDfWRbHxOO/jTHxxFLgUYCR+3AIcGlMrx5aRwiej8r5180MZqct3mYgk\nhLjL2i1/KqX8vRDiKiHEAWAY62w87jgTXwD/HQgD/zwySi1LKf/oYmjnI2foh7f9y5gbOUac4fdj\njxDiD8A2wAB+KqXcdQ7NPiuc4efifwK/OCXN78+klP3nyOSzihDi34HlQEQI0Qp8B7DzAXSzMvmo\nQoUKFcYR43s9qAoVKlT4mFER9QoVKlQYR1REvUKFChXGERVRr1ChQoVxREXUK1SoUGEcURH1ChUq\nVBhHVES9QoUKFcYRFVGvUKFChXHE/wdVa6Mb9Y1w/wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa6e14a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# your code\n",
    "L = 1.\n",
    "kappa = 2.1e-6\n",
    "dt = 25\n",
    "nt = 200\n",
    "nx= 100\n",
    "umax = 10\n",
    "sigma = 0.05\n",
    "dx = L/nx\n",
    "\n",
    "beta = (kappa*dt) / (dx**2)\n",
    "\n",
    "##########################################\n",
    "\n",
    "def tridiag(nx,beta):\n",
    "    A = np.zeros(shape=(nx,nx))\n",
    "    i,j = np.indices(A.shape)\n",
    "    A[i==j] = 1-2*beta\n",
    "    A[i==j-1] = beta\n",
    "    A[i==j+1] = beta\n",
    "    return  A\n",
    "\n",
    "\n",
    "#########################################\n",
    "\n",
    "u = np.ones((nx,nx))\n",
    "\n",
    "\n",
    "xvals = np.linspace(0,1,100)\n",
    "u =  umax * np.exp(-((xvals - 0.5)/0.05)**2)\n",
    "\n",
    "b = np.zeros(nx)\n",
    "b0=0\n",
    "b1=0\n",
    "b[0] = b0\n",
    "b[-1] = b1\n",
    "\n",
    "A = tridiag(nx,beta)\n",
    "\n",
    "for i in range(nt):\n",
    "    u = np.dot(A,u) + np.dot(beta,b)\n",
    "    #plt.plot(u)\n",
    "    if np.mod(i,10) == 0: plt.plot(xvals,u,color=plt.cm.copper_r(i/200.))\n",
    "\n",
    "print beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stability estimation\n",
    "\n",
    "Unfortunately, the previously described explicit scheme is not unconditionally stable. \n",
    "\n",
    "**Your task**: try out with your code above what happens when you change parameters such that $2 \\kappa t > x^2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# your code\n",
    "L = 1.\n",
    "kappa = 2.2e-6\n",
    "dt = 25\n",
    "nt = 200\n",
    "nx= 100\n",
    "umax = 10\n",
    "sigma = 0.05\n",
    "dx = L/nx\n",
    "\n",
    "beta = (kappa*dt) / (dx**2)\n",
    "\n",
    "##########################################\n",
    "\n",
    "def tridiag(nx,beta):\n",
    "    A = np.zeros(shape=(nx,nx))\n",
    "    i,j = np.indices(A.shape)\n",
    "    A[i==j] = 1-2*beta\n",
    "    A[i==j-1] = beta\n",
    "    A[i==j+1] = beta\n",
    "    return  A\n",
    "\n",
    "\n",
    "#########################################\n",
    "\n",
    "u = np.ones((nx,nx))\n",
    "\n",
    "\n",
    "xvals = np.linspace(0,1,100)\n",
    "u =  umax * np.exp(-((xvals - 0.5)/0.05)**2)\n",
    "\n",
    "b = np.zeros(nx)\n",
    "b0=0\n",
    "b1=0\n",
    "b[0] = b0\n",
    "b[-1] = b1\n",
    "\n",
    "A = tridiag(nx,beta)\n",
    "\n",
    "for i in range(nt):\n",
    "    u = np.dot(A,u) + np.dot(beta,b)\n",
    "    #plt.plot(u)\n",
    "    if np.mod(i,10) == 0: plt.plot(xvals,u,color=plt.cm.copper_r(i/200.))\n",
    "\n",
    "print beta\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens? And why? $Unstability occurs$\n",
    "\n",
    "Preview to lecture: The stability for the explicit scheme based on the Neumann stability analysis. The stability criterion for our scheme (CD in space, FD in time) is:\n",
    "\n",
    "$\\beta = \\frac{\\kappa \\Delta t}{\\Delta x^2} < 0.5$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implicit scheme to solve the diffusion problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will talk in the lecture about other possibilities to determine the spatial gradients. If we take the gradients at time step $k$, then we obtain the explicit scheme from above. Let's now have  a look at the case when we consider $k+1$. Applying the same Finite Difference schemes, but now for $k+1$ in the spatial gradients, we obtain:\n",
    "\n",
    "\n",
    "$$\\frac{u_i^{k+1} - u_i^k}{\\Delta t} = \\kappa \\frac{u_{i-1}^{k+1} - 2 u_i^{k+1} + u_{i+1}^{k+1}}{\\Delta x^2} $$\n",
    "\n",
    "If we now separate the variables for a new time step $k+1$ and the current time step $k$, we obtain:\n",
    "\n",
    "$$u_i^{k+1} - \\frac{\\kappa \\Delta t}{\\Delta x^2} \\left( u_{i-1}^{k+1} - 2 u_i^{k+1} + u_{i+1}^{k+1} \\right) = u_i^k$$\n",
    "\n",
    "For simplification, we substitute $\\frac{\\kappa \\Delta t}{\\Delta x^2} = \\beta$ and obtain:\n",
    "\n",
    "$$- \\beta u_{i-1}^{k+1} + (1 + 2 \\beta) u_i^{k+1} - \\beta u_{i+1}^{k+1} = u_i^k $$\n",
    "\n",
    "So, we now obtain an equation with several more unknowns on the left hand side! However, this equation should, again, look very similar to the solutions of the Laplace equation when we write it in matrix form! We separate the left hand side into a matrix $B$ and a column vector $u^{k+1}$:\n",
    "\n",
    "$$B = \\begin{bmatrix}\n",
    "1 + 2\\beta & -\\beta & 0 & 0 & \\cdots & 0\\\\\n",
    "-\\beta & 1 + 2\\beta  &  -\\beta &  0 & \\cdots & 0\\\\\n",
    "0 &  -\\beta & 1 + 2\\beta  &  -\\beta & \\cdots & 0\\\\\n",
    "\\vdots & & & & & \\vdots \\\\\n",
    "0 & \\cdots & 0 &  -\\beta & 1 + 2\\beta  & - \\beta \\\\\n",
    "0 & \\cdots & 0 & 0 &  - \\beta & 1 + 2\\beta \n",
    "\\end{bmatrix}$$\n",
    "\n",
    "The scheme, written in matrix form and, again, considering Dirichlet boundary conditions, is now:\n",
    "\n",
    "$$B u^{k+1} = u^k + \\beta b$$\n",
    "\n",
    "\n",
    "The implementation is now a bit more complicated as for the simple explicit formulation because we need to solve the entire system of algebraic equations to consider all dependencies (this is why it is an *implicit* method). \n",
    "\n",
    "However, we now obtain an equation that is in fact very similar to the Laplace problem and it is directly clear that we can solve the equation with the methods that we discussed previously for solving the Laplace equation!\n",
    "\n",
    "**Your task**: Solve the problem from above with the implicit method!\n",
    "\n",
    "Hint:\n",
    "\n",
    "1. write the equation in a form that you can solve\n",
    "2. use your favourite method to solve it\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# your code:\n",
    "# your code\n",
    "L = 1.\n",
    "kappa = 9e-6\n",
    "dt = 25\n",
    "nt = 200\n",
    "nx= 100\n",
    "umax = 10\n",
    "sigma = 0.05\n",
    "dx = L/nx\n",
    "\n",
    "beta = (kappa*dt) / (dx**2)\n",
    "\n",
    "##########################################\n",
    "\n",
    "def tridiag(nx,beta):\n",
    "    B = np.zeros(shape=(nx,nx))\n",
    "    i,j = np.indices(B.shape)\n",
    "    B[i==j] = 1+2*beta\n",
    "    B[i==j-1] = -beta\n",
    "    B[i==j+1] = -beta\n",
    "    return  B\n",
    "\n",
    "\n",
    "#########################################\n",
    "\n",
    "u = np.ones((nx,nx))\n",
    "\n",
    "\n",
    "xvals = np.linspace(0,1,100)\n",
    "u =  umax * np.exp(-((xvals - 0.5)/0.05)**2)\n",
    "\n",
    "b = np.zeros(nx)\n",
    "b0=0\n",
    "b1=0\n",
    "b[0] = b0\n",
    "b[-1] = b1\n",
    "\n",
    "B = tridiag(nx,beta)\n",
    "\n",
    "for i in range(nt):\n",
    "    u = np.dot(np.linalg.inv(B),u) + np.dot(np.linalg.inv(B),np.dot(beta,b))\n",
    "    #plt.plot(u)\n",
    "    if np.mod(i,10) == 0: plt.plot(xvals,u,color=plt.cm.copper_r(i/200.))\n",
    "\n",
    "print beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "u = np.dot(np.linalg.inv(B),u) + np.dot(np.linalg.inv(B),np.dot(beta,b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional questions and tasks\n",
    "\n",
    "1. Which solution is faster, explicit or implicit? Can you think of ways to speed-up the solution?\n",
    "2. We discussed the fact that the implicit solution is unconditionally stable! Convince yourself that this is the case!\n",
    "3. Implement the Crank-Nicolson method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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
}