homework4/notebooks/HW4-applied-transport-phenomena.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8097202d",
   "metadata": {},
   "source": [
    "# HW4: Applied Transport Phenomena  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9626a33b",
   "metadata": {},
   "source": [
    "What is our goal? \n",
    "\n",
    "Solve various forms of the convection-diffusion equation. Possible even add a reaction term; \n",
    "\n",
    "We expect to go through the following stages: \n",
    "1. extend 1D (only x) diffusion equation (HW1) to the convection-diffusion equation; \n",
    "2. extend linear convection (prescribed(u) times ux) to non-linear convection-diffusion (u times ux); \n",
    "3. extend from 1D (only x) to 2D (both x and y);\n",
    "4. extend from stationary (time independent) to transient (time dependent); \n",
    "5. extend from single component (only u(x) or u(x,y)) to (both u(x,y) and v(x,y)); (potential flow, Stokes, laminar Navier-Stokes); "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48616c9f",
   "metadata": {},
   "source": [
    "## Import Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c6f395c0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Succesfully imported numpy -- Version: 1.20.1\n",
      "Succesfully imported scipy -- Version: 1.6.2\n",
      "Succesfully imported matplotlib.pyplot\n",
      "Succesfully imported pandas -- Version: 1.2.4\n",
      "Succesfully imported sympy -- Version: 1.8\n",
      "Succesfully imported scipy.optimize\n",
      "Succesfully imported fsolve\n",
      "Succesfully imported find_peaks\n",
      "Succesfully imported orth\n",
      "Succesfully imported odeint\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "print(\"Succesfully imported %s -- Version: %s\"%(np.__name__,np.__version__))\n",
    "import scipy\n",
    "print(\"Succesfully imported %s -- Version: %s\"%(scipy.__name__,scipy.__version__))\n",
    "import matplotlib.pyplot as plt\n",
    "print(\"Succesfully imported %s\"%plt.__name__)\n",
    "import pandas as pd\n",
    "print(\"Succesfully imported %s -- Version: %s\"%(pd.__name__,pd.__version__))\n",
    "import sympy as sym \n",
    "print(\"Succesfully imported %s -- Version: %s\"%(sym.__name__,sym.__version__))\n",
    "from scipy import optimize\n",
    "print(\"Succesfully imported %s\"%optimize.__name__)\n",
    "from scipy.optimize import fsolve\n",
    "print(\"Succesfully imported %s\"%fsolve.__name__)\n",
    "from scipy.signal import find_peaks\n",
    "print(\"Succesfully imported %s\"%find_peaks.__name__)\n",
    "from scipy.linalg import orth\n",
    "print(\"Succesfully imported %s\"%orth.__name__)\n",
    "from scipy.integrate import odeint\n",
    "print(\"Succesfully imported %s\"%odeint.__name__)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84348980",
   "metadata": {},
   "source": [
    "## Section 1/: Introduction and Model Description"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39d86584",
   "metadata": {},
   "source": [
    "Describe the Brusselator model as ODE (see wiki), a PDE in one spatial dimensions [see here](https://rveltz.github.io/BifurcationKit.jl/dev/tutorials/tutorials3/#d-Brusselator-(automatic)) (possiby better references do exist) and as PDE in two spatial dimensions (see above, and book by Hairer and Wannes Part 2).  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0d94535",
   "metadata": {},
   "source": [
    "## Section 2/: Convection-Diffusion Equation in 1D"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ee62f53",
   "metadata": {},
   "source": [
    "### Section 2.1/: Problem statement \n",
    "\n",
    "We aim at solving the convection-diffusion equation on the interval $0 \\leq x \\leq 1$. Our objective is to numerically approximate the function $u(x)$ that is the solution of the following problem.\n",
    "\n",
    "\\begin{equation}\n",
    " -ε\\frac{\\text{d}^2u}{\\text{d}x^2}+\\frac{\\text{d}u}{\\text{d}x} = 0 \\tag{1}\n",
    "\\end{equation}\n",
    "\n",
    "where the first term will be referred to as the \"diffusion term\" and the second one will be referred to as the \"convective term\". \n",
    "\n",
    "Given the small positive number  $ε > 0$, find the function $u(x)$ that is the solution of the differential equation written above, supplied with the following non-homogeneous Dirichlet boundary condition in $x = 0$ and homogeneous Dirichlet boundary condition in $x = 1$:\n",
    "\n",
    "\\begin{equation}\n",
    " u(x=0)=1  \\text{          and         }  u(x=1)=0 \\tag{2}\n",
    "\\end{equation}\n",
    "\n",
    "Compared to the tutorial of this week, we now added a convective term and set the source term equal zero. The latter was done to simplify matters for the time being. The number $ε$ is the relative weight of the diffusion\n",
    "and the convective term in the differential equation. The number $ε$ plays an important role in the practice of a computational chemical engineer. It will reappear in later courses under the aliases such as the Peclet or\n",
    "Reynolds number. The solution $u(x)$ clearly depends on the value of $ε$.\n",
    "\n",
    "#### Analytical Reference Solution \n",
    "The differential equation $(1)$ supplied with the boundary conditions $(2)$\n",
    "can be solved analytically by first supposing that u(x) = exp(r x), solving the characteristic equation $-ε r^2 +r = 0$ for $r$, setting $u(x) = C_{1} exp(r_{1} x)+C_{2} exp(r_{2} x)$ and using the boundary conditions to find\n",
    "$C_{1}$ and $C_{2}$. Doing so, one finds the solution:\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    " u(x)=\\frac{\\exp (x/ε)-\\exp (1/ε)}{1-\\exp (1/ε)} \\tag{3}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c9f7b00",
   "metadata": {},
   "source": [
    "#### Assignment 1\n",
    "\n",
    "First, verify that the function $u(x)$ given by $(3)$ indeed satisfies the differential equation $(1)$ and the boundary conditions $(2)$. Please do so using pen and paper (touch screen?). The use of a symbolic\n",
    "manipulation toolbox is too elaborate for this purpose. Next, plot the function $u(x)$ as a function of $x$ for $0 < x < 1$ by avoiding for-loops in Python. Compare different values for $ε$ using for instance $ε = 1/10$,\n",
    "$ε = 1/50$ and $ε = 1/100$. Give a physical interpretation for the way in which $u(x)$ changes with $ε$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "375ec786",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U(x) changes into a sharper step function with decreasing e-values.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA2eklEQVR4nO3deZgU5bnw/+89vcy+z4DAsMouomFVg0ZjUNEoZgU1atTEJTHHeJI35uScaEzMq74/sxwTlxiToIlGjRqDiopxw2gMoEEQUEAQGdZhBpiNWXr6/v1RNUPP3gNd3dPT9+e65uruqqeq7uru6bue56l6SlQVY4wxqSst0QEYY4xJLEsExhiT4iwRGGNMirNEYIwxKc4SgTHGpDhLBMYYk+IsESQBERklIioi/sNc/iIRWRrruNx1XyMiu0WkVkSKvdhGh+2li8g6ETkqirL/ISK39VLmkyKy0Y3//JgFemj9XxWRf0S8rhWRMe7zTBF5WkQOiMhf3Gm3iMheEdkV61gSSUReFZGvJTqOViIyQUT+LSI1IvIfiY4n0SwRxJGIzBGRN91//CoReUNEZsZ4G52Shqo+pKpnxHI77rYCwM+BM1Q1R1UrY72NLlwJLFPVaH4o7wO+IiKDeijzY+DXbvxPxSLAnrjb2ey+/CIwGChW1S+JyHDgO8BkVe010cWa+70ZG+/tJsj3gFdVNVdV70x0MIlmiSBORCQPeAb4FVAEDANuBhoTGdcRGgxkAGv7uqA4Duf7dxXwx2gKqmoD8BxwSQ/FRnIY8QMcbg2tw7Y3qGoo4nWlqu45jFgO9/1MVYf9uQ9Iqmp/cfgDZgD7e5ifBvwPsBXYAzwI5LvzRgEK+N3XHwGfiVj2R8Cf3Ocfu2Vr3b8Tga8C/4gofxKwAjjgPp4UMe9V4CfAG0ANsBQo6SLe8UBdxLZejnLdP3XXfRAY28V6NXI6sAi4xX0+wl2u9X0IAquAb7mvfe66b4xY/iLglW7e8w+BsLvOWiAdGAosBqqATcDXO7zPjwN/AqqBr3WxzmJ3+Wpgufte/qPj/uEcBDQBze62r3LjCLuvF7nlTwDeBPYD7wKn9vR+AhOBF934PwC+3OG9vAt41v1s/wUc7c5b5sZW525/QYf9SndjmBIxrdTd7iCgEOdApwLY5z4v6xDr1zp+X7v5fucDvwN2AtuBWwCfO28s8BrO92sv8GgP/1Pn4fzY73e3P8md/jLQAjS4+zq+i2W7jWEg/iU8gFT5A/KASuABYB5Q2GH+5Tg/PGOAHOBJ4I/uvI7/KB/RfSJoV9ad9lXcHyOc2sg+4GLAD1zgvi5257+K8wM5Hsh0X9/WzT51jCuadX8MHOPOD3Sxzp4SwTnA2g7lp7jbmAT8N/BW5D8sMA2o6uFz6fhevgbcjVPTOR7nh+30iPe5GTgfJ3FndrG+R4DHgGw3tu10kQg6fm7u61OB8ojXw3C+M2e725vrvi7t5v3MB7YBl7mvp+H8WB4T8V5WAbPc+Q8Bj3T33nexb78Hfhrx+pvA8+7zYuALQBaQC/wFeCqi7KtEnwieAn7jvoeDcBLqVe68P7ufc5r7Gc3pJtbWA5W5QACnKWgTEOwYTzfLdxvDQPyzqmScqGo1MAfnC/9boEJEFovIYLfIRcDPVXWzqtYC/wUsjEHzQ0fnABtV9Y+qGlLVPwPvA+dGlPmDqm5Q1YM4P2rHx3Ddi1R1rTu/uY+xF+AcybZR1fdwjtb+CnwXuFhVWyKK1OD8QPbKbaOfA9ygqg2qugq4Hyextfqnqj6lqmH3/Ylc3ofzY3ijqta5sT3Qh/3r6CvAElVd4m7vRWAlTmJo1fZ+AmcBH6nqH9z39x3gCZy+iFZPqupyt/xDRP/ZAjyMk9xbXehOQ1UrVfUJVa1X1Rqcmsqn+ra74P4/zAO+7b6He4BfAAvdIs04zTpD3c/oH92sagHwrKq+6H7P7sA5sDkpBjEMOJYI4khV16vqV1W1DOdocSjwS3f2UJxmoVZbcY7aBhNbHbfTuq1hEa8jO2LrcWoosVr3tijX1ZV9OEebHT2Ac1S5RFU3dpiXi9OMEI2hOLWHyGTTl/hLcT6zyDId34++GAl8SUT2t/7hJKoh3cQzEpjdofxFQGTH8+F+tuA0qWSKyGwRGYmTRP4KICJZIvIbEdkqItU4TU0FbnLsi5E4R/A7I/bhNzhH5eAc2QuwXETWisjl3ayn3XdRVcM479Wwbsr3JYYBJ9ZHmyZKqvq+iCzCaRsG2IHzBWw1AggBu4GyDovX4VTBW0X+o/c2nGzH7bRu6/neo+5VNOvuLb56Ou9buft8NTBGRPx6qIMVnKacZ4AzRWROh6PESTht69HGXyQiuRHJYARO80408VfgfGbDcWpCrcsfrm04zYNf76FMZDzbgNdUde4RbLP7DamGReQxnFrBbuCZiPfpO8AEYLaq7hKR44F/4/xod9TT93cbzgkUJR0+49YYdgFfB+csPODvIrJMVTd1KLoDOLb1hYgIzueynd71GMNAZDWCOBGRiSLyHREpc18Px/mHesst8mfgehEZLSI5wP/F6Qjr6ou4CqfZKCAiM2hf9a/A6XAc000oS4DxInKhiPhFZAEwGeeH9EjFYt2rgAtFxCciZxHRvKCq5cBGnDZuAETkYmA6Tj/IfwAPuO9fq0/hnDnUK1XdhtMxe6uIZIjIVOAKnCaUaJZvwenb+ZF7hDwZuDSaZbvxJ+BcETnTfT8yROTU1u9QF57Bef8vdr8bARGZKSKTotzebrr/3rR6GKfZ5SL3eatcnI7j/SJSBNzUwzpWAaeIyAgRycdpBgVAVXfinKDwMxHJE5E0ETlaRD4FICJfitj/fTiJsIXOHgPOEZHT3dOcv4Pz4/5mL/vXawwDkSWC+KkBZgP/EpE6nATwHs4XFJyOuD/iVKm34JzR8K1u1vVD4Gicf4SbifiHVNV63DNJ3GrtCZELqnOu/2fd7VbiVLU/q6p7j3QHY7Tu63D6FPbj/Ng81WH+b3Db7EVkBE7T2iWqWquqD+O0of/CnZ+B057el3b6C3CamXbgNHvc5LbNR+tanOaWXTids3/ow7LtuIlpPvADnAS/Dfg/dPN/6x6dn4HTlr3DjeF2nDN+ovEjnES6X0S+3M02/oVzRD+U9gn2lzht8Htxvtvd1jDd9/NRnBre23Q+ULgE54ywdTjf8cc51Bw2E+d/qBbn7KzrVHVLF9v4AKeP5VduTOcC56pqU3dx9SGGAUdUe6upG9N/iEg6TpPD6e6RW09lvwUMV9XvxSU4Y5KUJQJjjElx1jRkjDEpzhKBMcakOEsExhiT4pLuOoKSkhIdNWpUosMwxpik8vbbb+9V1dKu5iVdIhg1ahQrV65MdBjGGJNURKTbq9ytacgYY1KcJQJjjElxlgiMMSbFWSIwxpgUZ4nAGGNSnGeJQER+LyJ7ROS9buaLiNwpIptEZLWITPMqFmOMMd3zskawCOeOSd2ZB4xz/64E7vEwFmOMMd3w7DoCVV0mIqN6KDIfeFCdUe/eEpECERnS24iSh+vVlX/lhdU9jAgsnZ4gXd5To6tlupvUzSvpPLXnLUXE1OM2O8/sqnx38XWKRror1X6JtIgdEzlUXtwApMO6xC3XqUyn6bStT3CeRMbnzHOWSRMQkbbX0vOOd7PfPU/pPDfKbUgU36UuF5OOE/pWvvegei/T6+Z7eZ/6tNvR7EEUMUWzUfGBLwBp/kOPaQHnuT8dAlngz+jzZ+cTHznBHHKDueQF88hPz6csp4yALxD1OhIhkReUDaP9bfbK3WmdEoGIXIlTa2DEiMO74dP6bct5Nm1z1OW1jz8kcdPXwWJtcFljEsonPkbljeKkYSdx+ojTmTZoWp8PVLzm6TDUbo3gGVWd0sW8Z4FbW28rKCIvAd9T1bd7WueMGTM01lcWR74HrU+1i/ntp7WW67xs64vI163lIqeFNdx5+x3XFT40L0y4/bx25cOdpoXb4m4NVttmdrsvGu4wLbJcxD64z8Pudpw/Z3pYoSWshMNhVNvPV1XCOPNVFQ1r2/xDZaFFnflhBQ2HO80Pq9KiSkiVUEuYUIvSEg7T0qI0hcO0tIRpCSvNLUqoRQnpoTLNLUpL+NDzULjFmdYSpimsNDW30BgK0xQK0xBqocl9HkmiyLB+n5ATTCMz3U920E9W0EdWup/sQBo56QHyMgPkZ/rJz/STlxmkIDNAflaAvAw/eZkBgmmRPxY9b08jP7QoaG/xd1hX1+V7iSniexJdRH2Lqevtd7+OtsVVQUMQaoRQQ4fHRmiqgYP7oL4Sra+Cmt2w/yM4sB3UvRlaVjGMOQ3GzXUefc4xdSgcora5lpqmGmqaaqhqqGLLgS2sq1zH8l3LaQ43M7VkKjeeeCMTiiZE88bEjIi8raozupqXyBpBOc49RFuV4dxVKe4is3PXibp/ZW8Tf6pKYyhMY3OYxlALDR0e65taqGsMUdsYoq4xRF1TCzUN7vPW6U0hKhpb+Kg6RE1DM/vqGzolmEg56X4KsgIMyk1nUG4Gg/PSGZSX4bzOc14Pzs2gICvQ744wB6RQE+zdANtXwpbXYcMLsPoxyBsGM6+A2ddARiGldDmcD3XNdTy35TnuXnU3Fzx7AbedfBtnjDojzjvRtUTWCM7Bua3f2Ti3cLxTVWd1LNeRFzUCYxJBVTnY3MK++mb21zexv76ZffVNzus653FffRN7ahrYU93I7uoGqhs638I6K+ijrDCT4YVZDC/Koqwwk7LCLIYXZTK6JJusYNINKZYcWpph41JYfh9sfhXyh8OZP4XJ83tcbF/DPq575TpWV6zmns/cw4lDT4xLuD3VCDxLBCLyZ+BUoATnptg3AQEAVb1XnEOYX+OcWVQPXKaqvf7CWyIwqayhuYU91Y3sqWlgt5scyvcdZNu+esr3HaS8qp6axvbJYlhBJmMH5TBuUA5jI/4KsoIJ2osBaOubsOR7sHsNHP8VOOdnEMjotnhdcx0XPnsh1U3VPHnekxRmFHoeYkISgVcsERjTPVXlwMFmtlUd5OOqejZX1LKpopaNu2v5sKKWxoimqCH5GUwZls/UYflMKcvn2GH5lOREe59700lLCF67HZb9Pxg5By58BNJzuy3+QdUHLHhmARdMvIAbZt3geXiWCIwxtISV7fsOsqmiho27a1m3s5o15QfYvLeurczQ/AyOH1HArFFFzBpdzISjcvGlWf9Dn6z+C/z1Khh5Elz0eI81gxvfuJFnNz/Ls59/lqOyj/I0LEsExphu1TQ0s3aHkxRWbz/AO1v3sX3/QQByM/zMHFXEiWOK+dSEUsYNyrGO6Wi8+yj89Ur4xMVw3q+6vbZhe+125j0xj6uOu4pvHv9NT0OyRGCM6ZPyffWs+KiK5Vuq+NfmqrZaw7CCTE4ZX8qpE0r55NgSctKtI7pbL/0EXr8DPv9bmPrlbotdufRKPq75mCWfX0KaeDfYgyUCY8wR2b7/IK99UMGrH+zhjU17qWtqIehP47QJpZwzdSinTxxEtiWF9sIt8Id5UPE+fHMF5A7ustjTHz7ND/7xAxadtYjpg6d7Fk5/vY7AGJMkhhVkcuHsEVw4ewRNoTArt1axdO1ulqzZyQtrd5MRSOPTEwfx2alD+cykwQT9NrAxaT44/x64aza89GM4/64ui316xKfxi5/Xy1/3NBH0xBKBMaZPgv40Tjq6hJOOLuHGz05mxUdVPLtmJ0vW7GLJml0UZQf5wrRhLJg5nLGDuj9rJiUUHw0nXANv/gpmXwVDpnYqkh3IZmrpVN7a+VYCAnRY2jbGHLa0NGH2mGJ+PH8K//rB6Sy6bCazRxfxhzc+4jM/X8aX7n2T59bspCWcXE3QMXXyd5zTSF+/o9siJww5gXWV6zjQeCCOgR1iicAYExO+NOHUCYO45yvT+ed/nc5/zZvIruoGrnnoHU694xX+8MYW6ho7Xxk94GUWwKwrYd1iqNjQZZHZQ2ajKCt3Jab/0xKBMSbmSnPTuepTR/Pqd0/j3q9MY1BuBjc/vY4Tb32JX720kdpUSwizr3aGuF5xf5ezp5RMwS9+3qvs8j5enrNEYIzxjC9NOGvKEJ645iSe/MZJzBpdzM9e3MDJt7/Mva99SH1TiiSEnFJnDKJ3/wxNdZ1mB31BxhSM4f2q9xMQnCUCY0ycTBtRyP2XzuBv3/wkU8sKuO259znl/73Koys+JpwKfQgzLofGalj/dJezJxZN5IOqD+IclMMSgTEmro4bXsADl8/i8atPZGRxFjc8sYbz736Dt7fuS3Ro3hpxIuSPgDWPdzl7QuEEKg5WsPfg3jgHZonAGJMgM0YV8fjVJ/LLBcez60ADX7jnTf7z0VXsrW1MdGjeEIEpn4fNr0BdZafZE4smAiSkVmCJwBiTMCLC+Z8YxsvfPZVrTj2ap1fvYO7PX+Nvq7aTbKMeRGXK5yEcgg+WdJo1rnAcAJv2b4p3VJYIjDGJl5Pu54azJvLsf5zMiOJsrntkFV9/8G0qagZY7eCoqc4dzTa+0GlWQXoBuYFcymvK4x6WJQJjTL8xfnAuT15zEv999iRe31jBvP9dxmsbKhIdVuyIOPc5/vAV59aX7WYJZbllbKvdFvewLBEYY/oVX5rw9VPGsPjaORRnp3Pp75fzf5es7/H+zkll3JnQVAsfv9lpVlluGdtrtsc9JEsExph+acJRufzt2k9y8QkjuW/ZZi787VvsqWlIdFhHbtQcEB9seb3TrOG5wymvLacl3BLXkCwRGGP6rYyAj5+cP4VfXfAJ3ttxgPN+9Qartu1PdFhHJiMPhh4PW9/oNKsst4xQOMSe+j1xDckSgTGm3zv3uKE8ec0n8fuEL//mn/xtVfybT2Jq1BwoXwlN9e0mD88dDkB5bXw7jC0RGGOSwuSheTx97Rw+MbyA6x5ZxW9e+zB5TzEdOQfCzVC+ot3kspwyALbVxLfD2BKBMSZpFGYHefCKWXx26hBufe59bn56XXIOcT18pvO4vf1oo4OznLuY7a7fHddw7MY0xpikku73cefCTzAkP4Pfvr6Fqromfv7l4/D7kui4NrMQisfC9nfaTQ74AhSmF7K3Pr7DTFgiMMYknbQ04b/PmUxRdjq3P/8+LWHllwuPJ5BMyWDYdNiyrNPk4sziuI83lETvmjHGtHfNqUfzP+dM4tk1O7n24XeS61qDYdOhZidU72g3uTSz1BKBMcb0xddOHsNN507mhbW7ue6RfydPn8HQac7j9rfbTS7NKqXiYHyvprZEYIxJepd9cjT/c84knntvF//z1JrkOJto8DGAwO617Sa3Ng3Fcx+sj8AYMyB87eQx7Ktv4q5XPqQoO8j/OXNiokPqWTALisZ0SgSlmaU0h5s50HiAgoyCuIRiicAYM2B894wJVNU1c9crH1Kcnc7lc0YnOqSeDT6my0QAsPfg3rglAmsaMsYMGCLCLedP4cxjBnPLs+t45YP4DtXQZ4OPgarN7e5jXJxZDBDXfgJLBMaYAcWXJvxiwfFMPCqPbz38bzbsrkl0SN0bfAygUHHopvWRNYJ48TQRiMhZIvKBiGwSke93MT9fRJ4WkXdFZK2IXOZlPMaY1JAV9HP/pTPICPi44oEVVNU19b5QIgya7DzuWd82qSSzBIDKg51vZ+kVzxKBiPiAu4B5wGTgAhGZ3KHYN4F1qnoccCrwMxEJehWTMSZ1DC3I5LeXTGd3dSPXPvxO/zyttGAkpAVg78a2SdmBbHzi40DTgbiF4WWNYBawSVU3q2oT8Agwv0MZBXJFRIAcoAoIeRiTMSaFfGJEIbfMn8KbH1Zy50sbe18g3nx+58yhykP3KRYR8tPzqW6sjlsYXiaCYUDkEHrl7rRIvwYmATuANcB1qtrp0kARuVJEVorIyoqKAXTbOmOM5740o4zPTxvGnS9v5B8b43vFblSKx7ZLBAB5wbwBUyOQLqZ1rJudCawChgLHA78WkbxOC6nep6ozVHVGaWlprOM0xgxgrWcSjS3N4duP/ps91f3sLmclY50zhyLuSpaXnseBxoGRCMqB4RGvy3CO/CNdBjypjk3AFqCfXwVijEk2WUE/d180jbrGFr796CrC/am/oHgstDTB/o/bJuUH86luGhhNQyuAcSIy2u0AXggs7lDmY+B0ABEZDEwANnsYkzEmRY0bnMuN507mzQ8r+eNbWxMdziHF45zHyg/bJg2YGoGqhoBrgReA9cBjqrpWRK4WkavdYj8BThKRNcBLwA2q2g8b8YwxA8HCmcM5dUIptz63ns0VtYkOx1F8tPNYdSgRDKQaAaq6RFXHq+rRqvpTd9q9qnqv+3yHqp6hqseq6hRV/ZOX8RhjUpuIcPsXppLu9/Gdv7xLqKUfDFudXQr+zPZNQ+n51DTV0BLRb+Alu7LYGJNSBudl8OP5x/Dvj/dz3+v9oCVaBApGwP5DzVV5QeecmZqm+FwVbYnAGJNyzjtuKPOmHMUv/76Rj/bW9b6A1wpHwr5DiSA/PR8gbs1DlgiMMSlHRPjReccQ9KXxw7+9l/j7F3SoEbQmgnh1GFsiMMakpMF5GXznjPG8vnEvz67ZmdhgCkZCwwE4uB841DRkNQJjjPHYxSeMZMqwPH789DpqGpoTF0jhSOfR7TDOS3cSgdUIjDHGY35fGj89/1gqahv52dINiQukYITz6CaC/KDbNBSnYSYsERhjUtpxwwu4aPYI/vjWVjbtSdC9C/LdRHDAGZ4tN5gL2FlDxhgTN9d/ZjxZAR+3Lnm/98JeyCoCXzpUO6PwBH1B/Gl+6pvr47J5SwTGmJRXnJPON04by0vv7+HNTQkY3EAE8oa2JQJw7ktQ1xyfU1stERhjDHDZJ0cxrCCTny5Zn5hB6fKGQs2hs5ey/FnUh6xGYIwxcZMR8PG9syawdkc1T/57e/wDyBsK1Ye2mx3ItqYhY4yJt/OOG8pxwwv42dIPaAzFZ5yfNrlDoHonuBe3ZQWyrGnIGGPiTUT43pkT2HmggUdXbOt9gVjKGwYtjVBfBUC2P9uahowxJhFOOrqYWaOKuOuVTTQ0x7FWkDfEeXSbh6xGYIwxCSIifHvuOHZXN/Ln5R/3vkCs5Lm3dHfPHLI+AmOMSaCTji5h9ugi7n71w/jVCnLdGkGNkwjsrCFjjEmw6+eOp6KmkYf+FadaQc4g57F2D2BNQ8YYk3AnjCnmxDHF3PtanGoFvgBkFrUlguxANs3hZppbvB8MzxKBMcZ049pPj6WippGn4nVdQc4gqDuUCIC4NA9ZIjDGmG6cdHQxxwzN477XN8fnauPs0kNNQ/4sgLg0D1kiMMaYbogIV54yhs0Vdbz8/h7vN5gzuF0fARCXM4csERhjTA/OPnYIwwoyuW9ZHG50nzMI6iqAQ01DdSGrERhjTEIFfGlcPmc0yz+q4t8f7/N2Y9ml0FQLTXXWNGSMMf3JwpnDycvwe18ryBnsPNbuaasRHGw+6O02sURgjDG9yk73c9EJI3lh7S7K93nYZt96LUFdRVsfgTUNGWNMP/GVE5wbzHs67ER2qfNYu9uahowxpr8ZVpDJpycO5tEV22gKhb3ZSMTVxW2dxZYIjDGm/7j4xJHsrW3i+bW7vNlAVonzWF9Fui8dQWhsafRmWxEsERhjTJROHlvCyOIs/vTPrd5swB+EYC7UVyIipPvSaQxZIjDGmH4jLU24aPYIln9UxQe7arzZSFYR1FcCkO5Pp6GlwZvtRPA0EYjIWSLygYhsEpHvd1PmVBFZJSJrReQ1L+Mxxpgj9aXpwwn60/jTWx7VCrKKDyUCXzoNoSROBCLiA+4C5gGTgQtEZHKHMgXA3cB5qnoM8CWv4jHGmFgozA7y2alDePKdcuoaQ7HfQFYxHHRuV5nhy0j6GsEsYJOqblbVJuARYH6HMhcCT6rqxwCqGofBPIwx5shcMGsEdU0tPP+eB53GkTUCf/L3EQwDIu/+XO5OizQeKBSRV0XkbRG5pKsViciVIrJSRFZWVFR4FK4xxkRnxshCRhZn8cQ75bFfeVZR2w3sM3wZSX/WkHQxreM4rn5gOnAOcCbwQxEZ32kh1ftUdYaqzigtLY19pMYY0wciwhemlfHmh5Wxv9I4q8gZbyjU6PQRJHnTUDkwPOJ1GbCjizLPq2qdqu4FlgHHeRiTMcbExOc+4TRw/PWdGN+0JqvYeayvGhBNQyuAcSIyWkSCwEJgcYcyfwNOFhG/iGQBs4H1HsZkjDExMbwoixPHFPP4O+WoxvCmNW2JoDL5O4tVNQRcC7yA8+P+mKquFZGrReRqt8x64HlgNbAcuF9V3/MqJmOMiaUvTC9ja2U9K7fGcHjqiESQ7ktP+j4CVHWJqo5X1aNV9afutHtV9d6IMv+fqk5W1Smq+ksv4zHGmFiaN+UosoI+nng7hp3GmUXOY30lGf6MpG8aMsaYAS073c+8KUN4ZvVOGppbYrPSDjWCpG4aMsaYVPC5TwyjtjHEK7G6p3FmofN4cP+AOH3UGGMGvBPGFFGcHeSZNTtjs0J/EAJZ0LDfOWuopTG2ndFdsERgjDFHwO9LY96xR/Hy+j3UN8VoyImMfCcR+NIBPK8VWCIwxpgjdM6xQznY3MLLsWoeyihoaxoCSwTGGNPvzRpdRGluOs+8G6PmocwCaDhAut+pEXg9AqklAmOMOUK+NOHsKUfxygd7qI3FiKRu01C/rBGISKGIHCMiY0TEkogxxrg+e9xQGkNhXlq/+8hXllHg1AjcPoKDoYNHvs4e9PpjLiL5IvIDEVkDvAX8BngM2CoifxGR0zyN0BhjksD0EYUMzkvnmdUxaB7KLICDB8jwx6dG4I+izOPAg8DJqro/coaITAcuFpExqvo7D+IzxpikkJYmnH3sEB5662OqG5rJywgc/soyCqDxABlpQaAfNA2p6lxV/WPHJODOe1tVv21JwBhj4Oxjh9DUEua1D47wvikZ+QCkh53+hn7TWSwiV3R47RORm2IfkjHGJKdpIwopyg7y9yPtJ8gsACCjuQnoBzWCCKeLyBIRGSIiU3D6C3I9issYY5KOL0349MRBvPL+Hppbwoe/oowCANLdmoDX4w1FnQhU9ULgAWANsAT4tqp+16vAjDEmGX1m0mCqG0Ks+Kjq8FfiNg1luCOPej0CaV+ahsYB1wFPAB/hdBJneRSXMcYkpVPGlxD0p/H3dUdwlbHbNJTe6NwGs9/UCICngR+q6lXAp4CNOHchM8YY48oK+pkztoQX1+86/MHiWpuGmp1E0J/6CGap6ksA6vgZcL4nURljTBL7zKTBbKs6yIbdtYe3gtamoQZn+YQ3DYnIHABVre44T1U3ikie23lsjDEGOH3SIIDDP3somA2SRlpzHX7x0xxujmF0nUVTI/iCiLwpIjeKyDkiMktEThGRy0Xkj8AzQKanURpjTBIZnJfBccMLeHHdYSYCEUjPhcYaAr4ATS1NsQ2wg16vLFbV60WkEPgi8CXgKOAgzg3p71XVNzyN0BhjktDcSYO4Y+kG9tQ0MCg3o+8rSM+DhmoCaYF+USNAVfcBQ4EtwD+BVUADcLpnkRljTBI7baLTPLRsw97DW0F6LjQ6iaAp7G2NoC+dxbURfyFgHjDKg5iMMSbpTToqj5KcdJZtOMzhJtLzoLGGoC9Ic4u3NYJoBp0DwD1LqI2I3AEsjnlExhgzAKSlCaeMK+HVDRWEw0pamvRtBem5UF9JIDerX9UIOsoCxsQqEGOMGWhOGV9KVV0Ta3d0Oumyd25ncdAXJBSO0b2QuxF1jcC9H0Hr1RE+oBT4sRdBGWPMQDBnXAkAyzZWcGxZft8Wbj1rKK0fnDUU4bMRz0PAblX1Nk0ZY0wSK8lJ55iheby2oYJvnja2bwtHdBb3i7OGAFR1a8TfdksCxhjTu1PGl/LO1n19v5dxeh401xNI83teI7D7DhtjjIdOHldCKKz888PKvi2YkQdAQNL6T43AGGNM380YWURW0Nf300jTndu9BBGrERhjTDIL+tM4cUwxyzYeXiIIgNUIjDEm2Z0yvpStlfVsrayLfqHWGoEmeSIQkbNE5AMR2SQi3++h3EwRaRGRL3oZjzHGJMInxzqnkfapnyDd7SNQ9fzKYs8SgYj4gLtwhqKYDFwgIpO7KXc78IJXsRhjTCIdXZpNaW46b23uSyJwm4bC4X59ZXFvZgGbVHWzqjYBjwDzuyj3LZzbXx7Bfd2MMab/EhFOGFPMW5uror9rWVuNoCWpm4aGAdsiXpe709qIyDDgc8C9Pa1IRK4UkZUisrKi4jAHcDLGmAQ6YUwRu6ob2FpZH90CrX0E4ZakPmuoqxGWOqbCXwI3qGpLTytS1ftUdYaqzigtLY1VfMYYEzcnjCkGiL55KJAFCIFwctcIyoHhEa/LgB0dyswAHhGRj3BufHO3iJzvYUzGGJMQY0r62E+QlgbBbIItLYTCIcIa9iy2vow11FcrgHEiMhrYDiwELowsoKqjW5+LyCLgGVV9ysOYjDEmITr2E4hEMSx1MJuAO/JoKBwi6At6EptnNQJ3LKJrcc4GWg88pqprReRqEbnaq+0aY0x/1ed+gmA2wRYnEXjZT+BljQBVXQIs6TCty45hVf2ql7EYY0yiRfYTjCrJ7n2BQDZ+9xoCL/sJ7MpiY4yJkz73EwSzCbo1AS9rBJYIjDEmTvp8PUEwm0DISQBWIzDGmAGiT/0EwWyCoUYAT68utkRgjDFxNHt0EQDLt1T1XjiYQ6DZSQRejjdkicAYY+JoTEkOBVkB3vl4X++Fg9kEQwcBaxoyxpgBIy1NmDaikLe3RpcI2moElgiMMWbgmD6ykI17ajlQ38uPezCHgJ01ZIwxA8+0EYUAvLOtl1pBMJuAe3aR1QiMMWYAOW54Pr404Z3emoeC2QRbE4F1FhtjzMCRFfQzaUhu7/0EwWwC7uUGdvqoMcYMMNNHFPLutv2EWnoYVTSYc6hGYE1DxhgzsEwbWUhdUwsf7K7pvlAwmwDWNGSMMQPS9JFuh3FPzUORfQRWIzDGmIFlWEEmg/PSe+4nCOYc6iOw00eNMWZgERGmjyzk7Z6uMLbTR40xZmCbNqKQbVUH2VPd0HWBYFZbIrAagTHGDEDTWvsJuqsVBHPwA4KdPmqMMQPSMUPzCPiEd8sPdF3AF0DSAvhJoyXc4lkclgiMMSZB0v0+JhyVy5ruEgFAIAs/Qsi9ib0XLBEYY0wCTS0rYHX5/u7vWBbIxC8QUksExhgzIE0dlk91Q6j7O5YFMvGr1QiMMWbAOrYsH4DV27tpHgpm4wNLBMYYM1CNH5xLuj+NNeX7uy4QyMSPJQJjjBmwAr40Jg/N6/7MoUAmflVa1M4aMsaYAWvqsHzWbj9AS7iLDuNAFn5VqxEYY8xANrWsgLqmFrbsre08060RWCIwxpgBbGprh3FXzUOBLHwattNHjTFmIBtTmkNW0NdNIsjEr2GrERhjzEDmSxOmDMtndVdnDgUy8YfDNsSEMcYMdFOH5bN2R3XnW1cGsvFpCyEPh6H2e7ZmQETOAv4X8AH3q+ptHeZfBNzgvqwFrlHVd/u6nebmZsrLy2lo6GYoV9NnGRkZlJWVEQgEEh2KMSnh2LJ8GkNhNuyuZfLQvEMzApkE1NthqD1LBCLiA+4C5gLlwAoRWayq6yKKbQE+par7RGQecB8wu6/bKi8vJzc3l1GjRiEisQg/pakqlZWVlJeXM3r06ESHY0xKOK6sAIA12/d3SATO6aMHk/TGNLOATaq6WVWbgEeA+ZEFVPVNVW0diPstoOxwNtTQ0EBxcbElgRgREYqLi62GZUwcjSjKIjvoY92O6vYzApn4SN6b1w8DtkW8LnendecK4LmuZojIlSKyUkRWVlRUdLmwJYHYsvfTmPhKSxMmDclj/c6a9jPariNIzhvTdPVL0uU4qyJyGk4iuKGr+ap6n6rOUNUZpaWlMQzRGGP6DycRVLcfktptGmpJ0tNHy4HhEa/LgB0dC4nIVOB+YL6qVnoYT1KprKzktNNOIycnh2uvvbbT/FtvvZWHHnqIZcuWMW3aNPx+P48//ni7Mg888ADjxo1j3LhxPPDAA/EK3RhzmCYNyaOmMUT5voOHJrpNQ15eR+DlWUMrgHEiMhrYDiwELowsICIjgCeBi1V1g4exJJ2MjAx+8pOf8N577/Hee+91mr906VIee+wx6urqWLRoEXfccUe7+VVVVdx8882sXLkSEWH69Omcd955FBYWxmsXjDF91NpJvG5nNcOLspyJwWwCqsl5+qiqhkTkWuAFnNNHf6+qa0Xkanf+vcCNQDFwt9smHVLVGUey3ZufXtu5s+UITR6ax03nHtNjmT/96U/ceeedNDU1MXv2bO6++258Pt9hbzM7O5s5c+awadOmTvOqq6tpamqitLSU1qaytLT2lbsXXniBuXPnUlRUBMDcuXN5/vnnueCCCw47JmOMtyYMziVNYN2Oas485ihnYluNwLsLyjy9jkBVlwBLOky7N+L514CveRlDPKxfv55HH32UN954g0AgwDe+8Q0eeughLrnkknblrr/+el555ZVOyy9cuJDvf//7UW/v73//O6effnqPZbZv387w4Yda5srKyti+fXvU2zDGxF9m0MfokmzW74w4mG0dfdTDYag9TQSJ0NuRuxdeeukl3n77bWbOnAnAwYMHGTRoUKdyv/jFL2Kyveeff57LLrusxzJd3f/UzgQypv+bNCSPdyOHmghk4lc8vR/BgEsEiaCqXHrppdx66609lotVjWD58uXcc889PZYpKyvj1VdfbXtdXl7OqaeeGvU2jDGJMWlIHs+s3kl1QzN5GQG3achqBP3e6aefzvz587n++usZNGgQVVVV1NTUMHLkyHblYlEjWLt2LRMnTuy1/+HMM8/kBz/4Afv2OdfrLV26tNdEZYxJvNYO4/d31jBrdBEEsggohDTcy5KHzxJBDEyePJlbbrmFM844g3A4TCAQ4K677uqUCPpq1KhRbR3DTz31FEuXLuW5557jrLPOaiuzYsUKPve5z7Fv3z6efvppbrrpJtauXUtRURE//OEP25qrbrzxxraOY2NM/zV5iHvm0I4DTiLwBfEDLZYI+r8FCxawYMGCmK7zo48+6jTtuuuu48EHH2x7PXPmTMrLy7tc/vLLL+fyyy+PaUzGGG8Nyk2nKDt46ApjEXxpfkIoqupJX58lgiTz4osvJjoEY4yHRITJQ/JYv+vQmUN+cZqCQxoiILEfEdjuR2CMMf3MpCG5vL+rpu3eBP4058ffq5vTWCIwxph+ZtKQPJpCYbbsrQPAn+Y03ng1zIQlAmOM6Wcih5qAQzUCSwTGGJMiji7NIeAT3t/ldBj7fG4iUEsExhiTEgK+NEaXZLNxt5MIrEaQoqIdhnrRokWUlpZy/PHHc/zxx3P//fe3lbFhqI1JXuMG57Jhdy0Afl8Q8G6YCTt9tJ+KdhjqZ599lgULFvDrX/+63XwbhtqY5DZuUA5L1uzkYFMLvrQAtHhXIxh4ieC578OuNbFd51HHwrzbeiySqGGou2PDUBuT3MYPzkUVPqyoJeAPepoIrGkoBiKHoV61ahU+n4+HHnqoU7nrr7++rQkn8u+223pOMh11HIb6iSeeYOrUqXzxi19k2zbnNtE2DLUxyW384BwANuyuwe9LB6xGEL1ejty9kMhhqM8991wuuOAC0tPTuffee7n00kt5+eWXbRhqY5LcyOJsAj5hw+5aZrt9BF6dNTTwEkECJHIY6uLi4rbpX//617nhhhsAG4bamGQX8KUxpiSHjbtr+GSh1Qj6vUQOQ71z506GDBkCwOLFi5k0aRJgw1AbMxCMHZzD6vL9+EszAO+GmLBEEAOJHIb6zjvvZPHixfj9foqKili0aBGADUNtzAAwflAuS9bsRMVtGgo1erIdSwQxkqhhqG+99dZuj/RtGGpjktv4wTmoQn2zO/poc70n27FEkGRsGGpjUse4wbkA7G9wTvQIhQ56sh07fdQYY/qpUcVZBH1pVFkiMMaY1OT3pTGmNJuKOicRtDQ3eLIdSwTGGNOPjR2Uw+5a57ogqxEYY0wKGj84lz11zk91KGQ1AmOMSTnjB+fQpO4FZZYIUku0w1AvW7aMadOm4ff7efzxx9uV6W4Y6i1btjB79mzGjRvHggULaGpq8nx/jDGHZ9zgXJq19ToCSwQppXUY6jvuuKPL+UuXLuWMM85gxIgRLFq0iAsvvLDd/NZhqP/1r3+xfPlybr755rarjG+44Qauv/56Nm7cSGFhIb/73e883x9jzOEZWZRFWDIBCLV4kwgG3HUEty+/nfer3o/pOicWTeSGWTf0WCZRw1C3DkWdltY+p3c3DPXChQt5+eWXefjhhwG49NJL+dGPfsQ111xz2LEaY7zj96VRnJ/Px0BLize1d6sRxECih6HuSnfDUFdWVlJQUIDf72833RjTfx3lHtDZEBNR6u3I3QuJHIa6O90NQ23DUxuTfMpKiqEKGpMxEYjIWcD/Aj7gflW9rcN8ceefDdQDX1XVd7yMyQuJHIa6O90NQ11SUsL+/fsJhUL4/X7Ky8sZOnRo1Ns2xsTfiEElUAV1Dd5cR+BZIhARH3AXMBcoB1aIyGJVXRdRbB4wzv2bDdzjPiaVRA5D3Z3uhqEWEU477TQef/xxFi5cyAMPPMD8+fOPOC5jjHdGH1UC70NNQ/KdNTQL2KSqm1W1CXgE6PiLMx94UB1vAQUiMsTDmDwROQz11KlTmTt3Ljt37jzi9Y4aNYr//M//ZNGiRZSVlbFu3bpOw1CvWLGCsrIy/vKXv3DVVVdxzDHHAO2HoZ45c2a7Yahvv/12fv7znzN27FgqKyu54oorjjhWY4x3jj6qCL8qB5uS76yhYcC2iNfldD7a76rMMKDdr6iIXAlcCTBixIiYBxoLiRqGeubMmZSXl3e5fHfDUI8ZM4bly5fHLE5jjLdy0v3MaMpjaOFYT9bvZSLoqgeyY09lNGVQ1fuA+wBmzJjRubczhdgw1Makpt9e+aZn6/ayaagcGB7xugzYcRhljDHGeMjLRLACGCcio0UkCCwEFncosxi4RBwnAAdU9bAa17s6LdIcPns/jUkdnjUNqWpIRK4FXsA5ffT3qrpWRK52598LLME5dXQTzumjPZ8c342MjAwqKyspLi62c+JjQFWprKwkIyMj0aEYY+JAku3Ib8aMGbpy5cp205qbmykvL6fBo1OrUlFGRgZlZWUEAoFEh2KMiQEReVtVZ3Q1b0BcWRwIBBg9enSiwzDGmKRkYw0ZY0yKs0RgjDEpzhKBMcakuKTrLBaRCmDrYS5eAuyNYTjJwPY5Ndg+p4Yj2eeRqlra1YykSwRHQkRWdtdrPlDZPqcG2+fU4NU+W9OQMcakOEsExhiT4lItEdyX6AASwPY5Ndg+pwZP9jml+giMMcZ0lmo1AmOMMR1YIjDGmBQ3IBOBiJwlIh+IyCYR6XRXeHfY6zvd+atFZFoi4oylKPb5IndfV4vImyJyXCLijKXe9jmi3EwRaRGRL8YzPi9Es88icqqIrBKRtSLyWrxjjLUovtv5IvK0iLzr7vNhjWLcX4jI70Vkj4i818382P9+qeqA+sMZ8vpDYAwQBN4FJncoczbwHM4d0k4A/pXouOOwzycBhe7zeamwzxHlXsYZ8vyLiY47Dp9zAbAOGOG+HpTouOOwzz8AbneflwJVQDDRsR/BPp8CTAPe62Z+zH+/BmKNYBawSVU3q2oT8Agwv0OZ+cCD6ngLKBCRIfEONIZ63WdVfVNV97kv38K5G1wyi+ZzBvgW8ASwJ57BeSSafb4QeFJVPwZQ1WTf72j2WYFccW5GkoOTCELxDTN2VHUZzj50J+a/XwMxEQwDtkW8Lnen9bVMMunr/lyBc0SRzHrdZxEZBnwOuDeOcXkpms95PFAoIq+KyNsickncovNGNPv8a2ASzm1u1wDXqWo4PuElRMx/vwbE/Qg66OoWZR3PkY2mTDKJen9E5DScRDDH04i8F80+/xK4QVVbBsid66LZZz8wHTgdyAT+KSJvqeoGr4PzSDT7fCawCvg0cDTwooi8rqrVHseWKDH//RqIiaAcGB7xugznSKGvZZJJVPsjIlOB+4F5qloZp9i8Es0+zwAecZNACXC2iIRU9am4RBh70X6396pqHVAnIsuA44BkTQTR7PNlwG3qNKBvEpEtwERgeXxCjLuY/34NxKahFcA4ERktIkFgIbC4Q5nFwCVu7/sJwAFV3RnvQGOo130WkRHAk8DFSXx0GKnXfVbV0ao6SlVHAY8D30jiJADRfbf/BpwsIn4RyQJmA+vjHGcsRbPPH+PUgBCRwcAEYHNco4yvmP9+DbgagaqGRORa4AWcMw5+r6prReRqd/69OGeQnA1sAupxjiiSVpT7fCNQDNztHiGHNIlHboxynweUaPZZVdeLyPPAaiAM3K+qXZ6GmAyi/Jx/AiwSkTU4zSY3qGrSDk8tIn8GTgVKRKQcuAkIgHe/XzbEhDHGpLiB2DRkjDGmDywRGGNMirNEYIwxKc4SgTHGpDhLBMYYk+IsERhjTIqzRGCMMSnOEoExR8i938FqEckQkWx3TPwpiY7LmGjZBWXGxICI3AJk4Az0Vq6qtyY4JGOiZonAmBhwx8FZATQAJ6lqS4JDMiZq1jRkTGwU4dwUJRenZmBM0rAagTExICKLce6eNRoYoqrXJjgkY6I24EYfNSbe3LuAhVT1YRHxAW+KyKdV9eVEx2ZMNKxGYIwxKc76CIwxJsVZIjDGmBRnicAYY1KcJQJjjElxlgiMMSbFWSIwxpgUZ4nAGGNS3P8PM1B9/xILsYwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First I define the function and the necessary variables\n",
    "def u(x,e):\n",
    "    F = (np.exp(x/e)-np.exp(1/e))/(1-np.exp(1/e))\n",
    "    return F\n",
    "e1 = 0.1\n",
    "e2 = 1/50\n",
    "e3 = 1/100\n",
    "x = np.linspace(0,1,300)\n",
    "\n",
    "# Then I plot the figures\n",
    "plt.figure()\n",
    "plt.plot(x, u(x,e1), label = 'e = 1/10')\n",
    "plt.plot(x, u(x,e2), label = 'e = 1/50')\n",
    "plt.plot(x, u(x,e3), label = 'e = 1/100')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('u(x)')\n",
    "plt.title('Solution for u(x) for different values of e')\n",
    "plt.legend(loc = 3)\n",
    "\n",
    "print('U(x) changes into a sharper step function with decreasing e-values.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e3e4d4a",
   "metadata": {},
   "source": [
    "### Section 1.2:/ Solve Symbolically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e80050bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "#..define symbol x (spatial coordinate), function f (source) and function u (unknown)\n",
    "x = sp.symbols('x')\n",
    "f = sp.Function(\"f\")\n",
    "u = sp.Function(\"u\") \n",
    "\n",
    "#..define the expression for the right-hand side function \n",
    "f = x**2 \n",
    "\n",
    "#..solve the differential equation - u''(x) = f(x) + BC symbolically\n",
    "computed = sp.dsolve(Derivative(u(x), x, x) - 100*Derivative(u(x), x),ics={u(0):0,u(1):1})\n",
    "\n",
    "#..print the expression for the symbolic solution \n",
    "print(\"The computed solution is\", computed.rhs)\n",
    "\n",
    "#..to switch from symbolic to numerical computation\n",
    "#..we need to perform a \"lambdification\" of the right-hand side of the solution\n",
    "#..the function created accepts vector arguments \n",
    "computed_sol = sp.lambdify(x, computed.rhs)\n",
    "\n",
    "#..sample the lambda expression for the computed solution \n",
    "nsamples = 100 \n",
    "xvec = np.linspace(0,1,nsamples)\n",
    "computed_sol_vec = computed_sol(xvec)\n",
    "\n",
    "#..plot the computed solution \n",
    "plt.plot(xvec,computed_sol_vec,label='computed u(x)')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('u(x)')\n",
    "plt.legend(loc='best')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc705146",
   "metadata": {},
   "source": [
    "### Section 1.3:/ Solve Numerically using in-build function \n",
    "Works for both the linear and non-linear convective term. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01860ec9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def righthandsidefunction(x, u):\n",
    "    return np.vstack((u[1], 100*u[0]*u[1]))\n",
    "\n",
    "def boundaryconditions(ua, ub):\n",
    "    return np.array([ua[0], ub[0]-1])\n",
    "\n",
    "nsamples = 100 \n",
    "xvec = np.linspace(0,1,nsamples)\n",
    "\n",
    "ustart = np.zeros((2, xvec.size))\n",
    "\n",
    "computed = scipy.integrate.solve_bvp(righthandsidefunction, boundaryconditions, xvec, ustart)\n",
    "\n",
    "plt.plot(xvec,computed.sol(xvec)[0],label='computed u(x)')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('u(x)')\n",
    "plt.legend(loc='best')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3be85752",
   "metadata": {},
   "source": [
    "### Section 1.4:/ Finite Difference Solution\n",
    "Next we proceed with the finite difference approximation.\n",
    "\n",
    "#### Mesh \n",
    "Assume that the interval $0 < x < 1$ is discretized by an uniform mesh consisting of $N$ elements with mesh width $h = 1=N$ and vertices $x_{i} = (i-1)h$, where $i$ runs from 1 to $N + 1$. This enumeration includes the end points of the interval, that is, $x_{1} = 0$ and $x_{N+1} = 1$. The grid nodes can then be denoted by\n",
    "\n",
    "\\begin{equation}\n",
    " G_h=\\Bigg({x_{i}|x_{i}=(i-1)h; h=\\frac{1}{N}, 1\\leq i \\leq N+1}\\Bigg) \\tag{4}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "Let $u^{h}$ denote the vector with $N+1$ components such its $i$-th component $u^{h}_{i}$ approximates the exact solution $u(x)$ in the $i$-th node.\n",
    "\n",
    "#### Discretization of the diffusion term in the interior nodes\n",
    "The differential equation holds in particular for all of the internal nodes, that is, we have that\n",
    "\n",
    "\\begin{equation}\n",
    " -ε\\frac{\\text{d}^2u}{\\text{d}x^2}\\Bigg|_{x=x_{i}}+\\frac{\\text{d}u}{\\text{d}x}\\Bigg|_{x=x_{i}} = 0 \\text{  for all  } 2\\leq i \\leq N    \\tag{5}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "To approximate the diffusion term, we proceed as in the tutorial and use a finite difference formula to discretize the second derivative $\\frac{\\text{d}^2u}{\\text{d}x^2}\\Bigg|_{x=x_{i}} = u''(x_{i})$ as follows\n",
    "\n",
    "\\begin{equation}\n",
    " u''(x_{i})\\approx\\frac{u'(x_{i}+h/2)-u'(x_{i}-h/2)}{h}    \\tag{6}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    " u''(x_{i})\\approx\\frac{u(x_{i+1})-2u(x_{i})+u(x_{i-1})}{h^2}    \n",
    "\\end{equation}\n",
    "\n",
    "(Using Taylor polynomials it can be shown that the local truncation error is of second order in $h$). The finite difference discretization thus leads to the following stencil for the approximation of $-ε\\frac{\\text{d}^2u}{\\text{d}x^2}\\Bigg|_{x=x_{i}}$ (beware of the minus-sign and the factor $-ε$)\n",
    "\n",
    "\\begin{equation}\n",
    " \\frac{1}{h^2}\\Bigg[\\begin{matrix}-ε &  2ε & -ε\\\\ x_{i-1} & x_{i} & x_{i+1} \\end{matrix}\\Bigg]       \\text{[central approximation for difussion]} \\tag{7}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "This stencil implies that each node $x_{i}$ is coupled to its left ($x_{i-1}$) and right neighbour ($x_{i+1}$) with a weight of $-\\frac{ε}{h^2}$.\n",
    "\n",
    "\n",
    "#### Discretization of the convection term in the interior nodes \n",
    "\n",
    "The convective term can be approximated in two ways. The first option is to use a central approximation on neighboring nodes\n",
    "\n",
    "\\begin{equation}\n",
    " u'(x_{i})\\approx\\frac{u(x_{i+1})-u(x_{i-1})}{2h}   \\tag{8} \n",
    "\\end{equation}\n",
    "\n",
    "leading to the stencil\n",
    "\n",
    "\\begin{equation}\n",
    " \\frac{1}{2h}\\Bigg[\\begin{matrix}-1 &  0 & 1\\\\ x_{i-1} & x_{i} & x_{i+1} \\end{matrix}\\Bigg]       \\text{[central approximation for convection]} \\tag{9}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "(The local approximation error of this scheme is of second order in $h$. This second order comes however at an expensive that we discover later in this assignment.) The second option is to approximate the first order\n",
    "derivative using the backward finite difference approximation\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    " u'(x_{i})\\approx\\frac{u(x_{i})-u(x_{i-1})}{h}   \\tag{10} \n",
    "\\end{equation}\n",
    "\n",
    "In this approximation information in $x_{i-1}$ is used to approximate the derivative in $x_{i}$. In this exercise the convective term propagates information from the left to the right. The node $x_{i-1}$  thus lies upstream from $x_{i}$. The above finite difference approximation is therefore referred to as the upwind approximation. This approximation leads to the stencil\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    " \\frac{1}{h}\\Bigg[\\begin{matrix}-1 & 1 & 0\\\\ x_{i-1} & x_{i} & x_{i+1} \\end{matrix}\\Bigg]       \\text{[upwind approximation for convection]} \\tag{11}\n",
    "\\end{equation}\n",
    "\n",
    "(The local approximation of this scheme is only of first order in h. Despite this lower order of accuracy, this scheme has a distinct advantage that we will discover later in this assignment.) The stencil of the diffusive and the convective term should be added to obtain the stencil for the coefficient matrix $A^h$ corresponding to the finite difference scheme.\n",
    "\n",
    "#### Discretization in the left and right boundary condition \n",
    "Verify that the Dirichlet boundary condition in the left and right end point can be enforced by requiring that\n",
    "\n",
    "\\begin{equation}\n",
    " u_{1}=1 \\text{          and that         }       u_{N+1}=0  \\tag{12}\n",
    "\\end{equation}\n",
    "\n",
    "The finite difference stencil in the left and right end point thus reduces to\n",
    "\n",
    "\\begin{equation}\n",
    " \\Bigg[\\begin{matrix}1 & 0 & 0\\\\ x_{1} & x_{2} & x_{3} \\end{matrix}\\Bigg]   \\text{          and         }    \\Bigg[\\begin{matrix}0 & 0 & 1\\\\ x_{N-1} & x_{N} & x_{N+1} \\end{matrix}\\Bigg] \\tag{13}\n",
    "\\end{equation}\n",
    "\n",
    "The non-homogeneous contribution in the left boundary condition (the term equal to 1) will be placed in the first component of the right-hand vector. The other components of this right-hand vector will remain equal to one."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da122b92",
   "metadata": {},
   "source": [
    "#### Assignment 2 \n",
    "\n",
    "Assume $h = 1/3$ (and thus $N=3$) and $ε = 0.5$. Assume that the central difference scheme $(8)$ is used to approximate the convective term in $(1)$. Determine the size of the global matrix $A^h$ and the global right-hand vector $f$. Give all the elements of this matrix and vector with pen (or pencil) and paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "02c22531",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The global matrix A is:\n",
      "[[ 1.  0.  0.  0.]\n",
      " [-6.  9. -3.  0.]\n",
      " [ 0. -6.  9. -3.]\n",
      " [ 0.  0.  0.  1.]]\n",
      "and its size is: 16\n",
      "\n",
      "The global right-hand vector is:\n",
      "[1. 0. 0. 0.]\n",
      "with size: 4\n"
     ]
    }
   ],
   "source": [
    "# First I define the necessary variables\n",
    "e = 1/2\n",
    "h = 1/3\n",
    "x = np.linspace(0, 1, 4)\n",
    "f = np.zeros(4)\n",
    "f[0] = 1\n",
    "\n",
    "# The convection term using central difference scheme\n",
    "c1 = 1/(2*h)*np.diag(np.ones(3), k = 1)\n",
    "c2 = -1/(2*h)*np.diag(np.ones(3), k = -1)\n",
    "\n",
    "# The diffusion term \n",
    "d1 = -e/(h**2)*np.diag(np.ones(3), k = 1)\n",
    "d2 = -e/(h**2)*np.diag(np.ones(3), k = -1)\n",
    "d3 = 2*e/(h**2)*np.diag(np.ones(4))\n",
    "\n",
    "# The global matrix\n",
    "A = c1 + c2 + d1 + d2 + d3\n",
    "A[0,0] = 1\n",
    "A[0,1] = 0\n",
    "A[-1,-1] = 1\n",
    "A[-1,-2] = 0\n",
    "\n",
    "# Lastly I print the results\n",
    "print('The global matrix A is:')\n",
    "print(A) \n",
    "print('and its size is:', np.size(A))\n",
    "print('')\n",
    "print('The global right-hand vector is:')\n",
    "print(f)\n",
    "print('with size:', np.size(f))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5ae9208",
   "metadata": {},
   "source": [
    "#### Assignment 3 \n",
    "\n",
    "As previous assignment, this time using the upwind scheme to approximate the convective term in $(1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a3147693",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The global matrix A is:\n",
      "[[ 1.   0.   0.   0. ]\n",
      " [-7.5 12.  -4.5  0. ]\n",
      " [ 0.  -7.5 12.  -4.5]\n",
      " [ 0.   0.   0.   1. ]]\n",
      "and its size is: 16\n",
      "\n",
      "The global right-hand vector is:\n",
      "[1. 0. 0. 0.]\n",
      "with size: 4\n"
     ]
    }
   ],
   "source": [
    "# First I define the necessary variables\n",
    "e = 1/2\n",
    "h = 1/3\n",
    "x = np.linspace(0, 1, 4)\n",
    "f = np.zeros(4)\n",
    "f[0] = 1\n",
    "\n",
    "# The convection term using upwind scheme\n",
    "c1 = 1/h*np.diag(np.ones(4), k = 0)\n",
    "c2 = -1/h*np.diag(np.ones(3), k = -1)\n",
    "\n",
    "# The diffusion term \n",
    "d1 = -e/(h**2)*np.diag(np.ones(3), k = 1)\n",
    "d2 = -e/(h**2)*np.diag(np.ones(3), k = -1)\n",
    "d3 = 2*e/(h**2)*np.diag(np.ones(4))\n",
    "\n",
    "# The global matrix\n",
    "A = c1 + c2 + d1 + d2 + d3\n",
    "A[0,0] = 1\n",
    "A[0,1] = 0\n",
    "A[-1,-1] = 1\n",
    "A[-1,-2] = 0\n",
    "\n",
    "# Print the results\n",
    "print('The global matrix A is:')\n",
    "print(A) \n",
    "print('and its size is:', np.size(A))\n",
    "print('')\n",
    "print('The global right-hand vector is:')\n",
    "print(f)\n",
    "print('with size:', np.size(f))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c69fc5e",
   "metadata": {},
   "source": [
    "#### Assignment 4\n",
    "\n",
    "Assume $ε = 1/100$ and that the central difference scheme $(8)$ is used to approximate the convective term in $(1)$. Assemble for $h = 1/10$, $h = 1/20$ and $h = 1/100$ the global matrix $A^h$ and the\n",
    "global right-hand vector $f^h$. Solve the linear system\n",
    "\n",
    "\\begin{equation}\n",
    " A^h u^h = f^h   \\tag{14} \n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "using the numpy.linalg.solve() function. Plot the various solution for $u^h$ found and compare this plot with the plot of the exact solution in Part 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a65bbb29",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The higher the N-value the more similar to the exact solution the approximation becomes\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAEWCAYAAABWqYxLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7LElEQVR4nO3deXxU9b3/8dcnk0wWWWQXTRAtuNReTDCEeisuWFxbbK3ValsVFH7eq3X50aq3i8ttvaWtS+u1XmRTe3vrbpX6c61LtXoVQom41ZYqSgTZVFADZPv8/jgnYRgmyYRkciaZ9/PxyCNz9s/3bPOZ7/mec8zdEREREYlSXtQBiIiIiCghERERkcgpIREREZHIKSERERGRyCkhERERkcgpIREREZHIZSwhMbPRZuZmlr+L03/TzB7v7rjSXPYxZvZAmuMuNrODMhBDj5W/q9uqtzCzq8zst50Y381sTPh5jpn9KGHYv5jZWjP7xMyGmNkXzOzvYfdXMhB+VkpcR9I5ZvZ9M5vfQ8saFe6bsQzMO+P7vpmtNLMvZmLevYmZvWZmR0YdR6Z0mJCY2WFm9oKZbTKzD8zseTOb0J1BpPpCdPf/cfdjunM5nfAfwOw0x70W+PfuDiCT5e+NB7eZPWNm50a1fHc/z91/HMZSAFwPHOPu/dx9I8E+cFPY/UBUcUrv4e7/4e4Z2aeTj3F3fzfcN5sysLhu3ffN7DYz+0nXw+rdUq0Hdz/I3Z/p5uVkzQ/SdhMSMxsAPAT8JzAY2Au4GtiW+dCiESZbA939xTQnWQQcZWYjMxiWdKCHD6YRQBHwWkK/vZO605YNJwLZNdp2gPb9nNat29Dd2/wDKoGP2hmeB/wQeAdYB/yG4MscYDTgQH7YvRL4YsK0VwG/DT+/G477Sfh3KHA28OeE8f8ZWAJsCv//c8KwZ4AfA88DHwOPA0PbiHmH+Yb9HBgTfr4CmJ+03A1AWdh9MPARcEDCOE8AZ7WxvNZytrFezgbeCuN+G/hmqjjDac4D/g58CPwasHBYDLgujPNt4ILEZSTF899AM7AlXNeXJsR0VrgtNgA/SNrOlwP/ADYCdwOD29kvTgJqgM3hNMeF/QcCC4A1wHvAT4BYYnkJapw+DMtxfDjsGqAJ2BrGfFPCOjk/XCdvh/1+BawKl70UmNTWtkgR9/fC2FYD05P2i9vCePcDPmX7/vpUWMbEdVqYRlmfB24APgiHFYZlfxdYC8wBisPxjwRqgVkEx9kaYFpC3MXh9n+H4Pj4c8K0nwdeINhnXwaObKf8l4Wxfgy8CRydsH99Pyznx+F6LUvYBin3y3D4dOCNcNhjwN5J+/S/htN+THAMfwb433D73Q3EE8b/EsF+9VFYpnHtlKWj/eBe4K5wuX8BDk4YvhL4N+D1MO5bgaKkbXEZ8D7B8VQI/JJgv1kdfi4E4mG830lYj88DV6Q4B44O18e0MO4Pw/U6AVgelvmmhBg/Q7DvbSQ4Xv8H2D2NY7zlvLMnwY+pD4AVwIyk9XM3wfn8Y4Jko7KN9Zxq3+9o3vcCvw23zblJ85sJNAD14fz+kLBNvhuui03htivaxX3jIIJz9gcEx9r3w/4pt2NHxyDBMfY+4fEd9vsqsDyd8ydwGNuP0VUE54f21sMXw3W8JWk+FeG+UNDRsZe0Pnb6/k3z2G099yasn0sT1s9XgBOAv4Xr+vttbZPW+bY7EAaEK/B24HhgUNLw6QQ73L5AP+B+4L+TDrB0EpIdxk38ggo/Dw5XyreBfOD0sHtIOPyZcGPvR3ByfgaY3UaZWuebtHJbvnjuAb6XNPwagoO/mOCAuCBp+I3A9W0sr7WcyWUFdiM4KPcPh40EDkoVZzjNQ8DuwChgPdu/6M8jOHmWAoOAPyavz6SYkrdFS0zzwjIeTFALdmA4/GLgxXD+hcAtwB1tzLuK4IQxheBA3IsweQMeCKfdDRgOLAb+T0J5G4AZBCfufyE4KbQkXc+w88nLCU4sg9n+BfwtYEi4fmcRnCiKUm2LpHkdR3By+lwY3+9IkZC0s78mr9OOytoIfCeMs5jg5LcoLEt/4A/ATxNOho0EVeMFBAd5HeHxSJAEPBOu6xhBEl0Ydm8Mx88Lt8lGYFiK8u9PcDLcM6GMnwk/fw94JRzHCPaPIQnboK398isE54cDw3L+EHghafstIjjPHESwzz1JcD4ZSLBPnxWOO57gRDcxLONZ4TovbGN7drQfNACnhOvzuwQn1YKEbfkqUBZuj+fZvu1btsXPwnVcHG6XF8PtPIzgy+XH4fifIzhXHQj8IBwvlhBH8jlwDkHt2zEECfgD4Xz3Cst/RDj+mHB7FobLfBb4ZRrHeMv5+E/AzeGyysPtdnRCXFsJ9psY8FPgxXa+J5KX1dG8Gwj2jTzC4zZpfre1rO+kZSwm+CIeTPBFeV5n9w2CY2sNwT5RFHZPDIe1tx1btntbx+A/gCkJy7kHuLyj8yfBMfMxwXdaAcE+W97Bevhi+Pkpdkz2fgHMSefYS5rnDvtGJ47d1nNvwvq5IizHjHC7/y5cxwcR7FP7trUfuXeQkIQLPjBcMbXhAhcBI8JhTwL/mnRSawgLsEMh6VpC8m1gcVJc/wucnfBl9cOEYf8KPNpGeVrnm7RyW754niDc0ROGFxD8ynoFeJSEX4Dh8GuAhW0sr7WcyWUl+LL6CPgaSQdmcpzhNIcldN/N9h3+KcIvu7D7i8nrs4MTSEtMpQn9FgPfCD+/QXhCCbtHtmznFPO+BbghRf8RBF84xQn9TgeeTijvioRhJWFMeyRs41QJyeQO9t8PCX/9Jm+LpPEWkpDEEiS3u5SQpFnWdxOGGUGty2cS+h3K9lqfIwl+DSUubx3BL7O8cNjBKcp0GeEPhIR+j5GiNo/gC25duO8UJA17EzipjfXW3n75CHBOwrA8gpP43gnTfiFh+FLgsoTu6wi/ZIH/IvxySIrriPa2fzv7wYtJca0hrEUJt+V5CcNPAP6RsC3q2fHX+T+AExK6jwVWJnTPAv4axjA2of9V7HwO3Cth+EbgtITu+4CL2yjfV4BlqfbH5H2WINFqAvonDP8pcFtCXH9MGPZZYEs767Z1WWnO+9kOttVtpP4i/lZC98/Z/uWb9r5BcBwua2O5bW5H2jkGw88/IfwOIPgC/pTt+3mb50+Cmrjfd3I9tKzrc4Gnws9G8IPi8HSOvaR5tu4bCf3SOXYnJwxvWT+xhHXghMmebz++v9Letu+wUau7v+HuZ7t7KUG2vyfBrznCz+8kjP5OuJJHdDTfTkpeTsuy9krofj/hcx1Bjc2u+JBgZbZy9waCneNzwHUert0E/QkSi05x90+B0whqONaY2f8zswPamaStMu5JsDO2SPzcGW3Nf2/g92b2kZl9RHCANZF6O5cRHNjJ9iZI7NYkzOcWgl8jOy3f3evCjx1txx3KamazzOyNsBH2RwS/tId2MA/YeR0m72+dkU5ZE5c1jCABW5ow/qNh/xYb3b0xobtl+wwl+KXX1jr/ess8w/keRnBC3IG7ryD4JXcVsM7M7jSzPcPBbW3TFu3tN79KWPYHBCfOxON2bcLnLSm6E+c1K6ksZQTbbSdp7Aet69/dmwl+cO2ZajjBvpA4bL27b03oTnUeTBz/doKT/sPu/vdU8SZIa32Y2fBwG71nZpsJLoGks5+3xPuBu3+cFHN759OiNNsKpDPvTJyf0t032tuXO9qObR2DENQEnGxmhcDJwF/cvWVe7Z0/Ozq22nMvcGh4nB5OkAA8l7DMjo699qQzffJ23OjbG01vCf+3dTyn1Knbft39r2z/YoagSn3vhFFGEdSirGVnnxKcdFvskTjrDhadvJyWZb3XwXSp7BCHme2RNHw5wa9jEsbZC7iS4FrydeFOl+hAguvzHS6PHcuNuz/m7lMIviT+SnDZpLPWEFQHtijrYPyO1neyVQTtOXZP+Cty91TrfxXB9e1U/bcRtO1pmccAd0/3lum2Ym7tb2aTCGoFTiWoSt2d4PKRpTH/Ney43kalGVcq6ZQ1sTwbCA7WgxLGH+ju6STVGwiqQtta5/+dtN12c/eUd5C5++/c/TCCY80JLku0zCfV/DuyiqDmLnH5xe7+wi7O65qkeZW4+x3JI6a5H5QljJ9HcPysTjWcYF9IHJa8L6Y6DyaOfzPBZa1jzeywdkuZvp+GcYxz9wEEl6gSy9feMb4aGGxmiT+8dvV8uivz7uj8syvnp7T2Ddrflzvajm1y99cJEpjjgTMIEpTEZbZ1/mwvnnbXg7t/RNBe8tRwmXck/FjuzLGXajnpTN/Z7dShju6yOSD8pVEadpcRVHm13IFyB3CJme1jZv0Ibpe9KymLbFEDfMPMCsyskuD6bYv1BA2j9m0jlIeB/czsDDPLN7PTCKoRH0qrlDt6GTjIzMrNrIjgF2Hyso5o6TAzI0jCFgDnEHxx/ThheCFwCMGlnlRqgMPD5wAMJKiia5l2hJlNNbPdCL7APiHInDvrbuAiM9vLzHYnOBm3Zy1tr+tU5gDXmNneAGY2zMxOamPcBcA0MzvazPLCmA5w9zUEB891ZjYgHPYZMzuijfnsSsz9CRLi9UC+mV1B0D4hHXcDZ5vZZ82shCAB3SWdLWv4C30ecIOZDYcgCTazY9NYVjPB5abrzWxPM4uZ2aHhfvlb4MtmdmzYv8jMjmw5nhOZ2f5mNjmcbitBgtSyL84HfmxmYy0wzsyGpLEq5gD/ZuFzesxsoJl9PY3pUpkHnGdmE8MYdjOzE5O++Fqksx8cYmYnh7/6LyY4/hLvrDvfzErNbDBBg9672ontDuCH4XExlOA6+m8BzOzbBOeHs4ELgdvDc2VX9Sc4X3wU/mD6XtLwNo8Xd19F0D7ip+E+MY7g3PY/XQ2qm+bd2fNTZ/aNh4A9zOxiMys0s/5mNjEc1uZ2TNPvCLbx4QRtSFq0d/78H+CLZnZq+N02xMzKw2HprIffAWcSXPZPTII6c+yl+v7tzmM3bR3VkHxM0FDoJTP7lOCAfZXgmigEJ8L/JmhQ9TbBiew7bczrRwSZ4IcEtw63rrywev4a4HkLqog+nzihB895+FK43I0ELXm/5O4b0ivmDvP6G0HDpD8StBD+c9LwvwCbEnbSCwmq1n4UZp/TCL5wJ4XDpwLPuHvKTNrdnyA4mS0nuIaWmETlhWVaTVAldgRB+5fOmkfwBbgcWEaQVDXSdnLzU4ID7yMz+24a8/8VQduhx83sY4L9YGKqEd19McE6uoHgV+mf2P6r40yCOw9a7l64lxSXD9qJ4RQz+9DMbmxjnMcIrn3+jeDXylbSrB5290cILkU+RdCY66k042pLZ8t6WbjcFy2ogv8jQZusdHyXoH3TEoL96GdAXvjlcBLBF+p6gnXxPVIf94UEz97ZQFA1PjycDoJnrtxNsI9tJkg6izsKyt1/H8ZyZ1imVwl+QXaau1cTNJS7iWB9riD4kk8lnf3gQYLLpR8StFE7Obw02+J3BOV9K/xr77kYPwGqCY6/Vwju2vmJmY0i2KfOdPdP3P134Xg3dFjgjl1N0JhzE/D/CG4oSNTRMX46wWWk1cDvgSvDc1V36Oq8FwCfDWN/oKORO7NvhJeSpgBfJtjP/w4cFQ5OuR07EfcdBG0pnkr6bmrz/Onu7xK0UZpFcOzWEDQah/TWwyJgLLDW3Vtr6Ttz7KX6/u3OY7czWu5gkARmdgxBY92vpDHuSwSNf17NeGBpMrPjCRp8JV/mEsl5ZnYVQWPlb7UxfCVBA+o/9mRcIrlOD6VJwd0fJ/h1lM64KWsKepKZFRNk+Y8T1OZcSfDLREREpFfQy/X6BiOowv2Q4JLNGwTXP0VERHoFXbIRERGRyKmGRERERCKnNiS9yNChQ3306NFRhyEi0qssXbp0g7sP63hMiZISkl5k9OjRVFdXRx2GiEivYmZdefKy9BBdshEREZHIKSERERGRyCkhERERkcgpIREREZHIKSERERGRyCkhERERkcgpIRERkZTqli1jwy1zqVu2LOpQJAfoOSQiIrKTumXLeHfadLy+HovHGXXrQkoqKqIOS/ow1ZBkgJktNLN1ZvZqO+McaWY1Zvaamf2pJ+MTEelI3eIleH09NDfjDQ3ULV4SdUjSxykhyYzbgOPaGmhmuwM3A1Pd/SDg6z0TlohIekqqJmDxOMRiWEEBJVUTog5J+jhdsskAd3/WzEa3M8oZwP3u/m44/roeCUxEJE0lFRWMunUhdYuXUFI1QZdrJOOUkERjP6DAzJ4B+gO/cvffpBrRzGYCMwFGjRrVYwGKiJRUVCgRkR6jSzbRyAcOAU4EjgV+ZGb7pRrR3ee6e6W7Vw4bppdViohI36QakmjUAhvc/VPgUzN7FjgY+Fu0YYmIiERDNSTReBCYZGb5ZlYCTATeiDgmERGRyKiGJAPM7A7gSGComdUCVwIFAO4+x93fMLNHgeVAMzDf3du8RVhERKSvU0KSAe5+ehrj/AL4RQ+EIyIikvV0yUZEREQip4REREREIqeERERERCKnhEREREQip4REREREIqeERERERCKnhEREREQip4REREREIqeERERERCKnhEREREQip4REREREIqeERERERCKnhEREREQip4REREQyqmZdDfNfmU/NupqoQ5Eslh91ACIi0nfVrKthxuMzqG+qJx6LM++YeZQPL486LMlCqiHJADNbaGbrzOzVDsabYGZNZnZKT8UmItKTqtdWU99UTzPNNDQ3UL22OuqQJEspIcmM24Dj2hvBzGLAz4DHeiIgEZEoVI6oJB6LE7MYBXkFVI6ojDokyVK6ZJMB7v6smY3uYLTvAPcBEzIfkYhINMqHlzPvmHlUr62mckSlLtdIm5SQRMDM9gK+CkxGCYmI9HHlw8uViEiHdMkmGr8ELnP3po5GNLOZZlZtZtXr16/PfGQiIiIRUA1JNCqBO80MYChwgpk1uvsDySO6+1xgLkBlZaX3ZJAiIiI9RQlJBNx9n5bPZnYb8FCqZERERCRXKCHJADO7AzgSGGpmtcCVQAGAu8+JMDQREZGspIQkA9z99E6Me3YGQxEREekV1KhVREREIqeERERERCKnhEREREQip4REREREIqeERERERCKnhEREREQip4REREREIqeERERERCKnhEREREQip4RERCRUt2wZG26ZS92yZVGHIpJz9Oh4ERGCZOTdadPx+nosHmfUrQspqaiIOiyRnKEaEhERoG7xEry+Hpqb8YYG6hYviTokkZyihEREBCipmoDF4xCLYQUFlFRNiDokkZyiSzYiIkBJRQWjbl1I3eIllFRN0OUakR6mhEREJFRSUaFERCQiumQjIiIikVNCkgFmttDM1pnZq20M/6aZLQ//XjCzg3s6RhERkWyihCQzbgOOa2f428AR7j4O+DEwtyeCEhERyVZqQ5IB7v6smY1uZ/gLCZ0vAqUZD0pERCSLqYYkeucAj7Q10Mxmmlm1mVWvX7++B8MSERHpOUpIImRmRxEkJJe1NY67z3X3SnevHDZsWM8FJyIi0oN0ySYiZjYOmA8c7+4bo45HREQkSqohiYCZjQLuB77t7n+LOh4REZGoqYYkA8zsDuBIYKiZ1QJXAgUA7j4HuAIYAtxsZgCN7l4ZTbQiIiLRU0KSAe5+egfDzwXO7aFwREREsp4u2YiIiEjklJCIiIhI5JSQiIiISOSUkIiIiEjklJCI5LC6ZcvYcMtc6pYtizoUEclxustGJEfVLVvGu9Om4/X1WDzOqFsXUlJREXVYIpKjVEMikqPqFi/B6+uhuRlvaKBu8ZKoQxKRHKaERCRHlVRNwOJxiMWwggJKqiZEHZJkSM26Gua/Mp+adTVRhyLSJl2yEclRJRUVjLp1IXWLl1BSNUGXa/qomnU1zHh8BvVN9cRjceYdM4/y4eVRhyWyEyUkIjmspKJCiUgfV722mvqmeppppqG5geq11UpIJCvpko2ISB9WOaKSeCxOzGIU5BVQOUKvzZLspBoSEZE+rHx4OfOOmUf12moqR1SqdkSylhISyTl1y5ap3YTklPLh5UpEJOspIZGcomdviIhkJ7UhkZyiZ29IT9MttyLpUQ1JBpjZQuBLwDp3/1yK4Qb8CjgBqAPOdve/9GyUuanl2Rve0KBnb0jG6ZZbkfQpIcmM24CbgN+0Mfx4YGz4NxH4r/B/Rrx8389Z/dwT7DlpCgdXHQkrn4PRk4KBqT6XVcGqxa3dNZtWUP3WY1TueyzlA8d0PH3S55q/3kd1URGVg/an/MM1nZqWlc9RM2gk1R++SeXWrZQf8LV2Y01n+lE/+x51zzxCyZHHUzK0AZ67rt04ulT+sipqXv1d9q6/NNZl2tNnYv0lln//r1K+rb7H198ufQ7XZfWyX1PftI1mnIambVS/+ftdKkOm4uvR5Yl0wNw96hj6JDMbDTzURg3JLcAz7n5H2P0mcKS7r2lvnpWVlV5dXd2pOF6+7+f4lbeS3wSNMfBjPubgAZ/SbDEA8rxph89NeQU8uff/5eh3rifW3MCyoiL+z4jBNJhR4M4taz+gYuvWNqdP/vxyYcEO0897fwP/VN+U1rQAr8RjzNhj6A7LP3hbQ8pYOzt9OmXoSvmb8gpYsPe3md/8GA1G1q2/9rZ7Nqy/5PLHLJ+5a9Yzfmtdj62/ljJUF8Wp2rI1rWkT1+WrBcaMPYbSaEb+Lqy/THxO57jZ1c+GYd6I5+VjgHkTxOLYcbNhy8bIkhMzW+ruut85y6mGJBp7AasSumvDfjslJGY2E5gJMGrUqE4vaPVzT1DWBDEHmmDVhjgVAz7Gmx3DyTN2+Nzc1EDJPx7CrIE8a2ZpYT4NZjSb0Qj8pTCfQ7Y2tzl98ue/FBbvMP3Sojjjtn6c1rSGs7QoafrCfCq21KWMtbPTp1OGrpS/uamBVRv/RMNgsnL9tbfds2H9JZcfb2JpYYzKHlx/y4virQnN/IH9uGXNesq3betwu7esy4rmZuauWc/S4kImbN3KP21tSHvZmfqcznHT1c9NTQ2AYwaNDVvxh/4veTjE4uQdPxvb8oFqTmQnatQaDUvRL2VVlbvPdfdKd68cNmxYpxe056QpNMagyYIakj2HNYLFiOUXkJdfuNPn/IJCDp96LvkFQXdVfRNxh5g7BQ4T6pvbnT7584T6ZuLurdNXbmtIe9q8/EIqtzUkLN+pqm9qM9bOTp9OHF0pf35BIV+vPDlr11972z0b1l9y+eN5+Uxs6Nn1V128PaFpMOMvu5Wktd0T1+X4JmfGJ9sor2/q1LIz9Tmd4ybdz3n5BVh+Id76OY5bDMsvgFicZmJYXh4xmonRjDduo+mhWTQ/+RO4fWpw2UgkpBqSaNQCZQndpcDqTCzo4K9dysvAqs62IRnxWVj5HOWjJzGvC20AyoF5XbiGX77yOeYltmH4YlIbhoRYd2n6Dj53ufxlVcwbPDp7118H6zLy9ZdY/l1oQ9LV9Vf51/uIv/8EDc1NFOQZlUf/B6RThg72y55sg5PJ+Cyd/sVD4NHL8aZ6LKzpysNpbqwnb+VzqiWRVmpDkiEdtCE5EbiA4C6bicCN7t7hUbkrbUhEpGtq1tVE9pTTPnOXTksD2uIh+KOX09y4jXrPZ/1hVzOqaEvGL9+oDUnvoBqSDDCzO4AjgaFmVgtcCRQAuPsc4GGCZGQFwW2/06KJVEQ6EuVTTvvMi/HKqloTDhvxWbb97Rn+8383ctHzV+LWhMXicNYi1ZbkOCUkGeDup3cw3IHzeygcEemlWl6M19Dc0HdejFdWRUlZFdO2/Qf5LzViNAeXc3T5JucpIZEeo3fIiHROX34x3vB/+iKN1b+isakeixUQa2l3IjlLCYn0CL1DRmTX9NkX45VVUf/NB/iv225jyH6TOVu1IzlPt/1Kj9A7ZEQkWclnDqVm9HR+s2pE1KFIFlBCIj2i5R0yxGJ6h4yItPrigSN4a8OnvLX+k6hDkYgpIZEeUVJRwahbFzLswgt1uUZEWk0+YDgAT/11XcSRSNTUhkR6TElFhRIREdlB2eAS9h/Rnz++sZZzJ+0bdTgSIdWQiESkbtkyNtwyl7ply6IORSRSRx84nMZ3XmLrU7/Q4+RzmGpIRCKgu45Etps6pJbv5F9D/NlGeOE6PSQtR6mGRCQCuutIZLv9trxMgTWSRzM0he8rkpyjhEQkArrrSGS7vH0m0WwFNJKHx+LbX9QnOUWXbNJgZoOAPYEtwEp3b444JOnlWu460pNrRYCyKl6cdCsvPvkgp5z0DfbV5ZqcpISkDWY2kOB9M6cDcWA9UASMMLMXgZvd/ekIQ5ReTncdiWw3YOwXuPlxOKToIHSvTW5SQtK2e4HfAJPc/aPEAWZ2CPBtM9vX3RdEEZyISF8ycmARAGs2bY04EomKEpI2uPuUdoYtBZb2YDgiIn3a0H6FxPKM95WQ5Cw1au2AmZ2T1B0zsyujikdEpC+K5Rkj+heqhiSHKSHp2NFm9rCZjTSzzwEvAv07msjMjjOzN81shZldnmL4QDP7g5m9bGavmdm0TAQvItJVNetqmP/KfGrW1WR0OXsMLOL9zVsyugzJXrpk0wF3P8PMTgNeAeqA0939+famMbMY8GtgClALLDGzRe7+esJo5wOvu/uXzWwY8KaZ/Y+712emJCIinVezroYZj8+gvqmeeCzOvGPmUT68PCPLGjmwmDfe35yReUv2Uw1JB8xsLHARcB+wkqAxa0kHk1UBK9z9rTDBuBM4KWkcB/qbmQH9gA+Axu6MXUSkq6rXVlPfVE8zzTQ0N1C9tjpjyxoxoIj3N23F3TO2DMleSkg69gfgR+7+f4AjgL8DHT1Wcy9gVUJ3bdgv0U3AgcBqgtqXi/R8ExHJNpUjKonH4sQsRkFeAZUjKjO2rJEDi6irb2LzVv02y0W6ZNOxKnffDOBB2n6dmS3qYBpL0S855T8WqAEmA58BnjCz51qW1Tojs5nATIBRo0Z1PnqRXqBu2TI9JC5LlQ8vZ94x86heW03liMqMXa6BoA0JwPubtjKwuCBjy5HspISkDWZ2mLv/OTlBAHD3v5vZAGCUu7+aYvJaoCyhu5SgJiTRNGB2mOSsMLO3gQOAHV516e5zgbkAlZWVqseUPkcvGsx+5cPLM5qItNj+LJIt7L9Hh/cOSB+jSzZt+5qZvWBmV5jZiWZWZWaHm9l0M/tv4CGguI1plwBjzWwfM4sD3wCSa1XeBY4GMLMRwP7AW5kpikj20osGpUVLDcnazbr1NxephqQN7n5J+A6bU4CvA3sQvMvmDWBOe3fauHujmV0APAbEgIXu/pqZnRcOnwP8GLjNzF4huMRzmbtvyGihRLJQy4sGvaFBLxrMccP7F2Gmp7XmKiUk7XD3D81sT+Dt8A+CtiBHA+3e+uvuDwMPJ/Wbk/B5NXBMtwYs0gvpRYPSIp6fx9B+hXpaa45SQtKxTxI+FwFfIqglEZFuohcNSouRA4tUQ5KjlJB0wN2vS+w2s2vZuT2IiIh0gz0GFPHOxrqow5AIqFFr55WA3o4tIpIJQQ2JHh+fi1RD0oGw0WnL7bYxYBjw79FFJCLSd+0xsJjNWxv5dFsjuxXqKyqXaGt37EsJnxuBte6uxwiKiGTAHgMLAXh/81Y+M6xfxNFIT1JC0gF3fyfqGEREcsUeA4LHO72/SQlJrlEbEhERyRrbn9aqO21yjRISERHJGtvfZ6OGrblGCYmIiGSNooIYg0oKeF+Pj885SkhERCSrHLnbSireWQirFnc8svQZSkhERCR7rFrM7E9+yFc+vA1un6qkJIcoIRERkeyx8jnyvYEYzdBUDyufizoi6SFKSEREJHuMnkRzXgGNnofH4jB6UtQRSQ9RQiIiItmjrIpnD13A9Y1fZ+1X7oayqqgjkh6iB6OJiEhWKdj789zclMdR/T7HHlEHIz1GNSQiIpJV+hUFv5U/3aa3dOQSJSQZYmbHmdmbZrbCzC5vY5wjzazGzF4zsz/1dIwiItmouCAGwNaGpogjkZ6kSzYZYGYx4NfAFKAWWGJmi9z99YRxdgduBo5z93fNbHgkwYqIZJmWhGSLEpKcohqSzKgCVrj7W+5eD9wJnJQ0zhnA/e7+LoC7r+vhGEVEslJxPExI6psjjkR6khKSzNgLWJXQXRv2S7QfMMjMnjGzpWZ2ZqoZmdlMM6s2s+r169dnKFwRkexRFNaQ1NWrDUkuUUKSGZainyd15wOHACcCxwI/MrP9dprIfa67V7p75bBhw7o/UhGRLKM2JLlJbUgyoxYoS+guBVanGGeDu38KfGpmzwIHA3/rmRBFRLJTPD+P/DxTG5IcoxqSzFgCjDWzfcwsDnwDWJQ0zoPAJDPLN7MSYCLwRg/HKSKSlYoLYmpDkmNUQ5IB7t5oZhcAjwExYKG7v2Zm54XD57j7G2b2KLAcaAbmu/ur0UUtIpI9iuIx1ZDkGCUkGeLuDwMPJ/Wbk9T9C+AXPRmXiGS/umXLqFu8hJKqCZRUVEQdTiSKC2JqQ5JjlJCIiGSRumXLeHfadLy+HovHGXXrwpxMSooLYrrLJseoDYmISBapW7wEr6+H5ma8oYG6xUuiDikSxfEYWxrUhiSXKCEREckiJVUTsHgcYjGsoICSqglRhxSJ4oIYW+t1ySaX6JKNiEgWKamoYNStC9WGJB5j/cfbog5DepASEhGRLFNSUZGziUiL4gLdZZNrdMlGRESyTlFBjC26ZJNTlJCIiEjWKY7nqYYkxyghERGRrFMSz1cNSY5RQiIiIlmnKGxD4p78XlLpq9SotZdraGigtraWrVu3Rh2KtKGoqIjS0lIKCgqiDkWk12h54++2xmaKws/Stykh6eVqa2vp378/o0ePxsyiDkeSuDsbN26ktraWffbZJ+pwRHqN4oKgAn9LfZMSkhyhSza93NatWxkyZIiSkSxlZgwZMkQ1WCKdVBwPkpA6NWzNGUpI+gAlI9lN20ek84rjQQW+GrbmDiUkIiKSdVrakOiNv7lDCYl0WSwWo7y8vPVv9uzZ3TbvmpoaHn744W6ZV79+/Toc55e//CV1dXWt3SeccAIfffRRtyxfRNLXkpDoWSS5QwlJhpjZcWb2ppmtMLPL2xlvgpk1mdkpPRlfdyouLqampqb17/LL2yxup3VnQpKO5ITk4YcfZvfdd++x5YtIoDi+vVGr5AbdZZMBZhYDfg1MAWqBJWa2yN1fTzHez4DHumO5V//hNV5fvbk7ZtXqs3sO4MovH9Tp6TZt2kRVVRWLFi1i//335/TTT2fy5MnMmDGDf/mXf2HJkiVs2bKFU045hauvvhqAJUuWcNFFF/Hpp59SWFjIE088wRVXXMGWLVv485//zL/9279x2mmntS7jtddeY9q0adTX19Pc3Mx9993H2LFjuf7661m4cCEA5557LhdffPEOsT3zzDNce+21PPTQQwBccMEFVFZWsnnzZlavXs1RRx3F0KFDefrppxk9ejTV1dUMHTo05XxXrlzJ8ccfz2GHHcYLL7zAXnvtxYMPPkhxcfGurG4RCRWphiTnKCHJjCpghbu/BWBmdwInAa8njfcd4D6gV79ffMuWLZSXl7d2tyQON910E2effTYXXXQRH374ITNmzADgmmuuYfDgwTQ1NXH00UezfPlyDjjgAE477TTuuusuJkyYwObNmykpKeHf//3fqa6u5qabbtppuXPmzOGiiy7im9/8JvX19TQ1NbF06VJuvfVWXnrpJdydiRMncsQRR1CRxovKLrzwQq6//nqefvpphg4dusOwtuY7aNAg/v73v3PHHXcwb948Tj31VO677z6+9a1vdW2liuS41ks2qiHJGUpIMmMvYFVCdy0wMXEEM9sL+CowmXYSEjObCcwEGDVqVLsL3ZWajO7Qcskm2ZQpU7jnnns4//zzefnll1v733333cydO5fGxkbWrFnD66+/jpkxcuRIJkwIVsWAAQM6XO6hhx7KNddcQ21tLSeffDJjx47lz3/+M1/96lfZbbfdADj55JN57rnn0kpI2tPWfKdOnco+++zTmpAdcsghrFy5skvLEpHg0fGgGpJcojYkmZHqPs/k5x//ErjM3ds92tx9rrtXunvlsGHDuiu+HtHc3Mwbb7xBcXExH3zwAQBvv/021157LU8++STLly/nxBNPZOvWrbh7p2+PPeOMM1i0aBHFxcUce+yxPPXUU2k9Zjo/P5/m5ubW7nSeEdLefAsLC1s/x2IxGhsbO5yfiLRPNSS5RwlJZtQCZQndpcDqpHEqgTvNbCVwCnCzmX2lR6LrITfccAMHHnggd9xxB9OnT6ehoYHNmzez2267MXDgQNauXcsjjzwCwAEHHMDq1atZsmQJAB9//DGNjY3079+fjz/+OOX833rrLfbdd18uvPBCpk6dyvLlyzn88MN54IEHqKur49NPP+X3v/89kyZN2mG6vffem9dff51t27axadMmnnzyydZhbS0vnfmKSPcpamnUqhqSnKFLNpmxBBhrZvsA7wHfAM5IHMHdW58jbma3AQ+5+wM9GGO3SW5DctxxxzF9+nTmz5/P4sWL6d+/P4cffjg/+clPuPrqq6moqOCggw5i33335Qtf+AIA8Xicu+66i+985zts2bKF4uJi/vjHP3LUUUcxe/ZsysvLd2rUetddd/Hb3/6WgoIC9thjD6644goGDx7M2WefTVVVFRA0Pk2+XFNWVsapp57KuHHjGDt27A7DZ86cyfHHH8/IkSN5+umnW/uPHz8+5Xx1eUYkM+KxPPJMzyHJJaY3KWaGmZ1AcFkmBix092vM7DwAd5+TNO5tBAnJve3Ns7Ky0qurq3fo98Ybb3DggQd2Y+SSCdpOIp130BWPcnrVKH74pc92aT5mttTdK7spLMkQ1ZBkiLs/DDyc1G9OG+Oe3RMxiYj0JsXxmN5lk0PUhkRERLJScTzGVjVqzRlKSEREJCsVF8TUqDWHKCEREZGspIQktyghERGRrFRUENNzSHKIEhIREclKxfGYbvvNIUpIpMvMjFmzZrV2X3vttVx11VVdnu8PfvADysrK6Nev3w79t23bxmmnncaYMWOYOHHiDs8Cuf322xk7dixjx47l9ttv73IMIhKdkniMOtWQ5AwlJLlo1WJ47rrgfzcoLCzk/vvvZ8OGDd0yvxZf/vKXWbx45xgXLFjAoEGDWLFiBZdccgmXXXYZAB988AFXX301L730EosXL+bqq6/mww8/7NaYRKTnFKkNSU5RQpJrVi2G26fCU9cE/7shKcnPz2fmzJnccMMN3RDgdp///OcZOXLkTv0ffPBBzjrrLABOOeUUnnzySdydxx57jClTpjB48GAGDRrElClTePTRR7s1JhHpOcUFumSTS/RgtFyz8jloqgdvCv6vfA7Kqro82/PPP59x48Zx6aWXtjnO008/zSWXXLJT/5KSEl544YW0l/Xee+9RVha8Kig/P5+BAweycePGHfoDlJaW8t5773WiFCKSTYrVqDWnKCHJNaMnQSweJCOxeNDdDQYMGMCZZ57JjTfeSHFxccpxjjrqKGpqarq8rFSvOzCzNvuLSO9UHA8u2ezK28Cl91FCkmvKquCsRUHNyOhJ3VI70uLiiy9m/PjxTJs2LeXw7qohKS0tZdWqVZSWltLY2MimTZsYPHgwpaWlPPPMM63j1dbWcuSRR3a2GCKSJYoKYjQ71Dc1U5gfizocyTAlJLmorKpbE5EWgwcP5tRTT2XBggVMnz59p+HdVUMydepUbr/9dg499FDuvfdeJk+ejJlx7LHH8v3vf7+1Ievjjz/OT3/60y4vT0SiURIPkpAt9U1KSHKAGrVKt5o1a1a33W1z6aWXUlpaSl1dHaWlpa23Ep9zzjls3LiRMWPGcP311zN79mwgSIh+9KMfMWHCBCZMmMAVV1zB4MGDuyUWEel5xQVhQqKGrTnBUl13l+xUWVnp1dXVO/TTa+17B20nkc57sOY9LrqzhqdmHcG+w/p1PEEbzGypu1d2Y2iSAaohERGRrFSkGpKcooQkQ8zsODN708xWmNnlKYZ/08yWh38vmNnBUcQpIpKtWi7Z6FkkuUEJSQaYWQz4NXA88FngdDP7bNJobwNHuPs44MfA3J6NUkQku7U0atXj43ODEpLMqAJWuPtb7l4P3AmclDiCu7/g7i3PNX8RKO3hGEVEslrrJRslJDlBCUlm7AWsSuiuDfu15RzgkYxGJCLSyxTH1YYkl+g5JJmR6pGCKW9nMrOjCBKSw9oYPhOYCTBq1Kjuik9EJOupDUluUQ1JZtQCZQndpcDq5JHMbBwwHzjJ3TemmpG7z3X3SnevHDZsWEaC7SozY9asWa3d1157beszQ3ZVXV0dJ554IgcccAAHHXQQl1++vV3wtm3bOO200xgzZgwTJ05k5cqVXVqWiGSnYl2yySlKSDJjCTDWzPYxszjwDWBR4ghmNgq4H/i2u/+tJ4OrWVfD/FfmU7OuplvmV1hYyP33399tD0Rr8d3vfpe//vWvLFu2jOeff55HHgmuai1YsIBBgwaxYsUKLrnkEi677LJuXa6IZIftl2yaI45EeoISkgxw90bgAuAx4A3gbnd/zczOM7PzwtGuAIYAN5tZjZlVtzG7blWzroYZj8/gP//yn8x4fEa3JCX5+fnMnDmTG264oesBhkpKSjjqqKMAiMfjjB8/ntraWgAefPBBzjrrLABOOeUUnnzyyZQv1hOR3q0wPw8z2FLfGHUo0gPUhiRD3P1h4OGkfnMSPp8LnNvTcVWvraa+qZ5mmmlobqB6bTXlw8u7PN/zzz+fcePGcemll7Y5zq6+XO+jjz7iD3/4AxdddBEA7733HmVlwRWx/Px8Bg4cyMaNGxk6dGgXSyEi2cTMKC6IqVFrjlBCkmMqR1QSj8VpaG6gIK+AyhHd8zTlAQMGcOaZZ3LjjTdSXFyccpxdebleY2Mjp59+OhdeeCH77rsvQMraEL2aXKRvUkKSO5SQ5Jjy4eXMO2Ye1WurqRxR2S21Iy0uvvhixo8fz7Rp01IO35UakpkzZzJ27Fguvvji1n6lpaWsWrWK0tJSGhsb2bRpk16iJ9JHFRXE2FKvNiS5QAlJDiofXt6tiUiLwYMHc+qpp7JgwQKmT5++0/DO1pD88Ic/ZNOmTcyfP3+H/lOnTuX222/n0EMP5d5772Xy5MmqIRHpo4rjMd32myPUqFW61axZs7rlbpva2lquueYaXn/9dcaPH095eXlrYnLOOeewceNGxowZw/XXX8/s2bO7vDwRyU4lcV2yyRWqIZEu++STT1o/jxgxgrq6ui7Ps7S0tM07Z4qKirjnnnu6vAwRyX5FBTHqdJdNTlANiYiIZK2gUavakOQCJSQiIpK1igtibNWTWnOCEhIREclaxWpDkjOUkIiISNYq0nNIcoYSEhERyVolcV2yyRVKSEREJGsVF8Soa2jS+6pygBIS6TIzY9asWa3d1157LVdddVWX5/uDH/yAsrIy+vXrt0P/bdu2cdpppzFmzBgmTpzIypUru7wsEclOxfEYTc1OQ5MSkr5OCUkOqlu2jA23zKVu2bJumV9hYSH3339/tzwQLdGXv/xlFi9evFP/BQsWMGjQIFasWMEll1zCZZdd1q3LFZHsUVQQA1A7khyghCTH1C1bxrvTprP+V7/i3WnTuyUpyc/PZ+bMmdxwww3dEOF2n//85xk5cuRO/R988EHOOussAE455RSefPJJVeeK9FHFYUKix8f3fXpSa46pW7wEr6+H5ma8oYG6xUsoqajo8nzPP/98xo0bx6WXXtrmOLvycr1U3nvvPcrKyoAgGRo4cCAbN25k6NChnQ9cRLJaSTysIVHD1j5PCUmOKamagMXjeEMDVlBASdWEbpnvgAEDOPPMM7nxxhspLi5OOU5nX67XllS1IXq5nkjf1HLJpk4JSZ+nhCRDzOw44FdADJjv7rOThls4/ASgDjjb3f+S6bhKKioYdevCoGakakK31I60uPjiixk/fjzTpk1LOby7akhKS0tZtWoVpaWlNDY2smnTJgYPHrzLcYtI9iqOqw1JrlBCkgFmFgN+DUwBaoElZrbI3V9PGO14YGz4NxH4r/B/xpVUVHRrItJi8ODBnHrqqSxYsIDp06fvNLy7akimTp3K7bffzqGHHsq9997L5MmTVUMi0kcVF8Q47tNFvPPz68iffBwHf63ty8LSu6lRa2ZUASvc/S13rwfuBE5KGuck4DceeBHY3cx2bsHZy8yaNavb7ra59NJLKS0tpa6ujtLS0tZbic855xw2btzImDFjuP7665k9e3b7MxKRXmvzn37NeU8/y35PrcGvvJWX7/t51CFJhqiGJDP2AlYldNeyc+1HqnH2AtYkjmRmM4GZAKNGjer2QLvDJ5980vp5xIgR1NXVdct8f/7zn/Pzn+988ikqKuKee+7plmWISHbb8pc/kd8EMQeaYNVzT6iWpI9SDUlmpLp+kNwSM51xcPe57l7p7pXDhg3rluBERHqLUUccQ2MMmgwaY7DnpClRhyQZohqSzKgFyhK6S4HVuzCOiEhOO/hrl/IyQc3InpOmqHakD1NCkhlLgLFmtg/wHvAN4IykcRYBF5jZnQSXcza5+xp2gburUWcW00PbRLrm4K9dqkQkByghyQB3bzSzC4DHCG77Xejur5nZeeHwOcDDBLf8riC47Tf1vbIdKCoqYuPGjQwZMkRJSRZydzZu3EhRUVHUoYiIZDXTr7feo7Ky0qurq3fo19DQQG1tLVu3bo0oKulIUVERpaWlFBQURB2KSE4ys6XuXhl1HNI+1ZD0cgUFBeyzzz5RhyEiItIlustGREREIqeERERERCKnhEREREQip0atvYiZrQfe2cXJhwLd80z33kNlzg0qc27oSpn3dnc9WTLLKSHJEWZWnWutzFXm3KAy54ZcLHOu0SUbERERiZwSEhEREYmcEpLcMTfqACKgMucGlTk35GKZc4rakIiIiEjkVEMiIiIikVNCIiIiIpFTQtLHmNlxZvamma0ws8tTDDczuzEcvtzMxkcRZ3dKo8zfDMu63MxeMLODo4izO3VU5oTxJphZk5md0pPxdbd0ymtmR5pZjZm9ZmZ/6ukYu1sa+/VAM/uDmb0clnmX3hieTcxsoZmtM7NX2xje585fksDd9ddH/oAY8A9gXyAOvAx8NmmcE4BHAAM+D7wUddw9UOZ/BgaFn4/PhTInjPcU8DBwStRxZ3gb7w68DowKu4dHHXcPlPn7wM/Cz8OAD4B41LF3sdyHA+OBV9sY3qfOX/rb8U81JH1LFbDC3d9y93rgTuCkpHFOAn7jgReB3c1sZE8H2o06LLO7v+DuH4adLwKlPRxjd0tnOwN8B7gPWNeTwWVAOuU9A7jf3d8FcPdcKLMD/c3MgH4ECUljz4bZvdz9WYJytKWvnb8kgRKSvmUvYFVCd23Yr7Pj9CadLc85BL+werMOy2xmewFfBeb0YFyZks423g8YZGbPmNlSMzuzx6LLjHTKfBNwILAaeAW4yN2beya8yPS185ckyI86AOlWlqJf8n3d6YzTm6RdHjM7iiAhOSyjEWVeOmX+JXCZuzcFP6B7tXTKmw8cAhwNFAP/a2YvuvvfMh1chqRT5mOBGmAy8BngCTN7zt03Zzi2KPW185ckUELSt9QCZQndpQS/njo7Tm+SVnnMbBwwHzje3Tf2UGyZkk6ZK4E7w2RkKHCCmTW6+wM9EmH3Sne/3uDunwKfmtmzwMFAb01I0inzNGC2uzuwwszeBg4AFvdMiJHoa+cvSaBLNn3LEmCsme1jZnHgG8CipHEWAWeGrdU/D2xy9zU9HWg36rDMZjYKuB/4di/+xZyowzK7+z7uPtrdRwP3Av/aS5MRSG+/fhCYZGb5ZlYCTATe6OE4u1M6ZX6XoEYIMxsB7A+81aNR9ry+dv6SBKoh6UPcvdHMLgAeI2ilv9DdXzOz88LhcwjuuDgBWAHUEfzK6rXSLPMVwBDg5rDGoNF78VtD0yxzn5FOed39DTN7FFgONAPz3T3lraO9QZrb+MfAbWb2CsGljMvcfUNkQXcDM7sDOBIYama1wJVAAfTN85fsSI+OFxERkcjpko2IiIhETgmJiIiIRE4JiYiIiEROCYmIiIhETgmJiIiIRE4JiYiIiEROCYmIiIhETgmJiOwyM5tgZsvNrMjMdjOz18zsc1HHJSK9jx6MJiJdYmY/AYoIXmpX6+4/jTgkEemFlJCISJeE71pZAmwF/tndmyIOSUR6IV2yEZGuGgz0A/oT1JSIiHSaakhEpEvMbBFwJ7APMNLdL4g4JBHphfS2XxHZZWZ2JsHbk39nZjHgBTOb7O5PRR2biPQuqiERERGRyKkNiYiIiEROCYmIiIhETgmJiIiIRE4JiYiIiEROCYmIiIhETgmJiIiIRE4JiYiIiETu/wMy+3Dc8IwqRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First I define the function I'll use and the necessary variables\n",
    "def u(x):\n",
    "    F = (np.exp(x/e)-np.exp(1/e))/(1-np.exp(1/e))\n",
    "    return F\n",
    "e = 1/100\n",
    "h = np.array([1/10, 1/20, 1/100])\n",
    "x1 = np.linspace(0, 1, 100)\n",
    "x2 = np.linspace(0, 1, 20)\n",
    "x3 = np.linspace(0, 1, 10)\n",
    "\n",
    "f1 = np.zeros(100)\n",
    "f1[0] = 1\n",
    "f2 = np.zeros(20)\n",
    "f2[0] = 1\n",
    "f3 = np.zeros(10)\n",
    "f3[0] = 1\n",
    "\n",
    "# The convection terms using the central difference scheme\n",
    "c1 = 1/(2*h[2])*np.diag(np.ones(99),k = 1)\n",
    "cc1 = -1/(2*h[2])*np.diag(np.ones(99), k = -1)\n",
    "\n",
    "c2 = 1/(2*h[1])*np.diag(np.ones(19), k = 1)\n",
    "cc2 = -1/(2*h[1])*np.diag(np.ones(19), k = -1)\n",
    "\n",
    "c3 = 1/(2*h[0])*np.diag(np.ones(9), k = 1)\n",
    "cc3 = -1/(2*h[0])*np.diag(np.ones(9), k = -1)\n",
    "\n",
    "# The diffusion terms used for the solution\n",
    "d1 = -e/(h[2]**2)*np.diag(np.ones(99), k = 1)\n",
    "dd1 = -e/(h[2]**2)*np.diag(np.ones(99), k = -1)\n",
    "ddd1 = 2*e/(h[2]**2)*np.diag(np.ones(100))\n",
    "\n",
    "d2 = -e/(h[1]**2)*np.diag(np.ones(19), k = 1)\n",
    "dd2 = -e/(h[1]**2)*np.diag(np.ones(19), k = -1)\n",
    "ddd2 = 2*e/(h[1]**2)*np.diag(np.ones(20))\n",
    "\n",
    "d3 = -e/(h[0]**2)*np.diag(np.ones(9), k = 1)\n",
    "dd3 = -e/(h[0]**2)*np.diag(np.ones(9), k = -1)\n",
    "ddd3 = 2*e/(h[0]**2)*np.diag(np.ones(10))\n",
    "\n",
    "# The global matrices\n",
    "A1 = c1 + cc1 + d1 + dd1 + ddd1\n",
    "A1[0,0] = 1\n",
    "A1[0,1] = 0\n",
    "A1[-1,-1] = 1\n",
    "A1[-1,-2] = 0\n",
    "\n",
    "A2 = c2 + cc2 + d2 + dd2 + ddd2\n",
    "A2[0,0] = 1\n",
    "A2[0,1] = 0\n",
    "A2[-1,-1] = 1\n",
    "A2[-1,-2] = 0\n",
    "\n",
    "A3 = c3 + cc3  + d3 + dd3 + ddd3\n",
    "A3[0,0] = 1\n",
    "A3[0,1] = 0\n",
    "A3[-1,-1] = 1\n",
    "A3[-1,-2] = 0\n",
    "\n",
    "# Lastly I plot the functions\n",
    "plt.figure()\n",
    "plt.plot(x1, u(x1), label = 'Exact solution')\n",
    "plt.plot(x1, np.linalg.solve(A1,f1), '.', label = 'N = 100')\n",
    "plt.plot(x2, np.linalg.solve(A2,f2), '.', label = 'N = 20')\n",
    "plt.plot(x3, np.linalg.solve(A3,f3), '.', label = 'N = 10')\n",
    "plt.title('Solution u(x) using the central difference scheme approximation for the convective term')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('u(x)')\n",
    "plt.legend(loc=3)\n",
    "\n",
    "print('The higher the N-value the more similar to the exact solution the approximation becomes')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0175130b",
   "metadata": {},
   "source": [
    "#### Assignment 5\n",
    "\n",
    "As in the previous assignment, but this time using the upwind scheme to approximate the convective term in $(1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9428b7e9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f80e9d832d0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAEWCAYAAACzATTWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA1yklEQVR4nO3deXxcdb3/8dcnk5k06b5RSpPSYgsUvJCWNKyFFixQuMBVuawqi9KLgoAXL3BV4KJyRX5sKiqyoyAgglC5ICCbVYQ2hbRIK1Kh0NBCS+lG0zbb5/fHOUkn05lkkk5mJpP38/HII3PmbJ+zfub7Pd9zjrk7IiIikh+Kch2AiIiIbKPELCIikkeUmEVERPKIErOIiEgeUWIWERHJI0rMIiIieaTHErOZjTMzN7Pibo5/upk9nem40pz3kWb2aJrDzjOzvXsghqwt/45uq2wzs7Fm9omZRbo5/jIz+0wG43nBzL6Sqen1JWY2zczezOL83jCz6T0w3VFm9icz22hm12d6+uE87jaz7/fEtHsTM7vFzC7PdRw9qdPEbGaHmNlLZrbezD42s7+Y2dRMBpEsMbj7fe5+ZCbn0wX/C1yT5rDXAd/NdAA9ufyZTkzZ5u7vufsAd2/OdSyyY9x9rrvv0RPTTpbI3H1vd3+hB2Y3G/gIGOTuF+/oxMzsTDP7846H1bslWw/ufq67f68H5pU358UOE7OZDQIeB34CDAPGAFcBW3s+tNwIf3QMdveX0xxlDjDDzEb3YFgieau31LT0sF2Bxd6NJzZp/fV+FshcDbS7p/wDqoB1HfQvAr4DvAusAn5JkNQAxgEOFIfdy4DPxI37P8C94ef3wmE/Cf8OBM4E/hw3/EHAfGB9+P+guH4vAN8D/gJsBJ4GRqSIud10w+8cmBB+vgK4PWG+HwEVYfe+wDpgz7hhngHOSDG/tuVMsV7OBN4O434HOD1ZnOE45wJvAWuBnwIW9osA14dxvgOcHz+PhHh+BbQAm8N1fUlcTGeE2+Ij4NsJ2/ky4J/AGuA3wLBurt+7gVvCdbYReBHYNex3FfCT8HMU2ARcG3aXAluAoUnWYYfbH/giwT66Bvg2CftiQqzHAIvD6bwPfDOu3wlALbAhXBdHpzn/A4CXCPabhcD0hH33+2H/T4DfA8OB+8L5zAfGxQ2/Z7juPgbeBE7q4Pg8C1gSxvQ28B9x/aYDdcC3wu29jHDf62w7xW3T8wj2x3fC784BloaxzQF2Cb//OfDbuHF/CDwLWGsccf2WAf8FLAq3/x3AKODJMI4/AkPjhn8I+IDgvPAnYO/w+9lAI9DQul4Tz0NACXATsCL8uwkoSVg/FxOc21YCZ6VYz3cnzOszaU770jD2XyVMbxLBvt4cTm9d3Hx+CvxfuC5eAT7VzX1jGHBXGNta4NG4fkm3Y0fnoXB51wGfjht2JMF5Zqew+18Jjp91BPv7PnHDVgCPAKsJjtObO1kP3w8/LwH+NW46xQT785TOjr3OzotpHrtXExz3m4EJ4fr5Wrh+NhKcFz4F/JXgeP4NEEu1Xdqm3WFPGBSupHuAWcQdEGH/s8MNuBswIFyxv0qRgJaROjG3GzbxBE+wE60lOMEWA6eG3cPjVtA/gd0JTuAvANd0M3E8BPxXQv+rgefCaS8Czk/o/2PghhTza1vOxGUF+ocba4+w32i2nVjaxRmO8zgwBBhLsAO3JoZzCZJJOUHi+mPi+kyIKXFbtMZ0W7iM+xLUikwK+18EvBxOvwT4BXB/N9fv3QQ77KHhtH4Ut50PB14PPx8UbtNX4votTLFvpdz+wF4EB1rr/G4AmkidmFcC08LPQ9l2gFcTnPxnEvxQGUP446yT+Y8hOIaOCcebGXaPjBt3KcHBOzjcjv8gOLkXE/zYvSsctj+wnCDhFgNTCE5Ce6dYlmPD6RpwGFAftzzTw/VwQ7heDiNIhHt0tp3itukzBMdmabh9PgpjKiGoZftTOGxZuExnAtPC4crj4khMzC8TJOMxBEnxVWByON3ngCsTzkED2ZYIa+P63U14Ak+27xNcgnoZ2IkgibwEfC9h/XyX4EfiMeH6G5piXbebV5rT/mEYd2max9HdBMmyOtz+9wEPdHPf+D/gQYJ9PAocFnecJd2OaZyH7gSujhv2POAP4ecp4bbcn6AgcUa4LUrC7oXAjeFy9AMO6WQ9tCbmK4D7Evb5v6dz7KVxXkzn2H0P2Dtc59Fw/cwhyJ17E5xHnyXIka3H9xnJ5t8ulk4HCH613E3wC68pnOmosN+zwNfiht2D4JdjMZlNzF8E5iXE9VfgzLgV9J24fl9r3SG6kTieAc5N6B8FFgCvA38gLKnG9b8auDPF/NqWM3FZCXbCdcDnSTg4E+MMxzkkrvs3wGXh5+doXxr6TOL67GQHbI2pPO67ecAp4eclwBFx/Ua3budurN+7CU8mYfcAgl/EFWwrFQ8nKKF/i2C/G0BQmv5xsv2lo+1PcODGz68/QckmVWJ+D/gPgmuF8d//ArgxxTgdzf9Sti8RPUV4cIbjxtdOXA88Gdd9HGGyAU4G5iaJ68rOjuNw2EeBC8PP0wmO5/4J+9TlnW2nuG16eFz/OwhrN+KGbyQs7RMkk48Jai5OjRtuOtsn5viS+8PAz+O6v05c6S5h+YaEcQ2OW4aOEvM/gWPi+h0FLIuLazPtz0mrgANSzLvdvNKYdgPQr4NtdSbJE1J8bd4xbEtCae8bBMdvC0l+ZKSxHTs6D30GeDuu31+AL4Wff074wySu/5sEPwgPJEjw6Z5P2tY1QSl1I1AWdt8HXJHOsZdkXm37RheO3e8m9Hfg4LjuBcClcd3XAzd1dqx2Wifu7kvc/Ux3Lwc+DexC8MuU8PO7cYO/S5BwRnU23S5KnE/rvMbEdX8Q97meYIfqjrUEv8DbuHsjwc7waeB6D9dwnIEECbZL3H0TwQF1LrDSzP7PzPbsYJRUy7gLwa/lVvGfuyLV9HcFfmdm68xsHUGibqb727ktPnf/hOCEvYu7bwZqCA7WQwmqT18CDg6/e7EbsbdbN+E6X9PBdD5PcMJ718xeNLMDw+8rCE62XZ3/rsC/t667cP0dQnBybPVh3OfNSbrjp7V/wrROB3ZOFpCZzTKzl8NGm+vC5RoRN8jacH20epdgfbVKup2S9SfhGA2HX0N4jLr7PILqdCM4mXckrfVhZhEzu8bM/mlmGwhOrCQsY0eSnb/il2+NuzfFdXflvNLZtFe7+5Y0pxWvo/0s3X2jAvjY3dd2FnfiduwkhueAUjPb38x2BSqB38XFd3FCfBXh/CqAdxPWdVrcfSnB+eg4MysDjgd+HTfPzo69jqQzfrJzbbrHc0pduljt7n9nW4KC4PrErnGDjCX4Ff4h29tEUKXVKn6HSUx0iRLn0zqv9zsZL5l2cZhZ4o67iKBKkrhhxgBXElyTud7MShLGmURQFdPp/Eg4UNz9KXefSbCx/05QndxVKwmqmVtVdDJ8Z+s70XJglrsPifvr5+7J1n9n67ddfGY2gKA6dEX41YsE1WmTCa6vvkhQ2qgmuIbYVSsT5ldGUCJPyt3nu/sJBFWQj7ItiSwnqBbuquUEv7rj111/d0+31X/itF5MmNYAd/9q4oDhPvowwV0Do9x9CPAEQWJsNdTM+sd1j2XbdoCOtxO034/aHaPhdIcTHqNmdh5BteUKgnYNmXAawXX/zxBUE45rnX2S+JJJdv5akWLYrups2p3F1p1jNK19Ixx2mJkNSdKvw+3YYcDuLQTHy6kE2+Zxd98YN8+rE+Irc/f7w35jUzSCS2c93B/O8wSCBnhL4+bZlWMvcV7pjN/V7ZSWzlpl72lmF5tZedhdQbACWlss3w98w8zGhwfu/wIPpvjlUwucYmZRM6sCTozrt5qgamW3FKE8AexuZqeZWbGZnUxw7fDxtJayvYXA3mZWaWb9CKqaE+d1WGuHmRnBj5E7gC8TnOi/F9e/BNiPoAo8mVrg0PDe28HAf8eNO8rMjg93/q0E10K7cwvQb4ALzWxMeLBd2snwH5J6XSdzC3B1+CsYMxtpZiekGLaz9QtwTHgbXoxgXb7i7q2/PF8EvkRwgDUQVBd9haCB0eouxNzqt8C/xs3vu6TY780sZsH944PDWpINbNsedwBnmdkRZlYUruuOajda3Uvwa/6osITXz8ymtx5TXfQ4wXHwxfA4iprZVDOblGTYGEEiXA00mdksINntd1eFyz2NoHHOQ3H9OtpOiX5NsH4qw2Pif8Phl5nZ7gQN3L5AcFnqEjOr7PLSb28gwXGzhuDH4P8m9O9sP78f+E64P48guOxxbwbiysS0PwTKw3WfjrT3DXdfSdCY7mdmNjQc9tCwd8rtmGYcvyaoBTydbSVXCAoc54alaTOz/mZ2rJkNJLhsthK4Jvy+n5kd3IX18ADBvv3VhHl29dhL3F8yeex2SWcl5o0EF+tfMbNNBAn5bwQtFSG42P8rgpLMOwTXB7+eYlqXE5Q41hJcL2xbge5eT9i6LawyOCB+RHdfQ3DSuJjgILyEoCXeR+ktZrtp/YPg5PxHgpZziffIvQqsN7P9w68uIKiyvTyswj6LYMedFvY/HnjB3ZP+0nb3ZwgaWSwiuN4Q/2OiKFymFQTVhIcRXJ/sqtsIWgIvAl4j+HHRROok/wOCk8Y6M/tmGtP/EUHbgqfNbCPBfrB/sgE7W7+hXxPUQHxM8KPm9Lh+LxFca24tHS8m2K+6U1rG3d8gaITya4KDfy3BdetUvggss6Bq9FyCZNJaFXsWQQOV9YStlNOY/3KCX/LfIkiSywlaHXf51oqw9HEkcArBPvMB2xoQJRv2AoIfbWsJSjBzEgb7IOy3guDa3LlhrVirjrZT4vyeJTjGHyZYz58i+CFeTHCC+6G7L3T3twjWxa+S1Dx11S8Jql3fJ9hPEm9xvAPYK9zPH00y/vcJLp0sImg/8mr4XSbs6LSfA94APjCzTs9zXdk3Ql8kuHb8d4Jr5xeF00m6HdMN2t1fIag124Ug+bd+X0PQ2vtmgn1uKcH1Yzx4HsFxBNeL3yM4Pk8OR+10PYQ/NP5K0GD0wbjvu3rstTsvZvLY7arW220kjpkdSdCo7d/SGPYV4Mvu/rceDyxNYenoFnfvNHFkm5ndTdDY5zu5jqUvs+DpV/eGbUeS9b8bbSeRnNCN7Um4+9MEJdB0hk1acswmMysFZhDEPIqglPO7DkcSEZG8pJdYFAYjuDywlqAqewnBNS0REellVJUtIiKSR1RiFhERySO6xtxFI0aM8HHjxuU6DBGRXmPBggUfufvIXMfRWygxd9G4ceOoqanJdRgiIr2GmSU+uVE6oKpsERGRPKLELCIikkeUmEVERPKIErOIiEgeUWIWERHJIwWbmM3sTjNbZWZJn2EdvuHkx2a21MwWmdmUbMcoIiKSqJBvl7qb4E0mv0zRfxYwMfzbH/g5Kd6YlAkLH76WFXOfYZdpM9m3ejosmwvjwhdUpfG59u8PU9OvH1VD96By7cqOh6+ohuXzko4/dVV/xtS8Ttn0WZTtNSHtOGqHjqZm7ZtUbdlC5Z6fTzvu1s+165dS8/ZTVO12FLu/D/UvPBnEcNRp28Xa2fiVg9OPu6P1V7946bY4Uq2LcF12af13cf2lE0day59ku3dl/I7iyMjyez1VVtbp+MniyOT+l9b+k7Aud3T/y/jnimqkcBX0IznNbBzBy7o/naTfLwhe13h/2P0mMD18hVhKVVVV3tX7mBc+fC1+5V0UN0NTBPzIjew7aBMtFgGgyJs7/LywJMp/jBpGoxlRd2774CP+paE56fDNRVGeH/efzFh2A5GWxnbjj1sB37m/hViTUxSBMYevo//wrZ3G8Xoswjk7j2ib/y8+/Jh9tzZ2Gnfr59f69WuLf1JdC1fc34I1g0Xgg9OP4eCme9ti7Wz81vlP3rJlh9bfhBXG+88NwcM4kq2L1nU5bMVPOHfUkLTWf1fX36Y1JZ3Gkc7yJ9vuXRm/ozi6sv+lXv6RNBQZsRbntg9Wpxw/WRxLd/GM7X/p7D+J63JH97/WdTCvtB9VWxq6PG6y+B759M8xg4oNC/hgaBVFZlRseJW1I6vZNGo/hpXFGDmwhJ0GljCkLErwWvncMbMF7l6V0yB6kUIuMXdmDMH7NVvVhd9tl5jNbDYwG2Ds2LFdntGKuc9Q0QwRB5ph+UcxJg/aiLc4hlNkdPj51ZJSGs1oMaMJWNAvxj5bNiYdvqW5kZK3HseKGimylnbjT3rPKW4Gw/AWZ8sHEQYOb+k0jgX9EuZfUszkzfWdxt36eUFJcdv4eywHWoAwhv6vvojtvS3WzsZvAl4tKWa/LZ3H3dH6K/+gBW/eFkeyddG6LhcMi6S9/ru6/rZ8EOk0jnSWP9l278r4HcXRlf0v9fKDA03W8fjJ4liwW1HG9r909p/WdbnlY9i6uox/jC+mcefu73+LSqLM3nkEDWFiv/2DVfzLlsa011+y+BpfvY/PR+YSpYnm8IpkhBYa/1HMVY1fZJh9wsstk3jVd2fMkFIOnjCcgyeM4IhJoxhQ0pdP+71DwV5jTkOyn5BJqw/c/VZ3r3L3qpEju/5UuV2mzaQpAs0WlJh3GdkEFiFSHKWouKTTz1MbWoi5E3En6lC1tTHl8MXREg474SsUR7cf/+9jwxK7OVYE/Ud7WnFUbW0k5oTzd6obmtOKu/VzdUNz2/j/qLBgrwtj2Pvzn2sXa2fjRx2mNrTs8PrrP7oFi2yLI9m6aF2X+zemv/67uv7SiSOd5U+23bsyfkdxdGX/63D5sU7HTxZHJve/dPaf4mgJU/eaxfsvDGP1ooHs+3gpe9V5t/e/Bf3709D6w8CMmtLSLq2/7T5HSzht/7H0K2qm2FqIWTMxCz73syZ+0O8e/iv6EA+VXcOvKhfzzbL/Y+Xf/sSFD9Ry9E1/YsG7a7t8DpPsUlV2FqqyQdeYdY1Z15h70zXmj56oZfWPfgQtLVBURNMph/FK5aZu7X+165dyzqvX0uhOtCjCbTvP7NYybPf5nuOhuQGKIoBBSxOYgbcEfxRBURG445EYrx/xS772p2JWrt/ChUdM5LwZE4gUZaeKW1XZXdOXE/OxwPnAMQSNvn7s7p22qOhuYhaR3qP+tdd476yz8cZGLBpl7F13UjZ5crenV7uqlpoPa6gaVUXlTpWZCTLZj7DS4fCHy4KEHZ+kLQKHf5sNUy/g8kf/xmO1Kzi1eiw/+Ny/ZCaWTigxd03BJmYzux+YDowAPgSuBKIA7n6LBa0hbgaOBuqBs9y904yrxCzSN9S/9hr18+ZTVj11h5JyJnQpsbcm7PgkHYnB0dfA5jUwbhpXvlrGr15+lz/+52HsNnJAj8evxNw1BZuYe4oSs4hkU+2qWs55+hwamhuIRWLcduRt6Ze6UyTptf/+Ww68dyOzPj2aG09Oc1o7QIm5a/py4y8RkbxX82ENDc0NtNBCY0sjNR92oWBQUQ3TLg5Kys0N4M3Q3MDQVa9wxoHjeKz2fZau+qTngpduUWIWEcljVaOqiEViRCxCtChK1ahuFDzHTQuqsy0S/B83jdmH7ka/aIQfP/tW5oOWHaIb2kRE8ljlTpXcduRtO9Z4rKIazpjTrgX/cOCMg8Zxy4v/5PzDJ7D7qIGZDl26SdeYu0jXmEWkUKzd1MAhP3yOoz69MzecVNlj89E15q5RVbaISB81tH+M6XvsRM0yPXQknygxi4j0NcvnwdzrYfk8Jo0eyHsf17NxS2Ouo5KQrjGLiPQly+dte2pYJMb+h90NwJsfbKRq3LDcxiaASswiIn3Lsrntbp3aY8tCAJas3JDjwKSVErOISF+ScOvUwD2mM6QsymIl5ryhqmwRkb4k4dYpq6hm0s4vs3jlxlxHJiElZhGRPNVjz+uuqA7+QpNGD+LX896lucWz9sYpSU2JWUQkD7W94aqhAYvFdvgNVx3Za5dBbGlsYdmaTXwqCy+1kI7pGrOISB6qnzcfb2iAlha8sZH6efN7bF6TRgdP/VIDsPygxCwikofKqqdisRhEIlg0Sln11B6b14SdBlBcZCxeocScD1SVLSKSh8omT2bsXXdm5Z3QJcURJuw0QCXmPKHELCKSp8omT+7RhBxv0uhB/PWfa7IyL+mYqrJFRIS9Rg/igw1b+HhTQ65D6fOUmEVE+rLwudnVxUsBNQDLB6rKFhHpq+Kem71PJMoUu4wlKydx8IQRuY6sT1NiFhHpq+Kem23NcETpW3o0Zx5QVbaISF+V8Nzsj0ZMZemqT3IdVZ+nErOISF+V8NzstX+Jsva9tbmOqs9TYhYR6cvinps9pOwN1tc35jggUVW2iIgAMKg0ysatTbS0eK5D6dOUmEVEBIDBpVHcYeOWplyH0qcpMYuIFLjaVbXc/vrt1K6q7XC4waVRANZvVnV2Lukas4hIAatdVcs5T59DQ3MDsUiM2468jcqdKpMO25qY121uYCxlWYxS4qnELCJSwGo+rKGhuYEWWmhsaaTmw5qUww4pU4k5Hygxi4gUsKpRVcQiMSIWIVoUpWpUVcphVZWdH1SVLSJSwCp3quS2I2+j5sMaqkZVpazGBiXmfKHELCJS4Cp3quwwIbdqu8ase5lzqmCrss3saDN708yWmtllSfoPNrPfm9lCM3vDzM7KRZwiIvmiXzRCrLiIDSox51RBJmYziwA/BWYBewGnmtleCYOdByx2932B6cD1ZhbLaqAiInlmcGlUVdk5VpCJGagGlrr72+7eADwAnJAwjAMDzcyAAcDHgO6qF5E+bYgSc84VamIeAyyP664Lv4t3MzAJWAG8Dlzo7i3JJmZms82sxsxqVq9e3RPxiojkBZWYc69QE7Ml+S7x4a9HAbXALkAlcLOZDUo2MXe/1d2r3L1q5MiRmYxTRCSvDC6NqvFXjhVqYq4DKuK6ywlKxvHOAh7xwFLgHWDPLMUnIpKXVGLOvUJNzPOBiWY2PmzQdQowJ2GY94AjAMxsFLAH8HZWoxQRyTODSqNqlZ1jBXkfs7s3mdn5wFNABLjT3d8ws3PD/rcA3wPuNrPXCaq+L3X3j3IWtIhIHhhSFrz6sbnFiRQluyooPa0gEzOAuz8BPJHw3S1xn1cAR2Y7LhGRfNb6kJENmxsZ2l93kOZCoVZli4hIN2x7w5Sqs3NFiVlERNroedm5V7BV2SIi0kXL57HHW08zxfqzfnN1rqPps5SYRUQEls+De45nTHMD98Ui1NTtDrsfm+uo+iRVZYuICCybC80NmDcTpYkBK1/OdUR9lhKziIjAuGkQieEWoZFi/tm/MtcR9VmqyhYREaiohjPmYMvmctbTxfxLdFKuI+qzlJhFRCRQUQ0V1bwz94+MVavsnFFVtoiItDOkNKbbpXJIiVlERNrRG6ZyS4lZRETaGaQ3TOWUErOIiLQzWG+YyiklZhERaUfvZM4tJWYREWlnSFmUTQ3NNDa35DqUPkmJWURE2tGLLHJLiVlERNpRYs4tJWYREWlHiTm3lJhFRKSdQUrMOaXELCIi7QwpCxKzbpnKDSVmERFpp7UqW0//yg0lZhERaUfXmHNLiVlERNqJRoooi0WUmHNEiVlERLajp3/ljhKziIhsR4k5d5SYRURkO4NLo6xX46+cUGIWEZHtqMScO0rMIiKynQH9ivlka1Ouw+iTlJhFRGQ7JcVFNOjtUjmhxCwiItuJRor02sccUWIWEZHtxCJFNDYpMedCwSZmMzvazN40s6VmdlmKYaabWa2ZvWFmL2Y7RhGRfBVVVXbOFOc6gJ5gZhHgp8BMoA6Yb2Zz3H1x3DBDgJ8BR7v7e2a2U06CFRHJQ0FVtuPumFmuw+lTCrXEXA0sdfe33b0BeAA4IWGY04BH3P09AHdfleUYRUTyViwSJOPGZs9xJH1PoSbmMcDyuO668Lt4uwNDzewFM1tgZl9KNTEzm21mNWZWs3r16h4IV0Qkv8SKg/SgBmDZV6iJOVm9S+LPvmJgP+BY4CjgcjPbPdnE3P1Wd69y96qRI0dmNlIRkTwUjQTpoUENwLKuIK8xE5SQK+K6y4EVSYb5yN03AZvM7E/AvsA/shOiiEj+ak3MKjFnX6GWmOcDE81svJnFgFOAOQnDPAZMM7NiMysD9geWZDlOEZG8FGstMSsxZ11BlpjdvcnMzgeeAiLAne7+hpmdG/a/xd2XmNkfgEVAC3C7u/8td1GLiOSPbdeY1fgr2woyMQO4+xPAEwnf3ZLQ/f+A/5fNuEREegNdY86dQq3KFhGRHRBtu11KiTnblJhFRGQ70WJdY84VJWYREdlOSWurbFVlZ50Ss4iIbEcl5txRYhYRke3oPubcUWIWEZHttDb+amjS7VLZpsQsIiLbKdGzsnOm19zHbGZDgV2AzcAyd9feIiLSQ3Qfc+7kdWI2s8HAecCpQAxYDfQDRpnZy8DP3P35HIYoIlKQdI05d/I6MQO/BX4JTHP3dfE9zGw/4Itmtpu735GL4ERECpUSc+7kdWJ295kd9FsALMhiOCIifUas7XYpNf7Ktl7R+MvMvpzQHTGzK3MVj4hIoYvpGnPO9IrEDBxhZk+Y2Wgz+zTwMjAw10GJiBSq6Mr5fC3yGMPX1uY6lD4nr6uyW7n7aWZ2MvA6UA+c6u5/yXFYIiKFafk8in/1b/xn8Vb89UehqgIqqnMdVZ/RK0rMZjYRuBB4GFhG0OirLKdBiYgUqmVzobmBYmuhqKUp6Jas6RWJGfg9cLm7/wdwGPAWMD+3IYmIFKhx0yASo8mLaLHioFuyprck5mp3fxbAA9cD/5bbkEREClRFNZwxh58XncKdE36kauwsy+vEbGaHALj7hsR+7v6WmQ0KG4OJiEgmVVRzX/RE3i7ZO9eR9Dn53vjr82Z2LfAHgnuWW5/8NQGYAewKXJy78ERECle02PSAkRzI68Ts7t8In5F9IvDvwM4Ez8peAtyiltkiIj0nGiliqxJz1uV1YgZw97VmtgvwTvgH4MARgBKziEgPiUWKaNQDRrIu7xNz6JO4z/2AfyUoNYuISA+JFRepKjsHekViDlthtzGz64A5OQpHRKRPiEaKaNSzsrMur1tld6AM2C3XQYiIFLJoxPSs7BzoFSVmM3ud4LoyQAQYCXw3dxGJiBS+aKSIjY1NuQ6jz+kViZngmnKrJuBDd9feIiLSg0qKi/hY15izrlckZnd/N9cxiIj0NcE1ZiXmbOut15hFRKSHRSNFusacA0rMIiKSlFpl54YSs4iIJBUrLqJBVdlZp8QsIiJJxSJ6VnYuFGxiNrOjzexNM1tqZpd1MNxUM2s2sxOzGZ+ISL7TNebcKMjEbGYR4KfALGAv4FQz2yvFcD8EnspuhCIi+S+qR3LmREEmZqAaWOrub7t7A/AAcEKS4b4OPAysymZwIiK9QSxs/OWuBmDZVKiJeQywPK67LvyujZmNAT4L3NLZxMxstpnVmFnN6tWrMxqoiEi+ihUHKUINwLKrUBOzJfku8SffTcCl7t7c2cTc/VZ3r3L3qpEjR2YiPhGRvBeNBKdS3TKVXb3iyV/dUAdUxHWXAysShqkCHjAzgBHAMWbW5O6PZiVCEZE8F40EZbfGphYoyXEwfUihJub5wEQzGw+8D5wCnBY/gLuPb/1sZncDjyspi4hs01qVrQZg2VWQidndm8zsfILW1hHgTnd/w8zODft3el1ZRKSvay0xb9UtU1lVkIkZwN2fAJ5I+C5pQnb3M7MRk4hIbxKLqMScC4Xa+EtERHZQ2zVmNf7KKiVmERFJSteYc0OJWUREkmq9XUrXmLNLiVlERJLSNebcUGIWEZGkoqrKzgklZhERSUol5txQYhYRkaRaW2Xr1Y/ZpcQsIiJJxYqDxl8Nul0qq5SYRUQkqVgkAoTPypasUWIWEZGkosWtb5dSYs4mJWYREUmq7RqzEnNWKTGLiEhSavyVG0rMIiKSVEmxnpWdC0rMIiKSVFT3MeeEErOIiCQVKTKKTFXZ2Vaw72POpsbGRurq6tiyZUuuQ5EU+vXrR3l5OdFoNNehiPQq0UiRSsxZpsScAXV1dQwcOJBx48ZhZrkORxK4O2vWrKGuro7x48fnOhyRXiVWXKRW2VmmquwM2LJlC8OHD1dSzlNmxvDhw1WjIdINMZWYs06JOUOUlPObto9I90QjRbrGnGVKzCIiklK02HS7VJYpMReISCRCZWVl298111yTsWnX1tbyxBNPZGRaAwYM6HSYm266ifr6+rbuY445hnXr1mVk/iLSNbGIrjFnmxp/FYjS0lJqa2t7ZNq1tbXU1NRwzDHH9Mj0E91000184QtfoKysDCBjPwpEpOuikSK9xCLLlJgz7Krfv8HiFRsyOs29dhnElcft3eXx1q9fT3V1NXPmzGGPPfbg1FNP5fDDD+ecc87hq1/9KvPnz2fz5s2ceOKJXHXVVQDMnz+fCy+8kE2bNlFSUsIzzzzDFVdcwebNm/nzn//Mf//3f3PyySe3zeONN97grLPOoqGhgZaWFh5++GEmTpzIDTfcwJ133gnAV77yFS666KJ2sb3wwgtcd911PP744wCcf/75VFVVsWHDBlasWMGMGTMYMWIEzz//POPGjaOmpoYRI0Ykne6yZcuYNWsWhxxyCC+99BJjxozhscceo7S0tDurW0TiqFV29ikxF4jNmzdTWVnZ1t2aQG+++WbOPPNMLrzwQtauXcs555wDwNVXX82wYcNobm7miCOOYNGiRey5556cfPLJPPjgg0ydOpUNGzZQVlbGd7/7XWpqarj55pu3m+8tt9zChRdeyOmnn05DQwPNzc0sWLCAu+66i1deeQV3Z//99+ewww5j8uTJnS7HBRdcwA033MDzzz/PiBEj2vVLNd2hQ4fy1ltvcf/993Pbbbdx0kkn8fDDD/OFL3xhx1aqiOg+5hxQYs6w7pRsMyFVVfbMmTN56KGHOO+881i4cGHb97/5zW+49dZbaWpqYuXKlSxevBgzY/To0UydOhWAQYMGdTrfAw88kKuvvpq6ujo+97nPMXHiRP785z/z2c9+lv79+wPwuc99jrlz56aVmDuSarrHH38848ePb/thst9++7Fs2bIdmpeIBGKRIhqb1Pgrm9T4q8C1tLSwZMkSSktL+fjjjwF45513uO6663j22WdZtGgRxx57LFu2bMHdu3xb0WmnncacOXMoLS3lqKOO4rnnnsO984O4uLiYlpZtv8LTuce4o+mWlJS0fY5EIjQ1NXU6PRHpXFRV2VmnxFzgbrzxRiZNmsT999/P2WefTWNjIxs2bKB///4MHjyYDz/8kCeffBKAPffckxUrVjB//nwANm7cSFNTEwMHDmTjxo1Jp//222+z2267ccEFF3D88cezaNEiDj30UB599FHq6+vZtGkTv/vd75g2bVq78XbddVcWL17M1q1bWb9+Pc8++2xbv1TzS2e6IpJZsYjpPuYsU1V2gUi8xnz00Udz9tlnc/vttzNv3jwGDhzIoYceyve//32uuuoqJk+ezN57781uu+3GwQcfDEAsFuPBBx/k61//Ops3b6a0tJQ//vGPzJgxg2uuuYbKysrtGn89+OCD3HvvvUSjUXbeeWeuuOIKhg0bxplnnkl1dTUQNNJKrMauqKjgpJNOYp999mHixInt+s+ePZtZs2YxevRonn/++bbvp0yZknS6qrYW6Tm6xpx9lk61o2xTVVXlNTU17b5bsmQJkyZNylFEki5tJ5Guu/CB11i4fB0v/NeMbk/DzBa4e1UGwypoqsoWEZGUghKzCnDZVLCJ2cyONrM3zWypmV2WpP/pZrYo/HvJzPbNRZwiIvksGiliq64xZ1VBJmYziwA/BWYBewGnmtleCYO9Axzm7vsA3wNuzW6UIiL5LxYxXWPOsoJMzEA1sNTd33b3BuAB4IT4Adz9JXdfG3a+DJRnOUYRkbwXK1bjr2wr1MQ8Blge110XfpfKl4EnU/U0s9lmVmNmNatXr85QiCIi+U+tsrOvUBNzsqdkJG29YGYzCBLzpakm5u63unuVu1eNHDkyQyGKiOS/1sZfLS1qAJYthZqY64CKuO5yYEXiQGa2D3A7cIK7r8lSbD3CzLj44ovbuq+77jr+53/+Z4en++1vf5uKiortXte4detWTj75ZCZMmMD+++/f7l7ie+65h4kTJzJx4kTuueeeHY5BRHInVhykicYWlZqzpVAT83xgopmNN7MYcAowJ34AMxsLPAJ80d3/kfUIl8+DudcH/zOgpKSERx55hI8++igj02t13HHHMW/e9jHecccdDB06lKVLl/KNb3yDSy8NKhw+/vhjrrrqKl555RXmzZvHVVddxdq1a7cbX0R6h1gkTMy6ZSprCjIxu3sTcD7wFLAE+I27v2Fm55rZueFgVwDDgZ+ZWa2Z1aSYXOYtnwf3HA/PXR38z0ByLi4uZvbs2dx4440ZCHCbAw44gNGjR2/3/WOPPcYZZ5wBwIknnsizzz6Lu/PUU08xc+ZMhg0bxtChQ5k5cyZ/+MMfMhqTiGRPNBJcGdRjObOnYB/J6e5PAE8kfHdL3OevAF/JdlwALJsLzQ3gzcH/ZXOhonqHJ3veeeexzz77cMkll6Qc5vnnn+cb3/jGdt+XlZXx0ksvpT2v999/n4qK4GpBcXExgwcPZs2aNe2+BygvL+f999/vwlKISD6JtlZlqwFY1hRsYs5r46ZBJBYk5Ugs6M6AQYMG8aUvfYkf//jHlJaWJh1mxowZSV8P2VXJHuVqZim/F5HeKRpWZavEnD1KzLlQUQ1nzAlKyuOmZaS03Oqiiy5iypQpnHXWWUn7Z6rEXF5ezvLlyykvL6epqYn169czbNgwysvLeeGFF9qGq6urY/r06V1dDBHJEyUqMWedEnOuVFRnNCG3GjZsGCeddBJ33HEHZ5999nb9M1ViPv7447nnnns48MAD+e1vf8vhhx+OmXHUUUfxrW99q63B19NPP80PfvCDHZ6fiORGW4lZiTlrCrLxV1938cUXZ6x19iWXXEJ5eTn19fWUl5e33YL15S9/mTVr1jBhwgRuuOEGrrnmGiD4YXD55ZczdepUpk6d2vYaSBHpnVoTc2OTWmVni1772EV67WPvpe0k0nUv/mM1Z9w5j4e/ehD77Tq0W9PQax+7RiVmERFJqfV2KV1jzh4lZhERSSmmVtlZp8QsIiIptV1jVok5a5SYRUQkpZhul8o6JWYREUlp2+1SaiicLUrMIiKSkq4xZ58Sc4Hoidc+1tfXc+yxx7Lnnnuy9957c9lll7X16+i1jyJSOKLFapWdbUrMOVK7qpbbX7+d2lW1GZleT7328Zvf/CZ///vfee211/jLX/7Ck08+CaR+7aOIFJaYGn9lnRJzDtSuquWcp8/hJ6/+hHOePicjybknXvtYVlbGjBkzAIjFYkyZMoW6ujog9WsfRaSwtL5dSlXZ2aPEnAM1H9bQ0NxACy00tjRS82FmXgV93nnncd9997F+/fqUwzz//PNUVlZu93fQQQd1OO1169bx+9//niOOOAJI/dpHESksMT0rO+v0EoscqBpVRSwSo7GlkWhRlKpRmXlSXU+99rGpqYlTTz2VCy64gN122w1I/dpHESkselZ29ikx50DlTpXcduRt1HxYQ9WoKip3qszYtHvitY+zZ89m4sSJXHTRRW3fpXrto4gUlkiRESkyXWPOIiXmHKncqTKjCblVpl/7+J3vfIf169dz++23t/s+1WsfRaTwRCNKzNmka8wFKFOvfayrq+Pqq69m8eLFTJkyhcrKyrYEneq1jyJSeKKRIraq8VfWqMRcID755JO2z6NGjaK+vn6Hp1leXp6ypXW/fv146KGHdngeIpL/YpEilZizSCVmERHpUKxYiTmblJhFRKRD0UgRjXpWdtYoMYuISIeiEdMDRrJIiVlERDoUjRTpASNZpMQsIiIdKtE15qxSYhYRkQ5F1So7q5SYC0RPvPYR4Nvf/jYVFRUMGDCg3fd67aNI3xGNFOkacxYpMedI/Wuv8dEvbqX+tdcyMr2eeu3jcccdx7x587b7Xq99FOk7osVFNKhVdtYoMedA/Wuv8d5ZZ7P6Rz/ivbPOzkhy7onXPgIccMABjB49ervv9dpHkb4jFimiUSXmrNGTv3Kgft58vKEBWlrwxkbq582nbPLkHZ7ueeedxz777MMll1yScpjuvMQimVSvfRwxYkTXAxeRvBYr1rOys0mJOQfKqqdisRje2IhFo5RVT83IdHvqtY/J6LWPIn2HbpfKroJNzGZ2NPAjIALc7u7XJPS3sP8xQD1wpru/mo3YyiZPZuxddwYl5eqpGSktt+qJ1z4mo9c+ivQdUVVlZ1VBJmYziwA/BWYCdcB8M5vj7ovjBpsFTAz/9gd+Hv7PirLJkzOakFtl+rWPqei1jyJ9x9il9zJxyUssfHgB+34+9aUyyYxCbfxVDSx197fdvQF4ADghYZgTgF964GVgiJlt38qpF8rUax8BLrnkEsrLy6mvr6e8vLztFiy99lGkb1j48LUcet/jHDz/Y/zKu1j48LW5DqngFWSJGRgDLI/rrmP70nCyYcYAKxMnZmazgdkAY8eOzWigmdITr30EuPbaa7n22u0PRL32UaRvWDH3GSqaIeJAMyyf+4xKzT2sUEvMyepUE1srpTNM8KX7re5e5e5VI0eO3OHgRER6i12mzaQpAs0GTZGgW3pWoZaY64CKuO5yYEU3hhER6dP2/fwlLCQoKe8ybaZKy1lQqIl5PjDRzMYD7wOnAKclDDMHON/MHiCo5l7v7ttVY6fL3dX4KY/p4Sci3bfv5y9RQs6igkzM7t5kZucDTxHcLnWnu79hZueG/W8BniC4VWopwe1Sye8vSkO/fv1Ys2YNw4cPV3LOQ+7OmjVr6NevX65DERHplKkk0TVVVVVeU1PT7rvGxkbq6urYsmVLjqKSzvTr14/y8nKi0WiuQxHpc8xsgbtX5TqO3qIgS8zZFo1GGT9+fK7DEBGRAlCorbJFRER6JSVmERGRPKLELCIikkfU+KuLzGw18G43Rx8BZOZZmb2Hlrnw9bXlBS1zV+3q7no6U5qUmLPIzGr6WstELXPh62vLC1pm6VmqyhYREckjSswiIiJ5RIk5u27NdQA5oGUufH1teUHLLD1I15hFRETyiErMIiIieUSJWUREJI8oMWeYmR1tZm+a2VIzuyxJfzOzH4f9F5nZlFzEmUlpLPPp4bIuMrOXzGzfXMSZSZ0tc9xwU82s2cxOzGZ8PSGdZTaz6WZWa2ZvmNmL2Y4x09LYtweb2e/NbGG4zN1+S10+MLM7zWyVmf0tRf+CO3/lJXfXX4b+CF4x+U9gNyAGLAT2ShjmGOBJwIADgFdyHXcWlvkgYGj4eVZfWOa44Z4jeMXoibmOOwvbeQiwGBgbdu+U67izsMzfAn4Yfh4JfAzEch37DizzocAU4G8p+hfU+Stf/1RizqxqYKm7v+3uDcADwAkJw5wA/NIDLwNDzGx0tgPNoE6X2d1fcve1YefLQHmWY8y0dLYzwNeBh4FV2Qyuh6SzzKcBj7j7ewDu3tuXO51ldmCgBS9iH0CQmJuyG2bmuPufCJYhlUI7f+UlJebMGgMsj+uuC7/r6jC9SVeX58sEv7h7s06X2czGAJ8FbsliXD0pne28OzDUzF4wswVm9qWsRdcz0lnmm4FJwArgdeBCd2/JTng5UWjnr7yk9zFnliX5LvF+tHSG6U3SXh4zm0GQmA/p0Yh6XjrLfBNwqbs3B4WpXi+dZS4G9gOOAEqBv5rZy+7+j54Oroeks8xHAbXA4cCngGfMbK67b+jh2HKl0M5feUmJObPqgIq47nKCX9JdHaY3SWt5zGwf4HZglruvyVJsPSWdZa4CHgiT8gjgGDNrcvdHsxJh5qW7b3/k7puATWb2J2BfoLcm5nSW+SzgGg8uwC41s3eAPYF52Qkx6wrt/JWXVJWdWfOBiWY23sxiwCnAnIRh5gBfCls3HgCsd/eV2Q40gzpdZjMbCzwCfLEXl57idbrM7j7e3ce5+zjgt8DXenFShvT27ceAaWZWbGZlwP7AkizHmUnpLPN7BDUEmNkoYA/g7axGmV2Fdv7KSyoxZ5C7N5nZ+cBTBC0673T3N8zs3LD/LQQtdI8BlgL1BL+4e600l/kKYDjws7AE2eS9+C01aS5zQUlnmd19iZn9AVgEtAC3u3vS2256gzS38/eAu83sdYJq3kvdvde+DtLM7gemAyPMrA64EohCYZ6/8pUeySkiIpJHVJUtIiKSR5SYRURE8ogSs4iISB5RYhYREckjSswiIiJ5RIlZREQkjygxi4iI5BElZpFeLHzf8yIz62dm/cN3An8613GJSPfpASMivZyZfR/oR/DiiDp3/0GOQxKRHaDELNLLhc9xng9sAQ5y9+YchyQiO0BV2SK93zBgADCQoOQsIr2YSswivZyZzQEeAMYDo939/ByHJCI7QG+XEunFzOxLBG/r+rWZRYCXzOxwd38u17GJSPeoxCwiIpJHdI1ZREQkjygxi4iI5BElZhERkTyixCwiIpJHlJhFRETyiBKziIhIHlFiFhERySP/Hzs3I5bxgs36AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First I again define the necessary variables and the function\n",
    "def u(x):\n",
    "    F = (np.exp(x/e)-np.exp(1/e))/(1-np.exp(1/e))\n",
    "    return F\n",
    "e = 1/100\n",
    "h = np.array([1/10, 1/20, 1/100])\n",
    "x1 = np.linspace(0, 1, 100)\n",
    "x2 = np.linspace(0, 1, 20)\n",
    "x3 = np.linspace(0, 1, 10)\n",
    "\n",
    "f1 = np.zeros(100)\n",
    "f1[0] = 1\n",
    "f2 = np.zeros(20)\n",
    "f2[0] = 1\n",
    "f3 = np.zeros(10)\n",
    "f3[0] = 1\n",
    "\n",
    "# The convection terms using the upwind scheme\n",
    "c1 = 1/h[2]*np.diag(np.ones(100), k = 0)\n",
    "cc1 = -1/h[2]*np.diag(np.ones(99), k = -1)\n",
    "\n",
    "c2 = 1/h[1]*np.diag(np.ones(20), k = 0)\n",
    "cc2 = -1/h[1]*np.diag(np.ones(19), k = -1)\n",
    "\n",
    "c3 = 1/h[0]*np.diag(np.ones(10), k = 0)\n",
    "cc3 = -1/h[0]*np.diag(np.ones(9), k = -1)\n",
    "\n",
    "# The diffusion terms for the solution\n",
    "d1 = -e/(h[2]**2)*np.diag(np.ones(99), k = 1)\n",
    "dd1 = -e/(h[2]**2)*np.diag(np.ones(99), k = -1)\n",
    "ddd1 = 2*e/(h[2]**2)*np.diag(np.ones(100))\n",
    "\n",
    "d2 = -e/(h[1]**2)*np.diag(np.ones(19), k = 1)\n",
    "dd2 = -e/(h[1]**2)*np.diag(np.ones(19), k = -1)\n",
    "ddd2 = 2*e/(h[1]**2)*np.diag(np.ones(20))\n",
    "\n",
    "d3 = -e/(h[0]**2)*np.diag(np.ones(9), k = 1)\n",
    "dd3 = -e/(h[0]**2)*np.diag(np.ones(9), k = -1)\n",
    "ddd3 = 2*e/(h[0]**2)*np.diag(np.ones(10))\n",
    "\n",
    "# The global matrices\n",
    "A1 = c1 + cc1 + d1 + dd1 + ddd1\n",
    "A1[0,0] = 1\n",
    "A1[0,1] = 0\n",
    "A1[-1,-1] = 1\n",
    "A1[-1,-2] = 0\n",
    "\n",
    "A2 = c2 + cc2 + d2 + dd2 + ddd2\n",
    "A2[0,0] = 1\n",
    "A2[0,1] = 0\n",
    "A2[-1,-1] = 1\n",
    "A2[-1,-2] = 0\n",
    "\n",
    "A3 = c3 + cc3  + d3 + dd3 + ddd3\n",
    "A3[0,0] = 1\n",
    "A3[0,1] = 0\n",
    "A3[-1,-1] = 1\n",
    "A3[-1,-2] = 0\n",
    "\n",
    "# Lastly I plot the figures\n",
    "plt.figure()\n",
    "plt.plot(x1, u(x1), label = 'Exact solution')\n",
    "plt.plot(x1, np.linalg.solve(A1,f1), '.', label = 'N = 100')\n",
    "plt.plot(x2, np.linalg.solve(A2,f2), '.', label = 'N = 20')\n",
    "plt.plot(x3, np.linalg.solve(A3,f3), '.', label = 'N = 10')\n",
    "plt.title('Solution u(x) using the upwind scheme approximation for the convective term')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('u(x)')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75044cd0",
   "metadata": {},
   "source": [
    "## Section 3/: Other Models "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ae834db",
   "metadata": {},
   "source": [
    "More later. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b81ec120",
   "metadata": {},
   "source": [
    "## References \n",
    "\n",
    "To be extended. \n",
    "\n",
    "1. [convection-diffusion equation on wiki](https://en.wikipedia.org/wiki/Convection–diffusion_equation)\n",
    "2. [CFD Python: 12 steps to Navier-Stokes - Lorena A. Barba](https://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/) \n",
    "3. [A CFD Tutorial in Julia: Introduction to Laminar Boundary-Layer Theory](https://www.mdpi.com/2311-5521/6/6/207) \n",
    "4. [2D Navier-Stokes Solved by pypde](https://pypde.readthedocs.io/en/latest/pages/example_code.html#navier-stokes-2d-parabolic)\n",
    "5. add references to [NGSolve](https://ngsolve.org), [FENICS](https://fenicsproject.org) and [FireDrake](https://www.firedrakeproject.org) among others; "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16e30ebc",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fb5a165",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f5604af",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}