{ "cells": [ { "cell_type": "markdown", "id": "81b1d2f0", "metadata": {}, "source": [ "---\n", "title: \"P3: Chapter 3\"\n", "subtitle: \"Problem Class 3: Finite Differences\"\n", "abstract: \"Exercises based on content from Lectures 12-15\"\n", "date: 2026-04-06\n", "format:\n", " html:\n", " other-links:\n", " - text: This notebook\n", " href: P3.ipynb\n", "---" ] }, { "cell_type": "code", "execution_count": 2, "id": "d3fe9815", "metadata": {}, "outputs": [], "source": [ "using Plots, LinearAlgebra, LaTeXStrings" ] }, { "cell_type": "markdown", "id": "13982f26", "metadata": {}, "source": [ "::: {#exr-}\n", "\n", "So far, we have considered homogeneous boundary conditions $u(0) = u(1) = 0$. How would you approximate the problem\n", "\n", "\\begin{align}\n", " &-u'' + p u' + q u = f \\qquad \\text{on }(0,1) \\nonumber\\\\\n", " &u(0) = \\alpha, \\qquad u(1) = \\beta\n", "\\end{align}\n", "\n", "with finite differences?\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "049ae0d3", "metadata": {}, "source": [ "::: {#exr-}\n", "\n", "So far, we have considered ADR\n", "\n", "\\begin{align}\n", " -u'' + p u' + q u = f \n", "\\end{align}\n", "\n", "with constant coefficients $p,q$. How would you write down a finite difference approximation to this equation if $p$ and $q$ are functions of $x$?\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "8ef63350", "metadata": {}, "source": [ "::: {#exr-}\n", "# FNC example 10.4.2\n", "\n", "Following on from the previous two questions, we try and solve the problem\n", "\n", "\\begin{align}\n", " &-u''(x) + a \\cos(ax) u'(x) - a^2\\sin(ax) u(x) = 0\\nonumber\\\\\n", " &u(0) = 1 \\nonumber\\\\\n", " &u(\\tfrac{3\\pi}2) = e^{-1}.\n", "\\end{align}\n", "\n", "Show that the exact solution is $u(x) = e^{\\sin ax}$.\n", "\n", "Spot and correct the error in the following code [we did this in class - instead tyr and understand the following code]:" ] }, { "cell_type": "code", "execution_count": 53, "id": "29825a8f", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdZ1gT2dsG8JNGQgsQpHdBRFSqgIAoioIFUREbFuy6Flxd267uyt/ee2+sDQt2QFFRaXYFFRFUUASVHpIQICQh836YfbNZbJTAAHl+lx+SycnMnWPgYdo5JAzDEAAAAKCoyEQHAAAAAIgEhRAAAIBCg0IImtzx48f9/f2lTzEMCwgI+Pvvv4lLJE979uzJzs6uY+PExMSoqKgGb6ugoMDV1fX58+cNXkNLcPfu3atXr9axcUFBwebNm5s0DwBQCEHTwjBs9erV5eXl0iVJSUlRUVEkEonAVPKyf/9+iURiaWlZx/Y9e/a8c+fOnTt3Gra5w4cPP3/+XEVFpWFvl4v3799zOJwGv/3Zs2fHjh0bPHhwHdvr6+ubmZmtXr26wVtsChiGffny5dOnT3CNRRuBAdCU4uLiqFRqamqqdMmYMWNcXV1ramoITCUXaWlpI0aMqO+7BAKBn58fm82u7xtramrMzMyWLFlS3zfK0YsXLygUSrdu3Rr29srKyt69e3O53Pq+cdasWUlJSQ3bqNzFx8cHBAQcOnTo8OHDI0eOTExMJDoRaCwohKBpjRgxYv78+dKnxcXFqqqqsnWxLnJyckQikbyjNYpEIvH29q7vB8GdO3du2rRp9X1XVFSUiYlJeXl5A7YoL7m5ucbGxuPGjWvY25cvX75jx44GvPHz589du3YVCAQN264cvXnzpl27doWFhfjTwsJCY2Pjp0+fEpsKNBIUQtCE8vPzTUxMOByOdMmGDRvmzZtX3/UcPHgwPz9frtEa68qVK35+fg17r1gstrS0zMrKqte7/P39L1y40LAttgTFxcUGBgZ8Pr9hb58wYcKRI0fkG6kBxo0bN3z4cNklM2fO9Pf3JyoPkAs4Rwia0NGjRzdu3KihoYE/xTAsOjo6LCysvutpgScUjxw5MnTo0Ia9l0KhDBw4MDw8vO5vycvLq6mpCQwMbNgWW4KTJ0/27NlTVVW1YW8fOnTokSNH5Bupvmpqaq5du2ZnZye70MLC4vbt25WVlUSlAo1HJToAaMVqamoiIyMFAkFwcLCSkpJ0YXZ2trW1NUJIIBCMHj1a2j4tLS00NFRTU7M5Q4pEosjIyJSUFGVlZQ8Pj/79+0vLKoZh169ff/r0KZlMtrCwGDFihPRTIISio6MLCwtFIpGBgQGJRHJ3d9fR0cFfEgqFt2/fXrt27deb43K5169fr6ioGD58uKam5qtXr+7evevs7Ozh4SHbzNPTc+PGjbLXgKSkpCQnJw8ePNjCwkK6MCsry9TUVElJ6e7duzt37pRjt/zU69evk5OTSSQSmUx2dnbmcrndu3e/fft2aWkpQmj8+PEIITabfe/evdLSUhaLFRAQ8PLly4SEBA0NDV9fX319/VorvHbtWkBAwNcbqqmpuXXrVlZWVv/+/a2srAoLCy9evKinpzd06FAy+d+/1D09PUeMGMFms1ksFr4kPz//woULDg4OPXr0kDYrKiqiUCja2tpy7xCE0MuXL9lsdq2LldTV1auqqp4+fdqzZ8+m2ChoDkTvkoJWbMWKFenp6QsXLgwNDZUu3L17d7t27eS7oUOHDjXs0OiXL1/s7e2XLFnC4/GEQmFERMTEiRPxl7hcrq+v7+7du4VCIYZhMTExrq6uHz58wF+dOXNmbGws/vjVq1csFuvVq1fS1T569EhdXV0ikdTaXEZGxrJly758+fL48WMbG5vIyMjDhw+/e/fO3Nz82rVrsi2zs7PJZHJlZSX+9MaNG+Hh4ampqXp6elVVVfjC0tJSOp0eHR3dgA/eSFeuXFm4cCF+XpbL5fr4+Kxdu7aiomLFihVmZmadO3fGm33+/HnNmjW6urpDhgw5ePDgmTNneDzehQsXWCxWRkaG7AolEom6uvrjx49rbaiiomLRokWpqaklJSWdO3eOior6/fffS0pKfH19Fy9eXKuxmZnZ9evX8cd5eXl//vknh8MxNDR88uSJtI2rq+uCBQvk2xtSsbGxCKG9e/fKLjxx4gRCKDIysok2CpoB7BGCBrp9+7ajo6Otre379++rq6uly2/duuXo6EhgMCmJRBIYGGhjY7N+/Xp8SURERHR09KZNm9q1azdhwgQDA4PZs2fjLw0cODAtLW3QoEEpKSl8Pv/KlSv79u3DX+rcufP06dNl1/z27VsdHZ1aB2zFYvGePXt27NhBJpMNDAw+f/588uTJy5cvHz58+NOnT+rq6rKNdXV1JRLJu3fv7OzsxGJxbGzs1q1bo6KiioqKysvLGQwGQighIaG6utrBwaGJ+ucHNm3atHXrViqVihBiMpmrV6++e/euiopKWFgYm82W3v5haGj4xx9/pKenP3nyJCQkZNiwYQihwMDA33777e+//5Z2O0IoPz+/vLxcV1e31obWrFmzaNEifFfbyspq6tSpHz58+Pz5c1xcnLe3d63Gurq6b9686d+/P0Jo165dYWFhQqGwsLDw/fv33bp1Qwix2eynT5/OnTv3mx/q/fv3kZGRdfn4AwYMqHX8E1dWVoYQqnV0F/+fxXeUQSsFhRA0UFJS0vLly6urq+/cufPXX3/hCyUSSVJS0oIFC4jNhrt+/frDhw8XLVokXTJ9+nQXFxcWi/Xw4cMrV66cP39etv3QoUOXLl164sSJESNGFBUVLV269Ndff8UP8QUHB8se6ysqKpIeoJNKSEgICgrCj+YVFxeXl5fjp/SmTp0aFBRU64CwmpqakpJSSUkJQujhw4e+vr4IIfz8k/QAbEJCgpWVlZGRkfy6pK7odPqSJUu2bNmC/03j6OiI/f8NczQarVZjGo2Wn58ve8ZUR0enVmEoKipCCNU6YllRUcFkMqWfNzs7e8CAAcrKylZWVqWlpV8fQtfS0sJXKxQK1dTUlJWV8dEJevfujTdISkqSSCS9evX65odq3779kiVL6t4JX5NIJF8vxP8eEovFjVkzIBYUQvCPsWPHPnz48KfN8ANcSkpK+DUv0dHRXC4X3xVA/38SRfqLqQHEYvG2bdtq/cZ58uRJTk5OrZ0qKyur4cOH/2BVCQkJCCFTU1PpksGDB+O3cuP7NLXOY+ElJy4uburUqWvXrv399983bNhgbW3t5+e3cOFC2V/iAoFATU2t1uZ8fHykjx89eoQQkv5G/uZpUXV1dfwKC/wUl0gkOn/+vOzeTHx8fL168sWLF/gFjT9tGRoaOm/evB80WLduna+vr5OTk7a2to+Pz/Tp02U/3ddMTU1l949JJFJFRYVsg6qqKhKJVGtfSlVVVVqZysvLMzMzFy5ciD/9Zo8xmUx8tUpKSn/++SdC6PTp0z169Kj1p4OJickPojZGu3btEEIikUh2oVAolL4EWikohOAfp06dasC7zpw5061bN3Nzc/xpQkKCmpqai4tLg2NQqVTZfTjc4cOH/f39v77+4sfw4Wyk16zK4vF46KudG/xKGS6XixBavHixr6/v5cuXk5KSDh48eOLEiZSUFOllLHQ6nc/n/2DTd+/ebd++vZmZ2Y/jKSsrS5/GxcWVlJQEBQXhT9lsdlpa2tKlS+vwQf9hb2+flZVV9/Y/4OLi8u7du4iIiOTk5Dt37pw7dy48PHzixInfa0+hUGotqVWPGQwGhmEVFRVf/wGBS0hIEIvFXx8OlVVeXi57oQqXy42NjZUdgC0+Pr5Pnz4/WEMj4RVXIBDILsSffn2EALQiUAhBozx48AC/gBCXkJDQo0cPGo325s2b9u3bf30YrTnhI58VFBR06NCh1kvt27dHCLHZbNmF+GE3S0vLT58+PX/+3N/fHz8/V1hYOGTIkF27dm3duhVvyWKxfjzM2N27d2UvZUxPT+/cubNsAz6fLxQKZfcyHzx4YGxsbGtriz/FB1Lp3bu3WCzOysqysbGp34dvnP3798+cOTM0NDQ0NFQsFi9dunTVqlU/KIQ/hdeJsrKy7xXC+Ph4ExMT6Z8Or1+/7tSpU62zsGw2W7bHUlJSBAJBv3798KccDufly5eLFy9GCL169apLly61NvHx48cDBw588/CmLBKJFBwc3LVr169f6tChA41Gww/zShUUFFAoFGdn5x+vFrRkUAhBw4nF4ry8vE6dOuFPa2pqEhMT8WNuFy5c+OOPPwhNh0aNGvXnn3/GxcV5eXlJF+KnkQIDAxctWnTv3j38ygvcgwcPSCRSSEgIl8s9cuSIdKBwPT29X375JTExUdrS2toaP70nKyoqis1mh4SEFBYWvnjxYsKECfhyNpsdGRlZqxAWFxeTSCTZCp2TkyOtguj/9yn19PSSk5MlEknzF8IJEybgu19UKnX58uXHjh1rzApNTEyUlZVLSkpkj1ump6dfvnx50aJFSkpKN27ckL3Gav/+/V/fLlJUVFSrx+h0unSg18TExJqaGg8Pj+rq6tjY2K8LoZmZ2TfveKk7NTU1Hx+fzMxM2YUZGRmenp6wR9iqwQ31oOGoVGqHDh2kl4zu2rVLIBBYWVlVVla2hGsHTExMdu3atXPnTul0DVVVVSdPnvTw8GjXrl14eHh4eHhubi7+EpfLXbVq1apVq/DjulFRUbJnTJ8/fy49D4oQsre3r6yslD0OiWHYyJEj9+zZgxA6dOhQ+/bt8d+MEolk06ZNoaGhtbI9efKkS5cusufMbGxspMfcXr9+HRMTY2VlhRCKjY2tdQ9iM6ioqJCtGampqdJrYbhcbq3Dwlwul8fjyR4L5XK5tW4wx29GfPr0qezChQsXrl69msvl3r59GyGkpaWFLw8PD5e9/RRXXFycl5eHXx2Ks7GxkUgk+Ck6Pp+/d+9eOp1uamp69erVH5/RbIyVK1fGxcUVFhbiT9lsdsPGiAAtCqkup9YB+J7k5OSwsDBvb+/q6mo/P79Xr15du3atffv2y5Ytk17C0HgNO0eIu3bt2vr1662trc3NzYVC4axZs6TrefDgwY4dO/T19Wk02vv378eOHYtf55mZmbl69WpPT08Oh0OhUMrKyqysrKZMmSK7Wh8fnwkTJoSEhEiX/P777xKJRFlZ2cnJydzcfMWKFZ6ennw+f9y4cXhJkxUaGkqn0zdt2iRdwuPxJk+ebGpqqqmpqampOWDAgJCQEB8fHycnJ9ka3Dz8/f1DQkJev36trKwsEon4fP4ff/zBZrNnz55dVVVFJpMxDFu/fj2DwViyZIlAICCTySQSae7cuVQqddu2bRiGSSQSVVXVsLAwJycnfJ2rV6/OysqSnX7rwoULMTExHTp0UFJSmjp16tSpU+3s7Gg0mpOTE34ZrayLFy+uXbu2Vin9888/s7Ozu3TpIhAI5syZM2nSpI4dOxoYGHx9mlmOLl26dODAgV9++YVCoezZs2fs2LHjxo1rus2BZgCFEMiBRCKRDgKCYZjcR0RrTCHEySas16vf+zinTp26cuXKuXPnftD4e++VSCQdO3a8fPlyreOlXydpis6sFzkGyMvL8/T0zMrKkh2+p9YmfrC5qVOn2tvbf/Mewab++n1NLBa/fPlSIpHY2toSOysWkAs4NArkQPZ3d1P8GnJ3d2/kOZgfVMEfv/q9jzNixIjMzMwPHz78oPH33nvp0iVnZ+evq+DXSQgfZFWOAUxMTPr373/mzJkfbOJ7myssLExKSpo8efI3X23qr9/XqFSqk5NTt27doAq2DbBHCEADJSUlnThx4uDBg/V6V01Nja+v74kTJwwNDZsoWItVVlY2fPjw69ev0+n0er1xwYIF3t7e3xyqFIDGgz1CABrIy8tLR0fn7Nmz9XrXX3/9NWPGDAWsggghLS2tRYsWSe+ar6PY2FihUAhVEDQdKIQANNzq1avT0tKys7Pr2P7u3bvm5uYjR45s0lQt2YABAzw9Pa9cuVLH9gUFBXFxcTt27GjSVEDBwaFRABqr7hdoEH7xSwsBPQZaFCiEAAAAFBocGgUAAKDQoBACAABQaFAIAQAAKDQohAAAABQaAYXw+vXr+LSlddQShm9uG6An5QV6Ul6gJ+UFerIxCCiEt27dun//ft3bV1VVNV0YhQI9KS/Qk/ICPSkv0JONAYdGAQAAKDQohAAAABQaFEIAAAAKDQohAAAAhQaFEAAAgEKDQggAAEChQSEELYtYLC4rK/veqxFnI32Hj/s9bA1cLA4AkBcohKAFuX0n3rizs41PUHefgQKBoNarcbdvz95y/JbDkm1vaNPnLyEkIQCg7YFCCFqQWYv/LJx5vWhW7Attz4gztWd+v5Vwn+M2GRnYVPeZl3y/HoMTAQDAD1CJDgDAv4QiEWKoIYQwOhNxS6rfp0vK2TVcNiaqJiurjTBT+fgsnMMS5X144+zW7eu3nzwTufR/6xBCK5YumBYyrrnTt1QCgSA5OZnoFN9WWVmpoqJCdIq2QBF60s3NTV1dvSnWDIUQNCsMw1JTU5WVlTt16lTrJYmgcvekwa9uDevUTt1Go0qNo8+L/kJWZ1E0WCQaXVxaYKlOW+xm8OnTCQttCosqKdw0i6ZvrmTagdGlO1XboKKiYkHYuuLQu4hEXraxb6D/AG1tbUI+Y0tz69atKVOm2NvbEx3kG2ACenlp8z35+vXrtWvXhoSENMXKoRCC5iORSPoMDkqvVEUC3mAny6N7tiKEJBW8qudJVWn3hTkZbpZd7ccHFyKGSZ+ByhpaX69BCyEnhBBCmEgoLswVfflQnZNRfvscWV1LZNalk7FRsZIKQoikY1lcXAyFECeRSDw8PC5fvkx0EAAabvLkyRKJpIlWDoUQNJ9Xr16lV6qUjD2KEIre4c3JfC5JuS1Ie8CwdVV1H6A9+U+SEgMhZFyHVZFoSjRjK5qxlYprPzRirjAnoyrt/hZrUVnKqKPlhm/EBdbW1k38aQAAbQQUQtB8GAwGScBDCHnxU5c6M6rObVPrMVhzyDSyKrNR6yWRlCxslSxs7fzGJx7bu+zTS1Z7g4r4i6oeA8gMVflEBwC0XVAIQfOxtrae7mJqnzDAkIFyO3Q3+GM1kutZDSU6ve/0+Qgh0afs8rvnC1ZNUu87Sq3nEBLl3++5RCIhk+FiaQDAv6AQgmZSU1bEvXZstiabGjBDxc3PU71xe4E/RDO2ZI1fIi7+zLm0v+LBdc1hMxidXN68eTNo9ERetdjKWP/WpTOqqrCzCABACAohaA4Yxr8Xzbt+Qs0rQGvZHBJduXk2S9Uxajd9leD1E86lA9SkqD+upWcH7EIm9tyEfTv3Hfx94fzmiQEAaOGgEAI543K5L1++7Nixo66uLkKopqyIfXobJhTo/rqNqmPU/HkYti56HR35dy/8z+gxUsq7iOyFWqaFpSnNn6RVY7PZqamp2dnZRkZGdnZ2JiYmRCcCQG7gZAmQpzdv3ti4eQ/dcL5Lr4F34xOqnicVbgmlt++sG7qFkCqII1Go6n1HZbkEznq/7/DjidZ3/jdrEtxuX1cCgWD+/PmGhoZ9+/adMWOGv7+/mZnZ8OHDCwsLG7lmDoezb9++DRs2lJWVrV69ms/nf69lcXHxli1bGrk5hNCePXtKSkoQQhcvXnz8+DG+8P3799u3b9+xYwdCiM1m79mzZ8OGDeXl5Y3fXKu2ZcuW4uJiea2ttLR08+bN8lqb3EEhBPK0cfehgoFr2UM3V44PZ5/YyLt1RmfWOmb/cagFXJ8ydOJ004XbHayN7gTamlGqiY7TOtTU1AwdOnT79u3VQiEyd0ZuY5B1T4xCu3jxopeXF15UGmzEiBHPnz/X0tKiUChnz579wUDqPB5Peh/kwIED371717AtrlmzpqioCCF0//59fCUCgcDT07OiokJTUxMhFBAQkJGRoaWlBVdUXbp0icfjNWYNY8aMefr0Kf6Yx+NdvHhRHrmaBBwaBfLEVFMhl3GMhQVHSjYVkWi6C3bIXrFJuA6dbFGnjYLMZ+xja9X9xqp5+hOdqKU7cuTIjRs3kIoWmnsJWXn8s7T4A9o99N27jOXLl+/fv7/BK3/06FF4eLixsTFCKC0tTfYlkUhUUlJiYGCAP7W0tExKSsIfp6Wl1SqZGIYVFRXp6up+b2iVqqqqyspK2QEWpHsnHz9+RAgtW7YMf/r48ePLly+3a9dO2rK4uFhLS4tK/c/XWCQSFRYWGhkZ4VssKirS1tamUCiybaqrqzkcjp6enuxCiURSWFjIYrHodLp0IZfLpdFo3xsgTSgUlpSU6Orq1sqAr0pPT69Wzebz+WKxGK/rsr58+cJisRgMBv6hNDU1aTTaN7eIqzUsn1gsLioq0tPTo1AoGIYVFBTo6OjIRiovL6+qqsJPiOAyMjKkO9YWFhb379+XXWFRURGLxar1oWpqaoqKivT19Zt5lBxF/6sHyNcfC+YGvt13+dW0O28+OC1c36KqoBTDxlln3taK5OiyczuxGjHRcVq08PBwhBAKWvdvFUQI6VigKX8jhE6dOvX1JCF1ZG1tzefze/To4ePjgxBiMpkFBQUIoYEDB86ZM8fBwaFXr17W1ta5ubkIoYyMDHNzc4TQnDlzCgoKBg0aZGlpeeHCBYTQ3r17TUxM/Pz8zM3Nr1279vWGli9fbmZm1qdPn0mTJmEYhi8MCQnZt29fenp63759S0pKLC0t582bZ2lpKRKJXFxcBg0ahBC6fft2hw4d+vbta2pqumnTJvyN3t7ec+bM6dSpk5eXV1VV1cWLF83NzX19fU1NTf/pK4Q6duy4atUqe3t7V1fX7t27S4vB2rVrDQwMBg4caGlpiUfNzs728PDo3r27ra3thAkThEJhrfArVqywtrYePHiwgYHBtm3b8IVhYWFDhw51c3MbMGCAsbFxQkICvpxOpy9YsMDNza1Lly6jRo2qrq5GCG3evNnf39/d3b13794xMTFPnjyxtrb29vY2NDRcsWIFQujt27e6urovXrxACKWnp+vq6r5+/RohZGZmlpGRgRAaO3bs1KlTHR0de/bsaWdnl5aW1qtXLx8fH2Nj45cvX+KbdnR0dHV17dOnj5WVFb4XuHz58oyMjLFjx1paWh49evTt27fS88qJiYkWFhZ9+vQxNDSU/kUyduzY6dOnOzo69unTx9zc/M2bNw34UjUc1uzmz5+/devWurfn8XhNF0ahNENP8u/FfPlzdNGTBIlE0tTbaiSJoLLkyP8Kty8Q89j1fW/r+k5evnx5yJAhDXvvP7spm3PRIWHtf5qG+K/OBgdjMBifP3/GH9Pp9C9fvmAYNmDAAA8Pj4qKCgzDQkJCFi1ahGFYenq6gYEB3tLY2PjFixf444SEBEtLy6KiIgzDXr58qa+vz+FwZDcRHx9vaGhYWFiIYdju3bulgYODg3ft2oW/y8jICG+Mj+DFZrMxDMP3S54+fYphGJvNtrS0fPToEYZhvXr18vDw4PP5GIa9e/dOX1//3bt3GIZ9/vzZwMAgOzsbwzBra+uRI0eKRCKxWOzt7b1v3z4MwyIjIy0sLPDPKxKJSktLJRKJi4vLoUOH8CVDhgzZsWNHrS7Kzc3Ff5QKCgr09fU/fPiAYdiKFSvU1NTev3+PYdj58+dNTEyqq6sxDFNSUpo2bRqGYQKBwMPDY/v27RiGbdq0SUVFJSMjA8MwoVBoZWV18OBBfIVGRkYxMTEYhoWHh1tbW+fn59va2oaHh+Ob1tfXl/aVvb09j8eTSCSDBw82NjbOycnBMOyvv/4KDg7GG3/8+BF/EBER4eHhgT+2t7e/c+cO/jgzM1NXVxfDsIqKCkNDw3PnzuHvYrFYycnJ+Fbs7Oy4XC6GYb/++uvUqVNrdcWkSZOOHj363S9T48AeIZATDOOc38NPjtKZt02nW8+WP/4via6sPelPRge74h0LxKX5+MKSkpKv/ypXWBKJBN+rQPRv3XPJUEMINcUMyePHj8cLsLe399u3b3/Q8uLFi507d46Pj4+MjMzMzFRWVn7+/Llsg9jY2KCgIPx43fTp0398MFDW7du3tbS03r9/HxkZGRcXZ2VlFR8fj780adIk/CbU6OhoCwuL1NTUyMjIe/fumZqa3rt3D28zdepUKpVKoVC8vLzwj3D58uVp06YZGhoihKhUKovF+vjx4/Pnz9XV1SMjIy9dumRqanr37t1aMTQ1NQ8dOjR//vy//vqLTCZLDyAPHDjQwsICITR8+HCRSJSeno4vnzVrFkKITqdPmTIlNjYWX9ivXz8bGxuEUFZWVmFh4ZQpUxBCenp6o0ePvn79OkJo4sSJ7u7udnZ2zs7OEydO/Lo3Ro8era6uTiKR3N3d3d3dzczMEEIeHh7Sk7UCgWDNmjVz5syJjY2t9V9QC74TOWLECISQqalpYGAgngHfCpPJRHX4f5e7lnjkCrQ+EknZmW3i0gKd0C1kRuuZC4ZEYg4MoWi0K965UGPSX74zFr0rFSBe4eHt6wcPHEB0OOKRyWRTU9MPHz6g3OfI2us/r1WWoZIcMpmM/06UL+lUO0pKSiKR6ActS0pKuFzus2fP8KcjR46sNdI6l8vV0dHBH9NotLqPolBSUiISiaRrdnBw6Nq1K/5YuomSkpKqqippG29vb7w4yX4EGo2GfwQ2m81isWQ3UVpaSqPRUlNT8acqKip+fn61YvTt29fJyWnQoEEMBuPx48fSC2s1NDSkbTQ1NTkcTq3lsgul2+VwOEwmU3pOUUtLCz/4iRDq1KnTsWPHevXq9c3eUFNTk36crz9aRkZGnz59li1b5uLiwuPxjh8//oPxmzgcjuz5Sy0trbKyslqd9tP/d7mDQggaC6sRs4+tw8TCdjPXkGhKRMepN1XPQWSmVsHepWSt9kUjNqEq7rw/BkMhxAUEBOzYsQNdXI5+i0W0/x8JAcPQ+T+QWOjRo4fsdSXNQ0lJSSz+58yura1tZWXl+vXrv9fYysoqMTERf5yTkyOtDT9la2srFApXrlyppPTdr7Stre3169fXrVtXl+Mftra2jx49mjFjhmw2sVg8Y8YMafmspbi4ODU19cGDB2Qyubq6GsaWOC8AACAASURBVD+HisNP6SGEOBzOx48fO3ToIF2Or+358+fSceel8fDDyPn5+fhVSCkpKQ4ODgihx48fb926NSoqavLkyR4eHl9PkfZjd+7c8fHxmTNnDkLo5s2b0uWy/1NSHTp0wP8j8HKYkpISEBBQr801BSiEoFGw6qqSI/+jqGmwQn5vmZfG1IVyV4900+67BPdCy1MSlTuLvvrpVVhLly6NiIgozn6AVndHfWYjXSvEK0LJR1FmPI1G27BhQ/NHsrOz27x5s7e3t7e39+zZsx0dHefPn+/n51dRUREbG7tlyxb88Bpu4sSJGzduXLNmjZOT044dO+o+dW3v3r27dOkyfPjw6dOnk0ik+/fvDx482N3dXbbNiBEjduzYMX78+ODgYJFIFB8fP3PmzI4dO35zhQsWLHByclq+fHnPnj1zcnLwy22WL18eEBCwbNkyTU3N169fa2try863x2KxdHV116xZ4+bmdujQIdlJiD58+LB48eK+ffvu3LlzyJAh+JW3CKF169ZhGFZaWrp3794bN27UyqCrqztx4sQxY8YsXLjw2bNnSUlJe/fu5XA4o0eP3r9/v7+//8qVK0eOHPn48WNl5XoM/9SlS5dVq1ZdvHhRLBZv375dutzOzm7Xrl0fPnxwd3eX/j1haWnp7+8fHBw8Z86cpKSkjIwM/KInYlHCwsKaeZM3btzQ0NCo9ZX6AaFQKHupMWgwefVkUlLSo0ePTExMaDWikgPLaLomrDELSP+9drzVMehk/+v2IxtpDz/fOTZx2lRXZ8cfNG5d38k3b948f/589OjRDXivmpqaj4/PjRs3uJ+y0Mtr6MFJlHIJleQwmcwTJ058fSivXgQCQe/evfGexDCsV69eSkpKQqHQzs5OX18fIVRTU8NisRwdHRFCKioqHh4eCCE/Pz8+n5+Xl9e+fXsLC4tJkya9efMmLi7uy5cvXl5ejo6OsgfllJWVhw0bFh8f//z584ULFxoZGXl5eeEHSG1tbfHLZKhUqpeXlzRS3759aTQaiUQaPXq0QCC4fft2RkaGtbV1nz59lJWVq6urHR0d8ZOOFAplwoQJRUVFcXFx2dnZDg4OPXr0oNFo1dXV3bt3x49S1tTUGBsb29jYqKurjx079tmzZ/Hx8RKJpG/fvhoaGr169bKwsLhz587Tp081NDQGDRqkpfXvNJxkMnngwIF37txJTU2dMGGCq6urvb29vr5+fHy8paWlra3t1atXXV1dN2zYgN+5sXr16rNnz0ZFReXm5m7YsMHNzQ0hJBaLDQ0NbW1t8XUOGjRILBbfvHmTwWAcOnTIyMgoMTGxa9euY8eORQh169ZNIBDQaDQTExMMwzw9PVVVVUUiUadOnfBrPmXXJpFIVFRUunfvbm5ubmZmdvXqVTabvXr1alVV1b59+5JIpD59+giFwry8PCMjIzMzMwaD4enpiRAaNmwYj8e7deuWpqbmkSNH8IMKsluRSCTq6uouLi6y35YrV66YmZnhXwa5I2H/fz1xs1mwYIGJicn8+XUd6bG8vFx67Bg0hlx68pcFS8+mfBLo2pi9vXJnXHdlyy6aQ6fLdxIJogiFwvSbV3TuR+pM/pNu1fUHLVvXd/LKlSvh4eGNmZiXz+efOHHi5s2bnz590tbW9vLymjp1aq075ECzCQsLKysrw4fCkUWn0z9+/Ij/DdH2TJ482cvLa9KkSU2x8tZ6LAsQ5fL1m2ULHjKQaFW7V/lC5NRWqiBCSElJydF/RLWtbenRVdqTl9PbdyE6UUuhpqb2yy+//PLLL0QHAaBJQCEE9cOg0ai8/H2cI6UV1e3cA9tMFZSit+/MGr+YHb5Ge/r/lExgmnvQ4nzvfNY/97qA+oP7CEH9HNu1+cizSdSspOesDt69exMdp0kwOjppjpxbeihMVPCR6CygWe3YsSM7O7uOjVNTU48ePdqkeUDzgEII6qdrQcrAnu5jT9zct20j0VmakHJXD42h00v2L6/hNGpcadCc5s6dGxER0Zg1/P333/joo98kFostLS0rKirwpwUFBdK7AEGrBoUQ1EP5rTPCvHfak5aTqHUdoaP1UnHyVus5pGT/MknVd6cHAo3E4/E+f/7802YikSgnJ6e0tLTW8pqampycHOkd2UVFRd+bMEEkEn348AEfsUy6UCgU5uTkfH2v2/dgGPb+/XvpGvr37y97xcr31sbj8T5+/Ci7XdDSQCEEdVWVmsi/F91ualizTTFPOPU+QfQO9qVHV8HY3HJXXl4+cuRIJyenYcOG2dnZ4QckJ0+ePGbMGIQQhmGBgYG//vorQujmzZvW1tbjx493cXHp37+/dFC3/fv3GxoaBgcHu7i47Nq169SpU7du3Vq3bl23bt3Wrl0ru63U1FRra+sJEyaMGjXK2dkZX7hv3z5jY+PRo0cbGxufO3euVrxJkyZJB7m+ffs2/i78ksWePXt269YtNTU1IiLC3/+fCUz27t0rXVtkZCS+sGPHjqGhoe7u7n5+fo6OjlwuV/79COShrhfLXL58OTo6+tOnT2ZmZvPmzZPelSL16tWrXbt2SZ/OmjXL3t5ebjEB0arfvyq7sEdn1nqKZnOPJEIszWEzS8NXcSJ3MYbMDD9+ks3hTg0ZV/eRulomrEYsePUQNeOtU3SrrmS1/0wMtGLFCiaT+fbtWzKZvHfv3tDQ0JiYmL1793bv3j08PJzL5ebm5p4+fRoh5OLi8vbtWxqNJpFIRowYceTIkTlz5jx8+HDZsmWPHz+2tLRE/z9w19WrV3v37j1z5sxaW9+3b9+sWbMWLVqEEMJHKXv9+vXixYtTU1OtrKwePnzYr1+/nj17/vTGg/Dw8FOnTiUmJuI3z+CzNCCE0tPTly5dmpKSYmVl9eDBAz8/v549e+K3l3A4nLS0NBKJ5Ovre/LkydmzZ8ulP4F81bUQ7tixY8SIEaNGjbp9+7aHh0daWpp0Tg1cXl7erVu3li5dij/9ejYs0HqJS/LZx9axxi+hGX57LKi2jExmjV9avHvx4cnD/8J6ClX1DvgMehwX3YruI/waJqisenmv2XZzSWQytZ1BrUJ45cqViRMn4qOK0Gi0pKQkDMMYDMapU6e8vb1JJNKDBw/we+2ZTGZkZOTDhw+rqqry8/PxUZtjYmKCgoLwKoh+9gunY8eOBw8eZDAY+BRICKG7d+/icwYhhLp3725jY5OcnBwUFNSwD3jnzh3p2tzd3Tt06HDv3r3AwECE0IQJE/Ab/N3d3et+GQ5oZnUthNJh0fv163fjxo3bt29/PUi5rq7u9OnT5RgOtASSSn7JgeXMgSGMjk5EZyEGSYmuOm6J97uJ/axcozS8KkrSU1JS8CnrWimyKpM1fgmxGdhsNpvNfv/+Pf502bJlGIaRSCR8dngtLS3prLxhYWH379//9ddf9fX1IyIiiouL0VdjN//YggULLCwsLl269Oeff/bq1evChQs8Hk92JDYNDY0fzMb+09N7P1ib7HDVdT8ZCZpZve8jrKqq+vTpU63dQVxeXl5ISAiLxQoKCsKH0gGtnkTCPrGe0dlN1c2X6ChEUtU1XPSMs0N5bzZJuyjnsbl5kwxvoVC6du1qY2MjOwg1QkgikYSEhMyePTsrKys0NPTw4cMIoVu3boWFhfXv3x8htG3bNnwHq0uXLsePH6+1zm+O8owQIpFIgYGBgYGBlZWVZmZmaWlpHTt2PHHiBF56BQLBy5cvV65cKfsWFotVVFSEP5ZOckShUKhU6tebsLGxiYiIkK4tLS0Nn/YItBb1LoS//vorPolwreW6urq//PKLlZVVZmbmwIED9+/fj5/0/lpGRsa5c+ekX2JlZeUzZ87IjrBXi3TmEdBIDehJwc1TNcJqeu+R0lm2FdZfGzbu2bTyWOWil7PHslisVtQhVVVVzT+S4k+tW7cuKCiorKzM1tY2Ly8vOzt769atK1euFIlEy5Ytq6qqcnFxOXny5Lhx4xwdHbds2SIWi+/du3fv3j18XNCQkJDdu3eHhIQMGzassLBQW1s7KCjI2dkZn1fWwcFBOnwoQmj58uV6enpWVlbp6ek0Gq19+/Zdu3ZdtWrVpEmTBg8efPz4cQcHB3wgU6kBAwaMGTOmQ4cOPB5PeksGmUx2cHBYuHChq6ur7JwJAQEBq1atmjx5sr+///Hjx52dnbt3794svahAMAwTCAQN+LljMBg/nYeyfoUwLCzswYMH8fHxX0874uzsLL0cq127dps2bfpeIbSwsLC1tQ0ODsafUigUExOTH89j0qrPx7QodenJ1Zu37zt6Qk1F5eySyfoZD3QX7CSrafz0XW1en97efXp7c68eNvuUzVBRaUXfSWVl5RY4T7Knp+e9e/ciIiJiYmIMDQ3Hjh0rEAjU1NROnz5NoVDU1NTOnz+fkJCAYdiWLVv27NkTFRXl4uJy4cKF/Px8hBCDwXj06NHhw4dv3rzJYrG8vb0RQnPnzrWwsMjKyqqsrJTd1sCBA69fv/7q1SsDA4P79+/jw2EnJSUdPHjw7t27/fv3x+eqRQjNmzcPP9XXt2/fI0eO3Lx509TUNCIiIjk5GW9w7dq1q1evlpSUiMViR0dH/CwmjUZLTk7+em0LFy40NTXFH/v4+NRKBeqFRCIxGIwm+rmrx6Db69evP378eHx8PD7y+g/ExsZOmzYtLy/vm6/CoNtEqUtPvnjxwmf6stIpkZ1KU8/lhVn/vkfJpEPzxGsdJJKSQ39hOiY6gTN+3rhlaPyg2wAQrkkH3a7rfYTbtm07evRoXFycbBUUCoW7du3C72bNzs7GTylXVVUdOHAAzhG2Unl5eUJjJyYmOFx2cENmBdXIkuhELQyZzBq3WPQyuepFMtFRAADyUadCiGHYggULPn/+3KVLFxaLxWKxtmzZghCqrKwMDQ0tLCxECG3atElPT8/JycnIyKisrGzr1q1NGxw0jR49emhlxmx7Mf9ukYinby07tRvAkVWZqsGLyiJ3iQpzic4CAJCDOp0jJJFIbDZbdgk+f7Gmpiafz8cnfd6/f39YWFhBQYGhoeFPj52CFktTU/PxpvlFiTHZfYOWDhlKdJwWimxgoTF4MvvvtboLdpJoSkTHAQA0Sl0vlvneVZ2yQ2zo6+u31TkhFYcw923N/ehOi7Z31Yb/yh9RdfOrfveCc2m/1shQorMAABoFDnyBf0mq+Oxj67RGhlKhCtaB1oi51dlplc/uEh0EANAoUAjB/8OwstNblbu6K9t5/LwxQIhEV9aeuIxzab+4+OfzJwAAWiyYoR78g590tYZTygr5g+ggrQnNwJzpG8w+sUFn3lYSpeX+NOXk5Bw8eJDoFAA03Js3b2QHSZCvlvujC5qTqCCXd+OU7q/bWvJv85ZJreeQ6qwXvKijUVUat5IeBvTrNXxYy7rIyN7e3s3N7dmzZ0QH+QaRSPTTUT9AXbT5nuzSpYurq2sTrRx+6wGEiUXsE+s1Bk+h6hgRnaVV0hq9IGvFhIhM+nXb6Vc27amsqh4fPIroUP8yNzc/cOAA0Sm+DYbLkBfoycaAc4QA8WL+prL0Vbv7ER2ktSKrqG3KQf+zxrQs7Lh+y89GxRKdCABQD1AIFV11dlplaoLW6F+JDtK6adrYx/DUNuVtV0455+EMU1ID0JpAIVRokio++9QmrZHzyKrMn7cG37dq2ZISdW2r4kc7O3IX/zqH6DgAgHqAc4SKqLCw8PXr1127diVFH1Du6smwdSE6UavHYDB2b90g+vKhw97fSTw2YsHgSgC0GrBHqHDu3b9v12vg8O0x80f483PeavjDHLNyQzO0UO8TxI7YhFre/H8AgO+BQqhwVm7dVzTmMGXA4gUORts/kWGoTPlS7z0cSTB+0lWigwAA6goKocJhqquSyovXfNl/WtKRp8wiOk6bQyJpjV3Iu3kahpsBoLWAc4QKZ3PY76xp4zqZK+9KqbgVfYHoOG0QVVuf6TuGHbFFd+5mBPNYAdDiwU+pwjHR1gzrbtL1t41pT+/BbCFNRM0rgESh8JOuEB0EAPBzUAgVTtm5naruAzQ7OhAdpE0jkbTG/Ma7dRYOkALQ8kEhVCyVz+6IS76o+44hOkjbR9XWZ/oFl53eCleQAtDCQSFUIFgFl3v5ECv4NxhZu3mo9RiMYYh/L4boIACAH4FCqECqYsJVXPvSjK2IDqIwSCTWmPm82JM1nGKiowAAvgsKoaIQZDyp+fSO6TeW6CCKhaprrOYVUHZuJ9FBAADfBYVQIWDCas75vcpDZpCUGERnUTjMfqNqeOzb+zZ29ejTb9iYrKwsohMBAP4DCqFC4MaE0626Ui3tiA6ikMiU8u7DdNLu5A/fH2c7x3/MZKIDAQD+Awph2yfMfVOVmqgRMJXoIIorvbQiSmiwjH8FWbiWVQokEgnRiQAA/4JC2NZJasrO7dQYOh0mWiKQm5vboQdp3cueeD9YY2VsQIbhZgBoSeAHso0rv3OewtRWcfImOohC09fXvxZ5/B5i7dZ6ef1MONFxAAD/AYWwLaspKyqPv6g5fBbRQQCytbVduH2/ficH9Oga0VkAAP8BhbANunn7Tnv77kadnJ5uWqzuHUjVhgFFWwrNoNn85ChRYS7RQQAA/4JC2NZgGBYya8GHSZfsp2yuLskvMLEnOhH4F0VDm9lvNCdyF4y7BkDLAYWwramurpbQVVVVVFcUhq/iWHz6kk90IvAfal5DsGpB5bO7RAcBAPwDCmFbw2Aw7NobL3kw+2E5LTszxc3NjehE4L/IZM2Rc7lXD0sEFURHAQAgBIWwTYravXa8TgXLqVtq0i0VFRWi44DalEysGZ1cymNPER0EAIAQFMI2CMN4F/boBE4fP32WhoYG0WnAt2kMnlzx9I6o4CPRQQAAUAjbnIrHtxCGqbr5ER0E/AhZTYPpF8y5sJfoIAAAKIRti0RQybt2TDNwJiKRiM4CfkLN019Sya9KTSQ6CACKDgphm8K7fpzR2U3JtCPRQUAdkMlaQbM4Vw5h1VVERwFAoUEhbDtEBbmVz+I1Bk0kOgioKyWLzvQOdrxbZ4gOAoBCg0LYdnAu7GX2HweDa7cuGgFTKx7Gios+ER0EAMVFJToAkI+q1EQJn6PmMYDoIKB+KOpa6n1HsSN3HxcZ5eYX/jJxbKdOnYgOBYBigT3CtgATCblRRzSHz0ZkCtFZQL2p9xyS/TotMTV3V1W33oHjPn/+THQiABQLFMK2oPzmaaX2nelWXYkOAhqETFn5nL1Mt0Cliw/XYVRSUhLRgQBQLFAIWz1xST7/fozG4ClEBwENl1tDT8F0f8k/rfo2ztbWlug4ACgWKIStHufSfnWfkRQNbaKDgIa7eurIreyPk4sv7p4xzM7Ojug4ACgWKIStmyDzmbjok1rPIUQHAY1iaWl58dJ5k2GT+5IKic4CgMKBQtgqVVRU3Lx5MyP9FefCXs3AmSQqjehEQA7UvAPFJZ8F6Y+IDgKAYoHbJ1qfsrIyp56+Ze17zyS/HGKj597JhehEQD5IFKrmsJmc83v0bJxJFPjZBKCZwB5h63PlatSXrqPoA34Lbq8+71Ym0XGAPDFsnGn6pvzEK0QHAUCBQCFsfTSY6krlBUsLjp9W61laAzcOtjUaQ2eUx52t4ZURHQQARQGFsPUJCAgIVvvUu/DW+fNHD21dR3QcIGfUdgaqLn15N04SHQQARQGFsPWhUCirPUzMg0Pfpaf26d2L6DhA/tT9ggVp90VfPhAdBACFAIWw9alKTZQIKjU9BxIdBDQVsrKauu9YzqUDRAcBQCFAIWxl/hlWdNhMRIb/u7ZMzWOghM8RpD8kOggAbR/8Mm1lyuMv0EysYVjRto9M1gycybl0ABOLiI4CQBsHhbA1qSkv48df0hg8ieggoDnQOzhQdU0qkqOIDgJAGweFsDXhRR1VdR9AbWdIdBDQTDSHTufFnZVU8IgOAkBbBoWw1RB9yhZkPlPvO4roIKD5UHWNVZx7866fIDoIAG0ZFMJWg3NpP3NgCJmhQnQQ0KyYfuOqXiR/fHY/IiIiIyOD6DgAtEFQCFuHytQEiaBS1bUf0UFAcyOrqJV19Hi2Y/nkqM9eo2ZcvHKV6EQAtDVQCFsBTCTkRR2FWyYU1sZ7WVrahj1cPUvHHdu09yjRcQBoa+o0wn11dfXly5fj4+PLy8sdHBxmzpyppqb2dbO4uLgTJ05QqdRp06Z1795d3lEVFz/xMs3YCm6ZUFiGerrrvrCWFxxNFA011NMhOg4AbU2d9jDS09P37t3buXPnQYMGxcTE9O/fH8OwWm0SExODgoK8vb2dnJz8/PzS09ObIK0ikvC55XfOawyeTHQQQJjffwut/JhS+On90pwdO9eFER0HgLamTnuEjo6OCQkJ+ON+/frp6Ojk5uaamZnJttm6devChQsnTZqEEHrz5s3u3bv37dsn97gKiBd7QqWbD1XHiOgggDDq6urJN66ICnO9di3W12ISHQeAtqZOe4QkEkn6OD8/n0ajsVisWm0ePXrUs2dP/HHPnj0fPoShoeRAVJhb+TyZ6RdMdBBAPJqeqYqDF+/GKaKDANDW1G8WbKFQOHPmzIULF6qrq8sul0gkxcXF2tra+NN27doVFBR8byXv3r2LioqKivpnvAwGg3Ho0CENDY3vta+oqJCtxIqgqKhowqwFObm54X3MHfyDKiUkxOc3frUK2JNNhKiepPQM5O2YR3byIWsbNP/WmwJ8J+UFevJ7GAwGlfqTSlePQigWi4ODg3V1dVeuXFnrJTKZTKfTq6ur8acCgUBF5bu3uxkYGBgZGQ0fPhx/SqFQ9PT0yN+/HrKmpuYHa2uTfl2+6n7HSe79jdWyVp7/UD7DVz4fXwF7sokQ1pMqKljvIGFcBGvSnwRsvQnAd1JeoCe/5wfFRaquhbCmpmbChAmVlZWXLl36ZnU1NjbOzc11cnJCCH38+NHY2Ph7q1JTUzMxMenXr663xJHJ5Lp8krYkK/sDydNred6fq5T6mWXlyOvjK2BPNhECe1K9d2DhumminNf09l0ICSBf8J2UF+jJxqhTx0kkksmTJxcXF1+8eJFOp0uX5+XlnT9/Hn8cFBR0/PhxhFBNTc2pU6eCgoKaIq6CmDgmaMKtKcJK7pMbpyaOGkZ0HNCCkChU5sAQzsV96KsrtwEADVOnQpiUlHT8+PGMjIzOnTtbWlpaWlo+ffoUIfTs2bPZs2fjbebPn//u3Ts3NzdnZ2e8cDZh6rZu0ezpf3WiFehZ3j59wM3Vleg4oGVRcexFotErU+KJDgJAG1GnQ6Ourq7Z2dmySwwNDRFC/fv3f/78Ob6kXbt2qampz549o1Kpjo6OsJPeGOW3z2nYdpsxfjHRQUCLRCJpDplW+vda5a4eJCX6z9sDAH6oToVQWVm5ffv2Xy9nMBgGBv9evUalUt3c3OQWTVHVcEv5SVf1fttNdBDQcimZd1Iy68hPvAyzkQDQeLDf1uLwrh1TdR9AYekSHQS0aBoBU8rvXqgpLyM6CACtHhTClkX0+b0g44l635FEBwEtHVXbQKWbTzncXw9Ao0EhbFm4Vw8z/caRGapEBwGtANNvbOXzZFHBR6KDANC6QSFsQQSvn4g5Jaru/YkOAloHsooas+9IbhRMzARAo0AhbDEkEu7Vw5oBUxCZQnQU0GqoegWIC3Or3z0nOggArVj9xhoFTafiYSxZXZPRGS67BfVAolA1Bk/OP7NzXgZSUVFet3yRubk50aEAaGVgj7BFwKqreDdOaQRMJToIaH2qzbqmZ31UMnE5y/T3GTqa6DgAtD5QCFuE8juRdGtHJZMORAcBrU9mZubmcvPFkmSGjVcFTaO0tJToRAC0MlAIiVfDLeUnR2sMnEB0ENAq2djYZGY8f4bpTcvYqiLkSGdDAwDUERRC4vGuHVP1GEjRgjvoQUNoampeizj8sKBkluThndOHiI4DQOsDhZBgoi8fBBlP1H1GEB0EtGKOjo5/nzyu13uo5uskorMA0PpAISRMdXV1VVUV9+phdd+xcAc9aDymb3Dl8yRRQS7RQQBoZaAQEuN/G7Ya23Uf3a9P/tvXcAc9kAuyipp6nyBeTDjRQQBoZaAQEoDD4ew+drbstwdzPe1WvizPLywiOhFoI9R6DRXl51S/hfvrAagHKIQEEAgEZBWNkdzbXKp6gkCjoqKC6ESgjfhn/vqrh2D+egDqDgohAfT19XvZGM//sGfX888OOjQrKyuiE4G2A+avB6C+oBAS4/C4Piwb+y0rF968dIZEIhEdB7QhJJLmkGncqKOYSEh0FABaByiEBKjhsflJV83G/ers7Ex0FtAGKZl3UjLtwE+8QnQQAFoHKIQE4F07rtrdj8rSIzoIaLM0Bk8pv3teUsEjOggArQAUwuYmKsgVvH6k3g8GRwZNiKpjpOLkzYs9SXQQAFoBKITNjXv5gLpvMNxBD5oas//4ytREuL8egJ+CQtisBBlPxexCVfcBRAcBbR/cXw9AHUEhbEYSCffqYY2AqSQKzIcMmgPcXw9AXUAhbD4Vj26QVdWVu3QnOghQFCQKlTloItxfD8CPQSFsJlh1FS/2pMaQaUQHAYpFxaFndQ1aFDTIpHO335atIDoOAC0RFMJmUn7nPN3aQcnEmuggQMGQSH8k54w0UysNvXP4YV50TAzRgQBocaAQNocabik/OUpjQAjRQYAiSvxQlKLWdQo7utzcK/NdNtFxAGhx4KqN5sCL+VvVYyCFBXPQAwIMGeC7/emHsx2e3n7yeejvZ4iOA0CLA3uETU70+b0g8ynMQQ+IsnPj6rXTAgqZBtdDh8EI7wB8DQphk+NePcz0Gwd30AOikEikYcOG+vy+WenjK1HBR6LjANDiQCFsKmVlZdHR0dm3Lok5xTAHPSAcWUVN3WcEN+oo0UEAaHGgEDaJ3Nzczu59xoc/yju3/7GaFSJTiE4EAFL1ChAXfRK8SSE6CAAt5WP4+wAAIABJREFUCxTCJhF+8kxBz98Gu7kUt+v8299RRMcBACGESBSqhv9E7tXDcH89ALKgEDaJdixNrfKP84rPrGIMYKrD2UHQUijbe5HpyhVP4ogOAkALAoWwSUyZOOF/1AePPxbwzi87umMj0XEA+JfG0Om8a8cwoYDoIAC0FHAfYZOgCcoDTVW05x/+pZ0+0VkA+A8l04709p3L715k+gUTnQWAFgH2CJsENypcrcdgOlRB0CIxB03iJ16u4bGJDgJAiwCFUP6EHzOrs1/CHfSgxaJq66u6D+DFHCM6CAAtAhRCecMwzuVDGoMmkZQYREcB4LuY/UYLMp8K894SHQQA4kEhlLPKlHhMVK3SrQ/RQQD4ERJdmTlgPOfifriVAgAohPKEiYTcmHDNwJmIRCI6CwA/oermh8Siqpf3iA4CAMGgEMpT+Z3zdLNO9PZdiA4CQB2QSBpDpnGvHsZEQqKjAEAkKIRyU1Nexk+6yhw8meggANQV3cqOZtS+8Pqpy5cvZ2ZmEh0HAGJAIZQbXtRRVff+VJYe0UEAqAeeU/+S2NO/nn7SY/Qvh/8+QXQcAAgAhVA+RJ+yBZnP1H1GEh0EgPoJj75zlmL3S2ed0mmXNu4+SHQcAAgAI8vIB+fSfubAEDJDheggANSPrrbW/9LFsYzHXWs6KmtqEB0HAALAHqEcVL1IkggqVF37ER0EgHqbNnlil/LUIy9y13zccnTHBqLjAEAA2CNsLKxGzI0K1xo5F5HhrwrQ+tDp9LvRF7CamqItc5gSHtFxACAA/O5uLH78RZqhOd3akeggADQciULRHDaDe+UgJhYRnQWA5gaFsIHevn1r5+nTxb7bl6hjTH+4ZQK0evQODlQ9M37SVaKDANDcoBA20JjpoWm+m8cPH3eWzbya9JjoOADIgeawGeVxZ2vKy4gOAkCzgkLYQCWl7A7a6v14j7fSe7/PzSM6DgByQG1noOLStzz2FNFBAGhWUAgbaHRgwNrU+Tt4JoykfUFDA4iOA4B8MH2Dq17eE+XnEB0EgOYDhbCBwgK9HQ2Y7k4Wj29esbCwIDoOAPJBVlFT9x3Dubif6CAANB8ohA2BiYTcq4eNJy6ZOXOmqakp0XEAkCc1z0GScrYg/SHRQQBoJlAIG6L89jkl0450KzuigwDQBMgUzcBZnEsHYFYKoCCgENZbDaeEn3RVY/AUooMA0FTo1g40Qwt+/EWigwDQHOo6sgyXy33y5MnLly8tLS2HDBnydYP3799HRkZKnw4fPtzKyko+GVsYzpWDaj2HUFi6RAcBoAlpDptZuHmOinMf+KqDNq+ue4SrVq1aunTpsWPHIiIivtngzZs3O3fuLPt/IlHbHJ+i+n26MCdTvXcQ0UEAaFoULV01rwBu9BGigwDQ5Oq6R7h582aE0MqVK9PT07/XxsTEZP369fLJ1TJJJJwLezSHzSAp0YmOAkCTU/cZWbh+enXWSzgdDto2eQ66XVRUtGLFChaLFRAQ0CbvKODfiyarMJXtPIkOAkBzINGUNAKmlEXufuEwtAZDvr6+VCoM0w/aILl9rdXV1b29vRkMRkpKyvLlyy9dutS3b99vtnz//n1cXFxycjL+lMFgbNmyhclkfm/NVVVVFApFXjkbDKvi825EqE8Jq6ysJDpLA7WQnmwDFKgnOzg/y1oTnX4ugqttvXFHfMxFEokkx9UrUE82MejJ71FSUvrpH3ByK4Q9evTo0aMH/tja2vqvv/76XiHU1tZmMpmDBg3Cn1IoFC0trR8EFQqFdDrxhyJ50UeUnXqpmHYgOkjDtZCebAMUpyc/ffq06gN1v6PgnNcfuSdnfvr0Sb4XwSlOTzY16MnvIddhgrwmOdDh4uKyb9++772qoaFhYmIyatSoOq6NQqEQ/peOKD9HkHZfb+kBMtFJGqMl9GTboDg9qaWllZ1fcMU7ZGHB8a3sPG1tbfl+cMXpyaYGPdkYjbqPUCwWR0dHV1RUIITKy8vxhRiGXbhwwc6uDZ1dxzDO+T3MAePJqt89fgtAm8RkMv/325wTZw8PKLi2eerwdu3aEZ0IAPmr6x7hxYsX165dm5+fX1VV1a1bt1GjRi1atIjP5w8ePDgjI8PGxmbmzJlv3741NTV98+aNQCCIjo5u0tzNqTLlLiYUqLoPIDoIAASYNXXirKkTKx7dMHx4E2EYkus5QgBagroWQi8vrwMHDkif6urqIoSYTGZaWlr79u0RQkeOHElNTS0oKDA0NHR0dFRSUmqKuM3s48ePSkiCoo5qT1wGP/9Akam6+lbcv1757K5Ktz5EZwFAzupaCHV0dHR0dGotJJPJXbp0wR8zGAx3d3d5RiNa4Lgp994V/GZSbWdh1N+8E9FxACAUiaQZOLP06CpGl+5khgrRaQCQJxhr9Ntev36d/JHHmrDD17xdaNxbPp9PdCIACKZkZsOwdeVdP050EADkDArhd5EQWv1l/0a9cRyhhOgsALQIGoMnV6UmiD5lEx0EAHmCQvhttra2v3WkqBW+uhOxZULgIDU1NaITAUA8soo6c+DEsnM7EIYRnQUAuYEBk75NUskfb0AS+i9NXdXRwMCA6DgAtBSqbr6VT+MqHsbCddSgzYA9wm/jxoQrO3iZuXlDFQTgP0gkzeGzedeOS/hcoqMAIB9QCL9BmPdOkPaAOWAC0UEAaIloBubKzr250eFEBwFAPqAQfgXDOOf3aAyeTFaB84IAfJvGgPGCzGfCD6+JDgKAHEAhrK3i/jUShaLSzYfoIAC0XCS6ssaQqWXndmI1YqKzANBYUAj/Q8Lncq+f0BwxF8aRAeDHVBx7UZisiqSrRAcBoLGgEP4H5+phVRcfmoE50UEAaAU0g2bz4s5yP398/fq1WAy7hqC1gkL4r+qstOqsF8z+44gOAkDrQNUxKjXuemFBSO85a62dPIqKiohOBEBDQCH8ByYSlp3drhU0l0RXJjoLAK3G9HPJNsYmHQNCP3afu3P/YaLjANAQUAj/UX7zNM3YkmHrQnQQAFoToQQt1xy59ss+ZUwEZ9ZBKwWFECGERAUf+fdjNIfOIDoIAK3MtlXLMyJWvC0qWV5yfN7MqUTHAaAhYIg1hDCs7OxOjUETKRraREcBoJXp6dUj62nS53cZauc3aQp5CNWerA2Alk+hC2FJScmTJ0+61JTQsRoYOBGAhlFVVbV26FZRNbHs7E7d+dvhAClodRT30Ojr16+7ePadcyS+7NrJVybd4acXgMZQ7d6fRFfm348hOggA9aa4hXDr/vBC//WLbdVOaA/848BpouMA0MqRSFoj5/Kun6zhlhIdBYD6UdxCqK3J9K1I6ST4sEtsr6XBJDoOAK0eVcdIrYc/58JeooMAUD+KWwiXzp66jpyw4WGW0fWlu9eHER0HgLZAve8oUWFe1ct7RAcBoB4U92IZ0v2rFl5+F7bOYjAYRGcBoI0gUWms4AWlR1fRrexh/hbQWijoHqEw901VaoJGwFSoggDIl5KZjbJ9D+6VQ0QHAaCuFLEQYmJR2eltGsNmkFXh1CAA8qfhP6k666Ug8xnRQQCoE0UshLzYk1RtfRXHXkQHAaBtIikxtEbN45zbiVVXEZ0FgJ9TuEIo+pRd+fim5sh5RAcBoC2jWzsoWXbhXT9BdBAAfk6xCiFWI2af2qQxZBqFqUV0FgDaOM2hMypT4oU5mUQHAeAnFKsQlt+MoLD0VJz7EB0EgLaPrMrUCJrzausik06O5l1dHz1+QnQiAL5NgQqh6PN7/r0YrRFziQ4CgKKIyyl5VskIGfvLxwnnQmb/RnQcAL5NYQqhpKbszDbNIdMomu2IjgKAoigqLlku6h7IiXehlVVU/l979x3QxNmAAfyyCISZhK0ouEBF0QoO1LaKrIo4WidOXKh1lIKjVala68JVZ8VRW60VRRTc1gVKqyIIrYqKKAhhCCGEsDK/P2j5LHUgRN+M5/fXXXIeD6/Ak0vu3sOJM6Ch9KUIxecP0U0tOB4DSAcB0COBgwIMbh2MKLTZkPHV1NGfko4D8HJ6MbOMTPBEkhhnE76NdBAA/WJtbZ1+/eLFixdtcuhz2+JWhaChdLkIlUrl0aMxeQLBiOr7FoMmMyzwewjwvnG53M8++0wp8S5cG2LUubdBSxfSiQDq0+W3RsdMmTV5f9LTjJzf/3ogcXQjHQdAf9FNzC0+nSk8GKmS1pDOAlCfLhfh1d9vOvjPDDbJnU/1v3T5Muk4AHrNyK2vgUO7spN7SQcBqE+Xi5Bvytn8dNUK64llWXccHR1JxwHQdxbDZ1X9+Xv1/WTSQQD+RZeLMHaav7Ao9+ZPS2cP6tWrVy/ScQD0Hd3QmDf6i9LDm5SVEtJZAP5PZ0+WqXn8p0nu/eE7jo7ELSYANAa7XVcj116i2J28oDDSWQD+pptHhMoqifDgOu7IObjREoCmMQ+cIn16vyotkXQQgL/pZhGKjmw1cu1l2N6DdBAAqI9mwOYFhYlidijEpaSzAFCUThZh5e1LMkGW+aDJpIMAwMsZOLY39vyk9JdISqUinQVA54pQISoWxe7iBoXRWAakswDAK5n5jFZWV0pvnCUdBEDHilClEh6MNO33qYFDO9JRAOC16AzeuPnlFw759Onr0NF9844o0oFAf+lUEZZfjqEUMtN+mNsXQAsweLbfJucv7NXu+awL3+46lJWVRToR6CndKUJZQbbkcgxv7AKKrjvfFIAOq6qqiitmZnBaLSg5rGzeOTc3l3Qi0FM6ch2hSiYV7l9lHjiFwbMmnQUAGoTD4Tg346/8SxzbLPlheZGHxyrSiUBP6UgRimK2sWxbcjy8SAcBgLdw4uC+8xcuZBdwI83uGEgrKSMj0olAH2n3u4ixcfFdPvRZPHa4JCPVYuRc0nEA4O2wWKxRI0cOn7vQtE+A8Oc1lFJJOhHoIy0uwqysrGmL15YFrP7MUjnjSjbdkEM6EQA0kplvEEWjl186QjoI6CMtLsJ79+6pnPvtKN0daT/x1vNqhUJBOhEANBaNxhsbLkk4Ln1yl3QU0DtaXIQeHh6LlVefVtGj7+U52VkxGAzSiQCg8RjmfO6o0JKf1igry0lnAf2ixUVoJngwrKNDusJgQRvRuZiDpOMAQFMZdvAw6ty79NBGTL0G75O2njUqL84vjdluF/Lt9uZtSGcBALUxD5z8fOuC8quxph8PI50F9IVWHhGqFHLhgTVmfkEstCCAbqExmPwJCyWXjtRk4cNCeE+0sgjLTkQxzC1N+gwiHQQA1I9hYcUdFSr8eY1SUkY6C+gF7SvCqvSk6rs3uaPmkQ4CAO+KYQcPjnt/4YG1+LAQ3gMtK0J5sUB0ZAtvwiK6kQnpLADwDpn7j6uplHz3mbd1e/cP+noLhULSiUBnaVMRqqQ1JXtXmPmPM2iBuywB6Do6Y9n9Gm8Hi47TtqZ3nBixaj3pQKCzGnrWqFgsTklJefTokaenZ8eOHV+6zb1796Kjo5lM5tixYx0dHdWW8R+l0ZtZzVobe36i9j0DgAbKfF42035SVN7mwbygfEEJ6Tigsxp6RNi/f/958+Z9/fXXV65ceekGaWlpvXr1UiqVIpGoW7du2dnZastIURRFSRLjZPnZ3BFz1LtbANBYiz6f8vToim3F3H2Fm+ZPG0c6Duishh4R/vHHH0wm08vrlbd3iIyMDAkJWb58OUVRxcXF27ZtW7t2bdPzlZSUJCcnd7M1l104ZDV3I41l0PR9AoBW8B7gdfNkq7S0tHbFdpys65Rnb9KJQDc19IiQyXxDZV6+fNnX17d22dfX9/Lly03KRVEURd24ebNDb+95208827m0tOcwJt+26fsEAC3i5OQ0ZMiQZhMWyApyJAknSMcB3aSemWWUSmVRUZGVlVXtqo2NjUAgeNXGOTk5SUlJd+/+fbWskZFRRESEiclLzgJduvb7opG7dkt/jpJ75x86f6g/LhxskurqahaLRTqFLsBIqkvDR9IkaH7Z9vkqvj2rdad3nUob4WfyVVgs1htnolZPEdJoNAaDofznXmJyufw1/yXGxsbNmzfv1q1b7SqDwTAyMnppUDMTY1pF6ffNRlx6WjzWJBfTajcRg8HAGKoFRlJdGj6SDEs78zHh4kOR3FlrGTy8OVQffiZfhUajvXEbtRWhra2tQCBwc3OjKEogENjb279qYz6f36VLlxkzZrxxt+u+WZQ8aHg628JJJv7u9DG83mkiFouFMVQLjKS6vNVIstp3U/kEifausP5iI93Q+J0G0zr4mWyKJl1HKBQKb926VbscEBAQExNTuxwTEzNw4MCmRqMoR0fHJ3/eunVsT2bajWbNmjV9hwCg1Uz6BBi26Szcvwr3sgc1amgRrlq1ytvb+86dO1u3bvX29k5MTKQoKiEhISAgoHaDL7/88vTp06NGjQoMDMzIyAgJCVFXRC6Xq65dAYC2Mx8WopLLyk7vJx0EdEdD3xodNGiQh4dH3aqzszNFUX379o2Li6t9xNHR8e7du+fOnWMymX5+fi89+QUAoIloDCZ/4tdFG+ay7Bw53fqRjgO6oKFF6Orq6urqWu9BPp/P5/PrVrlc7qhRo9QWDQDgZejGZvwpEc+3LVSacBMeC6ysrF58mQ7wtrT1xrwAoM9Ydo6cz2Y/2bxgYVF7QaFgRPfW29evIh0KtJU2TboNAFAnSVC2odxxRyuRauyOY6fPq3DDJmgsFCEAaCULC4u43MozZp5RT5Yb0ht0uRjAS6EIAUAreXp6Du5oFXVgpyjrTvy43riFLzQaihAAtNXerRsEd29NOnjBhk0Xnz1AOg5oKxQhAGgxBoNBYxlYTltemXoVs3JD46AIAUDr0Y3NLKd/W34xuiotkXQW0D4oQgDQBUy+rWXISlHM9pqsv0hnAS2DIgQAHcGyc+SNW1iy91tZXhbpLKBNUIQAoDvYbd0shoUUR0UohEWks4DWQBECgE7hfPCxqdfw5zsWKcSlpLOAdkARAoCuMekbyPEYkLcp9GMvP6fOPTZu+4F0ItBoKEIA0EFmPqOj03PCujoUh5xauT8uNTWVdCLQXChCANBN6++V/Wnm+lPud6rWno8ePSIdBzQXihAAdJN7544rU4tzxdXbaL991LcP6TiguVCEAKCbjv28e/PgNiW2Dp6evRgnf1DJZaQTgYbC/QgBQDexWKzJEydQFEUplcKDkSV7lvGDl9JYBqRzgcbBESEA6Do6nRcURjc2K9mzTCWTkk4DGgdFCAB6gE7njQmjm5ijC+G/UIQAoB9qu9DUAl0I9aAIAUBv0Om80V/STbklu79BF0IdFCEA6BM6nTc6lG5qURwVcf706ZlhX0UfjVHh7vb6DUUIAHqGTueNCRNIqqUxW35RdJ+65fgmzMGm31CEAKB/6PRVD6S/83rFsk8ZDwz/JTaedCAgCUUIAPrI3a3DpodVB829j+V+4+PWjnQcIAlFCAD6aN7M6ZPa0q5Eb0ySm87i5EuzM0gnAmIwswwA6CMmk7k18rva5eq7fxRHRfDGzjd06UY2FRCBI0IA0HeGHXvyJy0WHlhXdSeRdBYgAEUIAECxW3eymrFSFLtTcu0k6SzwvqEIAQAoiqJYzVpbzVlfcS1edHQbpVSQjgPvD4oQAOBvTL6t9RebFKLnxT8sqSwtjo2NPXfunFKpJJ0L3i0UIQDA/9HYRvzgpQz7Vn8uCoo4lDRyU/wnw8eSDgXvFooQAODf6PRMe7fdpZY/tch19ZuYmiUQi8WkM8E7hMsnAADqs7S0PP6oIM0zcuezNfssazgcDulE8A7hiBAAoD4nJ6e5owdlb5sYcjlzjodj+bHtOH1Gh6EIAQBeYnH4vMKHaTdvJ7deulshLCyOilBWV5AOBe8EihAA4HXohsaW01YwrZsXRc6W5WaSjgPqhyIEAHgTOt1iaIh54JTiXUskV4+TTgNqhiIEAGgQo86e1vM2V6ZeKdmzXFkpyc3NTUpKqqmpIZ0LmgpFCADQUAyetdXn6xhcq8ffTJgwMmjQigMdun9YWlpKOhc0CYoQAOAt0Jgsi2Ezlt8ujOxlN/Rjr9xOo389cpR0KGgSFCEAwFv7o1QVaLsksCzhB5ObZkz8IdVu+P8DAHhr29Yul0WNmXnhrrSitH9uojTnAelE0HiYWQYA4K35DPDKvXtbLBbzeLyq9KSS3d8Y9/Az9QuiMfBHVfvgiBAAoDGYTCaPx6Moyqizp838HbKiZ0XrZkmfPSSdC94aihAAoKnoJhb8SYvN/MaWREWUxe9VyWUURQkEAplMRjoavBmKEABAPYy69LWZv0P+PC9/7cwhXh93HTq1Rafu15N+J50L3gBFCACgNnQTC37wkgQFP8LVfPKgoSXBRz5fGEE6FLwBihAAQM0eMHj+Vf5ta3JOFq9tzZaTjgNvgCIEAFCz6cHj2dej5icXHbh9b90H3NJDG5QSEelQ8Eo40xcAQM1sbW0f3L6enp7u5ORkaW5afuloweoQM+9RJn0DKToOPzQOihAAQP0MDQ27d+9eu2zmN5bTrZ8oZnvFzfMWn85it+pINhvUg9cmAADvHNOqmWXISvNPxgsPrBEeXHf810NObj0dO/f49Wgs6WiAIgQAeF8MO/a0XfCDgm3cOmGfz/BZBVPj5yxZWVZWRjqXvkMRAgC8PzS2kaSb/9RHJh7Sx1dzvhzSoUVhfj7pUPoORQgA8F45OjpWSISfPzScU9gi0ExiFru+6k4ipVKRzqW/cLIMAMB7RaPRblw8fejXwxRl02/kV8y8h2Xxe8svHTEPCGa360I6nT56iyJ8/PhxXFyckZHRiBEjaqeafVFBQUFiYmLdau/eve3t7dWTEQBAt3A4nMnBk/5eadfVOvT7qrRrpUe2MLjW5oMmsZq3jT4ac/POnyMGD+zxz6mn8O40tAhv377t5eUVHBxcWFi4evXqlJSUel2YmpoaEhLi5eVVu+rk5IQiBABoEBrNqEtfw069Km+cL9mzPKtCtfMh/UqLT3+c/OXJqHW9evYknU/HNbQIV69ePWfOnOXLl1MU5e/vv2fPnvDw8HrbtG3bNjo6Ws0BAQD0A43BNPb8hNPDZ+WgARGdWs+n/7Z3wJBDx06iCN+1hp4sc/78+YCAgNrlgICACxcu/HebsrKyH3/8MS4uDmcDAwA0Do3BfMZt5S3x32g1cmLVtZmMzMpbv6kUmLD0HWrQEWFFRYVYLLa1ta1dtbOzEwgE9bZhMBi2trZJSUkZGRlTp049e/Zs165dX7o3gUCQmppatwdDQ8PQ0FAOh/Oqr15TU2NgYNCQnPB6GEl1wUiqC0bypXas+zZ4TlhGQuZZf++V4z+VJMSKTu/n9Bls2N2bZmD40n+CkXwVFotFf9O0dg0qQhqNRlGU6p+ze1UqVe0jL/Lx8fHx8aldDg0NDQ8P/+233166NwaDweFwLCws6lb/uzcAAL3F5/NPHNxXt8pu1VH27FHV1WMVl6KN3L0Me/iGr//h+JnzfK7Fz9s3uLi4EIyqGxpUhBwOx9zcvKCgoGXLlhRF5efn29nZvWZ7f3//w4cPv+pZGxsbd3f3L774ooERpVIpm81u4MbwGhhJdcFIqgtGsoHYbVxN2rgqykoqky8Wfh/mK5bmTI48K7EYNys0/dpFCiPZNA39jNDX1zc+Pr52OT4+3tfXl6IopVJ59+5dqVRau1y38dWrV9u0aaPuqAAAeo1hzjf1GnG2pdduqssY4bkb5WtHW8oVwiLSubReQ88aXbRoUb9+/crKygoKCh4/fhwcHExRlFgsdnV1vX//vouLy/Tp08vKylq2bJmRkXH9+vVTp069y9gAAHpq4MCByzYFnuo1060oP6J9y8L1n7NbudK7fER17UPRGaTTaSWaqsHz+uTk5MTHxxsbGw8dOtTc3JyiKLlcfubMmX79+pmYmOTk5CQmJhYWFtrb23t7e/P5/FftJzQ01MHBoeFvjZaXl5uamjZwY3gNjKS6YCTVBSPZOM+ePTt2PK5dm1Z+fn6UrKYy5ao46bRKWMDp8iGn28fHkh+ELl4uVyqnTQj6dvEC0mG1wFsUobqgCEnBSKoLRlJdMJLqUl5ezlFUV6Vdl9w8n5v56LTdwKN8n/yfvrx+cHO7du1Ip9N0mHQbAEAXMCysTD4aYjYrMvTPahmLsydn5a/uBjXX4uTCQtLRNB2KEABAd3A4HG7L1tuuJg/IdNydKbbnMIs2zC1cE1J2cp/0aQalUkml0vAly7sPCFi2OvLFkxz1Ge4+AQCgU47/su/y5cvl5eV+fkvZbDalUklzM6vv3hDFbJMLCx9KDQQFZvd9N9y7uMXCLGruzOmk85KHIgQA0Ck0Gq1///4vrhs4tDVwaGvmN1YuLNwza0pAC6slBQuTOzUvSrksK/Bj2bSg/j2riVQq1at5avDWKACAvmDybBjuA6bdkXoYh8U/LOjtaFmye5lgyeiSH1dKrp2UFeRQFBU09fPmXXrbtnOLjdeXq+BwRAgAoEcWfjHblLP7wrU9fp8G9hg/lqIohai45lFaTWa65HJMjUTcX65gTFz8B7PV3IWThg4aSDrv+4AiBADQI3Q6ffaMabNnTKt7hGFhyfHw4nh4URR1PubXK8evdK16NK3iuOXH/OfbFxq0cDZo4WzQ0plhzqco6lby7fEzQ8srKiaPHbVsURixb0OtUIQAAPC3Pv6Bs1d+H19mwiwWTf3QPaLfUGnOw4ob50RHtlB0ukGLdufjr7QeODuD32XHL+EBAz7y8PCotwehUJient6+fXsbGxsi30IjoAgBAOBvHA4n/ffLV65csbUd5+bmRlGUYXv32qcUwiLpsweyExeHKP5yyTll41GjPLFZ+LALy86RZdOcYWnP5Nvde/DQa1iQ1Lk/81FY9M4NH3/04X+/xPPnzw0MDGqnJ9MQKEIAAPg/Nptde1uFehg8ayOe9V2uy/rboipLD8fkPb8fXsauEMrzsyWZafLifIWTz8mKAAAHF0lEQVSwkJJR3w/snm3e/Gnniad2bPBsZccw5zNMuXVnpY4KnnEp5T4lrZwbHPR12Nx6X0KlUp05c+ZBZtawwYNqb3ZUT3V19YMHDxwcHHg8nhq/ZRQhAAA01OF9O2OPHxfkF45Ye9qq3pufSuX2BQuuVdm3MGrTsvx2T66q9PAmpViorCinm5gzLKwqKIZ7dZ71yCmFDF5a/ArJwL5sCx7dyITOMam94fDchUv2JwvKW/Ras33IH2djHR0dX9x9YWFhzwEBFbadaLnpP276zt/XR13fFIoQAAAaik6nfzps2KuemxseHuM/NOO6AUcumXT8V5vWrSmKUinkynKRQlT8PC05Oymbo5T2rLlr3ZxTce5ARU2lsrJcWSWhFAo6x3REabF/d+ca2r3qAZ2e71pq2s6ZxmDR2IY0OoPGNrr9+40xH/aV2XfaYhQRvmwSihAAADSOtbV1xu3rYrHYzMys7kEag8mwsGRYWHZo6Xx/b1z6+fNUZemo/p5B81bUbaOSy5SVkmmfjci28Te0bMa7v/4r9x6cLm4qhUxVU61SyFXSajnLQKyg1TBNFAo5k6HOG06hCAEAQJ1ebMEX0Wi0305E//XXX8bGxk5OTv96islimHEjt2wdOeVzQUlpj4G+H078vN4/7/uB9/wBAaI/H/OLV27du12NgVGEAADw/ri6ur7qKWdn5zuJF171LJfLvXcrMS8vz9rams1mqzESihAAALQDnU53cHBQ/27Vvkf1UqlU69atI51CFygUivXr15NOoQtqamo2b95MOoUukEgk27er8w0uvVVaWrpr1y7SKbSYphehQqFYtWoV6RS6QCKRbNy4kXQKXVBSUoI/32qRn58fFRVFOoUuyM7O/umnn0in0GKaXoQAAADvFIoQAAD0GooQAAD0Gk2lUr3nL+nv75+cnPyqC03+6+nTp/Um2oFGUKlUOTk5L52+D96KQqHIy8tr0aIF6SBaTy6XFxQUNG/enHQQrSeTyYqKipo1a0Y6iCYaM2bMihUrXr8NgSIsLy8vKChgNHhegJqaGvVeMqK3MJLqgpFUF4ykumAkX8XOzs7IyOj12xAoQgAAAM2BzwgBAECvoQgBAECvoQgBAECvoQgBAECvafSk23K5fN++fffv33d1dZ0wYULDTzSFejIzM1NSUkQiUVBQkLGxMek42koul1+5cuXatWtVVVU9e/YcMmQIjUYjHUoricXiY8eOZWRkqFSq7t27Dx06lE7Hi/ImycjISEhIGDx4sE29u8ZDA2j0D9+kSZP279/ftm3bXbt2zZgxg3QcbfXkyZMePXrs3Llz+vTpIpGIdBwtduHChbCwMIVCYW1tHR4ePn36dNKJtNXTp08vXbpkbW1taWm5cOHCWbNmkU6k3aRS6bhx42bPnv348WPSWbSS5l4+8fTpUxcXl9zcXEtLy/z8fCcnp6ysLHt7e9K5tI9SqaTT6eXl5WZmZrm5ubjqttGqq6sNDQ1rl1NSUnr06CEWi994iRK8XkJCQmBgIF6iNcXSpUvZbPb69etPnjzp6elJOo720dwjwsTERDc3N0tLS4qi7OzsXFxckpKSSIfSSnjTSV3qWpCiqOrqajabbWBgQDCPDlCpVNevX+/UqRPpIFosPT391KlT8+fPJx1Ei2nuZ4QFBQVWVlZ1qzY2NgKBgGAegDoymSw0NDQsLAyfWzdFmzZtiouLeTzepUuXSGfRVnK5fMqUKVu2bGGxWKSzaDHNPVZgMpkKhaJuVSaT4dU3aAKFQjF+/HgbG5uvv/6adBbtduvWrdu3b/v4+AwfPvzFX3ZouDVr1vTu3RtvhzaR5h4RNmvWLC8vr241Ly8PHxACcUqlcuLEiaWlpSdOnMBr8CbicrlcLnfTpk3GxsaZmZnOzs6kE2mfX375haIod3d3iqLEYnFwcPBXX301fvx40rm0jOYW4YABAyZNmpSRkeHi4pKWlpafn9+vXz/SoUCvqVSqGTNmZGdnnzlzBhMcN0VlZSWHw6ldTklJodPpeJnbONHR0dXV1bXL/fv3X7BggY+PD9lI2khzi5DH4y1evNjb29vHx+fs2bPLli0zNTUlHUpb9enTp6KigqKoTz75hMVi3bhxAx9uNcLx48d37drVoUOHjz76qPaR2NhYBwcHsqm00YYNG44dO+bq6lpaWpqQkLBp0yb8djdOx44d65YZDIazs7OtrS3BPFpKcy+fqHXnzp3aC+pxXllTpKamKpXKutUPPvgAV4I3glAofPLkyYuPuLq64tCwERQKRWpq6qNHj0xNTd3d3fG3Wy3u3LnTpk0bExMT0kG0j6YXIQAAwDuluWeNAgAAvAcoQgAA0GsoQgAA0GsoQgAA0GsoQgAA0GsoQgAA0GsoQgAA0GsoQgAA0GsoQgAA0GsoQgAA0GsoQgAA0Gv/A+DsT7iC5cPvAAAAAElFTkSuQmCC", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", " adr_1d(n, p, q, f; a=0,b=1, α=0,β=0)\n", "\n", "Approximate the solution to \n", " -u′′ + p u′ + q u = 0 on (a,b)\n", " u(a) = α, u(b) = β\n", "(linear Advection-Diffusion-Reaction equation with Dirichlet boundary conditions)\n", "using central finite differences\n", "\"\"\"\n", "function adr_1d(n, p, q, f; a=0.,b=1., α=0.,β=0.)\n", " # mesh size\n", " h = (b-a)/n \n", " # interior grid points \n", " x = [ a + j *(b-a)/n for j ∈ 1:n-1 ] \n", "\n", " # Finite Difference Coefficients\n", " # (using a second order central difference) \n", " d = (2/h^2) * ones(n-1) + q.(x)\n", " d₋ = (-1/h^2)*ones(n-2) - p.(x[2:end])/(2h)\n", " d₊ = (-1/h^2)*ones(n-2) + p.(x[1:end-1])/(2h)\n", " Lh = Tridiagonal(d₋, d, d₊)\n", " F = f.(x)\n", " F[1] += (1/h^2 + p(x[1])/(2h))*α\n", " F[end] += (1/h^2 - p(x[end])/(2h))*β\n", " U = Lh \\ F\n", "\n", " # add the boundary conditions\n", " U = [α, U..., β]\n", " x = [a, x..., b]\n", "\n", " return x, U\n", "end\n", "\n", "p(x) = cos(x)\n", "q(x) = -sin(x)\n", "f(x) = 0. \n", "\n", "x, u = adr_1d(50, p, q, f; a=0.,b=3π/2, α=1.,β=exp(-1))\n", "P1 = scatter(x, u, title=L\"$-u′′ + \\cos(x) u′ - \\sin(x) u = 0$\", \n", " label=\"finite difference approximation\", markersize=2)\n", "u_e(x) = exp(sin(x))\n", "plot!( u_e, label=\"exact solution\" )" ] }, { "cell_type": "markdown", "id": "6812ca8a", "metadata": {}, "source": [ ":::" ] }, { "cell_type": "markdown", "id": "1ed14937", "metadata": {}, "source": [ "* Exercises from the end of lecture 15" ] }, { "cell_type": "markdown", "id": "f0d03bbf", "metadata": {}, "source": [ "::: {#exr-convergence}\n", "\n", "*(i)* Show that $\\hat{U}_0 = \\frac{1}{\\sqrt{n}} \\sum_j U_j$. \n", "*(ii)* Write $-L_{\\rm per} U = \\bm f$ in terms of $\\hat{U}_k$ and $\\hat{F}_k$, \n", "*(iii)* Hence show that, if $F$ has zero mean $\\sum_{j=0}^{n-1} f(x_j) = 0$, then there exists $U$ solving $-L_{\\rm per} U = F$.\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "8a2a3cb1", "metadata": {}, "source": [ "::: {#exr-}\n", "\n", "Read the introduction to the book \n", "\n", "\"Brigham, E. O. (1988). The Fast Fourier Transform and Its Applications.\"\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "eba77266", "metadata": {}, "source": [ "* Trefethen (2000). Spectral Methods in MATLAB\n", "* Brigham, E. O. (1988). The Fast Fourier Transform and Its Applications." ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.11", "language": "julia", "name": "julia-1.11" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.11.6" } }, "nbformat": 4, "nbformat_minor": 5 }