Lectures/Lecture09_Solving_BVPs.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyOgehy2fMH9HiUCNWmuAHH1",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/hamidrezanorouzi/numericalMethods/blob/main/Lectures/Lecture09_Solving_BVPs.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# **Solving Boundary Value Problems**"
      ],
      "metadata": {
        "id": "jp9WfBVmiNm3"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### **Notes**❕\n",
        "This content forms a part of the instructional presentations for the **`numerical methods in chemical engineering`** course designed for undergraduate chemical engineering students at Amirkabir University of Technology.\n",
        "\n",
        "Feel free to utilize the information and source codes provided in this material, ensuring appropriate acknowledgment of the original document.\n",
        "\n",
        "The visual elements featured in this document have been sourced from the following sources, unless specified otherwise:\n",
        "\n",
        "* Steven C. Chapra, Applied Numerical Methods with Matlab for Engineers and Scientists, 3rd edition, McGraw-Hill (2012).\n",
        "\n",
        "\n",
        "<div align=\"center\">\n",
        "🟪 🟪 🟪\n",
        "</dive>\n",
        "\n",
        "---"
      ],
      "metadata": {
        "id": "HKYl2J2LijoI"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 🔵 1) What is a boundary value problem (ODE)?\n",
        "* In general form, it is:\n",
        "$$\n",
        "\\frac{d^2y}{dx^2} = f(x,y,\\frac{dy}{dx})\n",
        "\\tag{1-1}\n",
        "$$\n",
        "* where $x$ is independent variable $a\\le x\\le b$ and $y$ is dependent variable.\n",
        "* The counditions for solving this ODE are defined at $a$ and $b$. It can be a combination of the following:\n",
        " * **Dirichlet** boundary condition (constant value):\n",
        "  $$ y(a) = Y_a \\ \\ \\text{or} \\ \\ \\ y(b) = Y_b \\tag{1-2}$$\n",
        "   in a heat transfer problem, it is called constant temperature condition.\n",
        "\n",
        " * **Neumann** boundary condition (constant flux):\n",
        " $$\n",
        " \\frac{dy}{dx}|_{x=a} = D_a  \\ \\ \\text{or} \\ \\ \\ \\frac{dy}{dx}|_{x=b} = D_b \\tag{1-3}\n",
        " $$\n",
        "  in a heat transfer problem, it can be used for isolation or constant heat flux condition.\n",
        " * **Mixed** boundary condition:\n",
        " $$\n",
        " c_1\\frac{dy}{dx}|_{x=a} + c_2y(a) = C_a \\ \\text{or} \\\\\n",
        " c_3\\frac{dy}{dx}|_{x=b} + c_4y(b) = C_b \\tag{1-4}\n",
        " $$\n",
        " in a heat transfer problem, it is used when a solid surface has convective heat trasnfer with a fluid.\n",
        "\n"
      ],
      "metadata": {
        "id": "e13dfTc8u-aB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<div align=\"center\">\n",
        "🟦 🟦 🟦\n",
        "</div>\n",
        "\n",
        "\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "5AmzZgNDuw9o"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 2) Shooting method"
      ],
      "metadata": {
        "id": "BihjiLedvivG"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* Consider the following problem:\n",
        "$$\n",
        "\\frac{d^2y}{dx^2} = f(x,y,\\frac{dy}{dx}) \\ \\ \\ \\text{with} \\ \\ y(a)=Y_a, \\ \\ y(b) = Y_b\n",
        "$$\n",
        "* **Step1** (transforming to canonical form): We can convert this equation into two first-order differntial equations using the change of variables:\n",
        "$$\n",
        "\\begin{cases}\n",
        "\\frac{dy}{dx} = w & y(a) = Y_a \\\\\n",
        "\\frac{dw}{dx} = f(x,y,w) & \\text{no initial condition for }w\n",
        "\\end{cases}\n",
        "$$\n",
        "* **Step2** (first guess/estimate):  Solve until you get a solution at $x = b$, and you obtain $y_{b1}$, then compare.\n",
        "$$\n",
        "w(a) = W_1\n",
        "$$\n",
        "* **Step3** (second guess/estimate): Solve until you get a solution at $x=b$, and you obtain $y_{b2}$, then compare.\n",
        "$$\n",
        "w(a) = W_2\n",
        "$$\n",
        "* **Step4:** estimate a new value for w(a) based on two previous estimates:\n",
        "$$\n",
        "W_{new} = W_1+(Y_b - y_{b1})\\frac{W_2-W_1}{y_{b2}-y_{b1}}\n",
        "$$\n",
        "* **Step5:** solve again upto $x=b$ to obtain $y_{b3}$ and compare:\n",
        " * $\\epsilon_a = |\\frac{y_{b3}-Y_b}{Y_b}|$\n",
        " * if $\\epsilon_a < tol ⇒$ **finished**.\n",
        " * otherwise, switch the third guess with either of the first or second guess and goto **step4**.\n"
      ],
      "metadata": {
        "id": "_o8_GVQGv-1R"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<div align=\"center\">\n",
        "<img src=\"https://drive.google.com/uc?id=1tkn4PIOyx0t_69IcbB2C1P9_jRyX3dYI\" width=\"500\">\n",
        "</div>"
      ],
      "metadata": {
        "id": "65AUYLm6Bpfs"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### ❓ **Example 1:**\n",
        "Solve this equation using shooting method.\n",
        "\n",
        "$$\n",
        "\\frac{d^2T}{dx^2} + 0.05(200-T)=0  \\\\\n",
        "T(0) = 300 K  \\\\\n",
        "T(10) = 400 K\n",
        "$$"
      ],
      "metadata": {
        "id": "uAhVDZME25C7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<div align=\"center\">\n",
        "<img src=\"https://drive.google.com/uc?id=1mcIVusDQgJNdNPTEzIokMp-oMlf1-6cS\" width=\"500\">\n",
        "</div>"
      ],
      "metadata": {
        "id": "iu6ugg2YCRyO"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "💡 *Solution*\n",
        "\n",
        "**Step1** (transforming to canonical form)\n",
        "$$\n",
        "\\begin{cases}\n",
        "\\frac{dT}{dx} = z & T(0) = 300 \\\\\n",
        "\\frac{dz}{dx} = -0.05(200-T) & ❔\n",
        "\\end{cases}\n",
        "$$\n",
        "**Step2** (first guess)\n",
        "$$\n",
        "Z_1 = -5 ⇒ T_{b1} = 569.754\n",
        "$$\n",
        "**Step3** (second guess)\n",
        "$$\n",
        "Z_2 = -20 ⇒ T_{b2} = 259.513\n",
        "$$\n",
        "**Step4&5** (new estimate)\n",
        "$$\n",
        "Z_{new} = -5 + \\frac{-20 -(-5)}{259.513 - 569.75}(400 - 569.75) = -13.2075\n",
        "$$\n",
        "\n",
        "&nbsp;\n",
        "\n",
        "Solving with this new $Z(0) = -13.2075$ will give $T(10) = 400.000777$.\n",
        "We found the solution!"
      ],
      "metadata": {
        "id": "nXkeyF1n39al"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Eewt8diLiMpI",
        "outputId": "5fb6a1c1-fe4f-4c24-cd96-08c41ea1d21c"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "T at x = 10 is : 400.00077670448013\n"
          ]
        }
      ],
      "source": [
        "import scipy.integrate as Int\n",
        "\n",
        "def odeSystemExp1(t, y):\n",
        "  T , z = y\n",
        "  return [z, -0.05*(200-T)]\n",
        "\n",
        "z_geuss = -13.2075;\n",
        "sol = Int.solve_ivp(odeSystemExp1, [0,10], [300, z_geuss],rtol=1.0e-5)\n",
        "print(\"T at x = 10 is :\", sol.y[0,-1])"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import scipy.integrate as Int\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "sol1 = Int.solve_ivp(odeSystemExp1, [0,10], [300, -5],rtol=1.0e-5)\n",
        "sol2 = Int.solve_ivp(odeSystemExp1, [0,10], [300, -20],rtol=1.0e-5)\n",
        "sol3 = Int.solve_ivp(odeSystemExp1, [0,10], [300, -13.2075],rtol=1.0e-5)\n",
        "\n",
        "plt.scatter(sol1.t, sol1.y[0,:], c='r')\n",
        "plt.scatter(sol2.t, sol2.y[0,:], c='b')\n",
        "plt.scatter(sol3.t, sol3.y[0,:], c='g')\n",
        "plt.legend([\"solution with z(0)=-5\",\n",
        "           \"solution with z(0)=-20\",\n",
        "           \"solution with z(0)=-13.2075\"])\n",
        "plt.grid(True)\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 430
        },
        "id": "0XRniCB_94I5",
        "outputId": "bbe54c59-64f3-4cfa-9929-58214af307cc"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGdCAYAAADJ6dNTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbpklEQVR4nO3df3xMV/4/8NdkZPJLJkElk8gIiyIkWCpSElEkSNFN9AeKbZWtTTWk1ZRSQdWP6kpKqljFLkppdJemJbISQfyoit98ypfGj0noLxMJSczc7x/TuTImvyaZTG6S1/PxmEdyzz333HNPJpl3zjn3HpkgCAKIiIiIJMKuritAREREVBqDEyIiIpIUBidEREQkKQxOiIiISFIYnBAREZGkMDghIiIiSWFwQkRERJLC4ISIiIgkpUldV6A69Ho9bt26BVdXV8hksrquDhEREVWBIAjIz8+Ht7c37OzK7x+pl8HJrVu3oFar67oaREREVA3Xr1+Hj49PufvrZXDi6uoKwHBxSqXSqmWXlJRg7969CAsLg729vVXLpkfYzrbBdrYNtrNtsJ1tozbbWavVQq1Wi5/j5amXwYlxKEepVNZKcOLs7AylUsk3fy1iO9sG29k22M62wXa2DVu0c2VTMjghloiIiCSFwQkRERFJCoMTIiIikpR6OeekKgRBwMOHD6HT6Sw6rqSkBE2aNMGDBw8sPpaqju1sG/W5neVyOZo0acLHBRA1Qg0yOCkuLoZGo0FhYaHFxwqCAJVKhevXr/OPYi1iO9tGfW9nZ2dneHl5QaFQ1HVViMiGGlxwotfrcfXqVcjlcnh7e0OhUFj0R1mv1+PevXto2rRphQ+IoZphO9tGfW1nQRBQXFyMO3fu4OrVq+jQoUO9qj8R1UyDC06Ki4uh1+uhVqvh7Oxs8fF6vR7FxcVwdHTkH8NaxHa2jfrczk5OTrC3t8dPP/0kXgMRNQ7166+VBerbH2IiMsffY6LGib/5REREZKDTAQcPGr4/eNCwXQcYnBARERGQnAy0aQNERBi2IyIM28nJNq8Kg5MG4Nq1a5DJZMjOzq5xWTKZDF9//XWNy7GWql5baGgopk2bZvXzz5kzB5MnT65y/uLiYrRp0wbff/+91etCRFRrkpOBUaOAGzdM02/eNKTbOEBhcNJIxcfHo3v37mbpGo0GQ4cOtX2FyqFWq6HRaNC1a1cAQHp6OmQyGX7//fdaP3dubi4SExPx3nvvmaQnJSWhTZs2cHR0RGBgII4dOybuUygUePvttxEXF2f1+rRp0wYymczktXjxYqufh4gaGZ0OiIkBBMF8nzFt2jSbDvEwOCmPTgekpwNffGH4Ws8eYFVdKpUKDg4OdV0NkVwuh0qlQpMmtr+x7J///Ceefvpp+Pr6imnbtm1DbGws5s6dix9++AHdunVDeHg4bt++LeYZO3YsDh48iHPnzlm9TvPnz4dGoxFfU6dOtfo5iKiRycw07zEpTRCA69cN+WyEwUkZ7HftguxPfwIGDADGjDF8reVxtx07dsDf3x9OTk5o0aIFBg0ahIKCAgCG20Hnz58PHx8fODg4oHv37vjuu+/KLWvDhg1wd3c3Sfv666/F571s2LAB8+bNw6lTp8T/wDds2ADAfFjnzJkzeOaZZ8R6TZ48Gffu3RP3//Wvf8Vzzz2HZcuWwcvLCy1atEB0dDRKSkrKrNvdu3chl8vFYQ+9Xo/mzZujT58+Yp5NmzZBrVYDMB3WuXbtGgYMGAAAaNasGWQyGf7617+Kx+n1erzzzjto3rw5VCoV4uPjy2/wP6718VebNm3E/Vu3bsXw4cNNjvnHP/6BSZMm4ZVXXoGfnx8+++wzODs74/PPPxfzNGvWDH379sXWrVsrPH91uLq6QqVSiS8XFxern4OIGhmNxrr5rIDByeOSk+E8YYJNx900Gg1Gjx6NV199FRcuXEB6ejoiIyMh/NGdlpiYiI8//hjLli3D6dOnER4ejhEjRuDHH3+s1vlefPFFvPXWW+jSpYv4H/iLL75olq+goADh4eFo1qwZjh8/ju3bt2Pfvn144403TPLt378fV65cwf79+7Fx40Zs2LBBDHYe5+bmhu7duyMjIwOAIfiRyWQ4efKkGPRkZGSgf//+Zseq1Wp89dVXAIBLly5Bo9EgMTFR3L9x40a4uLjg6NGjWLp0KebPn4/U1NRy26F0D8Tly5fRvn17hISEAAB+/fVXnD9/Hr169RLzFxcX48SJExg0aJCYZmdnh0GDBiErK8uk7N69eyOz1H8ZmzdvRtOmTSt8ZVbhv5LFixejRYsW6NGjBz766CM8fPiw0mOIiCrk5WXdfFbQ4B7CViM6HWTTpwOCALNnygoCIJMZxt1GjgTkcqudVqPR4OHDh4iMjBSHEPz9/cX9y5YtQ1xcHF566SUAwJIlS7B//34kJCQgKSnJ4vM5OTmhadOmaNKkCVQqVbn5tmzZggcPHuBf//qX+B/6ypUrMXz4cCxZsgSenp4ADD0FK1euhFwuR6dOnRAREYG0tDRMmjSpzHJDQ0ORnp6OSZMmISMjA4MHD8bFixdx8OBBDBkyBOnp6XjnnXfMjpPL5WjevDkAwMPDw6x3KCAgAHPnzgUAdOjQAStXrkRaWhoGDx5cZj2M1y4IAqKiouDm5obVq1cDAHJyciAIAry9vcX8P//8M3Q6nXjdRp6enrh48aJJmre3N3766Sdxe8SIEQgMDCyzHkatWrWqcP+bb76JP//5z2jevDkOHz6MmTNnQqPR4B//+EeFxxERVSg4GPDxMfwTXta8E5nMsD842GZVYnBSWmYmZFUddwsNtdppu3XrhoEDB8Lf3x/h4eEICwvDqFGj0KxZM2i1Wty6dQt9+/Y1OaZv3744deqU1epQlgsXLqBbt24mQwd9+/aFXq/HpUuXxA/pLl26QF4qWPPy8sKZM2fKLbd///5Yt24ddDodMjIyEB4eDpVKhfT0dAQEBODy5csIrUb7BgQEmGx7eXmZzAUpz6xZs5CVlYXvv/8eTk5OAID79+8DQLWfSurk5GSytpOrqytcXV2rdOzrr7+OTZs2idvGHqXY2FgxLSAgAAqFAn/729+waNEiSc0TIqJ6Ri4HEhMNowOPL/di3E5IsOo/5ZXhsE5pdTTuJpfLkZqaim+//RZ+fn5YsWIFOnbsiKtXr1arPDs7O3FIyKi8OSDWYG9vb7Itk8mg1+vLzR8SEoL8/HycOnUKmZmZCA0NFXtTMjIy4O3tjQ4dOtR6PQDD/Jbly5dj586dJj0XTzzxBADgt99+M0mTy+XIy8szKSMvL8+sB+rXX39Fy5YtxW1LhnXmz5+P7Oxs8VWewMBAPHz4ENeuXavwGomIKhUZCezYATzeg+vjY0iPjLRpdRiclFaH424ymQx9+/bFvHnzcPLkSSgUCuzcuRNKpRLe3t44dOiQSf5Dhw7Bz8+vzLJatmyJ/Px8cUItALMPOYVCAV0ldyB17twZp06dMinn0KFDsLOzQ8eOHS28wkfc3d0REBCAtWvXwt7eHp06dUJISAhOnjyJ3bt3lznfpHS9AVRa96rIysrCa6+9htWrV5tMyAWAdu3aQalU4vz58ybn7tmzJ9LS0sQ0vV6PtLQ0BAUFmRx/9uxZ9OjRQ9weMWKEScBR1ss4v8XDwwPt27cXX+XJzs6GnZ0dPDw8atQOREQADAHItWvAN98Ytr/5Brh61eaBCcBhHVPBwRD+GHeT2XDc7ejRo0hLS0NYWBg8PDxw9OhR3LlzB507dwYAzJgxA3PnzkW7du3QvXt3rF+/HtnZ2di8eXOZ5QUGBsLZ2RmzZs3Cm2++iaNHj5pNUG3Tpg2uXr2K7Oxs+Pj4wNXV1WxoYOzYsZg7dy4mTJiA+Ph43LlzB1OnTsW4cePM5l1Yqn///li5ciWioqIAAM2bN0fnzp2xbdu2CufR+Pr6QiaTYffu3Rg2bJg4f8ZSubm5+Mtf/oKXXnoJ4eHhyM3NBWDoxWrZsqU40fXgwYN47rnnxONiY2MxYcIE9OrVC71790ZCQgIKCgrwyiuvmJSfmZmJBQsWiNuWDOuUJSsrC0ePHsWAAQPg6uqKrKwsTJ8+HS+//DKaNWtW7XKJiEzI5UC/fkBKiuGrDYdySmPPSWlyOYTlywEAgg3H3ZRKJQ4cOIBhw4bhySefxOzZs/Hxxx+LD0N78803ERsbi7feegv+/v747rvv8N///rfcoY/mzZtj06ZNSElJgb+/P7744guz22qjoqIwZMgQDBgwAC1btsQXX3xhVo6zszP27NmDX3/9FU899RRGjRqFgQMHYuXKlTW+5pCQEOh0OpNektDQUOh0ugrnm7Rq1Qrz5s3Du+++C09PT7M7h6rq4sWLyMvLw8aNG+Hl5SW+nnrqKTHPa6+9hq1bt5oMDb344otYtmwZ3n//fXTv3h3Z2dn47rvvTIK1rKws3L17F6NGjapW3cri4OCArVu3on///ujSpQsWLlyI6dOnY82aNVY7BxGRVMiExycn1ANarRZubm64e/culEqlyb4HDx7g6tWraNu2bbUmM+r1etzfvBnOs2aZTo5Vqw2BSR10bzVEer0eWq0WSqVSsivPCoKAwMBATJ8+HaNHj67ycS+++CK6deuGWbNm1WLtqqY+tHNFavr7bCslJSVISUnBsGHDzOY+kfWwnW2jNtu5os/v0jisU4aS4cMhvPQSZIcOGSa/enkZhnLqqHuL6oZMJsOaNWsqvPPoccXFxfD398f06dNrsWZERA0bg5PyyOVWvV2Y6qfu3buXuQZReRQKBWbPnl17FSIiagTqXz8vERERNWgMToiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJMWi4CQ+Ph4ymczk1alTJ3F/aGio2f7XX3/dpIycnBxERETA2dkZHh4emDFjBh4+fGidqyEiIqJ6z+Keky5dukCj0YivgwcPmuyfNGmSyf6lS5eK+3Q6HSIiIlBcXIzDhw9j48aN2LBhA95///2aX0kjdu3aNchksgpXsK0qmUyGr7/+usblWEtVry00NBTTpk2z+vnnzJmDyZMnVzl/cXEx2rRpg++//97qdSEiaiwsDk6aNGkClUolvoxLyxs5Ozub7C/9eNq9e/fi/Pnz2LRpE7p3746hQ4diwYIFSEpKQnFxcc2vxop0OiA9HfjiC8NXKyyCKynx8fFlPlxMo9GIa/pIgVqthkajQdeuXQEA6enpkMlk+P3332v93Lm5uUhMTMR7771nkp6UlIQ2bdrA0dERgYGBOHbsmLhPoVDg7bffRlxcnFXrcu3aNUycOBFt27aFk5MT2rVrh7lz55r93pw+fRrBwcFwdHSEWq3GRx99ZNV6EBHZgsVPiP3xxx/h7e0NR0dHBAUFYdGiRWjdurW4f/Pmzdi0aRNUKhWGDx+OOXPmwNnZGYBhQTR/f3+TRdLCw8MxZcoUnDt3zmSJ+dKKiopQVFQkbmu1WgCG5/+XlJSY5C0pKYEgCNDr9SYLtlWVIAjYtcses2bJUHppHR8fAcuXC5JcWsd4nZZcs3FJpcfze3h4lJlubcbzG39W5ZHJZCZ1quhaKyvLUmvXrkVQUBDUarVY7rZt2xAbG4tPP/0UgYGBSExMRHh4OC5cuCDWc/To0Xjrrbdw5swZdOnSxSp1OX/+PHQ6HVatWoX27dvj7Nmz+Nvf/oZ79+6JAYhWq0VYWBgGDhyITz/9FGfOnMFrr70GBwcHTJ06tdZ/prVBr9dDEASUlJRALuHlI4x/hx7/e0TWxXa2jdps56qWadHCf99++y3u3buHjh07QqPRYN68ebh58ybOnj0LV1dXrFmzBr6+vvD29sbp06cRFxeH3r17Izk5GQAwefJk/PTTT9izZ49YZmFhIVxcXJCSklLuf+zx8fGYN2+eWfqWLVvEwMfI2LOjVquhUCiqemmiXbvsMWGCMwyt8mhlYpnM0EwbNxZi+HDr/8D+85//YMmSJbh69SqcnJwQEBCAzZs3w8XFBXq9HsuWLcPGjRvx888/48knn8TcuXMxaNAgAIZ5PN26dcOBAwfg7++PLVu2YObMmfjpp5/E8r/55hu8/PLL+O2337BlyxZER0ebnD8pKQljxoxBs2bNsGnTJkRERAAAzp07h5kzZ+L48eNwcnLCiBEj8MEHH6Bp06YAgL///e+4e/cu+vTpI/aARUZGYtGiRWUuGHX37l386U9/wr59+9CjRw/o9Xq0a9cO7du3R2pqKgBDADB//nycO3fO5Nrc3NzQrVs3k/JGjx6NTz/9FM8++yy6dOkCBwcH/Pvf/4ZCocArr7yCd999t9w2b9asmVmaWq3G6dOnAQBBQUF49dVXMWnSJHH/oEGD0KNHDzEg0Ov16Nq1KyZNmmSyns6IESMQGBho1utiTZ988gk+//xzcchr3bp1+OCDD3Dp0iXxvR8fH4+UlBST3p36pLi4GNevX0dubi7nphE1AIWFhRgzZox1F/4rHTwEBAQgMDAQvr6++PLLLzFx4kSTsXl/f394eXlh4MCBuHLlCtq1a1eNyzCYOXMmYmNjxW2tVgu1Wo2wsLAyVyW+fv06mjZtavEqpjodMGuWzCwwAQBBkEEmE/Dee8546SXBqmsAajQavPbaa1iyZAmee+455Ofn4+DBg3B1dUXTpk2RkJCApKQkrFq1Cj169MD69esxZswYnDlzBh06dBADBRcXFyiVSjg6OkImk5m0jZOTEwBAqVRiwoQJuHLlCvbs2YO9e/cCANzc3MQ8Tk5OUCqVKCgowPPPP48+ffrg6NGjuH37NiZPnoz33nsP69evBwDY29vj4MGDUKvV+N///ofLly9j9OjReOqpp0w+1I2USiW6d++O48ePo0ePHrh27Rrs7Oxw+vRp2NnZoWnTpjh+/DhCQ0OhVCpNrq1z587Yvn07nn/+eVy4cAFKpVKsa5MmTbB161ZMnz4dR44cQVZWFl599VUMGDAAgwcPLrPdb968KX5fUFCAYcOGoU+fPlAqlfj1119x6dIl9OvXT2zH4uJiZGdnY9asWSZtO2jQIJw8edIkLSgoCEePHhXTNm/ejClTplT4Pvjmm28QHBxcYZ7SioqK8MQTT4jnyM7ORkhIiMlQ67PPPovExEQ8fPgQzZs3r3LZUvHgwQM4OTkhJCRE8qsSp6amYvDgwVwttxaxnW2jNtvZOPJRmRot/Ofu7o4nn3wSly9fLnN/YGAgAODy5cto164dVCqV2X9weXl5AACVSlXueRwcHODg4GCWbm9vb9ZwOp0OMpkMdnZ2Fi8Rf+AATIZyHicIMly/Dhw6JLPqmoB5eXl4+PAhoqKi4OvrCwAmPQQff/wx4uLiMGbMGADA0qVLkZ6ejk8++QRJSUnidRqvufS2Uek0FxcXuLq6okmTJvD29jarj7GMrVu34sGDB/j3v/8NFxcXAMDKlSsxfPhwLF26FJ6enpDJZGjWrBmSkpIgl8vh5+eHiIgI7N+/H3/729/KvN7Q0FBkZGRg8uTJOHDgAAYPHoyLFy/i8OHDGDJkCDIyMvDOO++YXYu9vb34watSqeDu7m5SbkBAAOLj4wEAHTt2xKeffor9+/cjPDy8zHoYr10QBERFRcHNzQ1r1qyBnZ0dbty4AUEQ4OPjI9bh119/hU6ng5eXl0nbqlQqXLp0ySStVatW+Omnn8S05557DkFBQWXWo/QxVX3PXr58GStXrsSyZcvEY/Ly8tC2bVuzuhn3PT4/rD6ws7ODTCYr83ddiupLPes7trNt1EY7V7W8Gj3n5N69e7hy5Qq8vLzK3G/sbjbuDwoKwpkzZ3D79m0xT2pqKpRKJfz8/GpSFavQaKybr6q6deuGgQMHwt/fH88//zzWrl2L3377DYAhyrx16xb69u1rckzfvn1x4cIF61bkMRcuXEC3bt3EwMR4Xr1ej0uXLolpXbp0MZkP4OXlZfIzflz//v1x6NAh6HQ6ZGRkIDQ0FKGhoUhPT8etW7dw+fJlhFYj+gsICDDZrqweRrNmzUJWVhb+85//iL1H9+/fB4Bq/7fu5OSEwsJCcdvV1RXt27ev8GU89+uvv46mTZuKr8fdvHkTQ4YMwfPPP19m7xQRUX1nUXDy9ttvIyMjA9euXcPhw4fxl7/8BXK5HKNHj8aVK1ewYMECnDhxAteuXcN///tfjB8/HiEhIeKHRlhYGPz8/DBu3DicOnUKe/bswezZsxEdHV1mz4itlRNjVTtfVcnlcqSmpuLbb7+Fn58fVqxYgY4dO+Lq1avVKs/Ozg6PTyWqzQlkj0fCMpmswsmXISEhyM/Px6lTp5CZmWkSnGRkZMDb2xsdOnSo9XoAwKZNm7B8+XLs3LkTrVq1EtONvQzGINGYJpfLxd4+o7y8PLOev19//RUtW7YUtzdv3mwScJT1yszMBADMnz8f2dnZ4qu0W7duYcCAAXj66aexZs0ak30qlarMuhn3ERHVFxYFJzdu3MDo0aPRsWNHvPDCC2jRogWOHDmCli1bQqFQYN++fQgLC0OnTp3w1ltvISoqCrt27RKPl8vl2L17N+RyOYKCgvDyyy9j/PjxmD9/vtUvrDqCgw135Rgnvz5OJgPUakM+a5PJZOjbty/mzZuHkydPQqFQYOfOnVAqlfD29sahQ4dM8h86dKjc3qaWLVsiPz8fBQUFYtrjH3IKhQK6Su6P7ty5M06dOmVSzqFDh2BnZ4eOHTtaeIWPuLu7IyAgAGvXroW9vT06deqEkJAQnDx5Ert370b//v3LPdY40bOyuldFVlYWXnvtNaxevRp9+vQx2deuXTsolUqcP3/e5Nw9e/ZEWlqamKbX65GWlmY2ZHP27FmTu89GjBhhEnCU9erVqxcAwx1TpXtUjG7evInQ0FD07NkT69evNxsCCgoKwoEDB0wC0X379qFDhw5lTv4lIpIqi+acbN26tdx9arUaGRkZlZbh6+uLlJQUS05rM3I5sHy5gBdeMEx+FYTSd+sYviYkwKqTYQHg6NGjSEtLQ1hYGDw8PHD06FHcuXMHnTt3BgDMmDEDc+fORbt27dC9e3esX78e2dnZ2Lx5c5nlBQYGwtnZGbNmzcKbb76Jo0ePYsOGDSZ52rRpg6tXryI7Oxs+Pj5wdXU1670aO3Ys5s6diwkTJiA+Ph537tzB1KlTMW7cOJPbwaujf//+WLlyJaKiogAAzZs3R+fOnbFt2zYkJSWVe5yvry9kMhl2796NYcOGwcnJqcyhj8rk5ubiL3/5C1566SWEh4cjNzcXgCGAbtmyJezs7DBo0CAcPHgQzz33nHhcbGwsJkyYgF69eqF3795ISEhAQUEBXnnlFZPyMzMzsWDBAnHb1dUVrq6uFtfTyBiY+Pr6YtmyZbhz5464z9grMmbMGMybNw8TJ05EXFwczp49i08++QQLFy6s9nmJiOqEUA/dvXtXACDcvXvXbN/9+/eF8+fPC/fv369W2TqdTvjXv+4JPj56ARDEl1otCF99VdOal+38+fNCeHi40LJlS8HBwUF48sknhRUrVpjUKT4+XmjVqpVgb28vdOvWTfj222/F/VevXhUACCdPnhTTdu7cKbRv315wcnISnn32WWHNmjVC6R/3gwcPhKioKMHd3V0AIKxfv14QBEEAIOzcuVPMd/r0aWHAgAGCo6Oj0Lx5c2HSpElCfn6+uH/ChAnCyJEjTa4nJiZG6N+/f4XX/NVXXwkAhKSkJJPjAAgXL16s8Nrmz58vqFQqQSaTCRMmTBAEQRD69+8vxMTEmJxj5MiR4v7H7d+/XwBg9vL19RXzpKSkCK1atRJ0Op3JsStWrBBat24tKBQKoXfv3sKRI0dM9h8+fFhwd3cXCgsLK2wDS6xfv77M+j7+K3zq1CmhX79+goODg9CqVSth0aJFwm+//WZ2DfVFTX+fbaW4uFj4+uuvheLi4rquSoPGdraN2mznij6/S7PoOSdSodVq4ebmVuZ90g8ePMDVq1fRtm3bak1m1Ov10Gq1cHFR4tAhO2g0hjkmwcHW7zFpzIztrFQqLb6rylYEQUBgYCCmT5+O0aNHV/m4F198Ed26dcOsWbNqsXZVUx/auSI1/X22lZKSEqSkpGDYsGG8i6QWsZ1tozbbuaLP79JqdCtxQyaXw6q3C1P9I5PJsGbNGpw5c6bKxxQXF8Pf39/kgWxERGQZBidEFejevXuZaxCVR6FQYPbs2bVXISKiRqD+9fMSERFRg8bghIiIiCSFwQkRERFJCoMTIiIikhQGJ0RERCQpDE6IiIhIUhicEBERkaQwOGkArl27BplMZra4X3XIZDJ8/fXXNS7HWqp6baGhoZg2bZrVzz9nzhxMnjy5yvmLi4vRpk0bfP/991avCxFRY8HgpBw6vQ7p19LxxZkvkH4tHTp9zVfBlZL4+PgyHy6m0WgwdOhQ21eoHGq1GhqNBl27dgUApKenQyaT4ffff6/1c+fm5iIxMRHvvfeeSXpSUhLatGkDR0dHBAYG4tixY+I+hUKBt99+G3FxcVavz8KFC/H000/D2dkZ7u7uZvt/+eUXDBkyBN7e3nBwcIBarcbUqVOh1WrLLfPatWuYOHEi2rZtCycnJ7Rr1w5z585FcXGxSb7Tp08jODgYjo6OUKvVWLp0qcn+0NBQyGQys1dERISY569//avZ/iFDhtSsUYioQeITYsuw6/IuzDowCzfyb4hpPkofJA5JRGTnyDqsWe0zrnArFXK5vM7q9M9//hNPP/00fH19xbRt27YhNjYWn332GQIDA5GQkIDw8HBcunQJHh4eAAyrOb/11ls4d+4cunTpYrX6FBcX4/nnn0dQUBDWrVtntt/Ozg4jR47EBx98gJYtW+Ly5cuIjo5GXl4evvzyyzLLvHjxIvR6PVavXo327dvj7NmzmDRpEgoKCrBs2TIAhrUwwsLCMGjQIHz22Wc4c+YMXn31Vbi7u4u9SsnJySYBzS+//IJu3brh+eefNznfkCFDsH79enH78ZWwiYgAcFXix20/u12QxcsExMPkJYuXCbJ4mfDV+dpZmnj79u1C165dxdV/Bw4cKNy7d08QBMOqxPPmzRNatWolKBSKSlclXr9+veDm5mZS/s6dO8UVbMta4bamqxJ/9NFHgkqlEpo3by78/e9/L3c1y99//12ws7MTjh49Kvz2229CSUmJ0KxZMyEwMFDM8+9//1vw8fExuzbj96VfpVclnjp1qjBjxgyhWbNmgqenpzB37twK2/zxsvDYqsRdunQRVq5caXJM7969hejoaHFbp9MJ3t7ewqJFi0zyDRgwQJg9e3aF56+usn6+5UlISBC8vb0tWpV46dKlQtu2bcXtTz/9VGjWrJlQVFQkpsXFxQkdO3Yst4zly5cLrq6u4ntYEMpewboyXJWYSmM724YUViXmsE4pOr0O0/dMhwDzhZqNadO+m2b1IR6NRoPRo0fj1VdfxYULF5Ceno7IyEgIfywYnZiYiI8//hjLli3D6dOnER4ejhEjRuDHH3+s1vlefPFFvPXWW+jSpQs0Gg00Gg1efPFFs3wFBQUIDw9Hs2bNcPz4cWzfvh379u3DG2+8YZJv//79uHLlCvbv34+NGzdiw4YN2LBhQ5nndnNzQ/fu3ZGRkQEAOHPmDGQyGU6ePIl79+4BADIyMtC/f3+zY9VqNb766isAwKVLl6DRaJCYmCju37hxI1xcXHD06FEsXboU8+fPR2pqarntYLx2jUaDy5cvo3379ggJCQEA/Prrrzh//jx69eol5i8uLsaJEycwaNAgMc3Ozg6DBg1CVlaWSdm9e/dGZmamuL1582Y0bdq0wlfp/NZw69Yt7Ny5E3379rXouLt376J58+bidlZWFkJCQqBQKMQ0Y2/Rb7/9VmYZ69atw0svvQQXFxeT9PT0dHh4eKBjx46YMmUKfvnlF4vqRkSNA4d1SsnMyTQZynmcAAHXtdeRmZOJ0DahVjuvRqPBw4cPERkZKQ4h+Pv7i/uXLVuGuLg4vPTSSwCAJUuWYP/+/UhISEBSUpLF53NyckLTpk3RpEmTCodMtmzZggcPHuBf//qX+CGzcuVKDB8+HEuWLIGnpycAoFmzZli5ciXkcjk6deqEiIgIpKWlYdKkSWWWGxoaivT0dEyaNAkZGRkYPHgwLl68iIMHD2LIkCFIT0/HO++8Y3acXC4XPzQ9PDzM5l0EBARg7ty5AIAOHTpg5cqVSEtLw+DBg8ush/HaBUFAVFQU3NzcsHr1agBATk4OBEGAt7e3mP/nn3+GTqcTr9vI09MTFy9eNEnz9vbGTz/9JG6PGDECgYGBZdbDqFWrVhXur6rRo0fjP//5D+7fv49nn30Wn3zySZWPvXz5MlasWCEO6QCGuTdt27Y1yWdsg9zcXDRr1sxk37Fjx3D27FmzoachQ4YgMjISbdu2xZUrVzBr1iwMHToUWVlZkMvlll4mETVgDE5K0eRrrJqvqrp164aBAwfC398f4eHhCAsLw6hRo9CsWTNotVrcunXL7L/fvn374tSpU1atx+MuXLiAbt26mfz327dvX+j1ely6dEn8gOrSpYvJh4uXlxfOnDlTbrn9+/fHunXroNPpkJGRgfDwcKhUKqSnpyMgIACXL19GaGioxfUNCAgw2fby8sLt27crPW7WrFnIysrC999/DycnJwDA/fv3AQCOjo4W1wMwBICFhYXitqurK1xdXat07Ouvv45NmzaJ28Yepapavnw55s6di//7v//DzJkz8d5772Ht2rWVHnfz5k0MGTIEzz//fLmBZVWsW7cO/v7+6N27t0m6MbgGDMF3QEAA2rVrh/T0dAwcOLDa5yOihofDOqV4uXpZNV9VyeVypKam4ttvv4Wfnx9WrFiBjh074urVq9Uqz87OThwSMiopKbFGVctkb29vsi2TyaDX68vNHxISgvz8fJw6dQqZmZkIDQ0Ve1MyMjLg7e2NDh061Ho9AGDTpk1Yvnw5du7cadJz8cQTTwCAybDFE088Ablcjry8PJMy8vLyzHqgfv31V7Rs2VLctmRYZ/78+cjOzhZfllKpVOjUqRNGjBiBVatW4fPPP4dGU3FAfevWLQwYMABPP/001qxZY1ZeWdds3FdaQUEBtm7diokTJ1Zazz/96U944okncPny5apcFhE1IgxOSgluHQwfVx/IICtzvwwyqJVqBLcOtvq5ZTIZ+vbti3nz5uHkyZNQKBTYuXMnlEolvL29cejQIZP8hw4dgp+fX5lltWzZEvn5+SgoKBDTHv+QUygU0OkqnjvTuXNnnDp1yqScQ4cOwc7ODh07drTwCh9xd3dHQEAA1q5dC3t7e3Tq1AkhISE4efIkdu/eXeZ8k9L1BlBp3asiKysLr732GlavXo0+ffqY7GvXrh2USiXOnz9vcu6ePXsiLS1NTNPr9UhLS0NQUJDJ8WfPnkWPHj3E7REjRpgEHGW9jPNbPDw80L59e/FVE8bgrKioqNw8N2/eRGhoKHr27In169fDzs70z0JQUBAOHDhgEuCmpqaiY8eOZkM627dvR1FREV5++eVK63bjxg388ssv8PKybrBPRPUfg5NS5HZyLA9fDgBmAYpxO2FIAuR21h0fP3r0KD788EN8//33yMnJQXJyMu7cuYPOnTsDAGbMmIElS5Zg27ZtuHTpEt59911kZ2cjJiamzPICAwPh7OyMWbNm4cqVK9iyZYvZBNU2bdrg6tWryM7Oxs8//1zmh9fYsWPh6OiICRMm4OzZs9i/fz+mTp2KcePGmc27sFT//v2xfft2cQJq8+bN0blzZ2zbtq3C4MTX1xcymQy7d+/GnTt3LB7yMMrNzcVf/vIXvPTSSwgPD0dubi5yc3Nx584dAI8muh48eNDkuNjYWKxduxYbN27EhQsXMGXKFBQUFOCVV14xyZeZmYmwsDBx29XV1STgKOtlHFIqT05ODrKzs5GTkwOdTicGNcY2SElJwfr163H27Flcu3YN33zzDf7+978jMDAQbdq0AWCYD9KpUyfcvHkTwKPApHXr1li2bBnu3LkjtoXRmDFjoFAoMHHiRJw7dw7btm1DYmIiYmNjzeq4bt06PPfcc2jRooVJ+r179zBjxgwcOXIE165dQ1paGkaOHIn27dsjPDy8wusmokbI6vcJ2UBt3kqs0+mEfx3/l+DzsY/JrcTqf6hr7Tbi8+fPC+Hh4ULLli0FBwcH4cknnxRWrFhhUqf4+HihVatWgr29faW3EguC4dbh9u3bC05OTsKzzz4rrFmzRij9437w4IEQFRUluLu7W+VW4tJiYmKE/v37V3jNX331lQBASEpKMjkOgHDx4sUKr23+/PmCSqUSZDKZya3EMTExJucYOXKkuP9x+/fvr/RW4pSUFKFVq1Zmt+GuWLFCaN26taBQKITevXsLR44cMdl/+PBhwd3dXSgsLKywDSw1YcKEMuu8f/9+QRAE4X//+58QFBQkuLm5CY6OjkKHDh2Ed955R7h27Zp4Dcbrvnr1qiAIZd9WbnyVdurUKaFfv36Cg4OD0KpVK2Hx4sVm9bt48aIAQNi7d6/ZvsLCQiEsLExo2bKlYG9vL/j6+gqTJk0ScnNzK7xm3kpMpbGdbUMKtxLLBOGxyQn1gFarhZubG+7evQulUmmy78GDB7h69Sratm1brcmMer0eWq0WLk1dcOjGIWjyNfBy9UJw62Cr95g0ZsZ2ViqVZsMIUiEIAgIDAzF9+nSMHj26yse9+OKL6NatG2bNmlWLtaua+tDOFanp77OtlJSUICUlBcOGDTOb+0TWw3a2jdps54o+v0vj3TrlkNvJrXq7MNU/MpkMa9asqfDOo8cVFxfD398f06dPr8WaERE1bAxOiCrQvXv3MtcgKo9CocDs2bNrr0JERI1A/evnJSIiogaNwQkRERFJSoMNTurhPF8iegx/j4kapwYXnBhnFpd+dDgR1U/G32PemUHUuDS4CbFyuRzu7u7imirOzs6Qycp+4mtZ9Ho9iouL8eDBg3p562V9wXa2jfrazoIgoLCwELdv34a7uzsXBiRqZBpccAI8Wu+jKou+PU4QBNy/fx9OTk4WBTVkGbazbdT3dnZ3d69w5WwiapgaZHAik8ng5eUFDw8Pixe8KykpwYEDBxASEsKu5FrEdraN+tzO9vb27DEhaqQaZHBiJJfLLf7jJpfL8fDhQzg6Ota7P+b1CdvZNtjORFQf1Z9BaCIiImoUGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFIuCk/j4eMhkMpNXp06dxP0PHjxAdHQ0WrRogaZNmyIqKgp5eXkmZeTk5CAiIgLOzs7w8PDAjBkz8PDhQ+tcDREREdV7Fj+ErUuXLti3b9+jApo8KmL69On45ptvsH37dri5ueGNN95AZGQkDh06BADQ6XSIiIiASqXC4cOHodFoMH78eNjb2+PDDz+0wuUQERFRfWdxcNKkSZMy17q4e/cu1q1bhy1btuCZZ54BAKxfvx6dO3fGkSNH0KdPH+zduxfnz5/Hvn374Onpie7du2PBggWIi4tDfHw8FApFza+IiIiI6jWL55z8+OOP8Pb2xp/+9CeMHTsWOTk5AIATJ06gpKQEgwYNEvN26tQJrVu3RlZWFgAgKysL/v7+8PT0FPOEh4dDq9Xi3LlzNb0WIiIiagAs6jkJDAzEhg0b0LFjR2g0GsybNw/BwcE4e/YscnNzoVAo4O7ubnKMp6cncnNzAQC5ubkmgYlxv3FfeYqKilBUVCRua7VaAIZFzSxd2K8yxvKsXS6ZYjvbBtvZNtjOtsF2to3abOeqlmlRcDJ06FDx+4CAAAQGBsLX1xdffvklnJycLKuhBRYtWoR58+aZpe/duxfOzs61cs7U1NRaKZdMsZ1tg+1sG2xn22A720ZttHNhYWGV8tVoVWJ3d3c8+eSTuHz5MgYPHozi4mL8/vvvJr0neXl54hwVlUqFY8eOmZRhvJunrHksRjNnzkRsbKy4rdVqoVarERYWBqVSWZNLMFNSUoLU1FQMHjyYq7jWIrazbbCdbYPtbBtsZ9uozXY2jnxUpkbByb1793DlyhWMGzcOPXv2hL29PdLS0hAVFQUAuHTpEnJychAUFAQACAoKwsKFC3H79m14eHgAMERmSqUSfn5+5Z7HwcEBDg4OZun29va19gatzbLpEbazbbCdbYPtbBtsZ9uojXauankWBSdvv/02hg8fDl9fX9y6dQtz586FXC7H6NGj4ebmhokTJyI2NhbNmzeHUqnE1KlTERQUhD59+gAAwsLC4Ofnh3HjxmHp0qXIzc3F7NmzER0dXWbwQURERI2PRcHJjRs3MHr0aPzyyy9o2bIl+vXrhyNHjqBly5YAgOXLl8POzg5RUVEoKipCeHg4Pv30U/F4uVyO3bt3Y8qUKQgKCoKLiwsmTJiA+fPnW/eqiIiIqN6yKDjZunVrhfsdHR2RlJSEpKSkcvP4+voiJSXFktMSERFRI8K1dYiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFAYnREREJCkMToiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFAYnREREJCkMToiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFAYnREREJCkMToiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFAYnREREJCkMToiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUlKk7quABERkVXodEBmJqDRAF5eQHAwIJfXda2oGhicEBFR/ZecDMTEADduPErz8QESE4HIyLqrF1ULh3WIiKh+S04GRo0yDUwA4OZNQ3pyct3Ui6qtRsHJ4sWLIZPJMG3aNDEtNDQUMpnM5PX666+bHJeTk4OIiAg4OzvDw8MDM2bMwMOHD2tSFSIiaox0OkOPiSCY7zOmTZtmyEf1RrWHdY4fP47Vq1cjICDAbN+kSZMwf/58cdvZ2Vn8XqfTISIiAiqVCocPH4ZGo8H48eNhb2+PDz/8sLrVISKixigz07zHpDRBAK5fN+QLDbVZtahmqtVzcu/ePYwdOxZr165Fs2bNzPY7OztDpVKJL6VSKe7bu3cvzp8/j02bNqF79+4YOnQoFixYgKSkJBQXF1f/SoiIqPHRaKybjyShWj0n0dHRiIiIwKBBg/DBBx+Y7d+8eTM2bdoElUqF4cOHY86cOWLvSVZWFvz9/eHp6SnmDw8Px5QpU3Du3Dn06NHDrLyioiIUFRWJ21qtFgBQUlKCkpKS6lxCuYzlWbtcMsV2tg22s22wnW2jzHZWqQAnp8oPVqkA/nyqpDbfz1Ut0+LgZOvWrfjhhx9w/PjxMvePGTMGvr6+8Pb2xunTpxEXF4dLly4h+Y8JSbm5uSaBCQBxOzc3t8wyFy1ahHnz5pml792712TIyJpSU1NrpVwyxXa2DbazbbCdbcOsnb/4ovKDtFogJaV2KtRA1cb7ubCwsEr5LApOrl+/jpiYGKSmpsLR0bHMPJMnTxa/9/f3h5eXFwYOHIgrV66gXbt2lpxONHPmTMTGxorbWq0WarUaYWFhJkNG1lBSUoLU1FQMHjwY9vb2Vi2bHmE72wbb2TbYzrZRbjvv2gWMG2f4vvTEWJnM8PXf/waGD7ddReu52nw/G0c+KmNRcHLixAncvn0bf/7zn8U0nU6HAwcOYOXKlSgqKoL8sQfeBAYGAgAuX76Mdu3aQaVS4dixYyZ58vLyAAAqlarM8zo4OMDBwcEs3d7evtb+ENRm2fQI29k22M62wXa2DbN2Nj7H5PHnnKjVQEICn3NSTbXxfq5qeRYFJwMHDsSZM2dM0l555RV06tQJcXFxZoEJAGRnZwMAvLy8AABBQUFYuHAhbt++DQ8PDwCGriOlUgk/Pz9LqkNERGQQGQmMHMknxDYQFgUnrq6u6Nq1q0mai4sLWrRoga5du+LKlSvYsmULhg0bhhYtWuD06dOYPn06QkJCxFuOw8LC4Ofnh3HjxmHp0qXIzc3F7NmzER0dXWbvCBERUZXI5bxduIGw6uPrFQoF9u3bh4SEBBQUFECtViMqKgqzZ88W88jlcuzevRtTpkxBUFAQXFxcMGHCBJPnohAREVHjVePgJD09XfxerVYjIyOj0mN8fX2RwlnTREREVAaurUNERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFAYnREREJCkMToiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFAYnREREJCkMToiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFAYnREREJCkMToiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFAYnREREJCkMToiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFAYnREREJCkMToiIiEhSGJwQERGRpDA4ISIiIklhcEJEREQAAJ1eh4M5BwEAB3MOQqfX1Uk9ahScLF68GDKZDNOmTRPTHjx4gOjoaLRo0QJNmzZFVFQU8vLyTI7LyclBREQEnJ2d4eHhgRkzZuDhw4c1qQoRERHVQPKFZLRJbIOILREAgIgtEWiT2AbJF5JtXpdqByfHjx/H6tWrERAQYJI+ffp07Nq1C9u3b0dGRgZu3bqFyMhIcb9Op0NERASKi4tx+PBhbNy4ERs2bMD7779f/asgIiKiaku+kIxRX47CDe0Nk/Sb2psY9eUomwco1QpO7t27h7Fjx2Lt2rVo1qyZmH737l2sW7cO//jHP/DMM8+gZ8+eWL9+PQ4fPowjR44AAPbu3Yvz589j06ZN6N69O4YOHYoFCxYgKSkJxcXF1rkqIiIiqhKdXoeY72IgQDDbZ0yb9t00mw7xNKnOQdHR0YiIiMCgQYPwwQcfiOknTpxASUkJBg0aJKZ16tQJrVu3RlZWFvr06YOsrCz4+/vD09NTzBMeHo4pU6bg3Llz6NGjh9n5ioqKUFRUJG5rtVoAQElJCUpKSqpzCeUylmftcskU29k22M62wXa2DbZz7TiYcxC/3PsFTnZOAGD2FQB+vvczDlw9gH6t+9XoXFX92VkcnGzduhU//PADjh8/brYvNzcXCoUC7u7uJumenp7Izc0V85QOTIz7jfvKsmjRIsybN88sfe/evXB2drb0EqokNTW1VsolU2xn22A72wbb2TbYztb3RcAXZmmfd/3cZFt7VouUsyk1Ok9hYWGV8lkUnFy/fh0xMTFITU2Fo6NjtSpWHTNnzkRsbKy4rdVqoVarERYWBqVSadVzlZSUIDU1FYMHD4a9vb1Vy6ZH2M62wXa2DbazbbCda8fBnIPiJFjA0GPyedfP8erZV3Fff19M/2bMNzXuOTGOfFTGouDkxIkTuH37Nv785z+LaTqdDgcOHMDKlSuxZ88eFBcX4/fffzfpPcnLy4NKpQIAqFQqHDt2zKRc4908xjyPc3BwgIODg1m6vb19rb1Ba7NseoTtbBtsZ9tgO9sG29m6QtqGoEXTFripvWky7+S+/j7u6+9DBhl8lD4IaRsCuZ28Rueq6s/NogmxAwcOxJkzZ5CdnS2+evXqhbFjx4rf29vbIy0tTTzm0qVLyMnJQVBQEAAgKCgIZ86cwe3bt8U8qampUCqV8PPzs6Q6REREVENyOzkShyQCAGSQmewzbicMSahxYGIJi3pOXF1d0bVrV5M0FxcXtGjRQkyfOHEiYmNj0bx5cyiVSkydOhVBQUHo06cPACAsLAx+fn4YN24cli5ditzcXMyePRvR0dFl9o4QERFR7YrsHIkdL+xAzHcx+OXeL2K6j9IHCUMSENk5soKjra9ad+tUZPny5bCzs0NUVBSKiooQHh6OTz/9VNwvl8uxe/duTJkyBUFBQXBxccGECRMwf/58a1eFiIiIYLhdODMnE5p8DbxcvRDcOtisJySycyRGdhyJA1cPQHtWi2/GfGOVoZzqqHFwkp6ebrLt6OiIpKQkJCUllXuMr68vUlJqNuOXiIiIKpd8IRkx38WYPGDNR+mDxCGJZj0icjs5+rXuh5SzKejXul+dBCYA19YhIiJqsKT25NeqYnBCRETUAEnxya9VxeCEiIioAcrMyTTrMSlNgIDr2uvIzMm0Ya2qhsEJERFRA6TJ11g1ny0xOCEiImqAvFy9rJrPlhicEBERNUDBrYPho/Qxe7CakQwyqJVqBLcOtnHNKsfghIiIqAGS4pNfq4rBCRERUQNlfPJrK2Urk3QfpQ92vLDD5k9+rSqrPyGWiIiIpMP45NfKnhArJQxOiIiIGji5nRyhbULruhpVxmEdIiIikhQGJ0RERCQpDE6IiIhIUhicEBERkaRwQiwRUV3R6YDMTECjAby8gOBgQC7dOyiIbIXBCRFRXUhOBmJigBulFmbz8QESE4FIaT57gshWOKxDRGRrycnAqFGmgQkA3LxpSE9Orpt6EUkEgxMiIlvS6Qw9JoJgvs+YNm2aIR9RI8XghIjIljIzzXtMShME4Pp1Qz6iRorBCRGRLWk01s1H1AAxOCEisiUvL+vmI2qAGJwQEdlScLDhrhyZrOz9MhmgVhvyETVSDE6IiGxJLjfcLgyYByjG7YQEPu+EGjUGJ0REthYZCezYAbRqZZru42NI53NOqJHjQ9iIiOpCZCQwciSfEEtUBgYnRER1RS4HQkPruhZEksNhHSIiIpIUBidEREQkKQxOiIiISFIYnBAREZGkMDghIiIiSWFwQkRERJLC4ISIiIgkhcEJERERSQqDEyIiIpIUBidEREQkKQxOiIiISFIYnBAREZGkMDghIiIiSWFwQkRERJLC4ISIiIgkxaLgZNWqVQgICIBSqYRSqURQUBC+/fZbcX9oaChkMpnJ6/XXXzcpIycnBxEREXB2doaHhwdmzJiBhw8fWudqiIiIqN5rYklmHx8fLF68GB06dIAgCNi4cSNGjhyJkydPokuXLgCASZMmYf78+eIxzs7O4vc6nQ4RERFQqVQ4fPgwNBoNxo8fD3t7e3z44YdWuiQiIiKqzywKToYPH26yvXDhQqxatQpHjhwRgxNnZ2eoVKoyj9+7dy/Onz+Pffv2wdPTE927d8eCBQsQFxeH+Ph4KBSKal4GERERNRQWBSel6XQ6bN++HQUFBQgKChLTN2/ejE2bNkGlUmH48OGYM2eO2HuSlZUFf39/eHp6ivnDw8MxZcoUnDt3Dj169CjzXEVFRSgqKhK3tVotAKCkpAQlJSXVvYQyGcuzdrlkiu1sG2xn22A72wbb2TZqs52rWqbFwcmZM2cQFBSEBw8eoGnTpti5cyf8/PwAAGPGjIGvry+8vb1x+vRpxMXF4dKlS0hOTgYA5ObmmgQmAMTt3Nzccs+5aNEizJs3zyx97969JsNG1pSamlor5ZIptrNtsJ1tg+1sG2xn26iNdi4sLKxSPouDk44dOyI7Oxt3797Fjh07MGHCBGRkZMDPzw+TJ08W8/n7+8PLywsDBw7ElStX0K5dO0tPJZo5cyZiY2PFba1WC7VajbCwMCiVymqXW5aSkhKkpqZi8ODBsLe3t2rZ9Ajb2TbYzrbBdrYNtrNt1GY7G0c+KmNxcKJQKNC+fXsAQM+ePXH8+HEkJiZi9erVZnkDAwMBAJcvX0a7du2gUqlw7Ngxkzx5eXkAUO48FQBwcHCAg4ODWbq9vX2tvUFrs2x6hO1sG2xn22A72wbb2TZqo52rWl6Nn3Oi1+tN5oOUlp2dDQDw8vICAAQFBeHMmTO4ffu2mCc1NRVKpVIcGiIiIqLGzaKek5kzZ2Lo0KFo3bo18vPzsWXLFqSnp2PPnj24cuUKtmzZgmHDhqFFixY4ffo0pk+fjpCQEAQEBAAAwsLC4Ofnh3HjxmHp0qXIzc3F7NmzER0dXWbPCBERETU+FgUnt2/fxvjx46HRaODm5oaAgADs2bMHgwcPxvXr17Fv3z4kJCSgoKAAarUaUVFRmD17tni8XC7H7t27MWXKFAQFBcHFxQUTJkwweS4KERERNW4WBSfr1q0rd59arUZGRkalZfj6+iIlJcWS0xIREVEjwrV1iIiISFIYnBAREZGkMDghIiIiSWFwQkRERJJS7bV1iIiIapNOr0NmTiY0+Rp4uXqhj1efuq4S2QiDEyIikpzkC8mI+S4GN7Q3xLT2bu2xrO2yOqwV2QqHdYiISFKSLyRj1JejTAITALiVfwsAsOvSrrqoFtkQgxMiIpIMnV6HmO9iIEAw22dMe3ffu9DpdbauGtkQgxMiqjmdDkhPB774wvBVxw8Oqp7MnEyzHpPH3ci/gcycTBvViOoC55wQUc0kJwMxMcCNUh8oPj5AYiIQGVl39aJ6SZOvsWo+qp/Yc0JE1ZecDIwaZRqYAMDNm4b05OS6qRfVW16uXlbNR/UTgxMiqh6dztBjIpjPDRDTpk3jEA9ZJLh1MHyUPpBBVm4eH1cfBLcOtmGtyNYYnBBR9WRmmveYlCYIwPXrhnxEVSS3kyNxSCIAmAUoxu3FgxZDbie3ed3IdhicEFH1aKo45l/VfER/iOwciR0v7EArZSuT9Fauhu3hHYfXRbXIhjghloiqx6uKY/5VzUdUSmTnSIzsONLsCbF7vttT11UjG2BwQkTVExxsuCvn5s2y553IZIb9wZwbQNUjt5MjtE2ouF1SUlJ3lSGb4rAOEVWPXG64XRgwBCKlGbcTEgz5iIgswOCEiKovMhLYsQNoZTo3AD4+hnQ+54SIqoHDOkRUM5GRwMiRhrtyNBrDHJPgYPaYEFG1MTghopqTy4HQ0LquBRE1EBzWISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKgxMiIiKSFAYnREREJCkMToiIiEhSGJwQERGRpDA4ISIiIklhcEJERESSwuCEiIiIJIXBCREREUkKF/4jsoROx9V3iYhqGYMToqpKTgZiYoAbNx6l+fgAiYlAZGTd1YuIqIHhsA5RVSQnA6NGmQYmAHDzpiE9Oblu6kVE1AAxOCGqjE5n6DERBPN9xrRp0wz5iIioxhicEFUmM9O8x6Q0QQCuXzfkIyKiGuOcE6LKaDTWzUeNlk6vQ2ZOJjT5Gni5eiG4dTDkdpxQTfQ4BidElfHysm4+apSSLyQj5rsY3NA+6oXzUfogcUgiIjtzQjVRaRzWIapMcLDhrhyZrOz9MhmgVhvyEZUh+UIyRn05yiQwAYCb2psY9eUoJF/ghGqi0iwKTlatWoWAgAAolUoolUoEBQXh22+/Ffc/ePAA0dHRaNGiBZo2bYqoqCjk5eWZlJGTk4OIiAg4OzvDw8MDM2bMwMOHD61zNTVU/KAYnyasBwC8/fZf8e/NM5F+JQ06PSc6NmpyueF2YcA8QDFuJyTweSdUJp1eh5jvYiDAfEK1MW3ad9P4d4aoFIuCEx8fHyxevBgnTpzA999/j2eeeQYjR47EuXPnAADTp0/Hrl27sH37dmRkZODWrVuILPX8B51Oh4iICBQXF+Pw4cPYuHEjNmzYgPfff9+6V1UN78xcDOdZbTDzwTQAwNondmL85cUYsGkQ2izyrP5/NjodkJ4OfPGF4WtDvqOjIV9rZCSwYwfQqpVpuo+PIZ3POaFyZOZkmvWYlCZAwHXtdWTmcEI1kZFFc06GDx9usr1w4UKsWrUKR44cgY+PD9atW4ctW7bgmWeeAQCsX78enTt3xpEjR9CnTx/s3bsX58+fx759++Dp6Ynu3btjwYIFiIuLQ3x8PBQKhfWuzALvzFyMjxxmAQ4CACez/TdKfkHUtih89eJXlo0NN6aHdjWGa42MBEaO5BNiySKa/KpNlK5qPqLGoNpzTnQ6HbZu3YqCggIEBQXhxIkTKCkpwaBBg8Q8nTp1QuvWrZGVlQUAyMrKgr+/Pzw9PcU84eHh0Gq1Yu+LrRU/KMY/ij4BIADlTCkwpMswOTmm6l2vjemhXY3pWuVyIDQUGD3a8JWBCVXCy7VqE6Wrmo+oMbD4bp0zZ84gKCgIDx48QNOmTbFz5074+fkhOzsbCoUC7u7uJvk9PT2Rm5sLAMjNzTUJTIz7jfvKU1RUhKKiInFbq9UCAEpKSlBSUmLpJZhY9em/oWj2O4w9Jk52pl9LK9T/gv2XD6B/234VF6rTAXFxgKNj2ftlMuDdd4Fhw+r/h1s1r9X4c6vpz48qxna2jYrauY9XH7R3a49b+bfKnHcigwytXFuhj1cf/pwqwfezbdRmO1e1TJkglPXYy/IVFxcjJycHd+/exY4dO/DPf/4TGRkZyM7OxiuvvGISRABA7969MWDAACxZsgSTJ0/GTz/9hD179oj7CwsL4eLigpSUFAwdOrTMc8bHx2PevHlm6Vu2bIGzs7Ml1SciIqI6UlhYiDFjxuDu3btQKpXl5rO450ShUKB9+/YAgJ49e+L48eNITEzEiy++iOLiYvz+++8mvSd5eXlQqVQAAJVKhWPHjpmUZ7ybx5inLDNnzkRsbKy4rdVqoVarERYWVuHFVcWnK/6FmQVTxW0nOyd83vVzvHr2VdzX3zfLP6PVN5g9rpKekx07gIkTKz/5unWGYY/6rJrXWlJSgtTUVAwePBj29va1WMHGje1sG1Vp512XdiFuXxxu5t8U03xcfbB40GIM7zi8zGPIFN/PtlGb7Wwc+ahMjR/CptfrUVRUhJ49e8Le3h5paWmIiooCAFy6dAk5OTkICgoCAAQFBWHhwoW4ffs2PDw8AACpqalQKpXw8/Mr9xwODg5wcHAwS7e3t69xw035+zi8NWsOdEqNyZyT+/r7psGJIAO0PgjpFwJ7+0qGYry8gPvmgQ0A6GRApi+gaQp4OeUhWG5Xv58QWcG1muUr42dljZ8hVY7tbBsVtXNk10iM9BvJJ8RaAd/PtlEb7VzV8iwKTmbOnImhQ4eidevWyM/Px5YtW5Ceno49e/bAzc0NEydORGxsLJo3bw6lUompU6ciKCgIffr0AQCEhYXBz88P48aNw9KlS5Gbm4vZs2cjOjq6zODDFhSOCsQ6vImPMKvshd0AGIeJWxxPQOhHVfhDYnxo182bJmUmdwZihgA33P5IODUdPlc/rt9PiCznWkUymWE/H1BGBLmdHKFtQuu6GkSSZ9HdOrdv38b48ePRsWNHDBw4EMePH8eePXswePBgAMDy5cvx7LPPIioqCiEhIVCpVEgudaeGXC7H7t27IZfLERQUhJdffhnjx4/H/PnzrXtVFlq66F3MKPoQdlrvsjNo1cCXO7BmemTV5q+W8dCu5M7AqBeAG4+NQtX7J0TyAWVERGRlFvWcrFu3rsL9jo6OSEpKQlJSUrl5fH19kZKSYslpbWLponfxwYNYRL+93ZCwbxHw2xNAfiv46IORuFxu2eM6jA/tiomB7uYNxAz5owPmsc9vAQJkkGHad9MwsuPI+tnFW+pazZ5zkpDQcJ5zQkRENsGF/0pROCrw6fIXkJKSgm/i/47cXPuaPWfrj4d2Ze5agRunppebrfQTIuttly8fUEZERFbC4KQc/fqVOX/TcnI5NO08gVOVZ/3P/zQIfdUK56wrxgeUERER1QCDExuo6pMfExZ4IdidoyDUcOj0Ot6dQkQWY3BiA8Gtg+Gj9MFN7c0ynxBpvE0ZOcGYNs0wOsLREKrvki8kI+a7GJNF73yUPvX77jQisolqr61DVSe3kyNxSKIhLBEenxH7x/Z3CYBejuvXDdM2iOqz5AvJGPXlKLPVeOv93WlEZBMMTmwksnMkpnntALStTHdofYAvdwAXHv0nqeHipFSP6fQ6xHwXU2YvoTFt2nfTqr6IJhE1OhzWsaGRT0Yi4fWRgG8m0FQD3PMCfgoGBNMxHC8uTkr1WGZOplmPSWkN4u40IqpVDE5sKDgY8Gklx82fQvkw1XJwAmX9p8mvWtdfVfMRUePD4MSGjA9THTXKEIiUDlD4MNXGNYGyIQdhVb07rar5iKjx4ZwTGzM+TLXVY1NPfHwM6Y31NuLGNIEy+UIy2iS2wYCNAzAmeQwGbByANoltGsw1Gu9Okz3+OOQ/yCCDWqlGcOtG3EVIRBVicFIHIiOBa9eA/fuBLVsMX69ebbyBSWOaQNkYgjDj3WkAzAIU43bCkIQG01NERNbH4KSOGB+mOnq04WtjHcoBLJtAWZ81piAssnMkdrywA62Upl2EPkof7HhhR4MbpiMi6+KcE6pzjWUCZWO7iyWycyRGdhzZYOfWEFHtYXBCda6xTKBsLEFYaXI7eYMItIjItjisQ3WusUygbCxBGBFRTTE4oTrXWCZQNpYgjIiophickCQ0hgmUjSUIIyKqKc45IcloDBMojUFYWQ+bSxiS0CCCMCKimmJwQpLSGCZQNoYgjIioJhicENWBxhCEERFVF+ecEBERkaQwOCEiIiJJYXBCREREksLghIiIiCSFwQkRERFJCu/WoXpFpwMyMwGNBvDyAoKDG/eKzkREDRGDE6o3kpOBmBjgRqmFfX18gMREIJLPLiMiajA4rEP1QnIyMGqUaWACADdvGtKTk+umXkREZH0MTkjydDpDj4kgmO8zpk2bZshHRET1H4MTkrzMTPMek9IEAbh+3ZCPiIjqPwYnJHkajXXzERFR2XQ64OBBw/cHD9ZdjzSDE5I8Ly/r5iMiInPJyUCbNkBEhGE7IsKwXRdz+hickOQFBxvuypHJyt4vkwFqtSEfERFZTmo3HTA4IcmTyw23CwPmAYpxOyGBzzshIipNpwPS04EvvjB8LW+IRoo3HTA4oXohMhLYsQNo1co03cfHkM7nnBARPWIcohkwABgzxvC1vCEaKd50wIewUb0RGQmMHMknxBIRVcQ4RPN4T4hxiObxf+ikeNMBgxOqV+RyIDS0rmtBRCRNlQ3RyGSGIZqRIx/9YyfFmw44rENERNRAVGeIRoo3HTA4ISIiaiCqM0QjxZsOGJwQNWDG2fU7dlQ8W5+IGobqDtFI7aYDBidEDVRyMuDvb/h+4sSKZ+sTUcNQkyGayEjg2jXgm28M2998A1y9Wjd3Q1oUnCxatAhPPfUUXF1d4eHhgeeeew6XLl0yyRMaGgqZTGbyev31103y5OTkICIiAs7OzvDw8MCMGTPw8OHDml8NEQF4NFv/5k3TdK7iTNSw1XSIRi4H+vUzfN+vX93dDWlRcJKRkYHo6GgcOXIEqampKCkpQVhYGAoKCkzyTZo0CRqNRnwtXbpU3KfT6RAREYHi4mIcPnwYGzduxIYNG/D+++9b54qIGjkpPlCJyFaq+uCxhkxqQzTVYdGtxN99953J9oYNG+Dh4YETJ04gJCRETHd2doZKpSqzjL179+L8+fPYt28fPD090b17dyxYsABxcXGIj4+HQqGoxmUQkZEls/V5WzY1JMnJhsC89Pvfx8fQk1AfPpCtqb4/F6pGzzm5e/cuAKB58+Ym6Zs3b8amTZugUqkwfPhwzJkzB87OzgCArKws+Pv7w9PTU8wfHh6OKVOm4Ny5c+jRo4fZeYqKilBUVCRua7VaAEBJSQlKSkpqcglmjOVZu1wyxXauPRoN4ORk+N7JqcTk6+P52PzWwfezbVTUzrt2AePGGYJv4/sfAH791ZAOAMOH26KW0tK376Pv9XrDqzK1+X6uapkyQSir87dyer0eI0aMwO+//46DxvWVAaxZswa+vr7w9vbG6dOnERcXh969eyP5j0HuyZMn46effsKePXvEYwoLC+Hi4oKUlBQMHTrU7Fzx8fGYN2+eWfqWLVvEoIeIiIikrbCwEGPGjMHdu3ehVCrLzVftnpPo6GicPXvWJDABDMGHkb+/P7y8vDBw4EBcuXIF7dq1q9a5Zs6cidjYWHFbq9VCrVYjLCyswourjpKSEqSmpmLw4MGwt7e3atn0CNu59uh0hrt0bt0CHB1L8PnnqXj11cG4f9/QzjKZYSz69On608Vbm3btAuLiTCcPt2oFLFlS9f+0pf5+tsY1SkF57XzwIBARUfnx33zzaLInla8238/GkY/KVCs4eeONN7B7924cOHAAPj4+FeYNDAwEAFy+fBnt2rWDSqXCsWPHTPLk5eUBQLnzVBwcHODg4GCWbm9vX2t/CGqzbHqE7Wx99vaGD51Rox6l3b9vj/v37cXZ+osXA46OdVM/KSlvDZIrV8peg6QyUnw/W/sapeDxds7NBe7fr/y43FzD7wdVTW28n6tankV36wiCgDfeeAM7d+7E//73P7Rt27bSY7KzswEAXn888SUoKAhnzpzB7du3xTypqalQKpXw8/OzpDpEVA7jbH1vb9P0+jRbv7Y1hruaGsM1AtJcG4ZqxqLgJDo6Gps2bcKWLVvg6uqK3Nxc5Obm4v4fIeuVK1ewYMECnDhxAteuXcN///tfjB8/HiEhIQgICAAAhIWFwc/PD+PGjcOpU6ewZ88ezJ49G9HR0WX2jhBR9URGAmfOGL5ftw7Yv7/uHqgkRVJcJt7aGsM1AtJcG4ZqxqLgZNWqVbh79y5CQ0Ph5eUlvrZt2wYAUCgU2LdvH8LCwtCpUye89dZbiIqKwq5du8Qy5HI5du/eDblcjqCgILz88ssYP3485s+fb90rIyJxTsmoUYbbhjnH5BEpLhNvbY3hGgFprg1DNWPRnJPKbuxRq9XIyMiotBxfX1+kpKRYcmoiIqtqDEMBjeEajYxDmWU95yQhgT2G9U2NnnNCRFRfGYcCbt4se06GTGbYX5+HAhrDNZZW3x88Ro9w4T8iapQaw1BAY7jGx8nlhiHM0aM5lFmfMTghokarIaxBUpnGcI3U8HBYh4gatcYwFNAYrpEaFgYnRNToGYcCGrLGcI3UcHBYh4iIiCSFwQkRERFJCoMTIiIikhQGJ0RERCQpDE6IiIhIUhicEBERkaQwOCEiIiJJYXBCREREksLghIiIiCSlXj4hVvhjeU2tVmv1sktKSlBYWAitVgt7e3url08GbGfbYDvbBtvZNtjOtlGb7Wz83BbKWia7lHoZnOTn5wMA1Gp1HdeEiIiILJWfnw83N7dy98uEysIXCdLr9bh16xZcXV0he3wd8BrSarVQq9W4fv06lEqlVcumR9jOtsF2tg22s22wnW2jNttZEATk5+fD29sbdnblzyyplz0ndnZ28PHxqdVzKJVKvvltgO1sG2xn22A72wbb2TZqq50r6jEx4oRYIiIikhQGJ0RERCQpDE4e4+DggLlz58LBwaGuq9KgsZ1tg+1sG2xn22A724YU2rleToglIiKihos9J0RERCQpDE6IiIhIUhicEBERkaQwOCEiIiJJYXBSSlJSEtq0aQNHR0cEBgbi2LFjdV2lBmXRokV46qmn4OrqCg8PDzz33HO4dOlSXVerwVu8eDFkMhmmTZtW11VpkG7evImXX34ZLVq0gJOTE/z9/fH999/XdbUaFJ1Ohzlz5qBt27ZwcnJCu3btsGDBgkrXZ6GKHThwAMOHD4e3tzdkMhm+/vprk/2CIOD999+Hl5cXnJycMGjQIPz44482qRuDkz9s27YNsbGxmDt3Ln744Qd069YN4eHhuH37dl1XrcHIyMhAdHQ0jhw5gtTUVJSUlCAsLAwFBQV1XbUG6/jx41i9ejUCAgLquioN0m+//Ya+ffvC3t4e3377Lc6fP4+PP/4YzZo1q+uqNShLlizBqlWrsHLlSly4cAFLlizB0qVLsWLFirquWr1WUFCAbt26ISkpqcz9S5cuxSeffILPPvsMR48ehYuLC8LDw/HgwYPar5xAgiAIQu/evYXo6GhxW6fTCd7e3sKiRYvqsFYN2+3btwUAQkZGRl1XpUHKz88XOnToIKSmpgr9+/cXYmJi6rpKDU5cXJzQr1+/uq5GgxcRESG8+uqrJmmRkZHC2LFj66hGDQ8AYefOneK2Xq8XVCqV8NFHH4lpv//+u+Dg4CB88cUXtV4f9pwAKC4uxokTJzBo0CAxzc7ODoMGDUJWVlYd1qxhu3v3LgCgefPmdVyThik6OhoREREm72uyrv/+97/o1asXnn/+eXh4eKBHjx5Yu3ZtXVerwXn66aeRlpaG//u//wMAnDp1CgcPHsTQoUPruGYN19WrV5Gbm2vy98PNzQ2BgYE2+Vyslwv/WdvPP/8MnU4HT09Pk3RPT09cvHixjmrVsOn1ekybNg19+/ZF165d67o6Dc7WrVvxww8/4Pjx43VdlQbt//2//4dVq1YhNjYWs2bNwvHjx/Hmm29CoVBgwoQJdV29BuPdd9+FVqtFp06dIJfLodPpsHDhQowdO7auq9Zg5ebmAkCZn4vGfbWJwQnViejoaJw9exYHDx6s66o0ONevX0dMTAxSU1Ph6OhY19Vp0PR6PXr16oUPP/wQANCjRw+cPXsWn332GYMTK/ryyy+xefNmbNmyBV26dEF2djamTZsGb29vtnMDxWEdAE888QTkcjny8vJM0vPy8qBSqeqoVg3XG2+8gd27d2P//v3w8fGp6+o0OCdOnMDt27fx5z//GU2aNEGTJk2QkZGBTz75BE2aNIFOp6vrKjYYXl5e8PPzM0nr3LkzcnJy6qhGDdOMGTPw7rvv4qWXXoK/vz/GjRuH6dOnY9GiRXVdtQbL+NlXV5+LDE4AKBQK9OzZE2lpaWKaXq9HWloagoKC6rBmDYsgCHjjjTewc+dO/O9//0Pbtm3rukoN0sCBA3HmzBlkZ2eLr169emHs2LHIzs6GXC6v6yo2GH379jW7Hf7//u//4OvrW0c1apgKCwthZ2f6cSWXy6HX6+uoRg1f27ZtoVKpTD4XtVotjh49apPPRQ7r/CE2NhYTJkxAr1690Lt3byQkJKCgoACvvPJKXVetwYiOjsaWLVvwn//8B66uruK4pZubG5ycnOq4dg2Hq6ur2TweFxcXtGjRgvN7rGz69Ol4+umn8eGHH+KFF17AsWPHsGbNGqxZs6auq9agDB8+HAsXLkTr1q3RpUsXnDx5Ev/4xz/w6quv1nXV6rV79+7h8uXL4vbVq1eRnZ2N5s2bo3Xr1pg2bRo++OADdOjQAW3btsWcOXPg7e2N5557rvYrV+v3A9UjK1asEFq3bi0oFAqhd+/ewpEjR+q6Sg0KgDJf69evr+uqNXi8lbj27Nq1S+jatavg4OAgdOrUSVizZk1dV6nB0Wq1QkxMjNC6dWvB0dFR+NOf/iS89957QlFRUV1XrV7bv39/mX+TJ0yYIAiC4XbiOXPmCJ6enoKDg4MwcOBA4dKlSzapm0wQ+Ig9IiIikg7OOSEiIiJJYXBCREREksLghIiIiCSFwQkRERFJCoMTIiIikhQGJ0RERCQpDE6IiIhIUhicEBERkaQwOCEiIiJJYXBCREREksLghIiIiCSFwQkRERFJyv8HaV41yhjQ1lEAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### ❓ **Example 2:**\n",
        "Solve this equation using shooting method.\n",
        "\n",
        "$$\n",
        "\\frac{d^2T}{dx^2} + 0.05(200-T)=0  \\\\\n",
        "\\frac{dT}{dx}(0) = \\frac{1}{200}(T(0)-200)    \\\\\n",
        "T(10) = 400 K\n",
        "$$"
      ],
      "metadata": {
        "id": "leOwQzAF_26P"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<div align=\"center\">\n",
        "<img src=\"https://drive.google.com/uc?id=1EfSF7gOP8lnCiRJ6pDxsE-TK2TGMwLWs\" width=\"500\">\n",
        "</div>"
      ],
      "metadata": {
        "id": "WP3j7B_TCisL"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "💡 *Solution*\n",
        "\n",
        "**Step1** (transforming to canonical form)\n",
        "$$\n",
        "\\begin{cases}\n",
        "\\frac{dT}{dx} = z & ❔ \\\\\n",
        "\\frac{dz}{dx} = -0.05(200-T) & z(0) = \\frac{1}{200}(T(0)-200)\n",
        "\\end{cases}\n",
        "$$\n",
        "\n",
        "**Step2** (first guess): the guess for T at $x = 0$\n",
        "$$\n",
        "T_1 = 300 ⇒ z(0) = \\frac{1}{200}(300-200) = 0.5\n",
        "$$\n",
        "If we integrate the above equations up to $x=10$ we will get $T(10) = T_{b1} = 683.51$.\n",
        "\n",
        "**Step3** (second guess): the guess for T at $x=0$:\n",
        "$$\n",
        "T_2 = 250 ⇒ z(0) = \\frac{1}{200}(250-200) = 0.25\n",
        "$$\n",
        "If we integrate the above equations up to $x=10$ we will get $T(10) = T_{b2} = 441.75$.\n",
        "\n",
        "**Step4&5** (new estimate):\n",
        "$$\n",
        "T_{new} = 300 + \\frac{250-300}{441.75-683.51}(400 - 683.51) = 241.365\n",
        "$$\n",
        "Solving with new $T(0) = 241.365$ upto x=10, we get $T(10) = 400.0038$.\n",
        "This is close enough to 400, so we stop the solution here."
      ],
      "metadata": {
        "id": "GbTsqVoUEZ9E"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import scipy.integrate as Int\n",
        "\n",
        "def odeSystemExp1(t, y):\n",
        "  T , z = y\n",
        "  return [z, -0.05*(200-T)]\n",
        "\n",
        "T_guess = 241.365;\n",
        "z_0 = 1.0/200.0*(T_guess - 200)\n",
        "sol = Int.solve_ivp(odeSystemExp1, [0,10], [T_guess, z_0],rtol=1.0e-5)\n",
        "print(\"T at x = 10 is :\", sol.y[0,-1])"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "gzO68aXW90uN",
        "outputId": "258c4653-c4c7-4d51-c19a-1ab3f92ac309"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "T at x = 10 is : 400.00385678339063\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 3) Finite difference method ODE (BVP)\n",
        "* Conside an ODE with boundary conditions that holds in $a\\le x \\le b$.\n",
        "* We divide this domain in equally-spaced intervals. These intervals are separated by $Δx$ and the ends are repesented by nodes.\n",
        "* Since the ODE holds on every point in domain, it also hold on nodes.\n",
        "* We use finite difference approximations for estimating the derivatives."
      ],
      "metadata": {
        "id": "TO4trNz5JA4Y"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* For example consider the 1D heat equation in a rod.\n",
        "$$\n",
        "\\frac{d^2T}{dx^2} + 0.05(200-T)=0  \\\\\n",
        "T(0) = 300 K  \\\\\n",
        "T(10) = 400 K\n",
        "\\tag{3-1}\n",
        "$$\n",
        "* We divide $x$ length  into $n$ intervals:\n",
        " * We will have $n+1$ nodes\n",
        " * Temperatures at the first node and the last node are are known (boundary conditions)\n",
        " * Temperatures at $n-1$ nodes are unknown.  \n",
        "<div align=\"center\">\n",
        "<img src=\"https://drive.google.com/uc?id=1dj-Rflc9y3i9lqPgfb_I_Wwqm0Fhu5Ea\" width=\"500\">\n",
        "</div>\n",
        "\n"
      ],
      "metadata": {
        "id": "SRqrkgPYVmFC"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* We use finite difference formulas for evaluating the derivatives for node $i$:\n",
        "$$\n",
        "\\frac{d^2T}{dx^2} ≃ \\frac{T_{i+1}-2T_i+T_{i-1}}{Δx^2}\n",
        "\\tag{3-2}\n",
        "$$\n",
        "\n",
        "* This is substituted into the ODE and we get:\n",
        "$$\n",
        "\\frac{T_{i+1}-2T_i+T_{i-1}}{Δx^2} + 0.05(200-T_i)=0\n",
        "\\tag{3-3}\n",
        "$$\n",
        "* re-arranging the above equation will give:\n",
        "$$\n",
        "-T_{i-1} +(2+0.05Δx^2)T_i - T_{i+1} = 10 Δx^2\n",
        "\\tag{3-4}\n",
        "$$\n",
        "* We can write the above equation for n-1 internal nodes and this is equal to the number of unknown temperatures.\n",
        "* The above equation forms a set of linear equations whose solution will give the temperature distribution in the rod."
      ],
      "metadata": {
        "id": "9jxrQ_WjZndv"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### ❓ **Example 3:**\n",
        "\n",
        "Solve the following equation using finite difference method with $\\Delta x = 2 \\ m$.\n",
        "$$\n",
        "\\frac{d^2T}{dx^2} + 0.05(200-T)=0  \\\\\n",
        "T(0) = 300 K  \\\\\n",
        "T(10) = 400 K\n",
        "$$\n"
      ],
      "metadata": {
        "id": "JumPGqfrJHIU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "💡 *Solution*\n",
        "\n",
        "Since $Δx = 2$, we have 6 nodes, indexed by $i=0,1, ...,5$.\n",
        "\n",
        "Nodes 0 and 5 are boundary nodes and their temperatures are known by boundary conditions.\n",
        "$$\n",
        "T_0 = 300 \\\\\n",
        "T_5 = 400\n",
        "$$\n",
        "\n",
        "using Eq. (3-4) and setting $Δx = 2$, we get:\n",
        "$$\n",
        "-T_{i-1} +2.2T_i - T_{i+1} = 40\n",
        "$$\n",
        "\n",
        "* For $i=1$:\n",
        "$$\n",
        "-T_0 + 2.2T_1-T_2 = 40 ⇒ 2.2T_1-T_2 = 340\n",
        "$$\n",
        "* For $i=4$:\n",
        "$$\n",
        "-T_3 + 2.2T_4-T_5 = 40 ⇒ -T_3 + 2.2T_4 = 440\n",
        "$$\n",
        "\n",
        "* For the rest of nodes:\n",
        "$$\n",
        "i = 2 ⇒ -T_1 +2.2T_2 - T_3 = 40   \\\\\n",
        "i = 3 ⇒ -T_2 +2.2T_3 - T_4 = 40\n",
        "$$\n",
        "\n"
      ],
      "metadata": {
        "id": "1NnIkOM5cJ_C"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* The above equations form a set that should be solved.\n",
        "$$\n",
        "\\begin{bmatrix}\n",
        " 2.2 & -1 & 0 & 0 \\\\\n",
        " -1  & 2.2& -1& 0 \\\\\n",
        " 0 & -1 & 2.2 & -1 \\\\\n",
        " 0 & 0 & -1 & 2.2\n",
        "\\end{bmatrix}\n",
        "\\begin{bmatrix}\n",
        "T_1 \\\\ T_2 \\\\ T_3 \\\\ T_4\n",
        "\\end{bmatrix}\n",
        "=\n",
        "\\begin{bmatrix}\n",
        "340 \\\\ 40 \\\\ 40\\\\ 440\n",
        "\\end{bmatrix}\n",
        "$$\n",
        "The solution of the above set is:\n",
        "$$\n",
        "T =\n",
        "\\begin{bmatrix}\n",
        "283.2660 \\\\\n",
        "283.1853 \\\\\n",
        "299.7416 \\\\\n",
        "336.2462\n",
        "\\end{bmatrix}\n",
        "$$\n",
        "\n"
      ],
      "metadata": {
        "id": "DfKzv6QqpRl7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### ❓ **Example 4 (derivative boundary condition):**\n",
        "Solve the following equation using finite difference method with $\\Delta x = 2 \\ m$.\n",
        "$$\n",
        "\\frac{d^2T}{dx^2} + 0.05(200-T)=0  \\\\\n",
        "\\frac{dT}{dx}|_{x=0} = -20  \\\\\n",
        "T(10) = 400 K\n",
        "$$\n",
        "\n",
        "* Generate the solution for $\\frac{dT}{dx}|_{x=0} = 0$ (thermal isolation)."
      ],
      "metadata": {
        "id": "TTJ74k_UqdGg"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "💡 *Solution*\n",
        "\n",
        "Since $Δx = 2$, we have 6 nodes, indexed by $i=0,1, ...,5$.\n",
        "\n",
        "Nodes 0 and 5 are boundary nodes. The temperature at node 5 is known because of constant temperature boundary condition. But the temperature at node 0 is not known (because the condition is derivative).\n",
        "$$\n",
        "T_5 = 400\n",
        "$$\n",
        "\n",
        "So, we have 5 unknown temperature from $T_0$ to $T_4$.\n",
        "\n",
        "Using Eq. (3-4) and setting $Δx = 2$, we get:\n",
        "$$\n",
        "-T_{i-1} +2.2T_i - T_{i+1} = 40\n",
        "$$\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "c4TmjN6ZxUQK"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* For $i=0$ (Node on $x = 0$ )\n",
        "\n",
        " For this node we have:\n",
        " $$-T_{-1} +2.2T_0 - T_1 = 40. \\tag{*}$$\n",
        " But node (-1) does not exist. But we can consider this node as a imaginary node and apply the boundary condition to this node.\n",
        "\n",
        " $$\n",
        " \\frac{dT}{dx}|_{x=0} = \\frac{T_1-T_{-1}}{2Δx} ⇒ \\\\\n",
        " T_{-1} = T_1 - 2Δx\\frac{dT}{dx} = T_1 - 2(2)(-20) = T_1 + 80 \\tag{**}\n",
        " $$\n",
        "\n",
        " Combining (\\*) and (**) to eliminate $T_{-1}$ yields:\n",
        " $$\n",
        " +2.2T_0 - 2T_1 = 120\n",
        " $$\n",
        "\n",
        " <div align=\"center\">\n",
        " <img src=\"https://drive.google.com/uc?id=1hjby9BMBeBHpwqn9bv0hUGaPmbM1ONOQ\" width = \"350\">\n",
        " </div>\n",
        "\n"
      ],
      "metadata": {
        "id": "qbdIknbt2a-l"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* The above equations form a set that should be solved.\n",
        "$$\n",
        "\\begin{bmatrix}\n",
        " 2.2 & -2 & 0 & 0 & 0 \\\\\n",
        " -1  & 2.2& -1& 0 & 0 \\\\\n",
        " 0 & -1 & 2.2 & -1& 0 \\\\\n",
        " 0 & 0 & -1 & 2.2& -1 \\\\\n",
        " 0 & 0& 0 & -1 & 2.2\n",
        "\\end{bmatrix}\n",
        "\\begin{bmatrix}\n",
        "T_0 \\\\T_1 \\\\ T_2 \\\\ T_3 \\\\ T_4\n",
        "\\end{bmatrix}\n",
        "=\n",
        "\\begin{bmatrix}\n",
        "120 \\\\ 40 \\\\ 40 \\\\ 40\\\\ 440\n",
        "\\end{bmatrix}\n",
        "$$\n",
        "\n",
        "&nbsp;\n",
        "\n",
        "The solution of the above set is:\n",
        "$$\n",
        "T =\n",
        "\\begin{bmatrix}\n",
        "328.2710 \\\\\n",
        "301.0981 \\\\\n",
        "294.1448 \\\\\n",
        "306.0204 \\\\\n",
        "339.1002\n",
        "\\end{bmatrix}\n",
        "$$\n",
        "\n",
        "&nbsp;\n",
        "&nbsp;"
      ],
      "metadata": {
        "id": "ZA_O6uI23oDH"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* For the second case $T^{'}(0) = 0$, we start with Eq. (**) and we have:\n",
        " $$\n",
        " T_{-1} = T_1 - 2Δx\\frac{dT}{dx} = T_1 - 2(2)(0) = T_1  \\tag{***}\n",
        " $$\n",
        "\n",
        " Combining (\\*) and (***) to eliminate $T_{-1}$ yields:\n",
        "$$\n",
        " +2.2T_0 - 2T_1 = 40\n",
        "$$\n",
        "\n",
        " &nbsp;\n",
        "\n",
        "$$\n",
        "\\begin{bmatrix}\n",
        " 2.2 & -2 & 0 & 0 & 0 \\\\\n",
        " -1  & 2.2& -1& 0 & 0 \\\\\n",
        " 0 & -1 & 2.2 & -1& 0 \\\\\n",
        " 0 & 0 & -1 & 2.2& -1 \\\\\n",
        " 0 & 0& 0 & -1 & 2.2\n",
        " \\end{bmatrix}\n",
        " \\begin{bmatrix}\n",
        " T_0 \\\\T_1 \\\\ T_2 \\\\ T_3 \\\\ T_4\n",
        " \\end{bmatrix}\n",
        " =\n",
        " \\begin{bmatrix}\n",
        " 40 \\\\ 40 \\\\ 40 \\\\ 40\\\\ 440\n",
        " \\end{bmatrix}\n",
        " $$\n",
        "\n",
        " &nbsp;\n",
        "\n",
        " The solution of the above set is:\n",
        " $$\n",
        " T =\n",
        " \\begin{bmatrix}\n",
        " 243.0278 \\\\\n",
        " 247.3306 \\\\\n",
        " 261.0994 \\\\\n",
        " 287.0882 \\\\\n",
        " 330.4946\n",
        " \\end{bmatrix}\n",
        " $$\n",
        "\n"
      ],
      "metadata": {
        "id": "n9szqaO51jWz"
      }
    }
  ]
}