{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Computing with Python\n",
    "# Lecture  5: Partial Differential Equations\n",
    "### Ramses van Zon\n",
    "#### 26 November 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style> h1, h2 { color: #2a7bbd ! important; } h1 { font-size: 2.5em !important; }  </style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%html\n",
    "<style> h1, h2 { color: #2a7bbd ! important; } h1 { font-size: 2.5em !important; }  </style>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrixprint(F):\n",
    "    for i in range(len(F)):\n",
    "        print ('|',end='')\n",
    "        for el in F[i]:\n",
    "            print(\"%010s \"%str(el),end='')\n",
    "        print('|')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Today’s class\n",
    "\n",
    "Today we will discuss the following topics:\n",
    "\n",
    "  - Basic approaches to solving PDEs.\n",
    "  - How to discretize equations.\n",
    "  - How to implement boundary conditions.\n",
    "  - Implicit versus explicit approaches."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Partial Differential Equations\n",
    "\n",
    "Partial differential equations (PDEs) are differential equations which\n",
    "contain **derivatives of more than one variable**, e.g., two spatial coordinates, or space and time coordinates.\n",
    "\n",
    "Second order PDEs are of the general form:\n",
    "\n",
    "$$\n",
    " A\\frac{\\partial^2T}{\\partial x^2} + B\\frac{\\partial^2 T}{\\partial x\\partial y}  + C\\frac{\\partial^2 T}{\\partial y^2}\n",
    "=f\\left(x,y,T,\\frac{\\partial T}{\\partial x}, \\frac{\\partial T}{\\partial y}\\right).\n",
    "$$\n",
    "\n",
    "For different combinations of A, B and C, three classes of PDEs show up repeatedly in\n",
    "physical systems.\n",
    "\n",
    "  - If $B^2 - 4AC < 0$, the equation is called elliptic.\n",
    "  - If $B^2 - 4AC = 0$, the equation is called parabolic (diffusive).\n",
    "  - If $B^2 - 4AC > 0$, the equation is called hyperbolic (wavelike)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How do we solve these problems?\n",
    "\n",
    "Let’s looking at a parabolic equation, in particular, let us\n",
    "look at the ****heat equation****.\n",
    "\n",
    "$$\\frac{\\partial T}{\\partial t}=k\\frac{\\partial^2 T}{\\partial x^2}$$\n",
    "\n",
    "where $T$ is the temperature and $k$ is the thermal diffusivity.  Temperature varies along the spatial direction ($x$) and changes in time ($t$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How do we solve this equation?\n",
    "\n",
    "By **discretizing** in both space and time,\n",
    "and marching an initial condition forward in time.\n",
    "\n",
    "*\"Be wise, discretize!\"* (Mark Kac)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating derivatives\n",
    "\n",
    "We need to calculate the second spatial derivatives. How best to do that?\n",
    "\n",
    "We **discretize** the $x$ domain, \n",
    "\n",
    "o--------o--------o--------o--------o \n",
    "\n",
    "$ j-2\\:\\,j-1\\:\\:\\: j\\:\\:\\:\\:\\:\\:  j+1 \\:\\:   j+2$\n",
    "\n",
    "and examine the **Taylor expansion** of the\n",
    "function $T$, centered around three different points:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "T(x_{j-1}) &= T(x_j) - \\Delta x \\frac{\\partial T(x_j)}{\\partial x}\n",
    "+ \\frac{(\\Delta x)^2}{2}\\frac{\\partial^2T(x_j)}{\\partial x^2}+\\mathcal O(\\Delta x^3)\\\\\n",
    "T(x_j) &= T(x_j) \\\\\n",
    "T(x_{j+1}) &= T(x_j) + \\Delta x \\frac{\\partial T(x_j)}{\\partial x}\n",
    "+ \\frac{(\\Delta x)^2}{2}\\frac{\\partial^2 T(x_j)}{\\partial x^2}+\\mathcal O(\\Delta x^3)\\end{aligned}\n",
    "$$\n",
    "\n",
    "Where $\\Delta x=x_j-x_{j-1}$. \n",
    "\n",
    "Writing $T'=\\partial T/\\partial x$ and $T''=\\partial^2 T/\\partial x^2$, we get\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "T(x_{j-1})\\\\\n",
    "T(x_{j})\\\\\n",
    "T(x_{j+1})\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "1 & -\\Delta x & \\Delta x^2 \\\\\n",
    "1 & 0  & 0 \\\\\n",
    "1 & +\\Delta x & \\Delta x^2 \\\\\n",
    "\\end{pmatrix}\n",
    "\\cdot\n",
    "\\begin{pmatrix}\n",
    "T(x_j)\\\\T'(x_j)\\\\T''(x_j)\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "To get the answer we **invert the matrix** (exactly):\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "T(x_j)\\\\T'(x_j)\\\\T''(x_j)\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "0 & 1 & 0 \\\\\n",
    "\\frac{-1}{2\\Delta x} & 0 & \\frac{-1}{2\\Delta x} \\\\\n",
    "\\frac{1}{\\Delta x^2} & \\frac{-2}{\\Delta x^2}  & \\frac{1}{\\Delta x^2} \n",
    "\\end{pmatrix}\n",
    "\\cdot\n",
    "\\begin{pmatrix}\n",
    "T(x_{j-1})\\\\\n",
    "T(x_{j})\\\\\n",
    "T(x_{j+1})\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "So we get the following for the **derivatives in terms of grid values**:\n",
    "\n",
    "$$\\displaystyle T'(x_j)=\\frac{\\partial T}{\\partial x}\n",
    "=\\frac{T(x_{j+1})-T(x_{j-1})}{2\\Delta x}$$\n",
    "\n",
    "$$\\displaystyle T''(x_j)=\\frac{\\partial^2T}{\\partial x^2}\n",
    "=\\frac{T(x_{j+1})-2T(x_j)+T(x_{j-1})}{\\Delta x^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discretized heat equation\n",
    "\n",
    "$$\\frac{\\partial T_{j}}{\\partial t} = k \\left[\\frac{T_{j+1}-2T_{j}+T_{j-1}}{\\Delta x^2} \\right]$$\n",
    "\n",
    "Note that by discretizing space, we have turned our PDE into a set of ****coupled ODEs****!\n",
    "\n",
    "If we write $T_{j}$ as a vector of values, we can rewrite our\n",
    "equation as a ****matrix**** operation:\n",
    "\n",
    "$$\\frac{\\partial T}{\\partial t}=F\\cdot T ,$$\n",
    "\n",
    "where $F$ is the matrix\n",
    "\n",
    "$$\n",
    "F= \\begin{pmatrix}\n",
    "&&&\\vdots&&&\\\\\n",
    "\\dots&k/\\Delta x^2&-2k/\\Delta x^2&k/\\Delta x^2&0&0&\\dots\\\\\n",
    "\\dots&0&k/\\Delta x^2&-2k/\\Delta x^2&k/\\Delta x^2&0&\\dots\\\\\n",
    "\\dots&0&0&k/\\Delta x^2&-2k/\\Delta x^2&k/\\Delta x^2&\\dots\\\\\n",
    "&&&\\vdots&&&\n",
    "\\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|   -2k/Δx²      k/Δx²          0          0          0          0 |\n",
      "|     k/Δx²    -2k/Δx²      k/Δx²          0          0          0 |\n",
      "|         0      k/Δx²    -2k/Δx²      k/Δx²          0          0 |\n",
      "|         0          0      k/Δx²    -2k/Δx²      k/Δx²          0 |\n",
      "|         0          0          0      k/Δx²    -2k/Δx²      k/Δx² |\n",
      "|         0          0          0          0      k/Δx²    -2k/Δx² |\n"
     ]
    }
   ],
   "source": [
    "# for example:\n",
    "N=6\n",
    "F=[[\"-2k/Δx²\" if i==j else \"k/Δx²\" if abs(i-j)<=1 else \"0\" for i in range(N)] for j in range(N)]\n",
    "matrixprint(F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What about the boundaries?\n",
    "\n",
    "The boundaries are a problem. Why? Well, consider the first point, $j = 0$:\n",
    "\n",
    "$$\\frac{\\partial T_{0}}{\\partial t}=\n",
    "k\\left[\\frac{T_{1}-2T_{0}+T_{-1}}{\\Delta x^2}\\right]\n",
    "$$\n",
    "\n",
    "There is no $j = -1$ point!\n",
    "\n",
    "The solution to the boundary problem is **not to use the above equation for the edge points**, but to use a different equation instead. These are known as **boundary conditions**.\n",
    "\n",
    "How these conditions are implemented, depends on the approach to solving\n",
    "the equation that is being used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example problem\n",
    "\n",
    "### Suppose:\n",
    "\n",
    "  * We have a rod of length 1.\n",
    "\n",
    "  * At the left boundary ($x = 0$), the temperature varies as $T (t, 0) = \\sin(10t)$.\n",
    "\n",
    "  * The right boundary ($x = 1$) is kept at a constant temperature of $T(t, 1) = 0$.\n",
    "\n",
    "  * The thermal diffusivity is $k = 0.2$.\n",
    "\n",
    "Show how the temperature evolves in time and space.\n",
    "\n",
    "### How should we solve this problem?\n",
    "\n",
    "As mentioned two lectures ago, there are two basic classes of approaches:\n",
    "\n",
    "  - explicit methods\n",
    "  - implicit methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the implicit method\n",
    "\n",
    "As mentioned two lectures ago, explicit methods can be unstable. So we will\n",
    "solve our probltm using an implicit method. We start by returning to our\n",
    "equation:\n",
    "\n",
    "$$\\frac{\\partial T_i}{\\partial t} = F\\cdot T_i$$\n",
    "\n",
    "Using **backward Euler**, we rewrite this approximately as\n",
    "\n",
    "$$\\frac{T_{i+1}-T_i}{\\Delta t} = F\\cdot T_{i+1}$$\n",
    "\n",
    "which, after some rearranging, gives:\n",
    "\n",
    "$$(\\mathbf{1}-\\Delta t F)\\cdot T_{i+1} = T_i$$\n",
    "\n",
    "where $\\bf 1$ is the identity matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementing boundary conditions\n",
    "\n",
    "The boundary conditions are implemented by modifying the operator $F$ in\n",
    "\n",
    "$$(\\mathbf 1 - \\Delta t F)\\cdot T_{i+1} = T_i$$\n",
    "\n",
    "The equation for the boundary condition\n",
    "($T_{i+1,j=0} = \\sin(10t_{i+1})$) could be implemented by:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "\\sin(10t_i)/\\sin(10t_{i+1})  & 0         & 0         & 0 & \\dots\\\\\n",
    "-\\alpha                      & 1+2\\alpha & -\\alpha   & 0 & \\dots\\\\\n",
    "0                            & -\\alpha   & 1+2\\alpha & -\\alpha & \\dots\\\\\n",
    "                             &           & \\vdots    &&\n",
    "\\end{pmatrix}\n",
    "\\cdot\n",
    "\\begin{pmatrix}\n",
    "T_{i+1,0}\\\\\n",
    "T_{i+1,1}\\\\\n",
    "T_{i+1,2}\\\\\n",
    "\\vdots\n",
    "\\end{pmatrix}\n",
    "=\\begin{pmatrix}\n",
    "T_{i,0}\\\\\n",
    "T_{i,1}\\\\\n",
    "T_{i,2}\\\\\n",
    "\\vdots\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Where we have defined a new constant: $\\alpha=\\Delta t k /\\Delta x^2$.\n",
    "\n",
    "This approach may cause singularities; better is to\n",
    "**replace the top-left element of the matrix with a $1$**,\n",
    "and to **modify $T_{i+1,j=0}$ after each step** to the correct value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|         1          0          0          0          0          0 |\n",
      "|  -Δtk/Δx² 1+2Δtk/Δx²   -Δtk/Δx²          0          0          0 |\n",
      "|         0   -Δtk/Δx² 1+2Δtk/Δx²   -Δtk/Δx²          0          0 |\n",
      "|         0          0   -Δtk/Δx² 1+2Δtk/Δx²   -Δtk/Δx²          0 |\n",
      "|         0          0          0   -Δtk/Δx² 1+2Δtk/Δx²   -Δtk/Δx² |\n",
      "|         0          0          0          0   -Δtk/Δx² 1+2Δtk/Δx² |\n"
     ]
    }
   ],
   "source": [
    "A=[[\"1+2Δt\"+F[i][j][2:] if i==j else \"-Δt\"+F[i][j] if F[i][j] != '0' else F[i][j] for i in range(N)] for j in range(N)]\n",
    "A[0]=[\"1\"]+[\"0\"]*(N-1)\n",
    "matrixprint(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementing the other boundary condition\n",
    "\n",
    "Assume that there are 100 points in $x$.\n",
    "\n",
    "The boundary condition\n",
    "equation at $x = 1$ ($T_{i+1,j=99}$ = 0) is implemented similarly, but\n",
    "slightly simpler: \n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "&&\\vdots&&\\\\\n",
    "\\dots & -\\alpha & 1+2\\alpha & -\\alpha    & 0 \\\\\n",
    "\\dots & 0      & -\\alpha    & 1+2\\alpha & -\\alpha\\\\\n",
    "\\dots & 0      & 0         & 0         & 1 \\\\\n",
    "\\end{pmatrix}\n",
    "\\cdot\n",
    "\\begin{pmatrix}\n",
    "\\vdots\\\\\n",
    "T_{i+1,97}\\\\\n",
    "T_{i+1,98}\\\\\n",
    "T_{i+1,99}\n",
    "\\end{pmatrix}\n",
    "=\\begin{pmatrix}\n",
    "\\vdots\\\\\n",
    "T_{i,97}\\\\\n",
    "T_{i,98}\\\\\n",
    "T_{i,99}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "This works automatically if we set $T_{i,99}=0$ initially, but\n",
    "alternatively, we can set $T_{i+1,99}=0$ after each step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|         1          0          0          0          0          0 |\n",
      "|  -Δtk/Δx² 1+2Δtk/Δx²   -Δtk/Δx²          0          0          0 |\n",
      "|         0   -Δtk/Δx² 1+2Δtk/Δx²   -Δtk/Δx²          0          0 |\n",
      "|         0          0   -Δtk/Δx² 1+2Δtk/Δx²   -Δtk/Δx²          0 |\n",
      "|         0          0          0   -Δtk/Δx² 1+2Δtk/Δx²   -Δtk/Δx² |\n",
      "|         0          0          0          0          0          1 |\n"
     ]
    }
   ],
   "source": [
    "A[N-1]=[\"0\"]*(N-1)+[\"1\"]\n",
    "matrixprint(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The time stepping algorithm\n",
    "\n",
    "So what is the actual process?\n",
    "\n",
    "$$(\\mathbf 1 - \\Delta t F)\\cdot T_{i+1} = T_i$$\n",
    "\n",
    "  * This equation needs to be solved for $T_{i+1}$, and is of the form $A\\cdot x = b$.\n",
    "  * Build the **matrix** operator $A=\\mathbf 1 - \\Delta t F$.\n",
    "  * Modify the operator $A$ to deal with the **boundary** conditions.\n",
    "  * **Initialize** $T$\n",
    "  * Then **loop**:\n",
    "      - Copy the current value of $T$ to a temporary variable $b$.\n",
    "      - **Solve** $A T = b$ (linear algebra!)\n",
    "      - **Correct boundary** values.\n",
    "  * And repeat **until time is up**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rzon/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:49: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAEKCAYAAACWrQcQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO29e7xtR1Xn+x1zrrX3Pq8kJzmQhBAIj+ADUNBcaOVeoQExvojaqMG2BYTORaXVpvUjtjbS2NwPtveD93aLYkQQ0asorXBUkFYeaoMgUd7xQUCURELI+5yzH2utOcf9o6rmrFnzsebae62919qnfp/PPHvNZ9WYVafmbzxqlKgqERERERH9kBx0BSIiIiJWCXHQjIiIiJgBcdCMiIiImAFx0IyIiIiYAXHQjIiIiJgBcdCMiIiImAELGzRF5HUicoeIfLzlvIjIfxORW0TkoyLyFYuqS0RERMS8sEim+avAtR3nvx642m43AL+4wLpEREREzAULGzRV9c+AuzsuuQ74NTV4P3CRiFy+qPpEREREzAODAyz7CuCz3v6t9tjnwgtF5AYMG0XW1r5y+MAH7ksFIyLOV4xuvfVOVX3Abu//un95TO+6O+t17V99dOcdqtqllS4VDnLQ7A1VvRG4EWD9yiv1Qf/hhw+4RhERhxuf+fc/8o97uf+uuzP+8h0P6XVtevknT+2lrP3GQQ6atwFXevsPtsfmA5nbk1YLu00lsEzva5HpEA5Kzv1K8bAk7ahATn7Q1VgIDjLk6DTwPdaL/i+A+1S1pppHRESsHhRlrFmvbdWwMKYpIr8JPAU4JSK3Aj8FDAFU9TXA24BvAG4BNoHnLaouERER+4/DyjQXNmiq6rOnnFfgBxZVfkRExMFBUbJDmnZyJRxBERERq4d83wy5+4vDOWguiTE8IuJ8hQJZHDQjIiIi+iMyzYjdoQ/rnVffOpx9NGIFocA42jQjIiIi+kHRqJ5H9EALq5zcey/3/cm7GX32Vka3/TM6HnPFS/8jw0surt03uede7v69t7L1d58EVY580dVc/K3XMTh5cmrxk3vu5e63vJWtv7f3Pqr/vRERc4VCdjjHzJhPcz8w/sKdnPvQR0iOHGH9EQ9rvS4fjbj951/D+PN3cOq7rucB3/1sxl+4k9t//jXkOzudZeSjEbf/wmsY32Hv/dfPZnznndz+6un3RkTMG2ZGUL9t1XD4Bs399pyLt7Vg4xEP5yGveBmXvvAFHHvcl7ded/Z9H2By11088AXP49iXP4ajX/YYHviC5zG55x7OvO/9ndU4+xf23u99Hsce+xiOPvYxPPD59t6/6L43ImL+ELKe26rhvFPPNcu4753v5uwH/pLs/vtZe9AVnPqu7yTf2eH2//4LXPGTL2Fw4YVzLVOSft+mzY9/gvWrHsrwAWX+guEll7D+sKvY/PgnuPBfPrn73od23PuU9nsjIuYN4whavQGxD86rQVOzjM//0msZ/fM/c/KbvoH0xAnu+p3f5Z63/RE6GXPif//q2oCpqpA3KBFhfxDpPTi2YXT75zn62EfXylm77FLOffij0+99zKNrx/vcu6uP/Sz3HFLbVkQ7TJxmHDRXHmf+1/vY/uQtXPaDP8DGw64CYOefPsv9f/rnAJz6rutr92zf8ik+/+rXTH32+iMfzuX/7vv3VL98c5PkyJHa8eToUfKtrYXdGxGxCOSRaa4+zrz3L9j4okcVAyZAcuQIur3NRdc+g/TYsdo961c+mMtf/EPVgw19IVlfbywz7DcF6ZLynExjYsvU9/bCSvswzmVipV2yzlLPZZJpnxCZ5iHA5P77Gd9xByee9FXBiQnJ0aNc8JSvabxP1tdZu+JBwcGmC/feQZIjRxpZYX6umUX2ureFgUZELBKKkB1CPzOcT4PmnXcBMHCxkYDmOWdv+isGpy4h2dhovG8/1fO1yy9l/Lnba8dHn/88w8sunX7v7Q333j793oiIRSCq56sOywTzzc3i0Jn3vo/x7Z9n7cFXtN5WUc87+oBTz3fTT9w9Rx7zaO556x8wvvMuhqcuAWB8193sfPoznPzmb+h8xpFHP5p7Tgf33n03O//wGU5+U/e9vbHX/wPCaqiqfeScxeQwS9HB8xYy7uzDWKYII00XX9AB4LwZNNcedDmyscG97/gTZOMI2X33cs/pP+Tolz+WzY/fzObNf8PG1Y8kGQ4r9yUbG6w/xK7KsYfOdu5DHwFg9NlbAdi6+W9Jjh8jPX6cjasfAcDxr34iZ/78vdzx2tdz8huvBYR73vZHDE5eVDErbN/yKW7/hV/i1PXfwfEnXAPAia96Imf+13u543Wv5+TXXwsi3PP2P2Jw0UWc+OrAJBERsWCY4Paonq80kvV1Hvi87+Hu33srX3jDG0mPH+OS7/x2jnzxFzG557Xc8cuv4yGv/C/tD5gyYE5jBHe+/o2V/bt/53cBo9ZfdvX3F3W89EUv5J7fPc0X3vibAGw86pFc/K3XlY4mLcOg1EuIkKyvc9n3v5C733KaL/yGvffq4N69YF7sZNnZ5qxyzkmeNmegf3wq65yFIe8DDqsjSHTFMpFMXY1yUe20x0FzpqK6mmS35zoLnNM1fbEIGfpiHgNPiGl17iHT1AgK9nfQ/MwP/chfqeo1/a6u41GPPaKvPn1Vr2uf8fC/3VNZ+43DxZ8P6MM2b7vT0tnP512fPT5PtL7NBbut1xzk2Rfsc7/KkV5bH4jItSLydyJyi4i8pOH8z4nIh+329yJyr3cu886d3qtc5416HhERsX8wjqD5DC8ikgKvBr4WuBX4oIicVtWbi/JU/713/b8DHu89YktVHzeXynDYmOai0PIxVJnCCpsoUU+KNDPbXBRbWSLW2/XK9sw45xEZsJvbZqjzgco3I5wjqM/WA08AblHVT6vqCPgt4LqO658N/ObepWhGHDQjIiIWgkyl19YDVwCf9fZvtcdqEJGHAg8D3uUd3hCRm0Tk/SLyLbuVxyGq59PQwTKbr5+RDrjrGx6o0vC4/fQ+79H5NU9b3Vye1dpmc3j2LrAbmUR3oYUcgHwzzgg6JSI3efs3quqNuyz6euDNqpp5xx6qqreJyMOBd4nIx1T1U7t8fhw0IyIiFoNcew+ad07xnt8GXOntP9gea8L1wA/4B1T1Nvv30yLyHoy9c9eDZlTPIyIi5g6TsCPptfXAB4GrReRhIrKGGRhrXnAR+WLgJPAX3rGTIrJuf58CngTcHN47CyLT7MKiVfOme/c73mgP0z77XDfTK+kR4D0NM6uvi37dc5BpTzioMDyE8ZymUarqREReBLwDSIHXqeonROTlwE2q6gbQ64Hf0mrw+ZcAvyQiOYYkvtL3uu8GcdCMiIiYO1Qh66+e93ievg14W3DspcH+yxruex/w2LlVhDhozoxGFtOHNvj3dV0eUKWZ2dq8MCvLbsFe638gsnvw5e2V93ROs39mxhKFhhn0D1xfNcRBMyIiYu5Q5ss0lwlx0GxDw0dyJpbZGfTuP3SGOu0X5sQy/fv2kzHuKiwnQNP9vTPtLxi95esZMrYoeWIS4oiIiIieUCQmIT6vsBeWeUCpxRaNg+j/c2VAje1XP9RHzr0w52XwmtfWrVoA4zRL+B7O4eVwShUREXHAkEObT/PwDJr73T7hZ3kv5Ydsc1aj3IIZ2fQ8jtNjTKeys/1m27tkmf61U5nZXmNUF9Cnu2Scp+1ZmWlG0Erh8AyaERERS4XDyjQX+inokTj0ISLybhH5kIh8VETmtALYHtCHgfRlmdKy9Sy3tfx9ROfspzAf2xwzAh+0h3oRaJVp1uNTC2p41LQ0ht5184CqkGvSa1s1LIxp9kkcCvwk8Nuq+osi8qWYiP+rFlWniIiI/YFxBMXVKGdFkTgUQERc4lB/0FTgAvv7QuCfF1ifiIiIfYPE4PZdoClx6BODa14G/E+bnv4Y8PSmB4nIDcANAOnJk3OvaFlQ/dDMqvksQcczTKdsxTyzebeEolSv6TlPcA96Xi8Hy14db+EjF2wCaZSpT9v1kXWPeU8XAeMIijbNReDZwK+q6oOBbwDeKCK1Oqnqjap6japekx47tu+VjIiImB1zTA23VFgk0+yTOPT5wLUAqvoXIrIBnALuWGC9emMhLDO8Xjv2Gys1YxnzxKJytPV9rDb8nuN0wlmxpxCdJXZ0zYMgHuYZQYsc5vskDv0n4GkAIvIlwAbwhQXWKSIiYp8wx4XVlgoLY5o9E4f+B+CXReTfY769zw0SiC4PuljmPAPblwRTWfZ+YlpQ/B6SDvdOKD2rfXmBr6t7BdQZrl0gVGGcr96A2AcLDW6fljjUhh89aZF1iIiI2H8Y9TwOmocbu2UrPRNBVLDLBStnec7UmY0d9tiDZJm1ovabsXUtrL5L2rZnmbrYdEdfPGiT4mGdERQHzYiIiLkjhhwtOxad2KAyVbCj3K5pkuF1bft7lKVpZuOeH1g71rL1uXeZMa2+85Bnn15JK4tu2haC+U6j7DEl+7ki8gUR+bDdXuCde46IfNJuz9mrZJFpRkRELATzWiOo55RsgDep6ouCey8Gfgq4BvPJ+it77z27rc/hYJp7RZctqA/L7MswO8qcB3ZNGtrsmbMmWd6NTF2seDfxm4QyNP9eqN22TaZW27T0z6jRhi7ZahWZ4dwuYbznaa+tB4op2ao6AtyU7D74OuCPVfVuO1D+MTY2fLeIg2ZERMTc4YLb+2zAKRG5ydtuCB7XNCX7ioZi/5XNlvZmEXETa/re2xtRPY+IiFgIZlDP71TVa/ZY3O8Dv6mqOyLyfwJvAJ66x2c2IjLNALtSzWvP0MZtvhWt7s70+Ja+PDVoel4mhSat/yBDjRbkDJmqmjep5LOo6H0vPQCHnPOe92Sa0zB1Sraq3qWqO3b3tcBX9r13VsRBMyIiYiGYo/d86pRsEbnc230m8Df29zuAZ4jISRE5CTzDHts1onq+FwdOjSR0f9FVFPG/rLNOoVwEYegzjXDWwP9dOnCmHt8L+k6XnPqcliD3mQPWp6Tqb6nXNGK2LNNfVYXJnGYE9ZyS/YMi8kxgAtwNPNfee7eI/DRm4AV4uarevZf6xEEzIiJiIZhncHuPKdk/Dvx4y72vA143r7rEQdPDVHtma/hK/695jW1OQZ/ktXMnE7tlmXMorhW7WZy7od1mCqnazXvtk2x43jNl5sGi95g4OsRhnhEUB82IiIiFIA6ahxGzTHnsyzKbnlkLvu7BNnt2uJlZ5jTP+TSW2UO+Xuhrz2zyLncI3fv/aZ92c8c6cnj0fXyBBQ4kU1n0PuIwJyE+vwfNiIiIhWFe0yiXDas/aM6pXRq/0lOYlva1/U1hLNWKTDk+12Qc9pG7STISnuub62JaQuHKfgctnppco/m2amXar50JfWWahXn5Mnb1w6Y27CvfLhh0X6jCJCYhjoiIiOiPqJ5HRERE9ES0aR5GzBKm0Ri6sr9hOfuCvnlDG++lVPdmCHBvVNvn+Z9tmgliVtm64LTpnjL5r7vXEkS96jpn08MeoHHQjIiIiOiP6Ag6xGjOP9jyF49l9gnJ0eBcyFj6sJg+CS5mRRcDawqzmgfagvJ3FbI0xVHShoPQEIKX3KjMzBJb3uUEarruAKAabZoRERERM0DIovf8EKGvPbMl7Vvt3B7DcRoKaaxOzWbWcG7X7GIaa5slqL1gkN1yNN+7e3ZST1YRHO9imV0aQmeh1b819uxVqs8SRHsiZ7My7930zRkQbZoRERERPRHnnh9idNr0gv2ZWaZ/TZtnua89s41lNu1PY4ttDKxNrq7n7ZKp9PEwh7LObPPrc36Xsu0mo9xcMK0Nu+7Zz5mVauyahxHn/aAZERGxGETv+WFHyLZqX/SW872fz0xe8n7p0masQxN6xqW23z9jPXokspgLK5uVTfdFRwLiLnvmzAjq21jsNC2pzW47RwbdBo2OoIiIiIjZENXzZcRuPuSVeMuWc/NmmXtFkz1zWpKIKXbNvrLvFU22y916mLu8y9pV73nYovtiN5EDu4zPLI/PoX8uwGlzXnrPReTbejxj26aij4iIiAAMy5znoCki1wL/L2aNoNeq6iuD8y8GXoBZI+gLwPeq6j/acxnwMXvpP6nqM/dSl2lM85eBt9L97foagrU7IiIiIuYVciQiKfBq4GuBW4EPishpVb3Zu+xDwDWquiki3wf8V+A77bktVX3cXCrD9EHz7ar6vV0XiMivz6syB4bOTOAzqD7u2r10FpWq+toVyN7vcbZuHRd1qOZt6x/Nss5R+bCG33t1AHVUQ9vU8S41vUcdmtqnbWrkLM/dFXY1lZR9CT+ao03zCcAtqvppABH5LeA6oBg0VfXd3vXvB757bqUH6HRvqerUgvtcExERcX5BEfI86bUBp0TkJm+7IXjcFcBnvf1b7bE2PB94u7e/YZ/7fhH5lr3K1ssRZOnxNwJX+feo6qv2WoF9RZMTqCPUqBfL7KIYfdlYG6NsKmYOX++K7Ltgme5cb7bZFkrV5izx75mVQdm/vbK0t4Xs7PYdt4Ua7bHNmhizCtMdQP7+rAx6DpjhMXeq6jXzKFNEvhu4Bniyd/ihqnqbiDwceJeIfExVP7XbMvoGUv0+ZvH1S4AT3tYJEblWRP5ORG4RkZe0XPMdInKziHxCRP6/nvWJiIhYZlhHUJ+tB24DrvT2H2yPVSAiTwd+Animqu4UVVG9zf79NPAe4PG7F6x/yNGDVfXLZnlwH+OtiFyNWeD9Sap6j4g8cJYy5o6mL/a0MI/O582YgcG/titEp+n6PtPoKnI1JK/oyTJ7wbfNQv+po32nhbbVr0nGWcOp2tjmtPbxrpmJtXWFiPl1aqhGZ58N79/vuMn5lfdB4GoReRhmsLwe+C7/AhF5PPBLwLWqeod3/CSwqao7InIKeBLGSbRr9GWabxeRZ8z47MJ4q6ojwBlvffxb4NWqeg+AL2xERMRqY15MU1UnwIuAdwB/A/y2qn5CRF4uIi586GeB48DviMiHReS0Pf4lwE0i8hHg3cArA6/7zOjLNN8P/J6IJMAY+91S1Qs67mky3j4xuOZRACLyXkz81ctU9Y/CB1nD8A0A6cmTPavcjtbAbu+3Ntn7imsaWFqlgL4VqR8S376pLYylbU5dD3bY6lEuzk+RzatLo12zI5i7KtsMjGwKEzP2vZaTU1imk3e30QC9gvfD32Hd2iAd7dWHaTedm9ZF5mjPzPP5xWnaWPC3Bcde6v1+est97wMeO7eK0H/QfBXwVcDHVOc6OWoAXA08BWOn+DMReayq3utfpKo3AjcCrF955SGdnBURcYigLGSW0TKgr3r+WeDjMw6YfYy3twKnVXWsqv8A/D1mEJ0/2uxiPgPxWWbTfeJRiz7spg/a3mirPbD/w3uxlKZogfD6vaLFFrtnb22LxlDxLPvXBTZbX96+NtwKWy5uLh7SMEW04frwvsZC/H3vljaveVd/nIY5e/uLx2i/bdXQl2l+GniPiLwd8L1SXSFHU423wFuAZwOvt0baR9myIiIiVh0rOCD2Qd9B8x/stma3qVDViYg4420KvM4Zb4GbVPW0PfcMEbkZyIAfVdW7ZhUiIiJi2dA7nGjl0GvQVNX/vJuH9zDeKvBiu+0vWtQxbVO/+2bK8YOylX5hR9b+U9PKQtVuj52wptr1DVcJMc2h4OrZ5gSiRTXfhbOkplmGsgQyLiJgv1bwnIPEgXanUItjq3Z7Rx7Q1v294pAyzU6bpoi8bNoD+lwTERFxnkFBc+m1rRqmMc0XiMj9HecFY6t82dxqtAi0BQb757vY1kxJO9jbF9YysnoexsB5Uj/VXp9ZQ6zCa8L69UBXQH6js2QausKOmrSCNhnD+/qU3+IJLtiz3zZNIUgtjyx+tPXJpps7HFx7mZAw7ymU9qnzfNjSoE9quGnTJX95TnWJiIg4TDik6nnnoLlbW+bSo8Vu2RjQ3joVr4PyFfZMOsJTvA2mhrK0FVkrP2RZQdUKhCFWfVh0Z9ltBdUZWfX6XZTVao+mzsRCFtbG6rS8rs2u2cXGmqaM7pVn1dL6hf027LsE12v1fLNds6W95oHzcdCMiIiI2BVaTBqHAefVoNnahtLCtho9zR2fz54JOupTD6sMrHO6Yc9im2232o9dd6HLZhuy5gbbbJPdb1c2P4LX1uQ178My/Xv62lapsufGpMQN9/gHG5WQkCmH59yNTQyzyy46RS6//nNjmbCSget9cF4NmhEREfuIFfSM90GvaZQi8igReaeIfNzuf5mI/ORiqzYnTPO4dnmWCa5pspeFz9nN8gNQYSnhtLwak2myfzahYFrN5+Ziy+wDn5H1YJm9qtPEvBrk0ZBdt7VbW72bim463sQyw3aaod18Ztm4ZIlMYZmzosEmu1e4vjttWzX0nXv+y5i8l2MAVf0oJtQoIiIioo7wA9+1rRj6qudHVfUvRSqftMkC6rMwNC7xUJyr2/pq9r8Wr3QFfe1i4HWa6kygklU2eJqbnhHWualO7pzHrqYm5/DPz2DQr9kzfYeAbzsL5Wioc/g6u+NR6/L0TkAyzU7rs2TtPl5hmW3osNPW6lX8brCv+9d0Gobb6lG9qTGCY9eQmfrNKqHvoHmniDwC+0pF5FnA5xZWq4iIiNXHCrLIPug7aP4AJp/lF4vIbZjkHXEVyoiIiHbkB12BxaBvwo5PA08XkWNAoqpnFlutHtgr829S4Ypzviob6F1t5fYMVxGV5msb1NnWgPdZZW/NDxrud3g+ulStUA3vqkIfh0OTjG3OntBREpofpjl7doPQ5OB3kS75eqrktVCqwLRijvdw4jmzQ83O4V3SZFuch1o9pT+sMvp6z/8vEblIVc+p6hkROSki/2XRlYuIiFhdzNN7Pm1lWxFZF5E32fMfEJGrvHM/bo//nYh83V7l6us9/3p/CQq7ENo37LXwhSMMuWlL2NDk8AlZZu2aoMX7MDYfgSOhMXVaxenQsA/tbKaNhdHh+CK4dgpRaGTNDU6ginPIPTZkPjPav1qnhTaFFvXFtPCjsMnD5m9qn/AZLeV2TXMN66cNzLMVs2hGDY7JPUF7blPgrWz79cCXAs8WkS8NLns+cI+qPhL4OeBn7L1fion0eTRwLfAL9nm7Rt9BMxWRdU+II8B6x/URERER80KflW2vA95gf78ZeJqYcJ/rgN9S1R27pM4t9nm7Rl9H0G8A7xSR19v953kVXC2Edq8ulhnajNrsn2CogruuyT5UY40NYSwEx8Nn9UQXC+tkmY2hR7YOPaaI1l5PA8s0x+lmYr4trrNAz+ZXKbZFnrCClXXjgzqF7RTcV5MraK/WorrkCqhrzZ5JS/vtEl3TXeeBGZ55SkRu8vZvtIspOvRZ2ba4xq4acR9wiT3+/uDeK3rXrAF9HUE/IyIfBZ5mD/20qr5jLwVHREQcYiizTKO8U1WvWWBt5orec89V9e3A2xdYl/liGstoui60ZVLdl5Y+0DcxQWfKsRb2OYtJqvrQZhbWyDLbvLAh82uDrbv73WanLa+tFl+cmsos1ZMrqH8fJt30vBmWt6jJ1cAyZy4qOFZbw32abO7htYfsAqGWs1fMj732WdnWXXOriAyAC4G7et47E/p6z79NRD4pIveJyP0icmZKRveIiIjzHHP0nhcr24rIGsaxczq45jTwHPv7WcC77Bpkp4HrrXf9YZglwv9yL3L1ZZr/FfhmVf2bvRR2EGhLGVZf3sF9xTW4VqsMs40a9WUtDWzFPauTZfax8bVhmsy1wrxjs7CFit1WKscLRtZi83PHZiFJYdt2Todtq2twLIwIqEcxhNSQukxNz5Yp8gVtEdoz58agPa0m7IuNNva9YE7P6bmy7a8AbxSRW4C7sbkx7HW/DdyMmfr9A6qa7aU+fQfNz6/igBkREXGAmKNzqcfKttvAt7fc+wrgFfOqS99B8yYReRPwFmDHq8zvzqsiERERhwezBK6vGvoOmhcAm8AzvGMKrMagGaho7aFGgbrqq+Yd4UaCdDuDGlTxLtWoEpoTPsfWTbAOqEAGLdQ5X64uc0SHfC6Mqk02T9UrrBr+39BZ4v9HCp8p5X3q7btzhZ/JM7F0qq9NYUe+LLM4t6bJ5R/37w/LbysqrKN7YM3h4/9uMTu4smd0clX6ZL+7puOQJiHuG3L0vEVXJCIi4nDhvGaaIrKBmab0aGDDHVfV711QveaCflPTPIYSMBUJHUO1AqZVwC+m/JL7RYdf+dq9/rHdfLjbWJjP1sLnq3dd4wqG1MOnfDbthx955zrDWXzWHNS/6beGbVZzkjTc7+8X5U1hZFPkqjHnNu2gqS7hpb6W4K4PtYSwP/rP7DtIhdN3w8cumSNo2dB3GuUbgcuArwP+FBPrdPCZjiIiIpYTWto15xBytFToO2g+UlX/E3BOVd8AfCP1aUzLgy5bkNsPvuZVu6bHMu15Ea1sjc/10cQaofKVDwPai2qFNtDwOV0dLZTHP05VphrLDp8R3huiia0E8jSyTG3YpsnVUbcaE6uxsRn/h/pyTJErlGfmAaGxrcp9DduruId6u7QdZwZNp0sbmBVN7dzW9iuEvo6gsf17r4g8BrgdeOBiqhQREXEYIOdzEmLgRhE5CfwkJsL+OPCfFlareSL4OtdZSfWvBOfaWKWgaJGkQ40H3dnHmtb3CRhXK3vBOx9Ancc3NCdKw98WO20jE+2y13bYNYvCAqYVMola0H4gW5tclfOigXza3p5OprDt1JPHldkml/tdYZvtcon/PhxEu19fg42201bbJFvloW0vsLpbZf1SfU0zEPLzFX0HzXfaHJp/BjwcwE5JioiIiGjGIR18+9o0/0fDsTfPsyLzRuPXvZNlUtoxE49leufcZu73WOiUSpT2pOmsrPU4LQygiYV457RBznrEQIMhrsU2VvGch8w4jPVrYGMVRuYv69FcVE2e5vq3ydRwf1uBgVxNcZldctU0DCdjh3zFJSGD7NSKvAc2aRJdcPUsfjfI57PQvSDoUofJEdTJNEXkizFhRheKyLd5py7ACz2KiIiIqGEFB8Q+mKaefxHwTcBFwDd7x88A/3ZRlZqKFtNN7VybdxhoZipVhllL1OGYg1DMABKxtk0omRbUbJnuEVNZWXivNAj42H0AACAASURBVPxuew8BK6ktvtXEOCnlLGY1tRnifITs12MpApC3sLHac2Q63ah4k0tZWuXz9zvq3nmuwsCknWXm7WzS2GqtfB3t12TPbNWKvPtqfUNpb7tpfVAxjpt5MU2vzMOGzkFTVd8KvFVEvkpV/2Kf6hQREbHiEDi03vO+Ns1vFZELRGQoIu8UkS+IyNR1z6etIOdd969EREVkZbI3R0REdOAQ2zT7DprPUNX7Mar6Z4BHAj/adUPPFeQQkRPADwEf6F/tGVFxiFBVd+xfKRw+nmpec46U+9JT/ZMW549ThZpU2YqKFP5ulE+nhKyE8pYy+yaIcL8VRX2kpup1OUq6/sNIKLevuuLL1WByCeVL/GubCqO63yQXQb2dXHn5t/Lohq1zQAidN179a301oXBOlrIF9/p/w6Ja1qmqmFP61nsWtLyXxn6+Qug7aA7t328EfkdV7+txT58V5AB+GrPc5nbPukRERKwCzvNB8/dF5G+Br8SsSvkApg9yTSvIVVaBE5GvAK5U1T/sepCI3CAiN4nITdm5c1Mr27SKX/13wDITrbHMRtI2tXSvEpXO4Rndc4+VhQ6TgMGEzqGQBWjIMpqcQIknb6JIonUm7RVWSVQSFuTVp84mSxkKWXIp2FrTf5ZGB1FNJmosuh5KFbAzj0mHExYq5cwgV1U+Kdpx5oGgxp47NAS/7Tq1n4b3F8hVla9ZM2h12O0C57V6rqovAb4auEZVx8A5mlljb4hIArwK+A89yr9RVa9R1WvSY8f2UmxERMR+4ZAyzWlxmk9V1Xf5MZpSNXp1JSGetgrcCeAxwHvsMy8DTovIM1XVXwO5P1o+kOp/zgKbmM8yxX3V3em2z6CKeWZbCI3/dbdfbtGGL3zPsBxBuz/+DcylZJlakbkz3Z3ac9Om4/lyhMw5YNeiwX0tdW/f92Sv2fxK9hzaaisURqUSJlaTx98P7JlNdugac64XVx4HFKlTqpAZtsnWYIcOUYaKNcjky9mQ/LqiGfiy7RXuuQuGiFwMvAm4CuNv+Q47e9G/5nHAL2LiyzPgFar6JnvuV4EnA87k+FxV/XBXmdPiNJ8MvItqjKaD0j1oFivIYQbL64HvKm42dtFTbl9E3gP8yK4HzIiIiOXC/rDIl2Cmeb/SRui8BPix4JpN4HtU9ZMi8iDgr0TkHap6rz3/o6rae4bjtDjNn7J/n9dbhPLePivILRZNNiLvryRasWUWdkxpoA+FYOWn3mcvrctduC94DuQeQ6nZMj2G01CcqO2DTQzF/93IVrByerI3yKZIN1OBIhogtPWFbEzy/mysSUMoZA3smrWt8CrXk62EArQuS6JQSarcxMa8v402z+IFdZP1yiu2LLqiISTldWHfbbQ1O7mwhXawzZq2E7DNio12Dtgne+V1wFPs7zcA7yEYNFX1773f/ywidwAPAO5lF5imnr+467yqvmrK+c4V5ILjT+l6VkRExIqh/6B5SkR8DfNGVb2x572Xqurn7O/bgUu7LhaRJwBrwKe8w68QkZcC7wReoqo7jTdbTFPPT9i/XwT8b5QLtH8ze1xwfVGoeZP93/7XvPCwlqyzYJme7cgxMjdVUg0fK6dOVgp3m/uSS9X25zFM/3fN9hdUvckm5u/V2IpjJE5Oj2VK4jHTWlmFdJa9TJOv/N3IMnOv2g0mvQqLbmovb19tu1VYtGfLrCSNpsqkKwx6is1WQhld3V275VLVEhpkC+VrkknDtnIaAlTYc6ghNNW5JpMG58N2Chhm4/m9YjbGeqeqtk5sEZE/wfg8QvxEpUhVlY4sOiJyOWYViueoqrO4/jhmsF0DbsSw1Jd3VXaaev6fbWF/BnyFqp6x+y8DOsOEIiIizl+0WBJ2BVV9ems5Ip8XkctV9XN2ULyj5boLMGPWT6jq+71nO5a6IyKvB35kWn36xmleCoy8/RFTaHBERMT5jX2K0zwNPMf+fg7w1lo9RNaA3wN+LXT42IEWMSE83wJ8fFqBfZMQ/xrwlyLye3b/W4Bf7Xnv/sNzgBQqT8U5YtVUoQgzEkCS3AZDW8dB4EzwVfI2e3tFZQunR+Z1da8pfKVSSAdaVbxARa+o5hLKZZ/V1nl9mZscCDV1L1DNO9S04h02yeksGwJqpxFqxeRAocK6sDGgdOoV9TfHnNmhNQOQr5LnFJMM8H4XoVWB/OGjivK1KlvYXsV+AppobRJCaHboct7VQqoqJqJms0PpvPNWEHDq+jwwr+d045XAb4vI84F/BL4DwOaxeKGqvsAe+xrgEhF5rr3PhRb9hp2sI8CHgRdOK7DvuuevEJG3A/+HPfQ8Vf1Qb7EiIiLOP+zDoKmqdwFPazh+E/AC+/vXgV9vuf+ps5bZl2miqn8N/PWsBRw4Kk4gn3lhmZeSeM4f8f7WH2Sog4aMxXP4hF94P5xDoJG5tBnf3SHfqVB1LmiFvVTYSvC7loikJp2U0tUykHuV8RmmdYo0sjP7l4biOsNxnM/Gd3747Nk/5zu4Kk6g4B1O+c/rszE3RbIWauQ7gKY5uPydFjkrbebL5iXnqIXEhfBZdJeMgXxF27h+6THqeTqC9inkaN/Re9CMiIiImAlx0FxyiLMZlV/uwr4XTLPzk1b4LDNx0ylbP5GGRzSGG4FHB7EMTErmUnzZqQUS12x/Uv4pivIznAesy4Xj1EKrxMlZMpfikjCUqisUx8nkMWm8qZMVmUI21vAq+7CxItVdJRUaRRhOoTH4U2ChRm9aGbSfUKUiH6WTwmNjfeRqhd8nA/lq/bSFZTZNDZ3KosP+5fe5BpY5V3smHNokxIdn0IyIiFgqRPV8xVBbP8Z9uQs2SRHInjjWmZhP4zTvsgQmJPWZomfPFDd1Mq8yMsc827ywlb6mNLMxn32J1tmKkzXRaoRAEOchlGyz7oF1m9TsfsVjguDoavIHW4a6NqnLJBrIW5GJgqWpSzziWKbHxHz5QlQYdItdsJDR1woKBiZ1W22Lvc5NDa1oCA2ylSzatVkpZ8gywyQyvlwVFu3OVt57VZupMGiv/cK1geZl04zqeURERMQsiIPmiqDJy1qxC4Ek2K94TpKoZ8ukZtP0WZi/H573mSaezVKCL3pxvIGR+SJUivHZZmAXq9gzC7ZJwTLxZZJ6arGCo2h12mRjGrXCJiaB11w8uyZVxuIzdK/eoQlVa+1WblVGRvFS/WTKifiFei/QJrFolCdg0IU8eGzayeT+4v2tmhgDG3RQXKD1OFttyaAptAPf5l7GDNcfODU9YaVPesd8O23oNW9h0bMiNMEeJhy+QTMiImIpIPnhHDUPx6BZfMW13K+xlsCmGbDMgm3Wnu2+9CatWGM6uCJOUwjjM6te89KWGXoqKwTJ2cUIGJnPNn253O8ixi+wZUqZqKMp/rQ1CUnAWHzmVfntsbEuxlLsBrbNSgyqeK8zcTKV7Ve11TYnWTGPDliYe6GVilBhYUVEQ1OEQ2jPbNIO/B8NLNqR9MJWW8hGsy26QTbX51x7NSZXdnI5Nu0vodKiGTT1yT0h2jQjIiIiZkNUzyMiIiJmQRw0lxtNKqyf/KASauSp5oV67jkT/MDv3DpI3L1lXk1XMFVVpKL+SLMzIVRfA2eJNqh4NdODaKG+ahEU7at3FA6gWpiVK9YrqDNhR+BUCAO/pWKC8NTXUIVNPCeQ31beb3WyFjLah3iOEidfOP01KA0vTUdgA2iRJ1BhC5myQH11TdFgAQitApV+GTjrStnqZoc22VwfVBTRppynYuvmOxylaiby+6RvemiZ9rpbRKYZERERMQvioLkCKMiEx8pc6IbHvnyWmSS5uUy0CG73kUDBNmsoLPsUDgRcEosgpEMyaf66u+e4+io2Q7nHyDz5iiLD9bC9oG/f2TV9eqhlmUVslMfICkYlFfbsJ3toYjCtwd8aOEzwnCLumJVNG5wkErLoFgYNkOdWtpD6FbIFDLpFK6iwzIYplO6VNfYP33FHKVdlamjYbpZl+rKBka901FntoPaCw61kmxUnUEYt0Uqbg2vXcO/sEOJwDZoRERFLgRinuSoQf3N2PryErhrYM3NSKVlYUrP5QZYbwqNFAe6kGNuSZSku3KgWlpOBH+IRft0rcESvZs8MZQzsYhKwlbSU000PbbNnmrAVy2J8u1+FhdHAnsXKVt8aw438kJ+A+NXbDhsAXtpp8WRzTMxn0LWQHCnlq7zfik3SY9CenMWxkGU2ME3HIgtbdE0uTxvwbOy+bJLmpZ22Qzso2aYLf9OiX6pfL2erbbJhhvbnRTFNoN1Qvto4XINmRETE0iAyzWWFHywMXuIKc64SCN3AMpMkLxhmErRyDmXAtA0khuADWqRKo8pUfJaZUfm6l17NBlnczwYvbDXhMJXphVIwFmef1WCKaBm473GXwFYWwDKxildZA5bSwDZrnnOh0bNcMZ8Wdr9yiiEJkCqklj0HLLOJQTvvslhjY1PylZrtL2DQNLDMJKvK5bdP6aGnqiVItd1q9kwnm8cym6b1FtX26LpJGtMgnBoGXQnY9/tl1tB23ruY20A3b9a6ROi7sFpERETETGgy3TSac/ZShsjFIvLHIvJJ+/dky3WZiHzYbqe94w8TkQ+IyC0i8ia7CFsnDteg6b7gULOFOfaVpnmFZaaJMkgz0iSvbGXsZnNR1XRwUnzdJfOYiv1bYZkt7Ey8L73vVQ/j/CqeczdtspDT2sQKOR2TViOv3dzv0B5o5BJqqe5yCbzLUqm3ZKVMNXtZzaMbsJmKHdP91UZbbcgyffmSmkwtNMe2l9MOanY/J4snVxK0YU1baCqqphkE9szUl61stzCyoy5Xs1j+VN5af5tmg27SFOaA/Rg0gZcA71TVq4F32v0mbKnq4+z2TO/4zwA/p6qPBO4Bnj+twMM1aEZERCwHFMMs+mx7w3XAG+zvN2BWyu0Fu2zvUwG3rG+v+1ffpukgvj3T+4o7FpKW9qKBZZgiyjDNijhNh1yFVJQMmr/ulqloLd1WyCjLr37is7EGxqVJYBsLiqvE9RWszHlflSTVOst0sgfxp6rSHBHgyvZsYr63tWDJniwhe6l4Yb1IAE2ox5264j2SVEnSUdhqLfNK6zZpbBs5uXIb4CpS2jUrcrkZM7lUGXRGJZbWl6vLc46VyZ1rks+3aZayedEASV7IlqZ5YX92DNPZnDP3SMFGBxiZtCHm1GkIPoOuxg1DbaaTe/68mGb/55wSkZu8/RtV9cae916qqp+zv28HLm25bsOWMQFeqapvAS4B7lXVib3mVuCKaQUenkEzIiJiudB/0LxTVa9pOykifwJc1nDqJyrFqaq0L/D1UFW9TUQeDrxLRD4G3Ne7hh7ioBkRETF3eAEte4aqPr21HJHPi8jlqvo5EbkcuKPlGbfZv58WkfcAjwf+B3CRiAws23wwcNu0+hwOm6YfQAxliwk20NtzAomSJkqa5AzTjFScul5ug8C54FCaYZxKZNU7p+Z5Kl2SlfuJp9I6FTfx1b8GFamQw8pXaGDiwnK0DMkppk0aGdM0L2QapLmV12wuvMqpgFUnkPtBqd7locoqddU8UNeNGqiV/JruNRavM1DLSycJZQKS1A/8LtXX1HPk+fIlzonX0k3UmRsqTi5KtdWXw9sv20+R3Gx4z6oNDp7zp6qaayGbeLL5qnli++MgzQpHZZLkXnu1OIN8tbwmVzV/Ziijay9/2zNUa89s2/aI08Bz7O/nAG8NLxCRkyKybn+fAp4E3KyqCrwbeFbX/SEOx6AZERGxfGiInGjc9oZXAl8rIp8Enm73EZFrROS19povAW4SkY9gBslXqurN9tyPAS8WkVswNs5fmVbgaqvnpV2/2A+TO1QCvcWEFw0ChgmlIyhXIadM1GHOlUUUgeB+Yo6CoUj5Za+xlNLwHnYWcUxEG5wlFSYWblbG1LAVxzKd8ydNynCVSv3zhBxzLvPL8sNxKunEvFR3DeyyCMlRrTIvy/Y1kcb/IH4G8zJJR+nM8x0laeprCqUjz2+3BCnkyitOIE8z0AbNoMkJFDjuJNcgVMrWU6VxiqjfZkWijrRk0FRkyyv9059sYZySTgvIyfLE9hnfiwi+gwufXXqyJIFc5q/WmHLbXIdZMbdA+Q6o6l3A0xqO3wS8wP5+H/DYlvs/DTxhljJXe9CMiIhYTigQ1whaYtTsmZQJHTw72CDNGKZ5hWUOJK98tROEXIRJwxfdBX5XE3X4TKzOUqTGNBviUpxx0dpMa2Etvo2sZhcr2YqzZSZJbuRssO9lYNh1nrics1Y4m4Akl+ly1WTUwj5YYRdJGUIVMujKekBQ2jPt1EJnz/RZppPN2aR9e7M4NpknZWrAinze5gd8B/bnOiPT0kYYtJnTDrpCxNRpPzaY3WkIScCgnY12mOYVBi0u3ChPKmniXIdXXzvwQt3KUKN6GzrNIHHt5mTw/w/NA4dzzFysTVNErhWRv7NTlGqR+iLyYhG5WUQ+KiLvFJGHLrI+ERER+wfnAJy2rRoWxjRFJAVeDXwtJmj0gyJy2jPAAnwIuEZVN0Xk+4D/CnznbAVZu5kN9K4s9+AYmMcyB9ZrPkhyUjH7PnLLVpznPGtIYFuwTOdddjYx+zeZeF/2ScuXndJ21MVYfEZWtWs6JlbaMweFN7kaAVBJypEnXpC0KTD0mpcRAXjTQuv2zMR5lNtstdi6ajUpSCVY37P7FWzMS9BRMDFPNsegK0wTKpMRwhR/1YB9L4mKzyonVZaZZNo4GcFNRCg0g6DNqklVaIgI0BqDdux5EDJo90qTnCxPCzunkcvTTnKpprrz7bQ1u6aWnvOwP8r8iOZhXcJ3kUzzCcAtqvppVR0Bv4WZ8lRAVd+tqpt29/2YOKmIiIhVR+gh79pWDIu0aV4BfNbbvxV4Ysf1zwfe3nRCRG4AbgBIT5ZJTMJkFr733H3JE5uQw2eZa0lWftVREsnJNaGwACY5eV793rrYzNDuR4ZlZiVbkdyylknzlx2v2sZNj2XMbukJn7Go530NbX51lukiA4ZpVokzze0Uw0RM4trKFEPAjzs1N9SZSuLky7BsWmtRAaIeO+76j+HsmYmz1eLZanPP5mc1A0+2QZLXPOfkiXUiKyq+XKZ8N+W1tPEFkQ5+PO0EkonWp4b6/cEyMkmMwI7MFrK5GM1Erdec0k47yEkHVdncdF5ftlyFxE0NxURx5IVctkwN+2JVG6i0VwbiyVW0m98hna12jzBddwVHxB5YCkeQiHw3cA3w5Kbzdh7qjQDrV155OFsiIuKwYe8ZjJYSixw0bwOu9PYbpyiJyNMxc0ifrKo7C6xPRETEPiIyzdnxQeBqEXkYZrC8Hvgu/wIReTzwS8C1qto4Z3QqnBPIBUTbcCPxptWFqvkwzRhIziDJPKN7ziRPyEmYkFi13XcmCHkuZRaiIthbSrV1IhXV3KlDlXyMXj9ykSKSNGuwhTPBqXkNzoTBwKp3g6yYFurkawqlUhUyq6YbuTxVzw+M9lTzZFIN1jeyqefgoqrC+iFgTo7A2eCmg2LNDqETKEkzBoOqar5u854moiTeG8udOp4nNqi9iiILUGNQe12dTZxzqym7kZMrqTuBKm1WrElPVa5B3eyw5hyToWwCE02MXDYbvTR1IHXmBqk4JEOzQ2L7YiWDky+Dy2c6D6yovbIPFjZoqupERF4EvANIgdep6idE5OXATap6GvhZ4DjwOzbV1T8FCUIjIiJWEnOaw76EWKhNU1XfBrwtOPZS73dr9pJe8BN1iOdE8BNX2K+4zzI30nGZkMP/HCYwyQ0r87/ouTpGZthY4QiqhHJIYXBPJh4bm3iMxf+ye2ysYK9NIUc+Y7GMzDkTBoNSxmGaMSySkFRDV5wTSFSYSFIjE06u2vRJj7EUDNM5ScLkFl7QvpkWGkwvLNrMMWcrl2NmnnPLtJ0WLHptMKnIFjq4RC2LriVY8ZxbRcC+L5N4jp+qbEW7hZMRxPxTyBac8xN0aKrF5gfrOwbttAM/BC6Uzc3pdc6gSpu5dvPXNZrWZpX+aMopnHdIGf42D0T1PCIiIqIn7AfqMGL1B00vQDoMah9Ye+YwzVhPJ6ylGWvJhLXEheOU2b8z/6ueJTXGYkKNEvNlz3z7kRhb5qT8spuvu5ZfevdVD0JyjAGW8ngoUxH4HTKW0p45TDPWBpmVLStCjQaSVRjLRFMmXtC+LxsVm59n0wwCvwvG4kKp7BRK31ZbyIYiiWcDdnLZvxpOm0wUBjkyMFnMS9kmNdkqYWJiWLRatgkUwfyqRvaqrdazYU6qmkEy0cIW7Qe2h7LZ5PC1xCp+e9W0g4EJo0oHVjbbZq5vDjx7Ziibs6+XdmhpsNNKNfStEm7kyeVCxcIpry5af57R7ZFpRkRERMyAwzlmrvCgWXz1nSHN85inJgDaMbCQZa6nEwbiUsKZr7rxLifFl95QhPKrXnjOnQ1pYhlmOA1vbFmm95WvJBrG2fww7CqXWnq7Mujb/k5tdZLS5ufsYmsD41V2Mg4kq7CWHClsY7lIhWX68hUM2reH2SmhIRNLJiUTK2219rlJyVQkpxYZoEkgXwo6qNozfdnWPCa2lkwCW21uWHSeMpGk6Ba+bIWNVgO2OanaoMWTzUUGVOQCRIzNTxLDZH0Gqo3tVmoHJqA9Kxl0mrE+mDBMMjbSCUkwpTdXE9yeJ4JmglRYZrM9s/SSe1rP2LNrev3RRARUZcvt/yPV+VBNyQ+nfr66g2ZERMTywoblHUas9qBZ2JHU2jQpUsENB4ZxrQ8mbKQT1lLLMpMJ6+kksB2ZxBw5Qp7Vv7K5Y5qZ2YpYv+BrXrDMcflVdzZNVKu2MZfaS6hN0aukggs8sG4K3mBQZZlORhd/6ttqJ3lqPbA5Y4/3VVLd+dMMPa9yKaNWPLHlch1aMGgji5aec6jYM7sS88qglG1tMGF9OGE9zTgyGLOWZAySrGKLBpjkaWn3a7D5GTZGYfMrWPREPIbptdm4jNFsTK5ip55WGLRnfw6nuxYM2kU5+NrBwPTDtSRjYzD27Ja2T6ow0YSJtqzh7ttq86psknl90spW2jbrqe6cbJKYdzaPzEOCxuD2iIiIiJkQB80lhMc0y4S1eRH7VtiMBmPWkgnrlmkOkoyUvEiKkJEYu1jgpcxVyFwiCN9zPhHDVqxtLHFf9nHwhQ9smkW1BXLrqZREUaRwXFZsmmndA5taz/K6ZZkbgzEb6YSjgxGDJGNYxPrlVgZj68tzKY452QoGXWFjFDa/kj37jFoDm2Y9PhPMuuO1uNOaPVONbIOqbFUGPS7s0U42h1SUcZ6SJ0KS1xfBK9mYjamdNDBoxzJ9Bp1pyTRzU1dftjJzsi3LsWj8qABMLOogJxmolcsw6I3BxDJoY4c+kpZM07RNzlgTM9vK65NlmyXWTu1kczPRpBJzKi7SobCx27/aLFuTDXpP2IdBU0QuBt4EXAV8BvgOVb0nuOZfAj/nHfpi4HpVfYuI/Com54Vbzve5qvrhrjLjwmoRERHzh7Np9tn2hpcA71TVq4F32v1qVUwKysep6uOApwKbwP/0LvlRd37agAlx0IyIiFgQJM97bXvEdcAb7O83AN8y5fpnAW/38vjOjJVVz9WFG4lzAmkRFO2SIKxbtfVI6tTzCUMxKrqvgo/VOEqSfNCgwkKeJ8axkCWlOpRBMi5VWKe6puNQHQrUcxtC5IzuxVQ88eQqHCQmHKd0AinDoQ34HkwK1XxjMOZIOvbU81I1z51aJ3kRllO8Q+pmBz8Up1TLnaPEyJWO1cvYrkX9VQQGJp+laRpPsCZniZXNqeZOtg1PtqODcWFSGYZTTPIBmQhJIJdrs6bA9mTsOUo8c0o6ViuXMTkkE8/skJequZOt7IdSOrfSwAk0MNNdExc+NcjYsA6go4NRzTHp4FRzEhgkGaM8rbaZCzeyzq1KWxWymTZLvb6ZZEZdr0x59WQrptHORavW/bJpXqqqn7O/bwcunXL99cCrgmOvEJGXYpnqtGxrKztoRkRELDGCONApOCUiN3n7N9ocugCIyJ8AlzXc9xOVIlVVGkMNiudcjlnK9x3e4R/HDLZrmJy9Pwa8vKuyqz1o+tMnvfVW1uzXfCM1bCVkmhvJuHjEWNPCruIHtYNrd+MsyTPfCeQM7Z6DZKykI8fGtHAsiBe6UoTcqNjwKALWQiVBR+EEGuQkw4zUBkZvWGfC0cGokO/YYIehZDWnwVhTxpKSBE4UtQ4Fn0EXqcTGUgTq+yyzYNEjLUON3H8MsbFTItaxUGWZFdlcGJXHxgaOZQ4Ny/RlK5lmVmk3F7ife7L5baZZlUGHAd8lG1P72zJMF0rl2gzXLaRcEdKx6JBBeyFHMjCrTg6HZtqkcwCZfmm0g6HkrCfj6vrtmpKKkmfDxlCqvJI0xrHNepulRbtpqfVkFOFvpj9alokgg/mEGxXor3nfqarXtJ3sSuwjIp8XkctV9XN2UOxKMfkdwO+pajEAeCx1R0ReD/zItMpGm2ZERMRCIKq9tj3iNPAc+/s5wFs7rn028JuVOpqBFjG5Kb8F+Pi0AleXafqJh13KrbScOnl0MOKotfUdS3dYTyZsJGPLxspPYJIb22GWJ0VAOJThRlmekGVJnbX49sxRyTLTkWOaLhVXXgZIp2XQdzn9sAxyL5mY2fLC5mcYy2CQsbE2LhiLk+9IOrKsxaZOs9Mnx3lKomqSdXgLvxjZvODvilx4jIWCsaQFg1GScW5sms4ei2WWqbUfKsV63bVpoc7uVwR+a5HAYt1j0EXbDXYKDSG17ZZpQmKnGWaJGBugJ5sL1q/Y/UIb9MiysZFhYo5Bm+B2z+Zn/8ktkxZn9yv6oRdK5cs2UJKhY9AZR4djjg5HVrYRFwx2WE/GJqmMlS1TF0pk7NGDRvu6TRzj2zMzEG+CRdFmI9teIy3CqGSi9r+PnYRgGXRjGr+9Yn9smq8EfltEng/8I4ZNIiLXAC9UTIlkKwAAFO9JREFU1RfY/aswK0n8aXD/b4jIAzBN/WHghdMKXN1BMyIiYnmhCtni51Gq6l3A0xqO3wS8wNv/DGaxx/C6p85a5moOmsXUNTsFL4E01UqCB9/WdzQdcTQZMZSswjRzTcz0wlzIbLKO1LMfObaZZwlMkkrwsGErdhuXLDMdqZ2Wl5fTDO0XV3MhT8VGEzfLVFmdcaAwNOnSBsOsEhh9fLjDBcNtjiQj1pMJR9NRxeaXaUKSaBkg7QVPO7ufYdAJWBZt5CrZWGoZtPmbF4wlHeUmqN39n0iMbKqQJXZN8bRZtjzVGhsbNsh2fDDiSDLiuLPVoiUbEyFVszzJRMqCHMvMbNsVtlrbbjUG7VimZWPpyDFohczYZ419WZDcrnJZTKuUBgZt7LROtnSQMRy6gPZxwTKPD0YcG+wUmo9rt0wTdtT8l8xIGGs1RWGeJ3aTStKYpMIyyzYrNJ9xTjK2dujMerVFTBspMLCaj6cdzAVxRlBERETEDIiD5pLBec6LRbiM9/WI/aIfG4y4YLDFkXTM8XSbDZmwnoxJyQvGkklClpuUcGPLWBxLyTF2v9xulemTY7FMhYJdDnYcYzFf9STLIfPWSbGsJQHyRJHcY5s2xtHZxvIBBWNxLNMxlmPDEceHO5wY7HAsLVm0ka1M+OCiAvJEGCQZZEOAwuOcOS/sxNrGxuU0vHQUsLGRkuxYxmJZNP76LwKkiXEmB+F5aqfolTGMxlbLQJGhSZe2MZxwdDguZLtgsMORdMTRdMTxdNvY/OzLypBCtkGSk+SlduBPM9Qg2qGJQSdjSHdKO20yypEst9EOWsZgpolNcC3VbOT+dNcUdOAz6Izh0Mh2ZGjY84nhDsdSwzJPpNusW6bpZCsiORIYasZIyv+eqv60V8Myi/SDtj+Wmo9tM8sy053c2Nbzsj8a1iwwSMx66em86KWrMNU+coiwuoNmRETEEsN5GQ8f4qAZERExfyj74gg6CKzkoFkEiacmFCdNzVrffijOscEOxwc7HE+3OZFsF0Z3k93IqLAjG/g9lrS6jrYmxsmQpUwmNrxjIshYTGjHyKqwIyUZwWBHjZo3MqqQCcnJS6M7dqraIIEkqYS0FPDVvMSo6KSlo8SpsE41v2CwxYWDrUI135BxJSRnW606nkslo1Ph3LIbWeBQ8NTywY51bu14ck1ykkkObs6wiJ2GZ4SQYRnUroUjpSHcaGhzZ64Z9fWoZ3a4YLBVmB2OJju1oPYUEya2ma15Tj3rBMpdYHvihVL5jrvSUTLYMSp5OsqRcU4yymzWJhuW48lmzCletHsQ0J4PlDwFHVqzwzBjfTguzA4nhjucGGxzPN3heGrU86FMarK5NvNXSjXZ6ZPCKamTUjY3EaE0Fbn+mJfmop3M9MdAtpwUyNEkrYaPRUdQJ1Zy0IyIiFgBxEFzuaBJuSZQmYPRBEZfMNzmwsFWwTIvSjdtaMfEC1tJSDRnW9cqgdO5Jkw0KQLb89x82Z0jIR05o7v3Vd92ITk56fbEGN0nOcViLiKQ2ABsMTk0/Wl4lSl4A8NWdJjDUG04jsdYLBO7YLDNhekWRxMXujIhRckQE0qVQy5JhckAVi6pBOz7ISvpyDhH0h3LXLZzyzQzZJwh4xzJsjLcKLUM2spG3iKbZZj5UNGhZdBrE8ugR1ywts2FQyuXZdBHkx2OJTuVELGRpmyrkuXCMMnM5ARKNpa7yQjeJATXZumIKsvcUQbbRqZknCOjiWkv64xBnHYgxbRKI2cTe6acFmq1gyNOtqGRzfXJo8mIE8kWa17421gHpJoz1kElNM5pBwqlU3LiO7hKrafCMndKB1AympT90cmWJCRinKAFA/Wnu+79f2gcNCMiIiJ6w4zyB12LhWA1B007hVJSYxdz9szjwx2ODUYcS3e4MDVf9QuSLU6kWwzJWHNBxCKMbRDxmkzYYVg8urD72eDvfJLA2DLNURnake44Rqak25m1IWUk2xOYODammNUZBew0Q3WZzj37kT8FL0/V2DMHSrKWsTaccGQ44fiaCWa/aLjJhYMtLky3uHhwlg0xLLOQDSObCY5OSWRQXf/c2mqzTMgnCTJOTBjVSOp22u3cMs2MZHtsmOYkNwZ+x6DzkmVqBxvL7eZCclyw/pHhmAvWti2D3ubkYJOTg3McS3bYkDEbyagIyRlpSqJDyGmwQ4tdD8mbjDBusGc6+7NlmclOhjgWnWWl8yJNKtoBw7QIoSlXnhQbQmXbb6jIWs5gOGFjbWzabG2bi9a2ODnY5ES6zfF0m2PJDieSLVLJC+1grBnb+ZBRMjB2Wy+2KVNvOq+1Z9b6o9MQRqbNkp2MZGQ2GU2qbZaY/qipmPWFtN5uc0FkmhERERF9sT/TKA8Cqzlo2kQdLomFm6bmAtoNE9vkRLrFBck2J5KtWhDxtv1dmVaJWQFwnKeMs5RskqDjpPQsjxzLNJ7lwbYy2Mot0zRfdNkeQ2a8laaugqYJMkgN+0kEcZ9yFwWQWBaWgg5Bh7llLBlH10ceyzSM5cLBJhelm1yUnrNM00QFGNkGjMgYSVqR2a0PnhVTQ019CrnGkOwYBu3sYoZpTkh2JiRbY5hkyCSres4HqVkPSAQZVCmKJl5gu7PVrimylrG2Nilku3C4zcVr57h4cI4LrVwF0xSTxSsjIcltRIC11ZbRAlJhY7nnWTaMrLTTOgY92DIsM9kaGybmZHPRDoPURDowQJOkypqkyqALW61NEr2+ZrQDZ6e9aLDJxYNznEi3OJEYpnk02am1WZLkbOuwcS2ncVa11YqbYOFszzvGTpva/pjsTEx/HGfIzrhiX9c0gTRFJon57a0m6pLG7BkKGuM0IyIiImZAnBG0REiARMukvDad2InBtrWLnSuY2DEZcTQZMyQnteub73hT8lIvU+o4T43NL08YTwzTZJwgo9Lmlzo2tu1sfhnp1tjYMkdjZDQ2aklmvdZpigxSVBWGKeKmynl2sdwysXxoYv0YWnvm2oTja6PCLnbRcJOLB2etbJtclGyyIWbpDmcb29aMRIesqYkWKDywmM3FnuYTOxVv7Mtm2bOzZW5NSLYnJDtjZGsHHMt0TDNNbYqxxKaFK6foaSJFvKnZlHzNeJfTYc6RtTHHrdfcyHXOtts5LknPckxGhWxgGFea5NaeOai0WxlXm1jZknKq647XZjtKumVYZrpt2LNs7iBjwzSLNgNEh5CmRkvI87qt1vOa55ZFs2YiAnwGfcnwHCeH52y7uf64w0YlPnPCOV0jzxMSjJ0zQa1cRuuZZKmxr0+Sii3T2ddNf7R22s0xycj2x/EExpNqm9kVVskGldhNvP44F0SbZkRERERPqEbv+TJBxSSwHQzyShKLk5axXJKe5ZL0LCeSbY7KhGNJjls5IhMltfFqI8wyECllAtiJpowy83XPx6mxjY2k9JhvGzaWbuUMNyekmxOSzRGyPYLxBN0ZQe55zgcD450EyNQwziIZhE3SkZqvez6A3Nr8hsOMI2tjTgy3Obm2ycXDc5wanuEBgzNclJ4z8afJiCFa5FoYK4UN08SkamnHzY1co0lKNkmNrXZUshbHxAo2tpmRbo6Q7TGyM4LtHXQyKY37nmwMUiRLO9mYsdUqyXrG2vqYY2sjLlzf4uK1TU4Nz3JqcIZLBme5KNnk4nSTDckK2TI1mdDIIUvMbKehTADf5mfbLEtr0Q6lDdqyzM2xiQbYNgxax8YO7XvOBcyqMZPEeM5d3/PstLnTDtZMXG26nlmv+YiL1re4eO0cp4ZnuXhwlgek93NBss3RZMyGZGzYl5UpbFsjYi5JJXYTsDZ2Y69VFxFgY08dw3Raz2DL9sdt2x9HY9Nmk0mpKg8GiLND57npj3jRALY/zgWRaUZERET0haJZNv2yFUQcNCMiIuaPmBpuyZBAMjCrTrokFiYcxxjcL07PcnG6yQmZcCwR1iUltUkkxpqTk5NJxpqWITnmXFKo5pNJAmMxTqBxVTUfOmeCU803t9HtHaMOjUalQyFJkLWhCfq26lDV6E6hlusQ8iHoWs5wLePI+ogT6zucXNviEk81f8Dgfi5KtjkhE04kwtDmAc1VGYuS5BkZwpAyDCnDTg1Vl/RBCoeCkcvKtqUMNj3VfHPHqK87I9ix6nmem6BvEWSYmXUMs2FhdgA/iYUUjpLcOkoGa0aFvXCtNDtcOryPi61J5aJkhwuTjA2pyratOSPJ2GDMmmReDk0j2yRPmExS8omYgH3nJNn2VNit3Kjmm2Nk28q2uWXNDlkpW2rKlSRB0tTo0ACJFI6SQrYhJoxq3TiBTqyPuHDNmB0uHd7PqcH9PHBwhgemZzmRjDkqMPRkG2tOqjl5nrFtHXdgwt/GmhROoPE4hbE/CaGcXDHYtqaic2PTH89umX7o1PPxuOyPmqOsIWlWVZ+tqcj1x7ngkIYcLXQ1ShG5VkT+TkRuEZGXNJxfF5E32fMfsIsfRURErDgU0Fx7bauGhTFNEUmBVwNfC9wKfFBETqvqzd5lzwfuUdVHisj1wM8A39n54ATkyIT1I2NOHt3i0iNnuHzjPi5bv4+r1u7kAen9XJxs84BUOSprrMuAVMy3YawZCRnndESOMCLlXL7GfdlR7pkc5e6dY9y7fYQzmxuMNtdIz6YMzgnDs8La/craWWV4Nmd4dsLwni3DxDa30TNn0dGIfDQ2TiDzAgxLSVNkkNug4hQdpuTrKZMNYXxEmBwVJkdhclTJjmckx8ecOL7FqaObXHr0fh529E4eOLyfywb38sD0DA8abHFCEtZljaPJWvFadnTMZj42DExhTMq5fJ37siPcNz7CvTtHuHfrCFuba+TnhqRnUobnhOFZGJ41sq3dN2FwbkJ6Zodkcxs9c458exsdjQxzccwkSZHhwMiXK5om6CAhH6ZkGwmTjcTKBuNjMDmWkx3LWT++wwXHtjl19BwPOXY3l63dz4PW7uGq4Z1cnG5yUTLhhCQcT44UTCzTnB0mZIxtG6Zs65Az2Qb3TY5w7/gI9+wc5czWOjtbQ/TcgMHZhOHZUrb1+3IG5zIGZ0ek925azWCbfGeEbm0Z25tLmTYYIIMBrA2LQHAdpuRrCdl6wmSjbLPxMSU7nsPxMRvHRlx4bIsHHbuv6I8PX79jan8ck7GtIyMrCdv5kLNWtvvGR7hvZ4OzW+uMt4aV/jg4C+v3l/1xcO/2nvpjtl7tj3tGsdTp4cMimeYTgFtU9dOqOgJ+C7guuOY64A3295uBp9n1hyMiIlYcmmW9tlWD6ILCAkTkWcC13rrD/wZ4oqq+yLvm4/aaW+3+p+w1dwbPugG4we4+hh4Lui8RTgF3Tr1qObBKdYXVqu8q1RXgi1T1xG5vFpE/wsjcB3eq6rW7LWu/sRKOIFW9EbgRQERuUtVrDrhKvbFK9V2lusJq1XeV6gqmvnu5f5UGwVmxSPX8NuBKb//B9ljjNSIyAC4E7lpgnSIiIiL2hEUOmh8ErhaRh4nIGnA9cDq45jTwHPv7WcC7dFH2goiIiIg5YGHquapORORFwDswMxhfp6qfEJGXAzep6mngV4A3isgtwN2YgXUablxUnReEVarvKtUVVqu+q1RXWL367hsW5giKiIiIOIxYaHB7RERExGFDHDQjIiIiZsDSDpqrNAWzR12fKyJfEJEP2+0FB1FPW5fXicgdNka26byIyH+zsnxURL5iv+sY1GdafZ8iIvd57/al+11Hry5Xisi7ReRmEfmEiPxQwzVL8X571nVp3u1SQVWXbsM4jj4FPByT1fAjwJcG13w/8Br7+3rgTUtc1+cCP3/Q79XW5WuArwA+3nL+G4C3Y1Yv+hfAB5a8vk8B/uCg36uty+XAV9jfJ4C/b+gLS/F+e9Z1ad7tMm3LyjRXaQpmn7ouDVT1zzCRCm24Dvg1NXg/cJGIXL4/taujR32XBqr6OVX9a/v7DPA3wBXBZUvxfnvWNaIByzpoXgF81tu/lXqDFteo6gS4D7hkX2rXUg+LproC/Curjr1ZRK5sOL8s6CvPMuGrROQjIvJ2EXn0QVcGwJqLHg98IDi1dO+3o66whO/2oLGsg+Zhw+8DV6nqlwF/TMmQI/aOvwYeqqpfDvx34C0HXB9E5DjwP4AfVtX7D7o+XZhS16V7t8uAZR00V2kK5tS6qupdqrpjd18LfOU+1W036PPulwaqer+qnrW/3wYMRaRvooi5Q0SGmEHoN1T1dxsuWZr3O62uy/ZulwXLOmiu0hTMqXUNbFbPxNiPlhWnge+xXt5/Adynqp876Eq1QUQuc7ZsEXkCpk8fSP4CW49fAf5GVV/VctlSvN8+dV2md7tMWMosR7q4KZgHVdcfFJFnAhNb1+ceRF0BROQ3MV7RUyJyK/BTwBBAVV8DvA3j4b0F2ASedzA1NehR32cB3yciE2ALuP6APp4ATwL+DfAxEfmwPfYfgYfA0r3fPnVdpne7NIjTKCMiIiJmwLKq5xERERFLiThoRkRERMyAOGhGREREzIA4aEZERETMgDhoRkRERMyAOGhGREREzIA4aEbUICIXicj3e/sPEpE3L6Ccl4nIbTamte2aR9i0ZGfnXX5ExG4Q4zQjarAJHP5AVR+z4HJeBpxV1f+7x7VnVfX4IusTEdEHkWlGNOGVgGN4PysiV7kkwGISKr9FRP5YRD4jIi8SkReLyIdE5P0icrG97hEi8kci8lci8uci8sXTChWRJ3sJbz8kIicWLGdExMxYymmUEQeOlwCPUdXHQcE8fTwGk0psAzMd8MdU9fEi8nPA9wD/D2Y1wxeq6idF5InALwBPnVLujwA/oKrvtdl3tuckT0TE3BAHzYjd4N02ce0ZEbkPk/oO4GPAl9kB76uB3/HyQq/3eO57gVeJyG8Av6uqt8653hERe0YcNCN2gx3vd+7t55g+lQD3OqbaF6r6ShH5Q0xCi/eKyNep6t/Oo8IREfNCtGlGNOEMZt2YXcEms/0HEfl2KBYT+/Jp94nII1T1Y6r6M5iUe1PtoBER+404aEbUoKp3YZjex0XkZ3f5mH8NPF9EPgJ8gn7rJv2wLfOjwBizAFlExFIhhhxFHBhiyFHEKiIyzYiDxFnghj7B7cDn969aERHtiEwzIiIiYgZEphkRERExA+KgGRERETED4qAZERERMQPioBkRERExA/5/u0y5nRs7GRcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "k = 0.2\n",
    "t1 = 0\n",
    "t2 = 2.5\n",
    "x1 = 0\n",
    "x2 = 1\n",
    "nx = 100\n",
    "x = np.linspace(x1,x2,nx+1)\n",
    "dx = x[1]-x[0]\n",
    "\n",
    "dt = 0.005\n",
    "nt = int((t2-t1)/dt)\n",
    "alpha = dt*k/dx**2\n",
    "\n",
    "# allocate the temparature array\n",
    "T = np.zeros(nx)\n",
    "\n",
    "# allocate an array to plot the image\n",
    "TasImage = np.zeros([nx,nt+1])\n",
    "\n",
    "# build matrix\n",
    "A = np.zeros([nx,nx])\n",
    "for j in range(1,nx-1):\n",
    "    A[j,j-1] = -alpha\n",
    "    A[j,j] = 1 + 2*alpha\n",
    "    A[j,j+1] = -alpha\n",
    "# first and last element are different because of the boundary conditions\n",
    "A[0,0] = 1\n",
    "A[nx-1,nx-1] = 1\n",
    "\n",
    "# add to the image\n",
    "TasImage[:,0] = np.flip(T).copy()\n",
    "\n",
    "for i in range(1,nt):\n",
    "    t = t1 + i*dt\n",
    "    b = T.copy()\n",
    "    T = la.solve(A,b)\n",
    "    # re-adjust boundary condition\n",
    "    T[0] = np.sin(10*t)\n",
    "    T[nx-1] = 0\n",
    "    # add to the image\n",
    "    TasImage[:,i] = np.flip(T).copy()\n",
    "       \n",
    "# visualize the result\n",
    "plt.imshow(TasImage, aspect='equal', extent=[t1,t2,x1,x2])\n",
    "plt.axes().set_aspect((t2-t1)/(x2-x1))\n",
    "plt.xlabel('time [s]')\n",
    "plt.ylabel('distance [m]')\n",
    "plt.colorbar()\n",
    "plt.text(t1+0.06*(t2-t1),0.9*(x2-x1)+x1,r\"$\\alpha=\"+str(alpha)+\"$\",size=16,color='k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes about the example\n",
    "\n",
    "Some things to note:\n",
    "\n",
    "  * When solving $(1 - \\Delta t F )\\cdot T = b$, **do NOT invert the operator A**.\n",
    "\n",
    "    Using the inverse matrix may seem to be the intuitive thing to do, and scipy.linalg\n",
    "    can do it, but it is incorrect.  There are better, faster, and\n",
    "    more accurate algorithms for solving this problem.\n",
    "\n",
    "    *[www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix](https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix)*\n",
    "\n",
    "  * Because we are using an implicit method, this solution will be\n",
    "    **numerically stable** for arbitrary timestep size (this can be\n",
    "    demonstrated mathematically).\n",
    "\n",
    "  * However, for **accuracy** of solution one\n",
    "    must have $\\alpha <1$.\n",
    "\n",
    "  * We’ve used the most general method to solve $Ax = b$. However, you\n",
    "    will notice that our operator is banded. In such situations, there are\n",
    "    **special algorithms** for solving $Ax = b$, the one to use in Python is\n",
    "    `linalg.solve_banded`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the explicit method\n",
    "\n",
    "As mentioned last lecture, explicit methods can be unstable. Let’s see how\n",
    "it performs in this case. We start by returning to our equation:\n",
    "\n",
    "$$\\frac{\\partial T_i}{\\partial t} = F \\cdot T_i$$\n",
    "\n",
    "We rewrite this as\\tightsep\n",
    "\n",
    "$$\\frac{T_{i+1}-T_i}{\\Delta t} = F\\cdot T_{i}$$\n",
    "\n",
    "which, after some rearranging, gives:\n",
    "\n",
    "$$T_{i+1} = (\\mathbf{1}+\\Delta t F)\\cdot T_i$$\n",
    "\n",
    "\n",
    "This is obviously much more direct: just **matrix-vector multiplication**.\n",
    "\n",
    "Boundary conditions, surprisingly, are dealt with in the same way as\n",
    "for the explicit method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explicit results\n",
    "\n",
    "Implementation:\n",
    "\n",
    "```python\n",
    "# heatexplicit.py\n",
    "...\n",
    "for i in range(1,nt):\n",
    "    t = t1 + i*dt\n",
    "    # explicit solve\n",
    "    b = T.copy()\n",
    "    T = A @ b\n",
    "    # re-adjust boundary condition\n",
    "    T[0] = np.sin(10*t)\n",
    "    T[nx-1] = 0\n",
    "...\n",
    "```\n",
    "\n",
    "![](xalpha10.0.png)\n",
    "\n",
    "There’s a problem here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "k = 0.2\n",
    "t1 = 0\n",
    "t2 = 2.5\n",
    "x1 = 0\n",
    "x2 = 1\n",
    "nx = 100\n",
    "x = np.linspace(x1,x2,nx+1)\n",
    "dx = x[1]-x[0]\n",
    "dt = 0.005\n",
    "\n",
    "nt = int((t2-t1)/dt)\n",
    "alpha = dt*k/dx**2\n",
    "\n",
    "# allocate the temparature array\n",
    "T = np.zeros(nx)\n",
    "\n",
    "# allocate an array to plot the image\n",
    "TasImage = np.zeros([nx,nt+1])\n",
    "\n",
    "# build matrix\n",
    "Atilde = np.zeros([nx,nx])\n",
    "for j in range(1,nx-1):\n",
    "    Atilde[j,j-1] = +alpha\n",
    "    Atilde[j,j] = 1 - 2*alpha\n",
    "    Atilde[j,j+1] = +alpha\n",
    "# first and last element are different because of the boundary conditions\n",
    "Atilde[0,0] = 1\n",
    "Atilde[nx-1,nx-1] = 1\n",
    "\n",
    "TasImage[:,0] = np.flip(T).copy()\n",
    "\n",
    "for i in range(1,nt):\n",
    "    t = t1 + i*dt\n",
    "    b = T.copy()\n",
    "    # explicit solve (skip to avoid nans)\n",
    "    if (t<0.99 or alpha<0.11 ):\n",
    "       T = Atilde @ b\n",
    "    # re-adjust boundary condition\n",
    "    T[0] = np.sin(10*t)\n",
    "    T[nx-1] = 0\n",
    "    TasImage[:,i] = np.flip(T).copy()\n",
    "      \n",
    "plt.imshow(TasImage, aspect='equal', extent=[t1,t2,x1,x2])\n",
    "plt.axes().set_aspect((t2-t1)/(x2-x1))\n",
    "plt.xlabel('time [s]')\n",
    "plt.ylabel('distance [m]')\n",
    "plt.colorbar()\n",
    "plt.text(t1+0.06*(t2-t1),0.9*(x2-x1)+x1,r\"$\\alpha=\"+str(alpha)+\"$\",size=16,color='k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Can we escape the instability?\n",
    "\n",
    "Yes, we can, we just need to pick a better value of α, i.e., a smaller timestep for the implicit method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "k = 0.2\n",
    "t1 = 0\n",
    "t2 = 2.5\n",
    "x1 = 0\n",
    "x2 = 1\n",
    "nx = 100\n",
    "x = np.linspace(x1,x2,nx+1)\n",
    "dx = x[1]-x[0]\n",
    "dt = 0.00005\n",
    "\n",
    "nt = int((t2-t1)/dt)\n",
    "alpha = dt*k/dx**2\n",
    "\n",
    "# allocate the temparature array\n",
    "T = np.zeros(nx)\n",
    "\n",
    "# allocate an array to plot the image\n",
    "TasImage = np.zeros([nx,nt+1])\n",
    "\n",
    "# build matrix\n",
    "Atilde = np.zeros([nx,nx])\n",
    "for j in range(1,nx-1):\n",
    "    Atilde[j,j-1] = +alpha\n",
    "    Atilde[j,j] = 1 - 2*alpha\n",
    "    Atilde[j,j+1] = +alpha\n",
    "# first and last element are different because of the boundary conditions\n",
    "Atilde[0,0] = 1\n",
    "Atilde[nx-1,nx-1] = 1\n",
    "\n",
    "TasImage[:,0] = np.flip(T).copy()\n",
    "\n",
    "for i in range(1,nt):\n",
    "    t = t1 + i*dt\n",
    "    b = T.copy()\n",
    "    # explicit solve (skip to avoid nans)\n",
    "    if (t<0.99 or alpha<0.11 ):\n",
    "       T = Atilde @ b\n",
    "    # re-adjust boundary condition\n",
    "    T[0] = np.sin(10*t)\n",
    "    T[nx-1] = 0\n",
    "    TasImage[:,i] = np.flip(T).copy()\n",
    "      \n",
    "plt.imshow(TasImage, aspect='equal', extent=[t1,t2,x1,x2])\n",
    "plt.axes().set_aspect((t2-t1)/(x2-x1))\n",
    "plt.xlabel('time [s]')\n",
    "plt.ylabel('distance [m]')\n",
    "plt.colorbar()\n",
    "plt.text(t1+0.06*(t2-t1),0.9*(x2-x1)+x1,r\"$\\alpha=\"+str(alpha)+\"$\",size=16,color='k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But a smaller time step slows down the time-to-solution substantially.\n",
    "\n",
    "****Implicit methods are more complicated, but worth it.****"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final notes on PDEs and linear algebra\n",
    "\n",
    "  * We only scratched the surface for PDEs. We've not discussed\n",
    "    finite-volume vs. finite-element, or spectral methods.  Nor have\n",
    "    we looked at hyperbolic and elliptic equations.\n",
    "    \n",
    "  * There is a lot more linear algebra in SciPy.    \n",
    "    *Decompositions, inverses, determinants, eigenvalues, special matrices*    \n",
    "    *[docs.scipy.org/doc/scipy/reference/linalg.html](https://docs.scipy.org/doc/scipy/reference/linalg.htmlhttps://docs.scipy.org/doc/scipy/reference/linalg.html)*\n",
    "\n",
    "  * SciPy's routines use low-level routines from BLAS and LAPACK,\n",
    "    which often use highly optimized codes (MKL, ESSL, ...).\n",
    "    \n",
    "  * Implicit is better (faster, more stable) than explicit.\n",
    "    \n",
    "  * Never invert a matrix numerically."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}