01 - FEM Basics 1D.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "mimetype": "text/x-python",
      "nbconvert_exporter": "python",
      "name": "python",
      "pygments_lexer": "ipython3",
      "version": "3.5.4",
      "file_extension": ".py",
      "codemirror_mode": {
        "version": 3,
        "name": "ipython"
      }
    },
    "colab": {
      "name": "01 - FEM Basics 1D.ipynb",
      "provenance": [],
      "include_colab_link": true
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "dfdaf8ee303f43bbbda0a76f09396220": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "VBoxModel",
          "state": {
            "_view_name": "VBoxView",
            "_dom_classes": [
              "widget-interact"
            ],
            "_model_name": "VBoxModel",
            "_view_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_view_count": null,
            "_view_module_version": "1.5.0",
            "box_style": "",
            "layout": "IPY_MODEL_cc3c2018633a4007b36230fe78352f4b",
            "_model_module": "@jupyter-widgets/controls",
            "children": [
              "IPY_MODEL_3ccbd2cf9c634c31936efc3c95c0e8cd",
              "IPY_MODEL_a50d77c2c715494a90dfc9a6e491ddfc"
            ]
          }
        },
        "cc3c2018633a4007b36230fe78352f4b": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "state": {
            "_view_name": "LayoutView",
            "grid_template_rows": null,
            "right": null,
            "justify_content": null,
            "_view_module": "@jupyter-widgets/base",
            "overflow": null,
            "_model_module_version": "1.2.0",
            "_view_count": null,
            "flex_flow": null,
            "width": null,
            "min_width": null,
            "border": null,
            "align_items": null,
            "bottom": null,
            "_model_module": "@jupyter-widgets/base",
            "top": null,
            "grid_column": null,
            "overflow_y": null,
            "overflow_x": null,
            "grid_auto_flow": null,
            "grid_area": null,
            "grid_template_columns": null,
            "flex": null,
            "_model_name": "LayoutModel",
            "justify_items": null,
            "grid_row": null,
            "max_height": null,
            "align_content": null,
            "visibility": null,
            "align_self": null,
            "height": null,
            "min_height": null,
            "padding": null,
            "grid_auto_rows": null,
            "grid_gap": null,
            "max_width": null,
            "order": null,
            "_view_module_version": "1.2.0",
            "grid_template_areas": null,
            "object_position": null,
            "object_fit": null,
            "grid_auto_columns": null,
            "margin": null,
            "display": null,
            "left": null
          }
        },
        "3ccbd2cf9c634c31936efc3c95c0e8cd": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "IntSliderModel",
          "state": {
            "_view_name": "IntSliderView",
            "style": "IPY_MODEL_663fdb6fa5b8424a9e5a8ce25a41da23",
            "_dom_classes": [],
            "description": "n",
            "step": 1,
            "_model_name": "IntSliderModel",
            "orientation": "horizontal",
            "max": 11,
            "_view_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "value": 5,
            "_view_count": null,
            "disabled": false,
            "_view_module_version": "1.5.0",
            "min": 1,
            "continuous_update": true,
            "readout_format": "d",
            "description_tooltip": null,
            "readout": true,
            "_model_module": "@jupyter-widgets/controls",
            "layout": "IPY_MODEL_de608d5fc21a4866b725cd74a04a6bd6"
          }
        },
        "a50d77c2c715494a90dfc9a6e491ddfc": {
          "model_module": "@jupyter-widgets/output",
          "model_name": "OutputModel",
          "state": {
            "_view_name": "OutputView",
            "msg_id": "",
            "_dom_classes": [],
            "_model_name": "OutputModel",
            "outputs": [
              {
                "output_type": "display_data",
                "metadata": {
                  "tags": [],
                  "needs_background": "light"
                },
                "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU5b348c8z2RNCEnYIJGELyCIiEAIioKCCCshOEhXbKr9q1VZtqbd28d6rt629dW9vRVxQEggKxKBsimzKGtYkkECIISQkhCUL2Zd5fn9MiAGyzGQmQ8n5vl+veTEz5znnex4YvnPmnOd8H6W1RgghRNtnutE7IIQQwjkk4QshhEFIwhdCCIOQhC+EEAYhCV8IIQxCEr4QQhiEQxK+UupDpVSeUiqpkeUTlVKFSqnDtY8/OiKuEEII67k6aDsfA+8CnzTRZqfW+kEHxRNCCGEjhxzha613AJccsS0hhBCtw1FH+NYYo5Q6ApwFfq21Tm6okVJqEbAIwMfHZ8TAgQOduItCCHFzO3DgwAWtdeeGljkr4R8EgrXWxUqp+4E4oH9DDbXWS4AlACNHjtQJCQlO2kUhhLj5KaVON7bMKaN0tNZFWuvi2ufrATelVCdnxBZCCGHhlISvlOqmlFK1z8Nq4150RmwhhBAWDjmlo5RaAUwEOimlsoA/AW4AWut/AXOAJ5VS1UAZsEBLmU4hhHAqhyR8rXVEM8vfxTJsUwghxA0id9oKIYRBSMIXQgiDkIQvhBAGIQlfCCEMQhK+EEIYhCR8IYQwCEn4QghhEJLwhRDCICThCyGEQUjCF0IIg5CEL4QQBiEJXwghDEISvhBCGIQkfCGEMAhJ+EIIYRCS8IUQwiAk4QshhEFIwhdCCIOQhC+EEAZhiIS/dOlSlFL4+vpSXl5+1bLc3FyUUvz9739vM3FvZGwj9lmIm4UhEv6hQ4fw8PCguLiYb7755rplAMOHD28zcW9kbCP2WYibhSES/uHDhwkLC2PYsGHExcVdtexKIrjtttvaTNwbGduIfRbiZtHmE77WmqNHj3Lbbbfx0EMPsW7dOsxmc93yw4cPExQURIcOHdpE3BsZ24h9FuJm0uYT/smTJykuLmb48OHMmDGDvLw8du3aVbf80KFDrfIz/0bFvZGxjdhnIW4mbT7hHz58GLD8lB8+fDhBQUF1P/cvX77MqVOnWu3Uxo2IeyNjG7HPQtxMDJHw3dzcGDx4MAAzZsyoSwRHjhxBa1135Hfq1CnGjRtHaGgow4cPJyEhwSlx//u//5vQ0FBMJtN1555bM3Z+fj4PPvggoaGhDBs2jHvvvZe0tLRWjwswf/58br31VoYPH05YWBhbtmyxp8s2xb7io48+QinlkL9zIW4GbT7hHzp0iFtuuQV3d3cAHnroIU6dOkViYuJ1Izd+/vOfs3DhQk6cOMFrr71GVFQUWutWj3vPPfewceNGxo8fb1dfbY2tlOJXv/oVJ06c4MiRIzz44IM8/vjjrR4X4L333uPo0aMcOnSI9957j7lz5151zr01YwNkZGTw/vvvEx4e3uKYQtxsHJLwlVIfKqXylFJJjSxXSqm3lVJpSqmjSqnbHRHXGocPH77qP/r48eMJCAggLi6OQ4cO0aFDB4KCgjh//jx79uzhscceAyxJWGvNgQMHWjUuQHh4OH369Gl5J1sY29/fn8mTJ9e1Gzt2LBkZGa0eF8Df37+uXWFhYYtjtiS22Wzm8ccf55133sHDw8Pu2ELcLFwdtJ2PgXeBTxpZPhXoX/sYDfxf7Z+t6ty5c+Tm5l517tbV1ZUHHniAuLg4tNZ1yzIzM+nevTtubm51bUNCQsjMzGTkyJGtFtfR7In95ptvMmPGDKtjRUdH89JLL5GZmUlgYKDNcZ977jm++OILCgsLWb16NSZTy44/bO3z66+/zh133MGIESNaFE+Im5VDjvC11juAS000mQF8oi32AP5Kqe6OiN2QnMIy8orKGx17PWPGDA4ePEhiYqJjR26Ya6C8yPlxgaLKIrTWLY79n//5n6Snp/PnP//ZqnjR0dEsWrSIs6dP4wlkZWUBcP78eavjvvHGG6SnpxMdHc3ixYuprKy0trsAVJRWAY2PsW8odlJSEqtXr+b3v/+9TbGEaAucdQ4/EDhT73VW7XsOV1xRzd3/u51/bjt11ciN+qZMmYKnpyfV1dV1y4KCgsjJyaGqqqquXUZGRt1pAKts/R94axiH9++xOq4jnC0+y6RVk1iVusqmPl/xyiuvsH79ejZs2IC3t7dVMV988UVKS0t5tXt3VgWHoGrfX7JkCZcu/fjdb02fp0yZQn5+PomJiVb2GHLTC/nghZ1kJF6wqc87d+4kIyOD/v37ExISwp49e1i0aBHvvvuu1bGFuGlprR3yAEKApEaWfQmMq/d6CzCykbaLgAQgISgoSLfEsysO6iF/3KiLy6tsWu/uu+/WS5Ys0VprvXnzZt2vXz9tNputW7miROs/B2n9p/Za73nP1l3WWms9YcIEvXbtWpvXez3hdT3k4yH6gTUP6BpzjU3rvvzyyzosLEwXFBQ027aiokKvXbtWz5gxQwO6q6urPho6QB8bMFCP9/HRQN2jX79+OjIyUr/55pt6165durS0tG47paWlOj09ve71rl27tL+/v7506ZLV+73x/UT97v/bouPeOGhTf6/V0r9zIf5dAQm6kTztqHP4zckGetV73bP2vetorZcASwBGjhzZoiEyj44J4YvDZ1l7KJuHw4OtXu9f//oXCxcu5G9/+xve3t5ER0ejlGp+RYCk1VBeAD6dYf/7EPYEWLnuyy+/zNKlSzl//jxJSUk8/fTT7Nmzh549eza7bkVNBWtOrqGDZwdOF51mz9k9jA0ca1Xc5ORkXn75Zfr27cuECRMAy7nv+sNRtdYcPHiQZcuWERMTw8WLF+natSu+vr7Mc/fABORXVxPpH8COkhK6dOnCc889x759+9i+fTsxMTF12x06dChhYWEMGjSIDz74gJqaGlxdXfHx8WH16tUEBARYtd8lhRWkHzyPl68bWSn55OeWENDNx6p1hTC0xr4JbH3Q9BH+A8AGQAHhwD5rtjlixIgWfcOZzWZ9/1s79L2vb7f+CN0eZrPW/zdO63+Ea30o2nKUf2pr68fVWsedjNNDPh6id5zZocevHK+f/uZph2z37Nmz+rXXXtODBw/WgPbw8NDz5s3TX331la6qqtLRH3+sd/brr98NDNRPdeyojw0YqEPbt9fLly+/ajvZ2dk6Li5O/+53v9OTJ0/Wfn5+db8CfHx89Pjx4/Wvf/1rvWrVKp2RkWHVv9feden63f+3RZ9NK9D//MW3evvKVIf0WYi2gNY+wldKrQAmAp2UUlnAnwC32i+UfwHrgfuBNKAU+Ikj4jaxPywcE8Li1UfZ+8Mlwvt0bM1wkLUfco/CA6/D4Fmw+few733oM7F14wIrU1bSx68P4wLHMbv/bJYmLiW7OJvAdrZfIikvL+eLL75g2bJlbNq0CbPZTHh4OP/3f//H/PnzrzoCf6BDR866uLDF3Z3vzp3j5x078d7MmUyMirpqmz169GDGjBl1o3/MZjNpaWns27ev7vH222/XXbDt0qULo0aNIiwsjLCwMEaNGkXHjj/++9XUmEnemY1rQDljJg1jYt9IysrGklmewMMLI1vyVyiEYTgk4WutI5pZroFfOCKWtaYN68Gr64/z6e7TrZ/w9y0Bj/Zw63xw84TbH4Xv34KCM+Dfq/n1WyjxfCJJF5P43ejfoZRi3oB5fJj0IbGpsTw/4nmrtqG1Zvfu3SxbtozY2FgKCwvp1asXL774Io8++igDBgxocL386Gjcg4NZvWE9ymQi+/nncf/ue8xlZZi8vBqNZzKZCA0NJTQ0lIcffhiAyspKEhMTr/oSWL9+fd1Nb3379q37EujfZQSlhdV8+O1fOX36NNvL4hjVfzIfvfkxylUTdc0XjhDiR232Tlsvdxfmj+rFxuRccgvLm1+hpYrzIDkObosEj3aW90b+1PLngY9aLy6wMnUl3q7eTOszDYBuPt24q9ddrD25lvLqpvucmZnJK6+8woABA7jjjjtYvnw506dP55tvviEjI4NXX3210WRflpxM2eHDBERGoGrHzgdERmIuKqLwyy9t7oe7uzsjRozgySef5KOPPiI5OZmCggK2bt3KX//6V2677Ta+//57nn/+eb6JPcCForMcSvsOgIy842SeP0F4/6m89NJLNscWwkjabMIHeHh0MGatidmX2XpBDiwDcxWMqleSwD8IQqdallVXtErYS+WX2PDDBqb3nU4793Z170cMjKCgooCNGRuvW6e4uJhPPvmESZMmERISwh/+8Ad69OjBhx9+SG5ubt2y5m6Ayo+JQXl54TdzZt17XiNG4DFgAPnRMS0uR1Ff+/btmThxIosXL+bzzz8nMzOTlMM/0L/HMHYmx6P1j2UYdiTH0aNDb9yrrLvoK4RRtemEH9TRm7sGdCFmbyaV1S2v09KommpI+BD63AWd+l+9LOxxKL1gOfpvBWtOrqHKXMWCgQuuen9Ut1H09evLipQVaK0xm81s3bqVxx57jG7durFw4UIyMjJ4+eWXSU9PZ9u2bfzkJz/B19fXqrg1BQUUffkVftOm4dK+fd37SikCIiOpSEmhrPZGKEfLSSrHxc1EVmnyVe8fSNtKSXkRU0YtaGRNIQS08YQP8MiYYC4UV7AxOdfxG0/9Ci6fhbBF1y/rPRE69rec33ewGnMNq1JXMbrbaPr69yU6OpqQkBBMJhO9e/em7+W+HLt4jKdffZo+ffpw9913s3btWiIiIti5cydpaWn88Y9/pHfv3jbHLli9Bl1RQUDU9RdI/aY9iMnXl/zl0Y7o5lUqSqtI3ZtL6Kiu/OHl3111g1hVTSX7T33NgG4jKc5vxdN3Qtzk2nzCn9C/M8EdvflkV4bjN77vffALgtD7rl9mMllO82QnQPZBh4bdnrWdnJIcFgxcUFfi4PTp02itOX36NG///G1qymqIPxPPgAEDiImJIScnh/fff59x48ZZf2/BNXRNDfkrV+I1cgSeDZzfN3l74z9rJkWbN1OVl2dvN6+SsjuX6kozQyf2JCoqiiVLlhAcbLnHwsXFhQceGQtKkbzzrEPjCtGWtPmEbzIpHgkPJuF0Psln7a/KWCfvOGTshFE/BZNLw21uiwA3H9i/1HFxgRUpK+jq3ZWJvSby0ksvUVpaetXy6tJqyvaX0fmOzsTExRAREWF1yYSmFO/cSdWZM3RoYiRMQEQEVFdT8Nlndse7Qps1iduz6NanPZ2DLKeeoqKiyMjI4Msvv7TcwOWtCRnSkeTvzlLTGqfvhGgD2nzCB5g7oheebiY+3X3acRvdvxRcPGD4o4238fSDYfMh8XMobaq2nPXSC9PZk7OHeQPm4WpyJTOz4QvSORtzqNJVrE1b65C4YLlY69q5M771Sipfyz0kBJ9x4yiIXYWuV5fIHmdSLlGYV8aQCdffeTx16lQGDx7Ma6+9xpAJgZQVVXLqkGN/XQjRVhgi4ft5u/HQbYHEHc6msNQBSai8CI6shCGzwaeZMf6jnoCaCjjYWOVo28SmxOJmcmN2/9ns3bu30dMz3dy6Mbr7aGJTY6k2V9sdt/L0aUp27MR//nxUvRLSDQmIiqQ6L4/Lds5idUXitmy8fN3od3uX65aZTCZ+85vfkJiYSHLmfvw6e5G4tcGqHUIYniESPlgu3pZXmfnswJnmGzfnyEqoLLaMxGlO10EQPA4SPrCUT7ZDSVUJ8afiuTfkXrZv2M7EiRPp2LEjnp6eV7Xz9vbm1VdfJWJABLkluWzP2m5XXID8FSvB1RX/eXObbdtu/HjcAgPJj46xO27RhTIyEi8waFwPXNwa/rhGRETQs2dP/vY3y1F+bnoh5zMv2x1biLbGMAl/cA8/RgYH8Ome05jNdowT19pSHC1whOVhjbAnoCATTm5ueVzgy1NfUlxVjD6omTNnDrfddhvJycksXbqU4OBglFIEBwezZMkSoqKimNBrAt18urEiZYVdcc1lZRSsWUP7e+/Brcv1R9nXUi4uBERGULp/P+WpJ+yKnbQjG6UUg+9svFSEu7s7zz33HNu2baPUIwdXdxOJ27PsiitEW2SYhA/w6NgQTl8sZfvJ8803bswP2+HCCcupGmsNfAB8u1tG9bSQ1poVKSvwKfbhtV+9xpw5c/j222/p3Llz3QVMs9lMRkZGXXkBV5Mr80LnsTdnL+kF6S2OXfjll5iLigiItL5Wjd+sWSgPD/JjWn6UX11Zw/Hvc+g9rBO+HTybbPvEE0/g5+fHG2//L6Fh3Ti57xzlJY65hiBEW2GohD9lcDc6tfOw7+LtvvfBuyMMntl82ytc3CzlFk5tgQtpLQq744cdnCo8RerKVH7729+ycuVKvJqoWXPFrP6zcDO5sTJ1ZYviaq3Jj47BY8AAvGyYEtA1IID2DzxAYXw8NUVFLYp9MiGP8pIqhk5svky0r68vTz31FKtXr8avdw3VVWaO78ppUVwh2ipDJXx3VxORo4PYmppH5sXS5le4VsEZSF1vKY7m1vQR53VuXwgmN8u5fBtlZWXx9AdPU11czf88+j/85S9/sXr+145eHZkSMoX4U/GUVJXYHLvs0CEqUlIIiIq0efx+QFQkuqyMwjjb7zbWWpO4LYuA7j4Ehvo3vwLw7LPP4u7uzpJP36F7Pz+Stmeh7Tl9J0QbY6iEDxAZFoRJKZbvbcFRfsKHlj+vFEezhW9XGDQDDkVDpfWJ99ChQ4TfE47uq5ncdTJPPvGkzaEXDFxASVUJ606ts3nd/OXRmHx98XvwQZvX9Ro8GK9hw8iPWYE22zY2/lxGEeczLzN0QqDVXzRXSkd8/PHHBA33pehCOaeTL9q830K0VYZL+N38PJkyuBux+89QVmnDqJmqcji4zFIUzd+GeW7rC3sCKgrh6Cqrmn/11VfceeedeI32wuRi4rf3/bZFYYd2GsrgjoPr6utYqyovj6LNm/GfNQtTC2/cCng4isqMDEp27bZpvaRt2bh5ujAgvJtN673wwgtUVlYS9+2nePu5k7RdhmgKcYXhEj5YhmgWllWx7ogNt+Efi4PSi5ak3VK9RkO3oZbrAM0k3n/84x9Mnz6d0FtC6X5fd8b3HE9P3+bPZTdEKcWCgQtIL0xnf+5+q9cr+OwzqK4mIKLlRcl877sPlw4dbLp4W1pUyckD5xgY3h13T9umbAgNDWXmzJn845/v0j+sE6eTL1J4vgWn74RogwyZ8Ef37sCArr4s251h/RHvvvctxdD6TGx5YKUso3vykiGz4SPempoann/+eZ5++mkeeOAB/vDJH8ivzL+uKqatpoRMwd/D3+ohmrqqioLYVfjceSfuISEtjmtyd8d/3lyKt26lMsu6o+3ju85irtYMnWj7rF0AixcvpqCggP3pX2NSSo7yhahlyISvlOKRMcEkny3iYGZB8ytkH7QUQbNhYvJGDZ1rKbnQQBXNkpIS5syZwxtvvMGzzz7L2rVrWfPDGoJ8gxjbw7qJyRvj6erJzP4z2XpmK7klzVcOvbxlC9V5eQRENjmZmVUC5s8HpSiIbX6kkLnGTNL2bHoODGjxxOSjR49mwoQJvP72a4QM68jxXTlU2XL6Tog2ypAJH2Dm8EB8PVz5ZHdG8433L7UUQRvmgHrr7t4w/BE4vg6Kfhw2mJuby8SJE4mPj+ett97irbfe4kTBCQ6fP8z8AfMxKfv/qeYPmI9Zm1mV2vw1hPzoGNx69qTd+PF2x3Xr3h3fSZMo+OxzzBVNTwiTkXiR4vwKq4ZiNmXx4sVkZWWRW3WcitJqTu4/Z9f2hGgLDJvwfTxcmT2iJ+sTczh/uYkkVHLRUvxs2ALLkbkjjPyppczCgY8BSE5OZvTo0Rw7doy4uDieffZZwDKFoZerFzP6zXBI2MB2gUzoOYHVJ1dTWVPZaLvy1BOU7t9PQMQClEsjlUBtFBAVaZk8Zf2GJtslbsuiXYAHIUPtm4d46tSpDBkyhNffe4UOgT4kbstyyExcQtzMDJvwwXLxtqpGE7u/iSkQD31qKX5mz8Xaa3XsC/0mw4GP2LJ5A2PHjqWyspIdO3YwbZplftrCikLWp6/n/t734+fhoC8aLFMgXiq/xObTjZd5yF8Rg/LwwG/WLIfF9R49Gve+fcmPbnxylEs5JWSl5DNkQiAmF/s+mkopfvOb35CUlITqWMCFM8XkprfsBjAh2gpDJ/y+ndtxZ/9ORO/NpLqmgXHi5hrY/wGE3AldbnFs8LBFUHyOD349k6CgIPbu3cuIeneyxqXFUV5TTsRA+8+h1xfeI5zg9sGsTGn4fHrN5csUxq+j/QMP4BrguDliLVMgRlCelETZ0aMNtknano3JVTHojh4OiRkREUGvXr14/7P/xd3LlcRtUl9HGJuhEz7Ao2NCyCks55vjDZzjPbkZCjOvnqDcAcxmMy99vJVTl8z8x6ROfPfddwQF/Ti236zNrExZye1dbmdAh+tnlrKHSZlYMGABR84f4djFY9ctL1wbhy4tbXAKQ3v5zXgIk49Pg0f5leXVpOzJof+Irnj5ujsknpubG8899xzfbvsG/95w6mAeJYWtM6m8EDcDwyf8uwd2IdDfi2W7Grjzdt8S8O1hKX7mIOXl5URGRvI/f/4LB12GM9TvMn5lV5ds/i77O7KKsxx+dH/F9H7T8XL1uu4oX5vN5MfE4DVsGF6DBzs8rks7H/xmzKBo/QaqL109IUzqnlyqymsY0sKhmI15/PHH8ff3J/67ZZhrNMe+kykQhXEZPuG7mBQPhwezO/0iJ8/Vq6F+IQ1OfQsjf2IpfuYAFy5cYPLkycTGxvKXv/yFOa/GgaunpdxyPStTVtLJqxOTgiY5JO612ru358E+D7L+h/UUlP84LLVk924qMzIIeLjxKQztFRAVaRnj/9nnde9prUncnk2XYF+6hrR3aLwrRdWiP/+IjsEeJO88S01Dp++EMADDJ3yA+aN64e5q4pP6VTT3L7UUO7t9oUNinDx5kjFjxpCQkMCqVav47W9/i/LuYBmXf3QVlFkS75miM3yX/R1zQ+fi5qAvmoYsGLiAipoK4tJ+LGyWHx2DS4cO+N7XwKTsDuLRty/e4eHkx65EV1tm4so+UUB+TglDJvRs8QTrTblSVC0h4xtKCir44fAFh8cQ4mYgCR/o4OPOtFt7sOZgFpfLq6CiGA7HWIqd+Xa1e/s7d+4kPDycgoICtm7dyty59WaNCnsCqkot8YDY1FhclAtzQufYHbcpoQGhjOg6gpWpK6kx11CZlU3xtm34z5uLyd0x59AbExAVSfXZHIq3bQMgaVsWnj5u9B/Z/OQqLdG1a1cee+wx/vHJX/H2dyNJJkcRBiUJv9ajY4IpqaxhzcFsSFxlKXIWtsju7cbExDB58mQ6d+7Mnj17GDNmzNUNug+z1NjZ/z5llSWsSVvDpOBJdPFuneRX34KBC8guzub7s9/X3QUbMH9+q8f1vesuXLt3Jz8mhuL8ctKPXOCWO7rj6u6YMf8NeeGFF6ioKCe38hjZJwq4mF3carGE+HflkISvlJqilEpVSqUppV5sYPljSqnzSqnDtQ/HDntxgGG9/BnWy59Pdv2A3ve+pchZr7AWb09rzSuvvEJUVBRjxoxh165d9O3bt+HGo56AS+ls2P8Glysvs2CAA+7otcKkoEl09upMbOJyCj77HN9Jk3Dr3r3V4ypXVwLmz6dk126OxCejtWbIeMderL1W//79mTVrFu9Ev4KLq9TXEcZkd8JXSrkA/wCmAoOACKXUoAaaxmqtb6t9LLU3bmt4NDyYjhcPoPKOWY7ubTifHB0dTUhICCaTieDgYCZMmMAf/vAHHnnkETZt2kSHDh0aX3nQDLRPF1akfUH/gP6M6Gr9zFL2cDO5MTd0Luqb76kpKGiVoZiN8Z87B+3uyfG95wkZ2on2nZqfvcteixcvJicvi0rvC6TszaWirLrVYwrx78QRR/hhQJrWOl1rXQmsBBxTC8DJHri1O497bqHE5AtDrD+HHh0dzaJFizh9+jRaazIzM9m5cyezZs1i2bJleHh4NL0BV3eODHmAFMpZ0OueVrlw2Zg5oXOYclBzuYcf3qNHOy2ua8eOFE16lAqzO4NH21dGwVphYWFMnDiRT756m+qKGlL3yBSIwlgckfADgfoDybNq37vWbKXUUaXU50qpXo1tTCm1SCmVoJRKOH/ejsnGW8CzLI9J7CWmcjzZpdYn3ZdeeonS0utrrh84cMDq5B3jXoOv2cyD55oo89AK2qXl0DdH88WwSsqqy5wa+0z7EXiVnqN96g6nxVy8eDEHjn2HqV0Fiduypb6OMBRnXbRdB4RorW8FvgaWNdZQa71Eaz1Saz2yc+fOTtq9Wgc+xqTNRNdMJsaGKRAzMxtO0o29f60LZRf4+ux3zHDrgveRFVDlvMSbHx2D9vZk8y2VrP9hvdPins+8TN65akLMJ8mPiXZa4p0yZQpDhgxh84FYCs6VkpWS75S4Qvw7cETCzwbqH7H3rH2vjtb6otb6yj3tSwHnnKS2RXUlHPgI1f8e+t9yKyv3naGi2roa6j17NlzKt365hKZ8fuJzqs3VzB/+CyjLh6TVVu+2PaovXaJo/XoCHppFULcBrExZ6bTEm7gtC1d3E0MeHEpl2ilK91k/E5c9lFIsXryYL3fEYHLXUl9HGIojEv5+oL9SqrdSyh1YAMTXb6CUqj/0Yzpw3AFxHStlHRSfg1FPsHBMCBdLKlmf2Pw5XrPZTPcGRrZ4e3vz6quvNrt+lbmKz058xtgeYwkZNBs632Ip6eCExFvw+Wp0VRUdoiKJGBhBan4qh/IOtXrc8pIqTuw/x4DR3ej00FRc/PyarKLpaAsWLKB7YDeSzn5HxtELFF107qksIW4UuxO+1roaeBrYhCWRr9JaJyul/kspNb222bNKqWSl1BHgWeAxe+M63L73ISAE+k3mjn4d6dPZp+H6OtdYvHgx+/btIzIykuDgYJRSBAcHs2TJEqKimi9RsDVzK3mleZahmErBqJ9BzhHISnBApxqna2rIX7kC7/BwPPr25f7e9xWJkNYAABy2SURBVOPr5ttoFU1HOv59DjVVZoZO7InJ0xO/ObO5vGULVbnNz8TlCG5ubjz//POs2PRPNJC8Q+rrCGNwyDl8rfV6rXWo1rqv1vrV2vf+qLWOr33+H1rrwVrrYVrru7TWKY6I6zC5iZY5Zkc9DiYTSikeDQ/m8JkCjmY1PgXiO++8w9///neeeeYZli9fTkZGBmazmYyMDKuSPVgmOenh04PxPWtnlhq2ANx9r6uv42jF27ZRfTanbiimt5s3D/V/iK9Pf8350ta7WG42a5J2ZNGjvz8dA9sBEBARAWYz+bGxrRb3Wo8//ji4VZFXfopj35+lukqmQBRtn9xpC5aje1cvuO3HJD1rRE+83V2urq9TT1xcHL/85S+ZMWMGb7zxRouGUp7MP8n+3P3MGzAPF1PtXaYevnBbBCSvheLWS7z50dG4du+O71131b03f8B8qnU1n5/8vIk17ZOZfJGiC+UMmfDjQC73nj1pN2ECBas+w1zZ+ExcjtSuXTueeuopYr/+F+XFVaQdyHNKXCFuJEn4ZfmQ+BkMnQPeP94c1d7TjVm3BxJ/5Cz5JVcnoT179hAREcGoUaOIiYnBpYXTAMamxuJucmdW/2tmlhr1BNRUwsFGBzPZpSI9nZJduwmYPx/l6lr3fnD7YO4IvIPPUj+jylzVKrETt2Xj7edOn+FXj8AKiIqi5uJFLm9qfCYuR3vmmWfIuJBMuS4kcZvceSvaPkn4h2MsxcsamMLw0TEhVFabiU348TaDU6dOMW3aNHr06MG6devw9vZuUdjLlZeJPxXPlN5TCPC8ZmapzqHQewIkfAg1jr8bND9mBcrNDf+5199cFjEggvNl5/k281uHxy3IKyUz+SKD7wzE5ZopDH3uGItbcBD5MTEOj9uYrl278pOf/IT1e6PJyyjiXIZMgSjaNmMnfLPZUga512hLEbNrhHb1JbxPB5bvOU2NWXPhwgWmTp2K1poNGzbQpUvLC5zFn4qnrLqMyIGNlDMIWwRF2XCi6Um/bVVTXEJhXBy+U6fg2vH6O1zHBY4jsF0gK1JWODQuQNKObEwmxeA7r5/CUJlMdIiMpOzQIcqPXT8TV2t54YUX2H18I2aqSZIhmqKNM3bCP/UtXEpvsirmo2NCyMovY9PRLKZPn05mZibx8fGEhoa2OKzWmpUpKxnaaSiDOzUys1ToFGjf0zJE04GK1sVjLi6mQ2TDXzQuJhfmD5jPgXMHOJF/wmFxqyprSNmVQ5/bO+Pj13CpCb+ZM1FeXlxy4lF+v379eGD6VPad/JoTCecoK3bONQQhbgRjJ/z974NPF7hleqNN7hnUlW7tPXjxg/Xs2bOH6Ohoxo4da1fYPTl7yCjKaHoKQxdXGPVT+GEHnE+1K94VWmvyY2LwHDwYz2HX/6K5Yma/mXi4eDh0iObJfeeoKK1m6MSGb1IDcGnfHr9p0yha9yU1BY2PjnK0xYsX883BzzBXa45/L/V1RNtl3ISfnwEnNsGIx8C18Qk/3FxMBFw6RpFPT37/17eZPXu23aFXpKwgwCOAe0Pubbrh7QvBxd0yisgBSvftp+JkGgGRkU2OKvL39Gdq76l8mf4lRZX2n9fWWnN0WxYdA9vRva9fk20DoiLRFRUUrFlrd1xrjRo1iluG9yHjwjESt2dhNkt9HdE2GTfh7/8AlMkyZ20T3nrrLTb/608obcZ14N12hz1bfJbtWduZHTobD5dmqmj6dILBs+DICii3P/Hmx8Tg4udH+wfub7ZtxMAIyqrLiE+Lb7Ztc3JPFXIxq5ihEwObHb7qOWAAXiNHkL9iBdrsvLlnFy9ezNcHYym+VMHpRJkCUbRNxkz4VWVw6FO45UFof/0FxCvWrl3Lc889x4x77+LBYYF8duAMpZX2jZpZlboKgHmh86xbIewJqCyGo/bdlFSVm8vlb77Bb85sTJ6ezbYf1HEQt3a+ldjUWMzavsSbuC0Ldy9XQsO6WdW+Q2QkVWfOULJzp11xbXHfffeh212muCKfo1vl4q1om4yZ8JNWW8bfj7p+KOYVu3fvJjIyktGjR7N8+XIeuyOEy+XVxB1q+W34FTUVrDm5hok9J9K9nZUzSwWOgB7DLad17KivU7BqFZjNlrtarRQxMIKMogz2nN3T4rglhRWcOnieW8Z2x83DuvsVfCdPxrVzZy45sb6OUorfLP41W4+sISsln/zcEqfFFsJZjJfwtbaMfOl8C4SMa7DJyZMnmTZtGoGBgcTHx+Pt7c3tQQEM6t6eT3ZntLii5KaMTeRX5LNgoA1TGCpl+WK6kGq5gNsCurKS/FWf0W7CBNwbqezZkHuD76WDZwdWpLZ8iOax785iNts2haFyd8d/3jxKdn5H5Wnry1Tba/78+Zy+fJQaXSNTIIo2yXgJPyvBUpws7PEGpzA8f/48999/P0opNmzYwJWa/EopFo4NJiX3MvszWlZDfWXKSkLahxDePdy2FYfMAq8OLa6vU7T5a2ouXCDAyvo+V7i7uDO7/2y2n9lOdrHtCbCmxkzyjmyCBnfAv6ttN6j5z5sHLi7kr2j9Ym5XuLm58fOnH+dA2laSvsuislymQBRti/ES/r4l4NEebr3+KLusrIzp06eTlZVFfHw8/fv3v2r59GGB+Hm5sWx3hs1hky4kkXghkQUDF9hed8fNC25/BFK+gkLbzy/nR0fjFhyEzx22DyedN2AeSiliU22/hvDD4QuUFFYydIL1vyqucOvaBd97JlOwZg3mMueVL/7Zz37G4cxvMVfBiX3nnBZXCGcwVsIvzoNjcTAsAjzaXbWopqaGqKgo9u7dS3R0NGPGjLludS93F+aN7MmmpFzOFZXbFHpFygq8Xb2Z0beF0/2O/JnldFTCRzatVn7sGGWHDtEhMhJlsv2fu5tPN+7udTdrT66lvNq2Piduy6J9J0+ChrRsztoOUVGYi4oo/PLLFq3fEu3atWP6/HvIPH+ChE2nZApE0aYYK+EfXGYpSjbq8esWvfDCC6xdu5Y33niDWbNmNbCyxcPhwdRoTcxe6+eezS/PZ+MPG5nWdxrt3Ns1v0JDAoItd98e+BiqK5ptfsWlmBiUlxd+M2e2LC6wYOACCioK2Jix0ep1LmYXc/ZkAYPHB2IytWxSdq8RI/AIDbVMw+jExPvMs8+w+8R6Si5Wc/ak824AE6K1GSfh11Rbjo77TLQUJ6vnzTff5K233uJXv/oVv/zlL5vcTHBHHyaGdiZmXyaV1dYNV1xzcg2V5krLJCf2CHsCSi/AsS+sal5TUEDRl1/hN20aLu3btzxstzD6+PVhRcoKqxNv4vZsXNxMDBrb+LDX5iilCIiKoiIlhbJDrT8T1xVdunRh6LggSsqL2L8xzWlxhWhtxkn4qestxciuqZuzevVqnn/+eWbPns3f//53qzb16JgQzl+uYFNy8zM01ZhrWJW6ilHdRtEvoF+Ldr1On7ugQ1+r6+sUrFmLLi+vm+SkpZRSLBi4gGMXj5F4IbHZ9hVl1aTuzaX/qK54tnOzK7bftAcx+fqSv9x5QzQBnnvhV+xJ3UjWsUKK863/RSXEvzPjJPz974NfL8tpkVq7du3i4YcfJjw8nE8//RSTlee4J4R2JqiDN582MjlKfTuydnC25GzTdXOsZTJZjvKz9sPZpo94tdlM/ooVeI0cgeeAAXaHnt53Oj5uPlbV10nZnUN1RQ23NlE3x1omb2/8Z82kaPNmqvKcN0lJv379aBdchdZwaEu60+IK0ZqMkfDzUixj2Ef+FGpnljpx4gTTp0+nZ8+exMfH4+XlZfXmTCbFI+HB7Mu4xPGcpkserEhZQRfvLtzV664m21ltWAS4ecO+pU02K9m5k6ozZxqtimkrHzcfpvWZxsaMjVwsu9hoO23WJG3Ppmvv9nQO8nVI7ICICKiupuCzzxyyPWs988LPOZa5jyNbM6mx8vSdEP/OjJHw9y+1FCG7/VEA8vLymDp1at1Y+06dOtm8ybkje+LpZmp0CkSAHwp/YHfObuaFzsPV5NpoO5t4+cOt8yHpcyi91GizSzExuHbujO/kyY6Ji+XO2ypzFWvTGi9slpWST8G50iarYtrKPSQEn3HjKIhdha5qnZm4GjJq1CgKXNJRNW6c2C9VNMXNr+0n/PIiS/GxIbPBpxOlpaVMnz6dnJwc1q1bR79+LTuv7u/tzoxhgcQdyqawrOEkFJsai6vJldmh9lfYvErYE1BdbqkH1IDKzExKduzEf948lHvjlUBt1ce/D6O7jSY2NZZqc8M3JR3dloWXrxv9bm/55DANCYiMpDovj8tbtjh0u81Z+PRc8gqz2LbmiFPjCtEa2n7CPxprKT426om6sfb79u0jJiaG8HAb73i9xiNjgimrquHzA9ffDFVaVcoXaV9wb/C9dPKy/RdEk7oOhuA7LBU/zTXXLc5fsRJcXCx3qzpYxMAIckty2Z61/bplRRfKOJ14gUHjeuDi5tiPVrsJ43ELDCQ/2nmTowDcd9+9nMpPwHzZk7zThU6NLYSjte2Er7Wl6FiP29GBt/Pcc88RFxfHW2+9xUMPPWT35ocE+jEiOIBPd2dcV0P9y/QvKa4qdszF2oaMehwKTsPJr69621xWRsHq1fjeMxm3ro49ygaY0GsC3Xy6NTgFYvJOS/mFwXdaXzfHWsrFhYCIBZTu3095quNm4mo2rlLctyCcyqpyNizf67S4QrSGtp3wf9hhKToW9gRvvPEG77zzDs8//zzPPPOMw0I8OiaYjIul7Ez7sYa61poVKSu4pcMtDOvc+MxSdrllGrTrdl19naKvvsJcVEQHG+vmWMvV5Mq80HnszdlLesGPo1eqq2o49l0OvW/rjG+H5ssvt4Tf7NkoDw/yVzj3KH9B1FyO5+yh4LSZ8hLnXUMQwtHadsLftwS8OrA61XIn7Zw5c/jb3/7m0BBTh3SnUzsPPtmVUfdewrkE0grSWlY3x1oubpbJW9K+gYunAMsXzaXoGDxCQ/EaMaJ14gKz+s/CzeTGytQfh2imJeRRXlLF0AmOP7q/wjUggPb3309h/DpqLl9utTjXcnNzY8jEnria3Nm4YrfT4grhaG034RecgdT1ZHWdRNTCn3LHHXfYNNbeWu6uJiLDevFtah5nLpUClqqY7d3bM7X3VIfGus6Ix8DkajmXD5QdOkzF8eMEREW13hcN0NGrI/eF3Ef8qXhKqix14xO3ZRHQzZvAAQGtFhcgICoKXVpK4dq4Vo1zrZ/9IorT549zcs8FtEyBKG5SbTfhH/gIDTz48hqCgoL44osv8LRipqeWiBwdjEkplu85zbmSc2zJ3MLMfjPxcrV+bH+L+HazTMB+aDlUlpAfHY3J1xe/aQ+2blws9XVKqkpYd2od534oIu/0ZYZO7NmqXzQAXkMG4znsVvJjYpw6BaKPjw9+fTXerv58v0lG7Iibk0MSvlJqilIqVSmVppR6sYHlHkqp2Nrle5VSIY6I25Do6GhC+wSTt/FvrEupJP1SNRs2bKBjx5ZVbLRGNz9P7hvcldiEM6xMWYVZm5k/YH6rxbtK2CKoKKR650cUbd6M/6yZmLxtqz3fErd2upVBHQexMmUliduycPN0YUC4dVMY2qtDVBSVGRmU7Hbu6ZWf/WoOhaUXWfnPTZhMJkJCQoh20qxc0dHRhISEOD3ujYwtfW6FuFprux6AC3AK6AO4A0eAQde0eQr4V+3zBUCsNdseMWKEtsXy5cu1t7e3jhrqpvWf2uvJfVy0p6enXr58uU3baYldaRd08G/j9Ojld+onv36y1ePVMZu1/udYnbdwuD42YKCu+OEHp4Vec2KNHvn+aP2Pp7bo7TEpTotbU1GhU8eM1ZlPPuW0mFpbPl8PjvqJfnvR17qjb3cNaG9v71b/fF35XAN1D2fEvZGxpc8tjwsk6EZyqtJ2lp1VSo0BXtZa31f7+j9qv0j+XK/Npto2u5VSrkAu0Fk3E3zkyJE6ISHB6n0JCQnh9OnTLH34f1C4UV6jr8RvtdM59VVU1wDVuKp2mLCvaJgtPMzluJeVYnZ3wdyt9Y/urzCjOYk7nUtCSLrldco9zzst9p3flzN6fyUn+zroDmYr1JhrqHHx43yP/8K16iwu1T+OzGrNE1lN/Sdp3RNoNy620fusdRmLlv8BgODgYDIyMqzejlLqgNZ6ZEPLHPG/JRA4U+91FjC6sTZa62qlVCHQEbhwTTuUUouARQBBQUE27UhmZibt3EGbumHGDdd6J6xqalr7YwIutSG0dsGZlVfKlDeV3u0xeQEVrd/P+jqb4HTnnWR6O2/uWYCvh2u6nwG/fOdNQ2j5D1mA2W0zxe2HganrNcuc70ZePpY+tx6TuaTueWam9XNvNKuxQ39rH8AcYGm9148A717TJgnoWe/1KaBTc9u29ZROcHBw3U8hVe9nUXBwsE3bEaIh9T9fOPHzdaPi3sjY0ueWx6WJUzqOuGibDfSq97pn7XsNtqk9peMHNF5ysYVeffVVvGsvWF75Jvb29ubVV191dChhQPU/X1c44/N1o+LeyNjS51aK29g3gbUPLKeF0oHe/HjRdvA1bX7B1RdtV1mzbVuP8LW2XPgIDg7WSikdHBzslIs8wjhu1OfrRn6upc83V59pzYu2AEqp+4E3sYzY+VBr/apS6r9qA8crpTyBT4HhwCVggda62VklbL1oK4QQRtfaF23RWq8H1l/z3h/rPS8H5joilhBCiJZpu3faCiGEuIokfCGEMAhJ+EIIYRCS8IUQwiAk4QshhEFIwhdCCIOQhC+EEAYhCV8IIQxCEr4QQhiEJHwhhDAISfhCCGEQkvCFEMIgJOELIYRBSMIXQgiDkIQvhBAGIQlfCCEMQhK+EEIYhCR8IYQwCEn4QghhEJLwhRDCICThCyGEQUjCF0IIg5CEL4QQBiEJXwghDEISvhBCGIQkfCGEMAhJ+EIIYRB2JXylVAel1NdKqZO1fwY00q5GKXW49hFvT0whhBAtY+8R/ovAFq11f2BL7euGlGmtb6t9TLczphBCiBawN+HPAJbVPl8GPGTn9oQQQrQSexN+V611Tu3zXKBrI+08lVIJSqk9SqkmvxSUUotq2yacP3/ezt0TQghxhWtzDZRS3wDdGlj0Uv0XWmutlNKNbCZYa52tlOoDfKuUStRan2qoodZ6CbAEYOTIkY1tTwghhI2aTfha68mNLVNKnVNKddda5yilugN5jWwju/bPdKXUNmA40GDCF0II0TrsPaUTDyysfb4Q+OLaBkqpAKWUR+3zTsAdwDE74wohhLCRvQn/L8A9SqmTwOTa1yilRiqllta2uQVIUEodAbYCf9FaS8IXQggna/aUTlO01heBSQ28nwA8Xvt8FzDUnjhCCCHsJ3faCiGEQUjCF0IIg5CEL4QQBiEJXwghDEISvhBCGIQkfCGEMAhJ+EIIYRCS8IUQwiAk4QshhEFIwhdCCIOQhC+EEAYhCV8IIQxCEr4QQhiEJHwhhDAISfhCCGEQkvCFEMIgJOELIYRBSMIXQgiDkIQvhBAGIQlfCCEMQhK+EEIYhCR8IYQwCEn4QghhEJLwhRDCICThCyGEQUjCF0IIg5CEL4QQBmFXwldKzVVKJSulzEqpkU20m6KUSlVKpSmlXrQnphBCiJax9wg/CZgF7GisgVLKBfgHMBUYBEQopQbZGVcIIYSNXO1ZWWt9HEAp1VSzMCBNa51e23YlMAM4Zk9sIYQQtnHGOfxA4Ey911m17wkhhHCiZo/wlVLfAN0aWPSS1voLR++QUmoRsAggKCjI0ZsXQgjDajbha60n2xkjG+hV73XP2vcai7cEWAIwcuRIbWdsIYQQtZxxSmc/0F8p1Vsp5Q4sAOKdEFcIIUQ99g7LnKmUygLGAF8ppTbVvt9DKbUeQGtdDTwNbAKOA6u01sn27bYQQghb2TtKZy2wtoH3zwL313u9HlhvTywhhBD2kTtthRDCICThCyGEQUjCF0IIg5CEL4QQBiEJXwghDEISvhBCGIQkfCGEMAhJ+EIIYRCS8IUQwiAk4QshhEFIwhdCCIOQhC+EEAYhCV8IIQxCEr4QQhiEJHwhhDAISfhCCGEQkvCFEMIgJOELIYRBSMIXQgiDkIQvhBAGIQlfCCEMQhK+EEIYhCR8IYQwCEn4QghhEJLwhRDCICThCyGEQUjCF0IIg5CEL4QQBmFXwldKzVVKJSulzEqpkU20y1BKJSqlDiulEuyJKYQQomVc7Vw/CZgFvGdF27u01hfsjCeEEKKF7Er4WuvjAEopx+yNEEKIVmPvEb61NLBZKaWB97TWSxprqJRaBCyqfVmslEptYcxOgNF+UUif2z6j9Rekz7YKbmxBswlfKfUN0K2BRS9prb+wcgfGaa2zlVJdgK+VUila6x0NNaz9Mmj0C8FaSqkErXWj1xXaIulz22e0/oL02ZGaTfha68n2BtFaZ9f+maeUWguEAQ0mfCGEEK2j1YdlKqV8lFK+V54D92K52CuEEMKJ7B2WOVMplQWMAb5SSm2qfb+HUmp9bbOuwHdKqSPAPuArrfVGe+Jaye7TQjch6XPbZ7T+gvTZYZTWujW2K4QQ4t+M3GkrhBAGIQlfCCEM4qZP+EqpKUqpVKVUmlLqxQaWeyilYmuX71VKhTh/Lx3Hiv4+r5Q6ppQ6qpTaopRqdEzuzaK5PtdrN1sppZsq83GzsKbPSql5tf/WyUqpGGfvo6NZ8dkOUkptVUodqv18338j9tNRlFIfKqXylFINDmJRFm/X/n0cVUrdbndQrfVN+wBcgFNAH8AdOAIMuqbNU8C/ap8vAGJv9H63cn/vArxrnz95M/fX2j7XtvPFMtR3DzDyRu+3E/6d+wOHgIDa111u9H47oc9LgCdrnw8CMm70ftvZ5/HA7UBSI8vvBzYACggH9tob82Y/wg8D0rTW6VrrSmAlMOOaNjOAZbXPPwcmqZu3FkSz/dVab9Val9a+3AP0dPI+Opo1/8YA/w38FSh35s61Emv6/ATwD611PljucXHyPjqaNX3WQPva537AWSfun8Npy82nl5poMgP4RFvsAfyVUt3tiXmzJ/xA4Ey911m17zXYRmtdDRQCHZ2yd45nTX/r+xmWI4SbWbN9rv2p20tr/ZUzd6wVWfPvHAqEKqW+V0rtUUpNcdretQ5r+vwy8HDtUPD1wDPO2bUbxtb/781yVi0d4WRKqYeBkcCEG70vrUkpZQJeBx67wbvibK5YTutMxPIrbodSaqjWuuCG7lXrigA+1lr/XSk1BvhUKTVEa22+0Tt2s7jZj/CzgV71Xvesfa/BNkopVyw/BS86Ze8cz5r+opSaDLwETNdaVzhp31pLc332BYYA25RSGVjOdcbf5Bdurfl3zgLitdZVWusfgBNYvgBuVtb0+WfAKgCt9W7AE0uRsbbKqv/vtrjZE/5+oL9SqrdSyh3LRdn4a9rEAwtrn88BvtW1V0RuQs32Vyk1HMv8BNPbwHldaKbPWutCrXUnrXWI1joEy3WL6Vrrm3miHWs+13FYju5RSnXCcoon3Zk76WDW9DkTmASglLoFS8I/79S9dK544NHa0TrhQKHWOseeDd7Up3S01tVKqaeBTViu8n+otU5WSv0XkKC1jgc+wPLTLw3LBZIFN26P7WNlf/8GtAM+q702nam1nn7DdtpOVva5TbGyz5uAe5VSx4Aa4Dda65v1l6u1fX4BeF8p9RyWC7iP3cQHbyilVmD50u5Ue13iT4AbgNb6X1iuU9wPpAGlwE/sjnkT/30JIYSwwc1+SkcIIYSVJOELIYRBSMIXQgiDkIQvhBAGIQlfCCEMQhK+EEIYhCR8IYQwiP8PyAR3IkVX1owAAAAASUVORK5CYII=\n",
                "text/plain": "<Figure size 432x288 with 1 Axes>"
              }
            ],
            "_view_module": "@jupyter-widgets/output",
            "_model_module_version": "1.0.0",
            "_view_count": null,
            "_view_module_version": "1.0.0",
            "layout": "IPY_MODEL_ff37659491e4482384481242650a91cf",
            "_model_module": "@jupyter-widgets/output"
          }
        },
        "663fdb6fa5b8424a9e5a8ce25a41da23": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "SliderStyleModel",
          "state": {
            "_view_name": "StyleView",
            "handle_color": null,
            "_model_name": "SliderStyleModel",
            "description_width": "",
            "_view_module": "@jupyter-widgets/base",
            "_model_module_version": "1.5.0",
            "_view_count": null,
            "_view_module_version": "1.2.0",
            "_model_module": "@jupyter-widgets/controls"
          }
        },
        "de608d5fc21a4866b725cd74a04a6bd6": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "state": {
            "_view_name": "LayoutView",
            "grid_template_rows": null,
            "right": null,
            "justify_content": null,
            "_view_module": "@jupyter-widgets/base",
            "overflow": null,
            "_model_module_version": "1.2.0",
            "_view_count": null,
            "flex_flow": null,
            "width": null,
            "min_width": null,
            "border": null,
            "align_items": null,
            "bottom": null,
            "_model_module": "@jupyter-widgets/base",
            "top": null,
            "grid_column": null,
            "overflow_y": null,
            "overflow_x": null,
            "grid_auto_flow": null,
            "grid_area": null,
            "grid_template_columns": null,
            "flex": null,
            "_model_name": "LayoutModel",
            "justify_items": null,
            "grid_row": null,
            "max_height": null,
            "align_content": null,
            "visibility": null,
            "align_self": null,
            "height": null,
            "min_height": null,
            "padding": null,
            "grid_auto_rows": null,
            "grid_gap": null,
            "max_width": null,
            "order": null,
            "_view_module_version": "1.2.0",
            "grid_template_areas": null,
            "object_position": null,
            "object_fit": null,
            "grid_auto_columns": null,
            "margin": null,
            "display": null,
            "left": null
          }
        },
        "ff37659491e4482384481242650a91cf": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "state": {
            "_view_name": "LayoutView",
            "grid_template_rows": null,
            "right": null,
            "justify_content": null,
            "_view_module": "@jupyter-widgets/base",
            "overflow": null,
            "_model_module_version": "1.2.0",
            "_view_count": null,
            "flex_flow": null,
            "width": null,
            "min_width": null,
            "border": null,
            "align_items": null,
            "bottom": null,
            "_model_module": "@jupyter-widgets/base",
            "top": null,
            "grid_column": null,
            "overflow_y": null,
            "overflow_x": null,
            "grid_auto_flow": null,
            "grid_area": null,
            "grid_template_columns": null,
            "flex": null,
            "_model_name": "LayoutModel",
            "justify_items": null,
            "grid_row": null,
            "max_height": null,
            "align_content": null,
            "visibility": null,
            "align_self": null,
            "height": null,
            "min_height": null,
            "padding": null,
            "grid_auto_rows": null,
            "grid_gap": null,
            "max_width": null,
            "order": null,
            "_view_module_version": "1.2.0",
            "grid_template_areas": null,
            "object_position": null,
            "object_fit": null,
            "grid_auto_columns": null,
            "margin": null,
            "display": null,
            "left": null
          }
        }
      }
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/chrisrichardson/FEMexamples/blob/main/01%20-%20FEM%20Basics%201D.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "KwmgNaUBQSvg"
      },
      "source": [
        "from numpy import *\n",
        "set_printoptions(suppress=True, formatter={'all':lambda x: '%5.2f'%x})\n",
        "import matplotlib.pyplot as plt\n",
        "from ipywidgets import interact\n",
        "%matplotlib inline"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qCh1pLRUQSvm"
      },
      "source": [
        "# Finite Element basics 1D\n",
        "Let us look at the one-dimensional heat equation. Probably, you will be familiar with the steady-state form:\n",
        "\n",
        "$$ {d\\over dx}\\left(k A {dT\\over dx}\\right) = 0 $$\n",
        "\n",
        "where $T$ is the temperature, $A$ is the cross-sectional area, and $k$ is the thermal conductivity.\n",
        "\n",
        "For a simple example, let us solve this equation in a rod of constant cross-sectional area $A$ which extends from $x=0$ to $x=1$, with a thermal conductivity $k$. We can fix the temperature $T=100$ at $x=0$ and $T=0$ at $x=1$. Obviously, this has a trivial, linear, solution, but it is instructive to see how to get a numerical solution with Finite Element, and we can introduce more complexity later."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BCC4aNeKQSvn"
      },
      "source": [
        "### Representation of a function by discrete variables\n",
        "\n",
        "First, let us consider how we can represent $T(x)$ computationally. Simply, we might choose to subdivide the domain $(0,1)$ into $n$ pieces (or \"elements\"), and represent $T(x)$ as having a continuous, piecewise linear values in each part. For example, suppose $T(x) = \\sin(5x)$, and $n=10$."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "RonMUqVbQSvo",
        "outputId": "6c5928f3-e301-4b60-ec9a-ce18ec40c6ff",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 299
        }
      },
      "source": [
        "# Create an array of n+1 values for position x, from 0 to 1 inclusive\n",
        "n = 10\n",
        "x = linspace(0.0, 1.0, n + 1)\n",
        "T = sin(5*x)\n",
        "\n",
        "print(x)\n",
        "print(T)\n",
        "plt.plot(x,T, marker='o');\n",
        "\n",
        "# Plot a more accurate representation\n",
        "n = 1000\n",
        "x = linspace(0.0, 1.0, n + 1)\n",
        "T = sin(5*x)\n",
        "plt.plot(x,T);\n"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[ 0.00  0.10  0.20  0.30  0.40  0.50  0.60  0.70  0.80  0.90  1.00]\n",
            "[ 0.00  0.48  0.84  1.00  0.91  0.60  0.14 -0.35 -0.76 -0.98 -0.96]\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3iT9f7G8fenSVoKtJRRVpnKFpRCGUc8ArIdgIgsEUQEF0cEReCox60o5+ceiAuULSKgqMgGUZAiIsuyV1lllNWZ5Pv7I4FTsYVC0z5N8nldV64mz0juh5G7zxZjDEoppYJXiNUBlFJKWUuLQCmlgpwWgVJKBTktAqWUCnJaBEopFeTsVge4EmXKlDHVqlWzOoZSSvmVtWvXHjXGRF843C+LoFq1asTHx1sdQyml/IqI7MluuG4aUkqpIKdFoJRSQU6LQCmlgpwWgVJKBTktAqWUCnI+KQIR+VREjojIxhzGi4i8LSLbReQPEWmUZVx/EdnmffT3RR5lrdnrEmkxZjHVR82jxZjFzF6XaHUkpdRF+GqNYALQ8SLjOwE1vY/BwAcAIlIKeAZoBjQFnhGRkj7KpCwwe10io2dtIDE5FQMkJqcyetYGLQOlCjGfnEdgjFkuItUuMkkX4HPjueb1KhGJEpEKQCtggTHmOICILMBTKFN9kUsVvLHzE8jMTKeyHKOyJFGe4xRzp3Fw3jeQXAlsDrCHQVgkRFaEiPIQVRWKRFodXamgVVAnlMUA+7K83u8dltPwvxGRwXjWJqhSpUr+pFRX5tgO2LUcEtfyScpSaoQlYhf3X6dxAssu8h5RVaF8A6gYC1e1ggoNweaX5zsq5Xf85n+aMWY8MB4gLi5O76ZjJWNg7y+weS5smw/HdwKQ7ijBYaqx0NWIPaYc+000B0xpzphwzhBOBnbqlytK6xolaFU1lGtLpGA/cwiO74BDGzyPP7+FxS9AWAm4qiXU7wY1O0BoUYsXWqnAVVBFkAhUzvK6kndYIp7NQ1mHLy2gTOpyJe+D3yfD+qlwYjfYwnBX/ye/le/F2O0xrD4VRdVSRTl4Mp0M1//WCMIdITzapiYiwtKEI7z/yxHeXmmICLNzQ83KtKrdiFYdHqZcZBE4e9SzdrFzCWydD1vmgqMY1L0NmtwHleJAxLo/A6UCkPjqVpXefQTfGmPqZzPuFmAIcDOeHcNvG2OaencWrwXOHUX0G9D43D6DnMTFxRm91lABSvwNfnkXNs0G44bqN+K6thdz0hvzxvJE9h1P5bpKJRjWrhYta0Uz5/cDjJ2fwIHkVCpGhTOiQ226xv5vi9/ptExWbj/G0oQjLE1I4tCpNADqVoikVe1oWtcuS6MqUdjFwJ6VsGEmbJwFGac9m46aPQD1u+umI6Uuk4isNcbE/W24L4pARKbi+c2+DHAYz5FADgBjzDgREeBdPDuCU4ABxph477z3Av/2vtVLxpjPLvV5WgQFZP9aWPw87FwKoRHQuD+uJoOZs9vG24u2sftYCvVjIhnWthY31SmLXMFv6sYYEg6fZsmfSSxNOMLaPSdwug0RRez8s2YZWtUuS6ta0ZQNy4Q/psOvH0HSn1CyOrR8Ahr00EJQKpfytQgKmhZBPkvaCoue82yvL1oaWjyKq1F/vk04w1uLtrEz6Sx1K0QyrG1N2tUrd0UFkJNTaZms3HaUpQlJLN16hMOn0gGoVyGS1nWiaVUrmkZpq7AtGwOH/oBSV0P7F6F2J91kpNQlaBGoS0s/A8tfg1/eA3s4XP8v3M0e5LttZ3hz4Ta2HzlD7XIRDGtXk/b1yhMSkr9fvMYYthw8zdKtnk1Ia/ecwOU2RHrXFvqU2Eizne9iP74VrmoNHV9hdmLkRTdLKRXMtAjUxW2eCz+MglOJ0LAv7jbPMn+3kzcXbiPh8Glqli3O0LY1ubl+hXwvgJycTM1k5faj5/ctHDmdjh0nI0r9RP/0qTjcKYxzdeGtjC5keLZMEu6w8Uq3BloGSqFFoHKSchy+exw2fgXl6mNu/i8LzlTnjYXb2HLwFFdFF2Nom5rcem1FbBYVQHaMMWw+eIqlCUksS0hi5949jLJNprttOQnuSjyROZj1pgYAMVHhrBx1k8WJlbKeFoH6u+0LYc4QOJuEaTmSJWXu4vXFO9mYeIpqpYsytG1NOl8XU6gKICcnUzO57rkfaRWyjpcdn1COE7zr6srbzm64sbFrzC1WR1TKcjkVgR5uEYxcmbDwWfjlXUx0HdZe/wEvrHWwfv/vVC4Vztju13J7bAx2m/9cnLZEuIOYqHCWJsfSIf01/mP/nKH2r2kesoVXwh+zOp5ShZoWQbA5dRBmDoC9v3Cg1t08euIOfp2TQkxUOK/e0YBujSrh8KMCyGpEh9qMnrWB05lFGeF8gJXu+rzo+JSprsdhawmo1cHqiEoVSv75P15dmV0r4MN/4jqwnjdKjOT6Pzqx/5Sbl26vz5LHW9GzSRW/LQGArrExvNKtATFR4QjwS/E29JEx7HGWxEzpCcvHei6PoZT6C10jCGCz1yWeP5TyoeLLeMz5MQdtFemfMoLTIVfzQpca9GhSmTC7zeqoPtM1NuYvRwjtPZZCvw/L8UTGe9y8+EU4vAm6vAehxSxMqVThokUQoM7dFyA9M5N/26cwyPkdi1yx/NsM5YFbG9K7aRWKOAKnAHJSpXRRJt7fkp7jQtnmrMYjmyYjx7ZDnxmey2ArpXTTUKAaOz8BMs8yzvEGg+zf8ZmzA4MyH8NWJJIBLaoHRQmcU7V0Mabe/w+mOG5naMho3Md2wiftISnB6mhKFQpaBAEqJfkIU0Nfok3Ibzyb2Y/nnP1xE8LBk2lWR7NE9TLFmDqoOatsjbnb/QzOjDRPGexdZXU0pSynRRCITh3kyyIvUlf28kDmMCa4/ncX0YpR4RYGs9ZV0cWZMqg5W0Ou5o6MZ8kIKwmfd4E/51kdTSlLaREEmuO7SP2wHeVNEgMyn2CB+3/njoQ7bIzoUNvCcNarUbY4Uwc1I1HKc/PZ/5BWui5Mv9tzZrVSQUqLIJAkJZD5cXvSz5zg6RIv0aVrr/OHUsZEhes1d7xqlI1g6qBmJBNBx2OPkVqhCXx1H6yfZnU0pSyhl5gIFEe34fq0E8kpmQyxP8Mb/+pD+RJFrE5VqCUcOk2fj1ZRPCSdH8p9QPj+ldD5HWh0t9XRlMoXOV1iQtcIAsGxHZgJt3E6NZP+rqd5csAdWgK5ULt8BJMHNeO0O4xORx4mtUpLmDsE1k60OppSBconRSAiHUUkQUS2i8iobMa/ISK/ex9bRSQ5yzhXlnFzfZEnqBzfhZl4G2dSUuiZPpqhvW+jfkwJq1P5jTrlI5l8XzNOOu10OvQAqVVawzdDPbfHVCpI5LkIRMQGvAd0AuoBvUWkXtZpjDHDjDENjTENgXeAWVlGp54bZ4zpnNc8QSV5H0zsTFrKaXqkjuLOTu1pV6+c1an8Tt0KkUy6rxnJThu3HL6ftIrN4Ov7IeF7q6MpVSB8sUbQFNhujNlpjMkApgFdLjJ9b2CqDz43uJ09Bl/cTsbZE3Q/O5LYpv9k4A3VrU7lt66pWIJJA5txLMPGbcf+RUZ0fZjR33O/ZqUCnC+KIAbYl+X1fu+wvxGRqkB1YHGWwUVEJF5EVolI15w+REQGe6eLT0pK8kFsP5ZxFqb0wH1iD/1Sh1OqRhOe63yNT+8dHIzqx3jK4HC6gy4nh5NZ8iqY2hv2rbE6mlL5qqB3FvcCZhpjXFmGVfXuxe4DvCkiV2c3ozFmvDEmzhgTFx0dXRBZCydXJszohznwG8Pcj3K0dGPe7dPIr68aWpg0qFSCSfc1Y39aEbqfHYmzaFmY2hOO7bA6mlL5xhffHolA5SyvK3mHZacXF2wWMsYken/uBJYCsT7IFJjcbpjzMGxfyH8dD/KTrSmf3dOEEuEOq5MFlGsrRTFpYDN2phSlX8ZIXMbApDvgTJCviaqA5YsiWAPUFJHqIhKK58v+b0f/iEgdoCTwS5ZhJUUkzPu8DNAC2OyDTIFp0bPwx3SmFu/PRyn/ZHy/xlQuVdTqVAHpuspRfD6wKX+klGaIGYU5fQim9PBsllMqwOS5CIwxTmAIMB/YAswwxmwSkedFJOtRQL2AaeavZ7DVBeJFZD2wBBhjjNEiyM5vn8PKt1hZsjOjj7ZnbPdraVy1lNWpAlpslZJMvLcpy1Oq8rR9GObg7zBzILicVkdTyqf0zGJ/sGsFfNGVvSXiaH3wYR5pW5ehbWtanSpoxO8+Tv9Pf2Vw+CKGpo+HJoPglv9aHUupy6ZnFvurYztgel9OF6vKrQcH0jm2Co+0qWF1qqASV60UE+5tyoepbZjh6AprPoI1n1gdSymf0SIozFKOw5QeZBJC1xOPULtaJcbc0UAPE7VAk2ql+OyeJjyX1pPV9jjM90941tSUCgBaBIWVKxO+7I854bmnQGZkVT68Oy6g7i/sb5pdVZpPBjRjSPrD7KUC7ul3w/FdVsdSKs+0CAqrH5+GXcsZG/YQa1y1+fSeJpQqFmp1qqDX/KrSvHXPjQzKfIyz6Zk4J/eCtFNWx1IqT7QICqM/ZsDqD/ghohvjTzbjg76NqVG2uNWplNf1V5fhmf63McQ5FDm2lYyZgzzneCjlp7QICptDGzBzH2FXsYYMSerKi13r06JGGatTqQu0qFGGQf3u5WVXP0K3/0Dq4tesjqTUFbNbHUBlkXIcpt1Fiq04dx67n4Eta9GraRWrU6kc3FCzDO6+TzJ38jZu/WkMj6x28M3pOlSMCmdEh9p6NzjlN3SNoLBwu2DWINynDnD36SE0vqYWIzvUsTqVuoQba5clockLbHVX4tmM16nAURKTUxk9awOz1+V0pRWlChctgsJi6RjYvpDnnf3JrBDHGz0bEhKih4n6g9mbTvJg5qM4cPF+6FuEkklqpoux8xOsjqZUrmgRFAbbFsLy1/g25CbmF+nEx/3jKBqqW+38xYHkVHaZCjyeeT8NQ3bwtP2L88OV8gdaBFY7dQAzazC7bNX4j3MAn9zTlHKRer9hf1IxKhyA+e6mjHPeyt32hXQN+en8cKUKOy0CK7mcmJkDyUg7y6DUhxnbuxn1KkZanUpdphEdahPu8JzoN9bZk9XuOrzk+ISnmuvlwZV/0CKw0vLXkL0/Mzr9Hvrc3I42dfV+w/6oa2wMr3RrQExUOG5sPGsfRgYOYn99DJzpVsdT6pJ0Q7RVdi7DLHuNr1w3Uqzp3QxoUc3qRCoPusbG/OVw0YmfnaD/nlEcmfUEZXu8ZWEypS5N1wiscOYIGTMGssNU4Meqj/PMbfX0QnIBplvv+5huu42ymyeQselbq+ModVFaBAXN7SZl2kDcacmMjRjJf/tej13vNxxwIoo4qHTnq2xwV8M560E4ud/qSErlyCffQCLSUUQSRGS7iIzKZvw9IpIkIr97H/dlGddfRLZ5H/19kacwS1n2FkX3L+f1kAE8dW8PIovoDsVA1aJODAuveQW3M4NTk+/RO5upQivPRSAiNuA9oBNQD+gtIvWymXS6Maah9/Gxd95SwDNAM6Ap8IyIlMxrpsIqI3E9octe5Ed3Ezr2H633Gw4Cg7u25/WwB4g8soaMpWOtjqNUtnyxRtAU2G6M2WmMyQCmAV1yOW8HYIEx5rgx5gSwAOjog0yFxux1ibQYs5g6o75m7/g+HDfFMbe9RSO933BQKBZmp32vR/ja1QL7irGwf63VkZT6G18UQQywL8vr/d5hF7pDRP4QkZkiUvky50VEBotIvIjEJyUl+SB2/pu9LpHRszaQmJzKSPtUash+RjofINUeZXU0VYCaX1WahNj/cNBEkTr9Xsg4a3Ukpf6ioPZSfgNUM8Zci+e3/omX+wbGmPHGmDhjTFx0dLTPA+aHsfMTSM10cWPIegbY5/OZswNLXNfqNWiC0NBbmzA2/FHCTu8h8/t/Wx1Hqb/wRREkApWzvK7kHXaeMeaYMebcmTUfA41zO68/O5CcSklO8V/Hh2x1xzDG2fv8cBVcwkNt9O3Vl4+dN+NYNwG2zrc6klLn+aII1gA1RaS6iIQCvYC5WScQkQpZXnYGtnifzwfai0hJ707i9t5hAaFiiSK87PiEKE7zaObDpOO51aRegyY4xVUrxfFmT7DFXYWMWQ/B2aNWR1IK8EERGGOcwBA8X+BbgBnGmE0i8ryIdPZO9oiIbBKR9cAjwD3eeY8DL+ApkzXA895hAeE/ldfRybaGsc6ebDbVAAh32BjRoba1wZRlHu3YgP8r/hikJZM5ewgYY3UkpRDjh/8Q4+LiTHx8vNUxLir96G6c7zZnM1cxNPQ5Dp7K0DtXKQDW7T3Bd+Of5En7ZOj8DjTqZ3UkFSREZK0xJu7C4XqtofxgDIe+GERpY3B2fpefGzeyOpEqRGKrlGRBiyH8/PM6mn43Cnv1llCyqtWxVBDTaxvkg0NLxlH15K98V/5B/qEloLIxtF1t3o8aTrrThVM3ESmLaRH4mOvYLkosf5ZVNKBN379dbUMpAMLsNkb2bM8rzruw71kO8Z9aHUkFMS0CX3K7OTJpEE4jnGz3OqUj9E5jKmcNKpWg1I33s8JVH+f8p+DEHqsjqSClReBDJ5Z/QIUTa/iy9AO0v76J1XGUHxhyU00+LjmcNKch8+uHwe22OpIKQloEPmKO76LosudZYa6jfd8Ren8BlSuh9hBG9mrHGOddOPaugLW6iUgVPC0CX3C7OTb5PtLdwoEbX6NSqWJWJ1J+pF7FSMq1up/lrgY4f3gaTuy2OpIKMloEPnD2p/cpcyyeCZH30711M6vjKD/0QOsaTCg9nDSXIWOWbiJSBUuLIK+O78Sx5HmWuhvS8a7HsIXoJiF1+Ry2EEb1bscYZ19C9/2E0aOIVAHSIsgLYzgx/SHS3CFsbfoitcpHWp1I+bFa5SKIafMAy10NcM1/GpL3Wh1JBQktgjxIW/M5JQ//wifhA+jf8Xqr46gAMOjGq/giejgZLhfpcx7VE81UgdAiuFKnD2HmP8lqdx1u7PUYYXab1YlUALDbQhjVuz2vu3oRtmsR5o8ZVkdSQUCL4AqdmPkoIc40Vl/zDI2rlbE6jgogV0cXp0K7f7HOXYOMeU/o5apVvtMiuAKZG+dScs/3fGrvwb1d21sdRwWge26owcToxwnJOE3qN09YHUcFOC2Cy5WaTPrcYWx2V6XOHf+meJhewFX5ni1EeLR3Zz50dyX8z68wW3+0OpIKYD4pAhHpKCIJIrJdRP52pTURGS4im703r18kIlWzjHOJyO/ex9wL5y1sTs4dTXj6MeZVH03repWsjqMCWLUyxYhsN5Jt7hhSZv0L0k9bHUkFqDwXgYjYgPeATkA9oLeI1LtgsnVAnPfm9TOB17KMSzXGNPQ+OlOIuXcup8SWKUySWxlwZzer46gg0LdFLT4v+xjhqYc58/0zVsdRAcoXawRNge3GmJ3GmAxgGtAl6wTGmCXGmBTvy1V4blLvXzJTOTPzYXa7yxF18zOUKR5mdSIVBEJChMF9ejOV9hT9/VPM3tVWR1IByBdFEAPsy/J6v3dYTgYC32d5XURE4kVklYh0zWkmERnsnS4+KSkpb4mvwOn5LxKZspep5YbTuUmNAv98FbwqlyqKrd2zHDSlODn9AXCmWx1JBZgC3VksIn2BOGBslsFVvffQ7AO8KSJXZzevMWa8MSbOGBMXHR1dAGmzfPaBdRSNf5+v3K3p27ufXllUFbieLeoyucwwos7u5OSPr1odRwUYXxRBIlA5y+tK3mF/ISJtgSeBzsaY87/SGGMSvT93AkuBWB9k8h2Xk5PTH+K4iSCl9XNULlXU6kQqCIkId919H/PMDRT79S3ch/+0OpIKIL4ogjVATRGpLiKhQC/gL0f/iEgs8CGeEjiSZXhJEQnzPi8DtAA2+yCTz5xd+QFRJzfzWeSD9Gl5rdVxVBCLiQrH2e4lzpgwjkx9SC8/oXwmz0VgjHECQ4D5wBZghjFmk4g8LyLnjgIaCxQHvrzgMNG6QLyIrAeWAGOMMYWnCE4mYlv6Esvc19G5z0N6ZVFluc4trmNW6cGUT17L0RV6hVLlG2L88LeKuLg4Ex8fn++fk/RJDyL2LmZS4xnc1/mmfP88pXLjUHIKB95sTU1JpOjwddgiCnafmfJfIrLWu0/2L/TM4hykbZxH9L75TArrxd03t7Q6jlLnlY8qyrHWrxLmTmHH5EetjqMCgBZBdjLOkjZ3OFvdMTTs+ZReWVQVOm1vbMmCqJ7UOvQtievmWx1H+Tktgmwc/uY5ojIOsazWk8RdXd7qOEr9jYjQpP/L7KMcrrmP0uqV76k+ah4txixm9rq/HbSn1EVpEVwgM/EPymz4iLkhbejVvYfVcZTKUdlSJVl09UiqmAN0PTsDAyQmpzJ61gYtA3VZtAiycrs5Ou0hkk0xStz2MhFFHFYnUuqiPkqsxhzX9Txom8tVcgCA1EwXY+cnWJxM+RMtgiyOLB1HhdMbmFfhYVrG1rE6jlKXdCA5jRcz+5JGKC/ZPwWMd3iqtcGUX9Ei8HKfOkSxFS/yK9fQqc8wq+MolSsVo8JJIooxzt78w7aZO0JWnB+uVG5pEXjtmfIodnc6x1u/SnRkEavjKJUrIzrUJtxhY5qrNfHuWjzpmER5+xlGdKhtdTTlR7QIgGPrv6f6oe+ZV6IXHW68weo4SuVa19gYXunWgIpRxXgq814iSOWl4l/SNfZiFwBW6q+C/j6LJiMF5zfD2W0q0KTvC3plUeV3usbGeL/4b2LNx3/SZv8E1q34hth/3mZ1NOUngn6NYNtXz1HOeYBNjZ6lctlSVsdRKk+uvetFDkg5Si4eSWpKyqVnUIogL4KTezdSPeFjloS1psOtes6A8n9h4RGcumkM1Uwia6borS1V7gRvERhD0rSHOWvCqNjjdey24P2jUIGlzj+7sT6yNc32fcaOP9dbHUf5gaD99kuY/yE1Un5n9dWPUPvqq6yOo5RPVb3rbTLEwemvHsHtclsdRxVyQVkEKcmHKbfqRTaE1KVlr8esjqOUz0WVq8KOBsNomPk7K+eMszqOKuSCsgi2fjGMYiYFc8vrFAnVy0iowHRd1+Fsd9Sm3voxJB05ZHUcVYj5pAhEpKOIJIjIdhEZlc34MBGZ7h2/WkSqZRk32js8QUQ6+CJPdmavS6TFmMX0HP1fGh6bx/cRd3Bt4+vz6+OUspzY7IR3e4cSnCZhsq75qpzluQhExAa8B3QC6gG9RaTeBZMNBE4YY2oAbwCveueth+cex9cAHYH3ve/nU7PXJTJ61gaOJJ/mRcen7HNH85+Tt+gVGlXAi6nbjA2VenPDyW9Z+9MPVsdRhZQv1giaAtuNMTuNMRnANKDLBdN0ASZ6n88E2ojnzK0uwDRjTLoxZhew3ft+PjV2fgKpmS4G276lZkgiTzvvITnToVdoVEHhmrvGcETKELVoBCmpejE69Xe+KIIYYF+W1/u9w7Kdxnuz+5NA6VzOC4CIDBaReBGJT0pKuqyA567EWEGO8a2rGUvdsX8ZrlQgCy0aSXKrl7ja7GXVlBesjqMKIb/ZWWyMGW+MiTPGxEVHX97Nus9difEp50CGZg7523ClAl2tlr3YFHkD/9j7EdsSNlkdRxUyviiCRKBylteVvMOynUZE7EAJ4Fgu582zc1doBHDh+RnusOkVGlVQqdTnXYyEkPzVUFx6boHKwhdFsAaoKSLVRSQUz87fuRdMMxfo733eHVhsjDHe4b28RxVVB2oCv/og01+cu0JjTFQ4AsREhfNKtwZ6hUYVVEqUr86Oax6hScYaVsz91Oo4qhARz/dxHt9E5GbgTcAGfGqMeUlEngfijTFzRaQI8AUQCxwHehljdnrnfRK4F3ACjxpjvr/U58XFxZn4+Pg851Yq2BhXJnvGNCM84zgyZA1lL3Mzq/JvIrLWGBP3t+G+KIKCpkWg1JU7uOknys24lWUlb6f1o59ZHUcVoJyKwG92FiulfKPCNTewKaY7N574ml9XLrQ6jioEtAiUCkK1+4wlOSSKyIUjOJuaZnUcZTEtAqWCUGjxkiTf+Dx1zE5WTnnF6jjKYloESgWpq1vdTUJEM67fO44/E/60Oo6ykBaBUsFKhAq938Mubo5/NQyX2/8OHFG+oUWgVBCLrFiTnfUe5vqMn1ky93Or4yiLaBEoFeTqdhvNfkdV6q17nkNHj1odR1lAi0CpICf2MByd36KiHOWPSf+2Oo6ygBaBUopyDVqzuXwXbjrxJb/8vMzqOKqAaREopQCoedfrnAkpTvEFj3MmLcPqOKoAaREopQBwRJThxA3P0MBsZdnUsVbHUQVIi0ApdV71mwayo3gjbtj9Llu2bbc6jiogWgRKqf8RoWyv9wiXDA5/+RhOvW9BUNAiUEr9RUSleuyucz+tMpay8NupVsdRBUCLQCn1NzXveJpD9hjq/fYcB44etzqOymdaBEqpvxFHOCG3vUkVOcxvk562Oo7KZ3kqAhEpJSILRGSb92fJbKZpKCK/iMgmEflDRHpmGTdBRHaJyO/eR8O85FFK+U7Z69qztdzNtD8xlZ9+WWl1HJWP8rpGMApYZIypCSzyvr5QCtDPGHMN0BF4U0SisowfYYxp6H38nsc8Sikfqn7Xm6SHhFP0xxGcScu0Oo7KJ3ktgi7ARO/ziUDXCycwxmw1xmzzPj8AHAH0RqlK+QFHZDmOX/8kjcwmFk59w+o4Kp/ktQjKGWMOep8fAspdbGIRaQqEAjuyDH7Ju8noDREJu8i8g0UkXkTik5KS8hhbKZVbVds8wJ5i19Jy91ts3rrN6jgqH1yyCERkoYhszObRJet0xhgD5HhBcxGpAHwBDDDGnDs4eTRQB2gClAJG5jS/MWa8MSbOGBMXHa0rFEoVmJAQSvceRzFJ4+jMYXpuQQC6ZBEYY9oaY+pn85gDHPZ+wZ/7oj+S3XuISCQwD3jSGLMqy3sfNB7pwGdAU18slFLKt4pXuoZd1wzhxowVLJnzmdVxlI/lddPQXKC/93l/YM6FE4hIKPA18LkxZuYF486ViODZv7Axj3mUUvmk1u1Pss9xFdetf4EDhw5ZHUf5UF6LYAzQTkS2AW29rxGROBH52DtND+BG4L6SfoAAABJRSURBVJ5sDhOdLCIbgA1AGeDFPOZRSuUTsYfiuON9SpPM9kmP4tkarAKB+ONfZlxcnImPj7c6hlJB6fdPH6Hh3on8+s/PaNqmm9Vx1GUQkbXGmLgLh+uZxUqpy1K/zyvsD6lIxRWjOH0q2eo4yge0CJRSl8VepBipHd+kEof56v8eovqoebQYs5jZ6xKtjqaukBaBUuqybXLUZ5KrLXfzHdfJdhKTUxk9a4OWgZ/SIlBKXbax8xMYk9mLw5TkVcd4QskkNdPF2PkJVkdTV0CLQCl12Q4kp3KGojyZOZDaIft5yD7n/HDlf7QIlFKXrWJUOABL3LF87WrBQ7Y51JU954cr/6JFoJS6bCM61CbcYQPgucx+nKQ4rzs+oF/TChYnU1dCi0Apddm6xsbwSrcGxESFc5IIXnM8SN2QvRRb9Tpn051Wx1OXyW51AKWUf+oaG0PX2Bjvq1s48vlWeu2YyVtTbmL4Pb3wXDlG+QNdI1BK+UTZO98gLawMnXe9wPRftlsdR10GLQKllG+ER1H0znHUDEnkzA/PsjHxpNWJVC5pESilfCakZhvSruvHvSHzGPfFJE7p7S39ghaBUsqnitz8MpnFYxiR+hZPz1itVyn1A1oESinfCosgrPuHVJXDNNr6JhN+3m11InUJWgRKKd+rdgOm2QP0ty9g+ffTWbf3hNWJ1EXkqQhEpJSILBCRbd6fJXOYzpXlpjRzswyvLiKrRWS7iEz33s1MKRUApO2zuErXZqx9HE9OXsaJsxlWR1I5yOsawShgkTGmJrDI+zo7qcaYht5H5yzDXwXeMMbUAE4AA/OYRylVWDjCsXX/mFIhZxme+g7Dp6/D7db9BYVRXougCzDR+3winvsO54r3PsU3AefuY3xZ8yul/ECFawlp+wxtQ+Ipv2M6Hy7faXUilY28FkE5Y8xB7/NDQLkcpisiIvEiskpEzn3ZlwaSjTHnzkffD8RkPzuIyGDve8QnJSXlMbZSqsA0fwhzVWueC53ErB8Xs3rnMasTqQtcsghEZKGIbMzm0SXrdMZzjFhO631VvffJ7AO8KSJXX25QY8x4Y0ycMSYuOjr6cmdXSlklJATp+gGOIsV4r8gHDJvyK0fPpFudSmVxySIwxrQ1xtTP5jEHOCwiFQC8P4/k8B6J3p87gaVALHAMiBKRc9c7qgTo7Y2UCkSRFZDO71DLvYMBGZMZOm0dLt1fUGjkddPQXKC/93l/YM6FE4hISREJ8z4vA7QANnvXIJYA3S82v1IqQNS9FRrfw30h3+LeuZy3F22zOpHyymsRjAHaicg2oK33NSISJyIfe6epC8SLyHo8X/xjjDGbveNGAsNFZDuefQaf5DGPUqow6/AylK7BuKLjmLw4nhXbdH9fYSD+ePp3XFyciY+PtzqGUupKHNqI+bgNv1GHB9yj+WZoK8qXKGJ1qqAgImu9+2v/Qs8sVkoVrPL1kU6v0tj5O32dX/Gvqb/hdLmtThXUtAiUUgWvUX9ocCePhMzEtncl//1xq9WJgpoWgVKq4InArW8ipa/mo6IfMHPZbyzactjqVEFLi0ApZY2w4nDnBIpzlo+Kf8jj039j3/EUq1MFJS0CpZR1ytdHOr1GrPN3BpmZDJnyGxlO3V9Q0LQIlFLWatQPruvNQzKT0geW8PJ3W6xOFHS0CJRS1hKBW9+ACtfxfvgHLPvlF77bcPDS8ymf0SJQSlnPEQ49JxEWVoTPi77JszNXs/voWatTBQ0tAqVU4RBVBblzApXcB3hZ3ufBSWtJy3RZnSooaBEopQqP6jci7V+gLatpnTSJ577ZZHWioKBFoJQqXJo/BA3uZIRjBkfjv2bWb/utThTwtAiUUoWLCNz2NlSM5Z2w95j09TdsO3za6lQBTYtAKVX4hBZFek/DUbwMH9he46kvFpCS4bz0fOqKaBEopQqniHLY7ppBGXs6T516jue/WoM/Xi3ZH2gRKKUKr/L1sfWYwDUhe2m9+Smm/brH6kQBSYtAKVW41WoPHV6mgy2ezHkj2ZSYbHWigJOnIhCRUiKyQES2eX+WzGaa1iLye5ZHmoh09Y6bICK7soxrmJc8SqnAFNL8AVIbP0C/kB/4ecK/OZWWaXWkgGK/9CQXNQpYZIwZIyKjvK9HZp3AGLMEaAie4gC2Az9mmWSEMWZmHnMopQKZCOG3vMKxE4cYtHMy775bkqnOmziQnEbFqHBGdKhN19gYq1P6rbxuGuoCTPQ+nwh0vcT03YHvjTF6rVml1OUJCaH0XR+zqVgzHjz9HvVPrcAAicmpjJ61gdnrEq1O6LfyWgTljDHnrg51CCh3iel7AVMvGPaSiPwhIm+ISFhOM4rIYBGJF5H4pCS94bVSQcnm4BHnUNabq3nb8S7NQzYDkJrpYuz8BIvD+a9LFoGILBSRjdk8umSdzniO68rx2C4RqQA0AOZnGTwaqAM0AUpxwWalC95/vDEmzhgTFx0dfanYSqkAtfMk3Jsxgj2mLJ84xhInfwJwIDnV4mT+65JFYIxpa4ypn81jDnDY+wV/7ov+yEXeqgfwtTHm/F4eY8xB45EOfAY0zdviKKUCXcWocJKJ4K6MJzlkSjEh9DUayVbKFM9xg4K6hLxuGpoL9Pc+7w/Muci0vblgs1CWEhE8+xc25jGPUirAjehQm3CHjSSi6J3xFEdMFBNDX6XS2U28s2gbLreedHa58loEY4B2IrINaOt9jYjEicjH5yYSkWpAZWDZBfNPFpENwAagDPBiHvMopQJc19gYXunWgJiocJIoybDwF6FYGaaEv8rChd/Re/wqEnUz0WURfzxlOy4uzsTHx1sdQylVWJzcDxNuIfPUEQY7H2et1Oflbg249dqKVicrVERkrTEm7sLhemaxUsr/lagEA37AUaoKn9pfpWfkRoZMWceIL9dzNl0vVncpWgRKqcAQWQEGfI+Ur8+/T7/E+w22MfO3/dzy9grW79PLUlyMFoFSKnAULQX95iDVWnDztmdY1nwdGZku7vjgZ95ful13JOdAi0ApFVjCIqDPl1D/DqqsG8vS2rO4+ZrSvPZDAn0/Xs3Bk7oj+UJaBEqpwOMoAnd8Ai1HErphMm9lvsAbnauyfn8yHd9cwQ8bD176PYKIFoFSKjCJQOt/w+3jkX2ruT3+bn7sU4aqpYvywKTfGD3rD73rmZcWgVIqsF3XE/p/AxkpVJp5K7Ou38uDra5m2pp93PrOT2xMPGl1QstpESilAl+V5nD/cohpjH3ug4zMHMfUe64jJd3F7e+v5KPlO3EH8Y5kLQKlVHCIKAf95kCLR2HtZzRf0I0fe0ZyU52yvPTdFvp/9itHTqVZndISWgRKqeBhs0O756DvLEg/TeTkjoyrvIhXutRlze7jdHxrBQs3H7Y6ZYHTIlBKBZ8abeChn6FeV2Tpy/Re348FdxalQoki3Pd5PE/N3kBqhsvqlAVGi0ApFZzCS0L3T6DH53D2GJVndWZu1Rk88o9STFq1l87v/sSWg6esTlkg9KJzSimVfhqWjoFVH0BYBLvqDObujQ05khrCqE51GNCiGp6r5Vtj9rpExs5P4EByap7u0ZzTRee0CJRS6pzDm2Dhs7DtR9zFyjGlSA9eSGxM81ox/PfO64iOKPib38xel8joWRtIzXThwEkmdsIdNl7p1uCyy0CLQCmlcmv3Slj0HOxbTVpoScantmVOaCfaNKrLvA2H8vSbudttOJPh5HSak1Opmf/7mZ7JqVQnp9MyOZXm/ZnqZOGWw1Ry7aOf7Uc6236hffqrJFGSmKhwVo666bI+O6cisF/WuyilVDCo1gLunQ97VlLk53d4ZOuXPOD6mgWrGrHN1ZJDXEticiojv/qDPw+dokFMFKfSMs9/eV/4Ze4Z5/l5Jt3JpX7/DrOHUKFIJu3tv/Gp/EiLsE2kGzvfupsTiudsaF/eozlPRSAidwLPAnWBpsaYbH9NF5GOwFuADfjYGHPuTmbVgWlAaWAtcLcxJiMvmZRSyidEoNoNnseRP/n6w+dpa5Zzi+1XjppIlrobstjVkMnLznCaoudnCxGIKOIgooidSO/PyqWKnn8dWcROZHjW8Q4iw+1EFHEQlXGIiAM/Y98xH7b+CGnpHLRF81pmT6a5WnOcyPOfUzEq3GeLmtc1go1AN+DDnCYQERvwHtAO2A+sEZG5xpjNwKvAG8aYaSIyDhgIfJDHTEop5Vtl6zDq7F3Y6UmrkN+51baKtiFr6W5bjssIztK1cVZsjK1SI8LK1kBKVofIGM95C9nJOAunD8GJrXBoo2ffRGI8HN/pGV+8HDS+B+rfweqjFfns602kuv53OGu4w8aIDrV9tnh5KgJjzBbgUnvTmwLbjTE7vdNOA7qIyBbgJqCPd7qJeNYutAiUUoVOxahwEpNTWeCOY4E7DhsuYmUbnYomMLDUMcJ2fAcbJ2eZQzyXxA4tDqFFwe0EZwZknIH0Cw5LjYyBCtdB08FQvSWUretZIwG6VgEkxCdHDeWkIPYRxAD7srzeDzTDszko2RjjzDI8xyUTkcHAYIAqVarkT1KllMrBiA61zx+9A+DCxib7NfS9tRfExoAxcHIfnNjteZzc7zksNf0UZKSAzeF5OIpBRHmIrOi5xWbZep4b6lxE19gYn37xX+iSRSAiC4Hy2Yx60hgzx/eRsmeMGQ+MB89RQwX1uUopBZz/Is7xN3MRiKrieVS/0cKkl++SRWCMaZvHz0gEKmd5Xck77BgQJSJ271rBueFKKVUo5fdv5lYpiEtMrAFqikh1EQkFegFzjecEhiVAd+90/YECW8NQSinlkaciEJHbRWQ/8A9gnojM9w6vKCLfAXh/2x8CzAe2ADOMMZu8bzESGC4i2/HsM/gkL3mUUkpdPj2zWCmlgkROZxbr1UeVUirIaREopVSQ0yJQSqkgp0WglFJBzi93FotIErDnCmcvAxz1YRx/oMscHHSZA19el7eqMSb6woF+WQR5ISLx2e01D2S6zMFBlznw5dfy6qYhpZQKcloESikV5IKxCMZbHcACuszBQZc58OXL8gbdPgKllFJ/FYxrBEoppbLQIlBKqSAXsEUgIh1FJEFEtovIqGzGh4nIdO/41SJSreBT+lYulnm4iGwWkT9EZJGIVLUipy9dapmzTHeHiBgR8etDDXOzvCLSw/v3vElEphR0Rl/Lxb/rKiKyRETWef9t32xFTl8SkU9F5IiIbMxhvIjI294/kz9EpFGePtAYE3APwAbsAK4CQoH1QL0LpnkIGOd93guYbnXuAljm1kBR7/MHg2GZvdNFAMuBVUCc1bnz+e+4JrAOKOl9Xdbq3AWwzOOBB73P6wG7rc7tg+W+EWgEbMxh/M3A94AAzYHVefm8QF0jaApsN8bsNMZkANOALhdM0wWY6H0+E2gj4r1btH+65DIbY5YYY1K8L1fhuSucP8vN3zPAC8CrQFpBhssHuVneQcB7xpgTAMaYIwWc0ddys8wGiPQ+LwEcKMB8+cIYsxw4fpFJugCfG49VeO72WOFKPy9QiyAG2Jfl9X7vsGynMZ6b55zEc3Mcf5WbZc5qIJ7fKPzZJZfZu8pc2RgzryCD5ZPc/B3XAmqJyEoRWSUiHQssXf7IzTI/C/T13iTrO+BfBRPNUpf7//2iLnnPYhV4RKQvEAe0tDpLfhKREOB14B6LoxQkO57NQ63wrPEtF5EGxphkS1Plr97ABGPM/4nIP4AvRKS+McZtdTB/EahrBIlA5SyvK3mHZTuNiNjxrFIeK5B0+SM3y4yItAWeBDobY9ILKFt+udQyRwD1gaUishvPttS5frzDODd/x/vx3BM80xizC9iKpxj8VW6WeSAwA8AY8wtQBM/F2QJZrv6/51agFsEaoKaIVBeRUDw7g+deMM1coL/3eXdgsfHuhfFTl1xmEYkFPsRTAv6+7RgusczGmJPGmDLGmGrGmGp49ot0Nsb4631Oc/PvejaetQFEpAyeTUU7CzKkj+VmmfcCbQBEpC6eIkgq0JQFby7Qz3v0UHPgpDHm4JW+WUBuGjLGOEVkCDAfz1EHnxpjNonI80C8MWYu8AmeVcjteHbK9LIucd7lcpnHAsWBL737xfcaYzpbFjqPcrnMASOXyzsfaC8imwEXMMIY47drurlc5seAj0RkGJ4dx/f4+S91iMhUPIVexrvv4xnAAWCMGYdnX8jNwHYgBRiQp8/z8z8vpZRSeRSom4aUUkrlkhaBUkoFOS0CpZQKcloESikV5LQIlFIqyGkRKKVUkNMiUEqpIPf/Sl7aEd0mLTEAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gtnNE48UQSvu"
      },
      "source": [
        "We can build up the \"piecewise linear\" version of $T(x)$ from the following formula:\n",
        "\n",
        "$$ T_{piecewise}(x) = \\sum_i T_i N_i(x) $$ where $T_i$ are the discrete values of $T$ at each \"node\" and $N_i$ are functions that are \"local\" to the node. The functions $N_i(x)$ are called shape functions, and in this case, are just linear between neighbouring nodes: a triangular or tent-like shape. $N_i$ has value 1 at node $i$ and value 0 at all other nodes. The interactive plot below lets you build up $T$ by adding together the contributions from shape functions from 1 to 11"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "ea96GzElQSvv",
        "outputId": "dde18530-560d-4378-b690-63bf2a9ad93a",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 301,
          "referenced_widgets": [
            "dfdaf8ee303f43bbbda0a76f09396220",
            "cc3c2018633a4007b36230fe78352f4b",
            "3ccbd2cf9c634c31936efc3c95c0e8cd",
            "a50d77c2c715494a90dfc9a6e491ddfc",
            "663fdb6fa5b8424a9e5a8ce25a41da23",
            "de608d5fc21a4866b725cd74a04a6bd6",
            "ff37659491e4482384481242650a91cf"
          ]
        }
      },
      "source": [
        "def plot(n):\n",
        "    x = linspace(0, 1.0, 11)\n",
        "    T = sin(5*x)\n",
        "    # Plot T as given by the first 'n' shape functions\n",
        "    T[n:] = 0.0\n",
        "    plt.plot(x, T, marker='o', color='k')\n",
        "    \n",
        "    # Plot shape function n \n",
        "    for i in range(1, n+1):\n",
        "        Ni = zeros(11)\n",
        "        Ni[i-1] = 1.0\n",
        "        plt.plot(x, Ni)\n",
        "        plt.text(0.1*(i-1), 1, '$N_{%d}$'%(i-1), size=16)\n",
        "    plt.ylim(-1.5, 1.5)\n",
        "\n",
        "interact(plot, n=(1,11,1));"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "application/vnd.jupyter.widget-view+json": {
              "model_id": "dfdaf8ee303f43bbbda0a76f09396220",
              "version_minor": 0,
              "version_major": 2
            },
            "text/plain": [
              "interactive(children=(IntSlider(value=6, description='n', max=11, min=1), Output()), _dom_classes=('widget-int…"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4YtnOTG-QSv0"
      },
      "source": [
        "We can also calculate the derivatives, $$ {d\\over dx}T_{piecewise}(x) = \\sum_i T_i {dN_i\\over dx} $$ \n",
        "\n",
        "For these simple shape functions $N_i$, the derivatives are piecewise constant (and hence the second derivative is zero)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "a3a3pAUCQSv1",
        "outputId": "e4ad3f0b-d869-49c9-cc41-6b1296d91215"
      },
      "source": [
        "# For example, dN5/dx looks like this:\n",
        "x = [0.0, 0.4, 0.4, 0.5, 0.5, 0.6, 0.6, 1.0]\n",
        "y = [0.0, 0.0, 0.1, 0.1, -0.1, -0.1, 0.0, 0.0]\n",
        "plt.plot(x, y)\n",
        "plt.title('$dN_5/dx$')\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAGa9JREFUeJzt3X+w3XV95/HnSwJ0uyvIj2DTBBoc4rQRa6x3QaeVskUwskqcFSpINezConZZZ9fdncI4gou6lbq7OI6uEiQWnEWg2DaxxmYCSOt0lOYiUQHLEKjChaxcCSKjRY2+94/zjR6u5+Z+knPuvSR5PmbOnO/38/l8P9/Ph4T7up/v95x8U1VIkjST58z3ACRJewcDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTCk3ZDk0iQf2VfOI+0OA0PaPcuBryc5LEkl+VJ/ZZKrklw5UydJfjXJxEznGXKs0kgZGNLueRHwNWAF8P+A5UkW9dWvALY09HM68NcN55GeNQwMaRpJnpPkkiQPJXk0ydnAccDd9IJhHNgEnNG1PwB4MXBXQ/enAxtmOk+SP0nyF31j+mCSW5McOMq5Si0WzPcApGexS4FTgVcC36X3A35bVT2V5KX0VhL3AecCVwG/DhwAfGNXnXY/7E8CVjec5wrggSQrgJcDK4Hfqaofj3KiUgsDQxogyULgvwIvqapvdWWfA07omqwA1gG3AR9P8tyu7O6q+nGSpcBm4J6u/VlVNdltnwR8tQuEXZ6nqh5P8iHgOuBQemHx5OzNXJqegSENdgrwjap6oK/s+fRueB8M/AawpaqeSPL3wGuAnauOnf6mqs4c0PfPLkft6jx9+3cBlwHnVtXDw0xKGob3MKTBjgQe27nTXUZ6Pb0b0ccD/wQ82FX/ZVf3Up55/+K3k3wxyf9Ikr7y04HPNZyHJC8GPgZcC/y7UU1O2hMGhjTYfcDvJHlhkkPp/dA+ht5v/i+ld0lp58Nk1tELgf4VxjZ6N65PAo4C/g1AkmOBg6vqH2Y6T5LFwGeBtwF/CLw4ycmzN2Vp1wwMaYCq2gTcQO+TUJuBSeBp4H6mfHS2u/fwTeB5wFe7sh9W1fe7UPlz4CVd83/Nzy9H7eo827p2/7uq1lfVD4APAu+fnRlLM4tP3JNGL8khVfW9bvuP6d2nuC7JBuAjVbVh1z1Izz6uMKTZ8btJ7kzyRWAxcH1XfjvwhXkblTQEVxiSpCauMCRJTQwMSVITA0OS1GSf+qb3kUceWUuXLp3vYUjSXuXOO+/8TlUtnKndPhUYS5cuZXx8fL6HIUl7lSTfamnnJSlJUhMDQ5LUxMCQJDUxMCRJTQwMSVKTkQRGkpVJ7kuyNcnFA+pPSvKVJDuSnDmlbnWS+7vX6r7ylyX5etfnh6c8T0CSNMeGDozuwfcfpffEseXAOUmWT2n2EHAeP/8H2HYeezi9J4mdSO+RlJclOayr/hhwIbCse60cdqySpD03iu9hnABsraoHAZLcAKwC7t3ZoKq+2dX9dMqxrwY2VdX2rn4TsDLJ7cAhVfWlrvw6ek8h+/wIxqv91PV3PMS6LY/My7lXrVjMm048Zl7OLY3KKC5JLQb6nzM80ZUNc+zibnvGPpNcmGQ8yfjk5GTzoLX/WbflEe7d9r05P++92743b0EljdIoVhiD7i20/pvp0x3b3GdVrQHWAIyNjflvtWuXli86hBvf+oo5Pecbr/rSnJ5Pmi2jWGFMAEf37S8BHh3y2Ilue0/6lCTNglEExmZgWZJjkxwEnA2sbzx2I3BaksO6m92nARurahvwVJKXd5+OeguwbgRjlSTtoaEDo6p2ABfR++H/DeCmqronyeVJzgBI8i+TTABnAVcluac7djvwXnqhsxm4fOcNcODtwCeArcADeMNbkubVSP612u6B9humlF3at72ZZ15i6m+3Flg7oHwcOH4U45MkDc9vekuSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqMpLASLIyyX1Jtia5eED9wUlu7OrvSLK0Kz83yZa+10+TrOjqbu/63Fl31CjGKknaM0MHRpIDgI8CrwGWA+ckWT6l2fnAE1V1HHAlcAVAVf3fqlpRVSuANwPfrKotfcedu7O+qh4bdqySpD03ihXGCcDWqnqwqn4E3ACsmtJmFXBtt30zcEqSTGlzDvDpEYxHkjQLRhEYi4GH+/YnurKBbapqB/AkcMSUNm/kFwPjk93lqHcPCBgAklyYZDzJ+OTk5J7OQZI0g1EExqAf5LU7bZKcCPygqu7uqz+3ql4MvLJ7vXnQyatqTVWNVdXYwoULd2/kkqRmowiMCeDovv0lwKPTtUmyADgU2N5XfzZTVhdV9Uj3/hRwPb1LX5KkeTKKwNgMLEtybJKD6P3wXz+lzXpgdbd9JnBbVRVAkucAZ9G790FXtiDJkd32gcBrgbuRJM2bBcN2UFU7klwEbAQOANZW1T1JLgfGq2o9cA3wqSRb6a0szu7r4iRgoqoe7Cs7GNjYhcUBwC3A1cOOVZK054YODICq2gBsmFJ2ad/20/RWEYOOvR14+ZSy7wMvG8XYJEmj4Te9JUlNDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNTEwJElNDAxJUhMDQ5LUxMCQJDUZSWAkWZnkviRbk1w8oP7gJDd29XckWdqVL03yT0m2dK+P9x3zsiRf7475cJKMYqySpD0zdGAkOQD4KPAaYDlwTpLlU5qdDzxRVccBVwJX9NU9UFUrutfb+so/BlwILOteK4cdqyRpz41ihXECsLWqHqyqHwE3AKumtFkFXNtt3wycsqsVQ5JFwCFV9aWqKuA64PUjGKskaQ+NIjAWAw/37U90ZQPbVNUO4EngiK7u2CR3JfmbJK/saz8xQ58AJLkwyXiS8cnJyeFmIkma1igCY9BKoRrbbAOOqaqXAu8Erk9ySGOfvcKqNVU1VlVjCxcu3I1hS5J2xygCYwI4um9/CfDodG2SLAAOBbZX1Q+r6nGAqroTeAB4Ydd+yQx9SpLm0CgCYzOwLMmxSQ4CzgbWT2mzHljdbZ8J3FZVlWRhd9OcJC+gd3P7waraBjyV5OXdvY63AOtGMFZJ0h5aMGwHVbUjyUXARuAAYG1V3ZPkcmC8qtYD1wCfSrIV2E4vVABOAi5PsgP4CfC2qtre1b0d+FPgnwGf716SpHkydGAAVNUGYMOUskv7tp8Gzhpw3GeAz0zT5zhw/CjGJ0kant/0liQ1MTAkSU0MDElSEwNDktTEwJAkNTEwJElNDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSk5EERpKVSe5LsjXJxQPqD05yY1d/R5KlXfmpSe5M8vXu/ff6jrm963NL9zpqFGOVJO2ZoR+g1D1i9aPAqfSexb05yfqqurev2fnAE1V1XJKzgSuANwLfAV5XVY8mOZ7eU/sW9x13bvcgJUnSPBvFCuMEYGtVPVhVPwJuAFZNabMKuLbbvhk4JUmq6q6qerQrvwf4pSQHj2BMkqQRG0VgLAYe7tuf4JmrhGe0qaodwJPAEVPavAG4q6p+2Ff2ye5y1LuTZNDJk1yYZDzJ+OTk5DDzkCTtwigCY9AP8tqdNkleRO8y1Vv76s+tqhcDr+xebx508qpaU1VjVTW2cOHC3Rq4JKndKAJjAji6b38J8Oh0bZIsAA4Ftnf7S4C/AN5SVQ/sPKCqHunenwKup3fpS5I0T0YRGJuBZUmOTXIQcDawfkqb9cDqbvtM4LaqqiTPAz4HXFJVf7ezcZIFSY7stg8EXgvcPYKxSpL20NCB0d2TuIjeJ5y+AdxUVfckuTzJGV2za4AjkmwF3gns/OjtRcBxwLunfHz2YGBjkq8BW4BHgKuHHaskac8N/bFagKraAGyYUnZp3/bTwFkDjnsf8L5pun3ZKMYmSRoNv+ktSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqclIAiPJyiT3Jdma5OIB9QcnubGrvyPJ0r66S7ry+5K8urVPSdLcGjowkhwAfBR4DbAcOCfJ8inNzgeeqKrjgCuBK7pjl9N7BviLgJXA/0lyQGOfkqQ5NIpHtJ4AbK2qBwGS3ACsAu7ta7MKeE+3fTPwkSTpym+oqh8C/9g98/uErt1MfY7Mf//sPdz76Pdmo2s9i9y77XssX3TIfA9jTl1/x0Os2/LIfA9Dc2D5rx7CZa970ayeYxSXpBYDD/ftT3RlA9tU1Q7gSeCIXRzb0icASS5MMp5kfHJycohpaF+3fNEhrFox8K/RPmvdlke4d5u/DGk0RrHCyICyamwzXfmgIJvaZ6+wag2wBmBsbGxgm5nMdipL82n5okO48a2vmO9haB8wihXGBHB03/4S4NHp2iRZABwKbN/FsS19SpLm0CgCYzOwLMmxSQ6idxN7/ZQ264HV3faZwG1VVV352d2nqI4FlgF/39inJGkODX1Jqqp2JLkI2AgcAKytqnuSXA6MV9V64BrgU91N7e30AoCu3U30bmbvAP5DVf0EYFCfw45VkrTnRnEPg6raAGyYUnZp3/bTwFnTHPt+4P0tfUqS5o/f9JYkNTEwJElNDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNTEwJElNDAxJUhMDQ5LUZKjASHJ4kk1J7u/eD5um3equzf1JVndlv5zkc0n+Ick9ST7Q1/68JJNJtnSvC4YZpyRpeMOuMC4Gbq2qZcCt3f4zJDkcuAw4ETgBuKwvWP5nVf068FLgt5O8pu/QG6tqRff6xJDjlCQNadjAWAVc221fC7x+QJtXA5uqantVPQFsAlZW1Q+q6gsAVfUj4CvAkiHHI0maJcMGxvOrahtA937UgDaLgYf79ie6sp9J8jzgdfRWKTu9IcnXktyc5OjpBpDkwiTjScYnJyf3dB6SpBnMGBhJbkly94DXqsZzZEBZ9fW/APg08OGqerAr/iywtKp+E7iFn69ifrGjqjVVNVZVYwsXLmwckiRpdy2YqUFVvWq6uiTfTrKoqrYlWQQ8NqDZBHBy3/4S4Pa+/TXA/VX1ob5zPt5XfzVwxUzjlCTNrmEvSa0HVnfbq4F1A9psBE5Lclh3s/u0rowk7wMOBf5T/wFd+Ox0BvCNIccpSRrSsIHxAeDUJPcDp3b7JBlL8gmAqtoOvBfY3L0ur6rtSZYA7wKWA1+Z8vHZd3Qftf0q8A7gvCHHKUka0oyXpHalu3R0yoDyceCCvv21wNopbSYYfH+DqroEuGSYsUmSRstvekuSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqclQgZHk8CSbktzfvR82TbvVXZv7k6zuK789yX3dw5O2JDmqKz84yY1Jtia5I8nSYcYpSRresCuMi4Fbq2oZcGu3/wxJDgcuA04ETgAumxIs51bViu6185ng5wNPVNVxwJX4TG9JmnfDBsYq4Npu+1rg9QPavBrYVFXbq+oJYBOwcjf6vRk4JcnAp/NJkubGsIHx/KraBtC9HzWgzWLg4b79ia5sp092l6Pe3RcKPzumqnYATwJHDBpAkguTjCcZn5ycHG42kqRpzfhM7yS3AL8yoOpdjecYtDKo7v3cqnokyXOBzwBvBq6b4ZhnFlatAdYAjI2NDWwjSRrejIFRVa+ari7Jt5MsqqptSRYBjw1oNgGc3Le/BLi96/uR7v2pJNfTu8dxXXfM0cBEkgXAocD2lglJkmbHsJek1gM7P/W0Glg3oM1G4LQkh3U3u08DNiZZkORIgCQHAq8F7h7Q75nAbVXl6kGS5tGMK4wZfAC4Kcn5wEPAWQBJxoC3VdUFVbU9yXuBzd0xl3dl/5xecBwIHADcAlzdtbkG+FSSrfRWFmcPOU5J0pCGCoyqehw4ZUD5OHBB3/5aYO2UNt8HXjZNv0/ThY8k6dnBb3pLkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJajJUYCQ5PMmmJPd374dN02511+b+JKu7sucm2dL3+k6SD3V15yWZ7Ku7YFC/kqS5M+wK42Lg1qpaBtza7T9DksOBy4ATgROAy5IcVlVPVdWKnS/gW8Cf9x16Y1/9J4YcpyRpSMMGxirg2m77WuD1A9q8GthUVdur6glgE7Cyv0GSZcBRwBeHHI8kaZYMGxjPr6ptAN37UQPaLAYe7tuf6Mr6nUNvRVF9ZW9I8rUkNyc5eroBJLkwyXiS8cnJyT2bhSRpRjMGRpJbktw94LWq8RwZUFZT9s8GPt23/1lgaVX9JnALP1/F/GJHVWuqaqyqxhYuXNg4JEnS7lowU4OqetV0dUm+nWRRVW1Lsgh4bECzCeDkvv0lwO19fbwEWFBVd/ad8/G+9lcDV8w0TknS7Br2ktR6YHW3vRpYN6DNRuC0JId1n6I6rSvb6RyeubqgC5+dzgC+MeQ4JUlDmnGFMYMPADclOR94CDgLIMkY8LaquqCqtid5L7C5O+byqtre18fvA6dP6fcdSc4AdgDbgfOGHKckaUhDBUZ36eiUAeXjwAV9+2uBtdP08YIBZZcAlwwzNknSaPlNb0lSEwNDktTEwJAkNTEwJElNDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNTEwJElNhgqMJIcn2ZTk/u79sGna/XWS7yb5qynlxya5ozv+xiQHdeUHd/tbu/qlw4xTkjS8YVcYFwO3VtUy4NZuf5APAm8eUH4FcGV3/BPA+V35+cATVXUccGXXTpI0j4Z9pvcq4ORu+1rgduCPpjaqqluTnNxfliTA7wFv6jv+PcDHun7f05XfDHwkSaqqhhyvNC/u3fY93njVl+blvMsXHTLn59W+adjAeH5VbQOoqm1JjtqNY48AvltVO7r9CWBxt70YeLjrd0eSJ7v235naSZILgQsBjjnmmD2ahDSbVq1YPHOjWbJ80SHzen7tW2YMjCS3AL8yoOpdQ547A8qqoe6ZhVVrgDUAY2NjrkD0rPOmE4/hTSf6y4z2fjMGRlW9arq6JN9OsqhbXSwCHtuNc38HeF6SBd0qYwnwaFc3ARwNTCRZABwKbN+NviVJIzbsTe/1wOpuezWwrvXA7n7EF4AzBxzf3++ZwG3ev5Ck+TVsYHwAODXJ/cCp3T5JxpJ8YmejJF8E/gw4JclEkld3VX8EvDPJVnr3KK7pyq8BjujK38n0n76SJM2R7Eu/uI+NjdX4+Ph8D0OS9ipJ7qyqsZna+U1vSVITA0OS1MTAkCQ1MTAkSU32qZveSSaBb+3h4Ucy4Jvk+zjnvH9wzvuHYeb8a1W1cKZG+1RgDCPJeMunBPYlznn/4Jz3D3MxZy9JSZKaGBiSpCYGxs+tme8BzAPnvH9wzvuHWZ+z9zAkSU1cYUiSmhgYkqQm+11gJFmZ5L4kW5P8wr+Cm+TgJDd29XckWTr3oxythjm/M8m9Sb6W5NYkvzYf4xylmebc1+7MJJVkr/8IZsuck/x+92d9T5Lr53qMo9bwd/uYJF9Iclf39/v0+RjnqCRZm+SxJHdPU58kH+7+e3wtyW+NdABVtd+8gAOAB4AXAAcBXwWWT2nzh8DHu+2zgRvne9xzMOd/Bfxyt/32/WHOXbvnAn8LfBkYm+9xz8Gf8zLgLuCwbv+o+R73HMx5DfD2bns58M35HveQcz4J+C3g7mnqTwc+T++ppS8H7hjl+fe3FcYJwNaqerCqfgTcAKya0mYVcG23fTO9Z3gMemTs3mLGOVfVF6rqB93ul+k9/XBv1vLnDPBe4E+Ap+dycLOkZc7/HvhoVT0BUFW784TMZ6OWORdwSLd9KD9/qudeqar+ll0/fXQVcF31fJneU00Xjer8+1tgLAYe7tuf6MoGtqneo2OfpPdwp71Vy5z7nU/vN5S92YxzTvJS4Oiq+qu5HNgsavlzfiHwwiR/l+TLSVbO2ehmR8uc3wP8QZIJYAPwH+dmaPNmd/9/3y0zPtN7HzNopTD1c8UtbfYmzfNJ8gfAGPC7szqi2bfLOSd5DnAlcN5cDWgOtPw5L6B3WepkeqvILyY5vqq+O8tjmy0tcz4H+NOq+l9JXgF8qpvzT2d/ePNiVn9+7W8rjAng6L79JfziEvVnbZIsoLeM3dUS8NmuZc4keRXwLuCMqvrhHI1ttsw05+cCxwO3J/kmvWu96/fyG9+tf7fXVdWPq+ofgfvoBcjeqmXO5wM3AVTVl4BfoveP9O2rmv5/31P7W2BsBpYlOTbJQfRuaq+f0mY9sLrbPhO4rbq7SXupGefcXZ65il5Y7O3XtWGGOVfVk1V1ZFUtraql9O7bnFFVe/PzfVv+bv8lvQ84kORIepeoHpzTUY5Wy5wfAk4BSPIb9AJjck5HObfWA2/pPi31cuDJqto2qs73q0tSVbUjyUXARnqfsFhbVfckuRwYr6r1wDX0lq1b6a0szp6/EQ+vcc4fBP4F8Gfd/f2HquqMeRv0kBrnvE9pnPNG4LQk9wI/Af5bVT0+f6MeTuOc/wtwdZL/TO/SzHl78y+AST5N75Likd19mcuAAwGq6uP07tOcDmwFfgD825Gefy/+bydJmkP72yUpSdIeMjAkSU0MDElSEwNDktTEwJAkNTEwJElNDAxJUpP/DxcWPZ7IImKHAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f671c3c6358>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q40eE2n6QSv5"
      },
      "source": [
        "### Variational method\n",
        "\n",
        "Returning to the mathematical description of the heat equation,\n",
        "assuming $k$ and $A$ to be constant, it can also be written as follows,\n",
        "\n",
        "$$ \\int_0^1 w(x) {d\\over dx}\\left({dT\\over dx}\\right) dx = 0 $$\n",
        "\n",
        "simply by multiplying by an arbitrary function $w(x)$ and integrating. At first, there may seem to be little gained by writing it in this way, but by manipulating this equation, we can first simplify, then localise the equation, and make it possible to solve numerically.\n",
        "\n",
        "Integration by parts gives:\n",
        "\n",
        "$$ - \\int_0^1 {dw\\over dx} {dT\\over dx} dx + \\left[w {dT\\over dx}\\right]_0^1= 0$$ \n",
        "\n",
        "For now, let's consider the integral, and worry about the boundary terms later."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true,
        "id": "HsGejC14QSv6"
      },
      "source": [
        "### Finite Element method\n",
        "\n",
        "For a numerical approximation, we will replace both $w(x)$ and $T(x)$ with piecewise approximations $w_{piecewise}(x), T_{piecewise}(x)$ to get:\n",
        "\n",
        "$$ \\int_0^1 \\sum_{i,j} w_i {dN_i\\over dx} T_j {dN_j\\over dx} dx = 0 $$\n",
        "\n",
        "Note that we need separate iterators $i$ and $j$ for $w_i$ and $T_j$. If we can solve this equation to get all the $T_j$ values, we will have a numerical solution to the problem.\n",
        "\n",
        "The derivatives $dN_i\\over dx$ are just $+1/\\delta x$, $-1/\\delta x$, or 0, depending on which segment, or \"element\" you are in, where $\\delta x$ is the spacing between points (0.1 in this case). \n",
        "\n",
        "For example, if we take the \"element\" between 0.4 and 0.5, we see that ${dN_4\\over dx} = -{1\\over \\delta x}$ and ${dN_5\\over dx} = {1\\over\\delta x}$. \n",
        "\n",
        "All other derivatives ${dN_i\\over dx} = 0$. \n",
        "\n",
        "As a result, we can break down the integral into small pieces, over each \"element\", which are easy to evaluate:\n",
        "\n",
        "$$ \\int_{0.4}^{0.5} \\sum_{i,j} w_i  {dN_i\\over dx} T_j {dN_j\\over dx} dx = {w_4 T_4 - w_4 T_5 - w_5 T_4 + w_5 T_5 \\over \\delta x^2} \\delta x$$ \n",
        "\n",
        "$$ \\int_{0.5}^{0.6} \\sum_{i,j} w_i  {dN_i\\over dx} T_j {dN_j\\over dx} dx = {w_5 T_5 - w_5 T_6 - w_6 T_5 + w_6 T_6 \\over \\delta x^2} \\delta x$$ \n",
        "\n",
        "etc., and we could continue for the whole range from $x=0$ to $x=1$, integrating each element individually, and summing the contributions from each of them. \n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true,
        "id": "vtYhocYwQSv7"
      },
      "source": [
        "We can also write out each integral in a matrix form, like this:\n",
        "\n",
        "$$ \\int_{0.4}^{0.5} \\sum_{i,j} w_i  {dN_i\\over dx} T_j {dN_j\\over dx} dx = {1\\over\\delta x}(\n",
        "\\begin{array}{cc}\n",
        "w_4 & w_5 \\\\\n",
        "\\end{array})\n",
        "\\left(\n",
        "\\begin{array}{cc} \n",
        "1 & -1 \\\\\n",
        "-1 & 1 \\\\\n",
        "\\end{array}\n",
        "\\right)\n",
        "\\left(\n",
        "\\begin{array}{c}\n",
        "T_4 \\\\\n",
        "T_5 \\\\\n",
        "\\end{array}\n",
        "\\right) $$\n",
        "\n",
        "Summing all the integrals gives the following:\n",
        "\n",
        "$$ \\int_0^1 \\sum_{i,j} w_i {dN_i\\over dx} T_j {dN_j\\over dx} dx = {1\\over\\delta x} \n",
        "[\\begin{array}{cccc} \n",
        "w_0 & w_1 & \\cdots & w_n \n",
        "\\end{array}]\n",
        "\\left(\n",
        "\\begin{array}{cccccc}\n",
        "1 & -1 & 0 & 0 & 0 & \\cdots & 0 \\\\\n",
        "-1 & 2 & -1 & 0 & 0 & \\cdots & 0 \\\\\n",
        "0 & -1 & 2 & -1 & 0 & \\cdots & 0 \\\\\n",
        "\\cdots & & & & \\\\\n",
        "\\end{array}\\right)\n",
        "\\left(\n",
        "\\begin{array}{c}\n",
        "T_0 \\\\ T_1 \\\\ \\cdots \\\\ \\cdots \\\\ T_n\n",
        "\\end{array}\n",
        "\\right) = \\left(\\begin{array}{c}0\\\\0\\\\ \\cdot\\\\ \\cdot\\\\ 0 \\end{array}\\right)\n",
        "$$\n",
        "\n",
        "Because of the arbitrariness of $w(x)$ and hence $w_i$, we can consider each $w_i$ in turn, and assume all the other $w_{j\\neq i} = 0$. Each row of the matrix, multiplied by the vector $T$ must be zero in turn, and we can discard $w$ altogether. We have thus reduced the problem to a matrix equation $A.T = 0$ where $A$ in this case is tridiagonal. Note that we could have got the same matrix equation from a finite-difference approximation."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true,
        "id": "uGs5gmrvQSv9"
      },
      "source": [
        "#### Boundary Conditions\n",
        "Above, we ignored the $[w dT/dx]_0^1$ term, and so far, we haven't specified any boundary condition. Since $w$ is arbitrary, we can just let it be zero, and ignore this term. This is equivalent to a \"natural\" or \"Neumann\" boundary condition. If we wish to impose a hard or \"Dirichlet\" boundary condition, we will have to do something to our matrix above. For example, if we wish $T_0$ to be zero, we can clear the first row of the matrix, and set the diagonal to 1. Let's fully illustrate the example by solving the problem with the Dirichlet BCs described above."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "e-trvaHNQSv-",
        "outputId": "c5672609-dab3-4048-d4c9-e36406064c3c"
      },
      "source": [
        "# Set up matrix problem A.T = b (where T is the unknown temperature)\n",
        "\n",
        "# Full matrix, as shown above\n",
        "A = zeros((11, 11))\n",
        "# Element matrix, representing laplace operator on each element\n",
        "Ke = array([[1.0, -1.0], \n",
        "            [-1.0, 1.0]])\n",
        "# Insert element matrix into full matrix 10 times\n",
        "for i in range(10):\n",
        "    A[i:i+2, i:i+2] += Ke\n",
        "print('A = \\n', A, '\\n')\n",
        "\n",
        "# Right hand side vector (zero)\n",
        "b = zeros(11)\n",
        "\n",
        "# Introduce a source term on the RHS(!)\n",
        "b += 0\n",
        "\n",
        "# Clear a row and set diagonal to one, for Dirichlet BC\n",
        "# Set BC T=100 at x=0\n",
        "A[0, :] = 0.0\n",
        "A[0, 0] = 1.0\n",
        "b[0] = 100.0\n",
        "\n",
        "# Set BC T=0 at x=1\n",
        "A[10, :] = 0.0\n",
        "A[10, 10] = 1.0\n",
        "b[10] = 0.0\n",
        "\n",
        "print('Modified A = \\n', A, '\\n')\n",
        "print('RHS b = ', b)\n",
        "\n",
        "from numpy.linalg import solve\n",
        "T = solve(A, b)\n",
        "\n",
        "print('Solution T =', T)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "A = \n",
            " [[ 1.00 -1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00]\n",
            " [-1.00  2.00 -1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00]\n",
            " [ 0.00 -1.00  2.00 -1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00]\n",
            " [ 0.00  0.00 -1.00  2.00 -1.00  0.00  0.00  0.00  0.00  0.00  0.00]\n",
            " [ 0.00  0.00  0.00 -1.00  2.00 -1.00  0.00  0.00  0.00  0.00  0.00]\n",
            " [ 0.00  0.00  0.00  0.00 -1.00  2.00 -1.00  0.00  0.00  0.00  0.00]\n",
            " [ 0.00  0.00  0.00  0.00  0.00 -1.00  2.00 -1.00  0.00  0.00  0.00]\n",
            " [ 0.00  0.00  0.00  0.00  0.00  0.00 -1.00  2.00 -1.00  0.00  0.00]\n",
            " [ 0.00  0.00  0.00  0.00  0.00  0.00  0.00 -1.00  2.00 -1.00  0.00]\n",
            " [ 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 -1.00  2.00 -1.00]\n",
            " [ 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 -1.00  1.00]] \n",
            "\n",
            "Modified A = \n",
            " [[ 1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00]\n",
            " [-1.00  2.00 -1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00]\n",
            " [ 0.00 -1.00  2.00 -1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00]\n",
            " [ 0.00  0.00 -1.00  2.00 -1.00  0.00  0.00  0.00  0.00  0.00  0.00]\n",
            " [ 0.00  0.00  0.00 -1.00  2.00 -1.00  0.00  0.00  0.00  0.00  0.00]\n",
            " [ 0.00  0.00  0.00  0.00 -1.00  2.00 -1.00  0.00  0.00  0.00  0.00]\n",
            " [ 0.00  0.00  0.00  0.00  0.00 -1.00  2.00 -1.00  0.00  0.00  0.00]\n",
            " [ 0.00  0.00  0.00  0.00  0.00  0.00 -1.00  2.00 -1.00  0.00  0.00]\n",
            " [ 0.00  0.00  0.00  0.00  0.00  0.00  0.00 -1.00  2.00 -1.00  0.00]\n",
            " [ 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 -1.00  2.00 -1.00]\n",
            " [ 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  1.00]] \n",
            "\n",
            "RHS b =  [100.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00]\n",
            "Solution T = [100.00 90.00 80.00 70.00 60.00 50.00 40.00 30.00 20.00 10.00  0.00]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_CILghlYQSwD"
      },
      "source": [
        "As predicted, the solution is linear."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "G9xk42hxQSwE",
        "outputId": "18bc0919-fd9c-4fd2-a70b-ed43b50882da"
      },
      "source": [
        "x = linspace(0.0, 1.0, 11)\n",
        "plt.plot(x,T)\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd4VHX6/vH3k0rvAZEWOgIiJSA1sdBRUAS7YkVBWuLa9rtrWXd/tjUUBSmKYhexgAjS1IQOAaT33oQgvbfP74+Me7EukpDMZDIz9+u6uDIzOTPnPiS5c3Jm5jnmnENERIJXmL8DiIiIb6noRUSCnIpeRCTIqehFRIKcil5EJMip6EVEgpyKXkQkyKnoRUSCnIpeRCTIRfg7AECpUqVcbGysv2OIiASURYsW7XPOxWS2XJ4o+tjYWNLS0vwdQ0QkoJjZ1qwsp0M3IiJBTkUvIhLkVPQiIkFORS8iEuRU9CIiQS7Tojez0Wa218xWXHBbCTObZmbrPR+Le243MxtiZhvMbJmZNfRleBERyVxW9ug/ANr/4bZngRnOuerADM91gA5Adc+/nsA73okpIiLZlWnRO+dSgf1/uLkLMMZzeQxwywW3f+gyzAOKmVlZb4X9o43pR3lz6lpOnjnnq1WIiAS87B6jL+Oc2w3g+Vjac3s5YPsFy+3w3PY/zKynmaWZWVp6enq2QkxbtYe3ftxApyEzWbT1j7+LREQEvP9krF3ktouefdw5N9I5F+eci4uJyfQdvBf1eEJVxjzUhJNnztNt+FxenLCSY6fOZuuxRESCVXaLfs/vh2Q8H/d6bt8BVLhgufLAruzHy1xCjRimJMZzf9NKjJm7hbYDU0ldl72/EEREglF2i34C0MNzuQcw/oLb7/e8+qYpcOj3Qzy+VCg6gpe61GXsY82Ijgzj/tEL+MuXSzl4/LSvVy0ikudl5eWVnwFzgZpmtsPMHgZeBdqY2Xqgjec6wCRgE7ABGAX09knqP9E4tgST+rWi93VV+WbJTlonpzJ5uc9/z4iI5Gnm3EUPoeequLg45+3plSt2HuLpcctYtfswHepewUtd6lC6cD6vrkNExJ/MbJFzLi6z5YL2nbF1yxVlfJ8WPN2+JjPW7KVNcipfpm0nL/xiExHJTUFb9ACR4WH0vq4ak/u3okaZQjw1bhn3j17A9v3H/R1NRCTXBHXR/65qTCG+6NmMf3Spw+KtB2g3KJUPZm/m/Hnt3YtI8AuJogcICzPubxbLlMR44mJL8OJ3q7h9xFw27D3q72giIj4VMkX/u/LFCzDmwca82f0a1u89SsfBMxn60wbOnDvv72giIj4RckUPYGbc1qg805MSaF27NG9MWUuXt2ezYuchf0cTEfG6kCz638UUjmbYPY0Yfm8j0o+eosvQ2bz2wxoNSRORoBLSRf+79nWvYHpiArc1LMc7P2+k4+CZLNyiIWkiEhxU9B5FC0Tyerdr+Pjhazl97jzdh8/l+fErOKohaSIS4FT0f9CyeimmDIjnwRaxfDRvK+0GpvLz2r2Z31FEJI9S0V9EwegIXri5DuMeb07+qHAeeH8hSWN/4cAxDUkTkcCjor+ERpWK832/lvS7oRoTftlFm4EpTFq+W2MURCSgqOgzER0RTlLbmkzo05KyRfPT+5PFPP7xIvYePunvaCIiWaKiz6LaVxbhm97Nea5DLX5em86NySmMXaghaSKS96noL0NEeBiPJVTlhwHxXFW2CE9/tYz73tOQNBHJ21T02VC5VEE+f7Qp/7ylLr9sP0jbgamMnrWZcxqSJiJ5kIo+m8LCjHubVmJqYjxNq5TgHxNX0X34HNbvOeLvaCIi/0VFn0NXFsvP6AcaM+iO+mzed4xOQ2bx1oz1nD6rIWkikjeo6L3AzLilQTmmJSXQru4VvDltHZ3fnsWyHQf9HU1EREXvTaUKRfPWXQ0YdX8cB46f5pahs3ll0moNSRMRv1LR+0Cb2mWYmpjAHY0rMCJ1E+0HpTJv02/+jiUiIUpF7yNF80fyStd6fPrItZx3cOfIefzfN8s5cvKMv6OJSIhR0ftY82oZQ9IeaVmZzxZso+3AVH5cs8ffsUQkhKjoc0H+qHD+dlNtvurVnML5InjogzQGfL6E/RqSJiK5QEWfixpULM7Evq3of2N1vl++m9bJKUxYuktjFETEp1T0uSwqIozENjX4rm9LKhTPT7/PlvDoh4v49ZCGpImIb6jo/aTWFUX4uncL/tbpKmZtSKdNcgqfLdimvXsR8ToVvR+FhxmPtKrClAHx1C1XlOe+Xs7do+az9bdj/o4mIkFERZ8HVCpZkE8fvZZXul7Nip2HaDcolXdnbtKQNBHxChV9HmFm3NWkItOSEmhZrRT//H41Xd+Zw9pfNSRNRHImR0VvZolmttLMVpjZZ2aWz8wqm9l8M1tvZl+YWZS3woaCK4rmY9T9cQy5qwHb9x/nprdmMnDaOg1JE5Fsy3bRm1k5oB8Q55yrC4QDdwKvAQOdc9WBA8DD3ggaSsyMztdcyfSkBDpdXZbBM9Zz01sz+WW7hqSJyOXL6aGbCCC/mUUABYDdwA3AOM/nxwC35HAdIatEwSgG3dmA0Q/EceTkWboOm80/J67ixGkNSRORrMt20TvndgL/BraRUfCHgEXAQefcWc9iO4ByOQ0Z6m6oVYapifHc1aQi787aTLtBqczZsM/fsUQkQOTk0E1xoAtQGbgSKAh0uMiiF33piJn1NLM0M0tLT0/PboyQUThfJP+69Wo+79mUMIO7353Pc18v47CGpIlIJnJy6KY1sNk5l+6cOwN8DTQHinkO5QCUB3Zd7M7OuZHOuTjnXFxMTEwOYoSWplVK8sOAeB5LqMIXC7fTJjmF6as0JE1E/lxOin4b0NTMCpiZATcCq4CfgG6eZXoA43MWUf4oX2Q4z3W4im+faEHxAlE88mEafT9bwr6jp/wdTUTyoJwco59PxpOui4HlnscaCTwDJJnZBqAk8J4XcspF1CtfjAl9WvJkmxpMWfErbZJT+HbJTo1REJH/YnmhFOLi4lxaWpq/YwS09XuO8PRXy1iy7SA31CrNP2+py5XF8vs7loj4kJktcs7FZbac3hkbJKqXKcy4x5vz/E21mbvxN9oOTOXjeVs5rzEKIiFPRR9EwsOMh1pWZsqAeK6pUJS/fbuCu0bNY/M+DUkTCWUq+iBUsWQBPn74Wl6/rR6rdh+m/aBURqRs5Ow5jVEQCUUq+iBlZtzeuALTkxJIqBHDK5PX0PWdOazefdjf0UQkl6nog1yZIvkYcV8jht7dkF0HT3DzW7NInrqWU2c1RkEkVKjoQ4CZ0aleWaYlJtC5/pUM+XEDnYbMYtHWA/6OJiK5QEUfQooXjCL59vq8/2Bjjp86S7fhc3jpu5UcP3028zuLSMBS0Yeg62uWZmpSAvc1rcT7s7fQdmAqs9ZrSJpIsFLRh6hC0RH8o0tdxj7WjMjwMO59bz5Pj1vKoRMakiYSbFT0Ia5J5RJM7t+KXtdV5avFO2mTnMKUlb/6O5aIeJGKXsgXGc4z7Wsx/okWlCoUzWMfLeKJTxaTfkRD0kSCgYpe/qNuuaKM79OCp9rVZNqqPbQZmMLXi3doSJpIgFPRy3+JDA/jieurMal/K6rGFCJp7FIeeH8hOw+e8Hc0EckmFb1cVLXShfjysWa8eHNtFm7ZT9vkFD6cu0VD0kQCkIpe/lRYmPFAi4whaQ0rFef58Su5Y+RcNqYf9Xc0EbkMKnrJVIUSBfjwoSa80a0ea389QofBMxn28wYNSRMJECp6yRIzo3tcBaY/mcCNtUrz+g9ruWXYbFbuOuTvaCKSCRW9XJbShfPxzr2NeOeehvx66BSd357NG1PWcPKMhqSJ5FUqesmWDleXZXpSPLc2KMfQnzbScchM0rbs93csEbkIFb1kW7ECUfy7+zV8+FATTp05T/cRc3lxwkqOndKQNJG8REUvORZfI4apifH0aBbLmLkZQ9JS16X7O5aIeKjoxSsKRkfwYuc6fPlYM6Ijw7h/9AL+8uVSDh4/7e9oIiFPRS9eFRdbgkn9WtHn+mp8s2QnrZNTmbx8t79jiYQ0Fb14Xb7IcP7SriYT+rSgTJFoen2ymF4fL2LvkZP+jiYSklT04jN1rizK+Cda8Ez7WsxYs5fWb6bwZdp2DUkTyWUqevGpiPAwel1Xlcn9W1HzisI8NW4Z949ewPb9x/0dTSRkqOglV1SNKcQXPZvxcpc6LN56gHaDUvlg9mYNSRPJBSp6yTVhYcZ9zWKZkhhP49gSvPjdKrqPmMuGvUf8HU0kqKnoJdeVL16ADx5sTPLt17Ax/SgdB89i6E8bOKMhaSI+oaIXvzAzujYsz7TEBNrUKcMbU9bS+e3ZrNipIWki3qaiF7+KKRzN0LsbMuK+Ruw7eoouQ2fz6mQNSRPxphwVvZkVM7NxZrbGzFabWTMzK2Fm08xsvedjcW+FleDVrs4VTE9MoFvD8gxP2UjHwTNZsFlD0kS8Iad79IOBH5xztYBrgNXAs8AM51x1YIbnukimihaI5LVu9fj44Ws5fe48t4+Yy9+/XcFRDUkTyRHL7ptXzKwIsBSo4i54EDNbC1znnNttZmWBn51zNS/1WHFxcS4tLS1bOSQ4HT99ln9PWcf7czZTtkg+/tX1aq6vWdrfsUTyFDNb5JyLy2y5nOzRVwHSgffNbImZvWtmBYEyzrndAJ6PF/3pNLOeZpZmZmnp6Zp0KP+tQFQEz99cm3GPN6dgdAQPvr+QpC9+4cAxDUkTuVw5KfoIoCHwjnOuAXCMyzhM45wb6ZyLc87FxcTE5CCGBLNGlYozsV9L+t1QjQlLd9FmYArfL9utMQoilyEnRb8D2OGcm++5Po6M4t/jOWSD5+PenEWUUBcdEU5S25p817clZYvm54lPF/PYR4vYc1hD0kSyIttF75z7FdhuZr8ff78RWAVMAHp4busBjM9RQhGPq8oW4ZvezXmuQy1S1qXTOjmFLxZu0969SCay/WQsgJnVB94FooBNwINk/PIYC1QEtgHdnXOXfJ2cnoyVy7V53zGe+WoZCzbvp0W1krxyaz0qlizg71giuSqrT8bmqOi9RUUv2XH+vOPTBdt4dfIazp13/KVdTR5oHkt4mPk7mkiuyI1X3Yj4VViYcW/TSkxNjKdplRK8PHEVt70zh3V7NCRN5EIqegl4VxbLz+gHGjP4zvps/e0YnYbMZMiM9Zw+qyFpIqCilyBhZnSpX47pSQm0r1uW5Gnr6Pz2LJZuP+jvaCJ+p6KXoFKyUDRv3dWAUffHceD4aW4dNptXJq3mxGkNSZPQpaKXoNSmdhmmJSVwR+MKjEjdRIfBqczd+Ju/Y4n4hYpeglaRfJG80rUenz5yLecd3DVqHn/9ZjmHT57xdzSRXKWil6DXvFoppgyI55GWlfl8wTbaJqfy45o9/o4lkmtU9BIS8keF87ebavNVr+YUyR/BQx+k0f/zJfx29JS/o4n4nIpeQkqDisWZ2LcVA1pXZ9Ly3bQZmMr4X3ZqjIIENRW9hJyoiDAGtK7BxL6tqFCiAP0//4VHxqSx+9AJf0cT8QkVvYSsmlcU5utezflbp6uYvXEfbZNT+XT+Ns6f1969BBcVvYS08DDjkVZVmDIgnrrlivLXb5Zz97vz2LLvmL+jiXiNil4EqFSyIJ8+ei2vdL2alTsP035wKqNSN3FOe/cSBFT0Ih5mxl1NKjItKYGW1Urxr0mr6TpsNmt+PezvaCI5oqIX+YMriuZj1P1xDLmrAdsPnOCmIbNInraOU2c1RkECk4pe5CLMjM7XXMn0pARuqleWITPWc/Nbs1iy7YC/o4lcNhW9yCWUKBjFoDsbMPqBOI6cPEvXd+bw8sRVHD991t/RRLJMRS+SBTfUKsPUxHjublKR92Ztpv2gmczZsM/fsUSyREUvkkWF80Xyr1uv5vOeTQkzuPvd+Tz71TIOndCQNMnbVPQil6lplZJM7h/PY/FVGJu2nbYDU5i2SkPSJO9S0YtkQ/6ocJ7reBXfPtGC4gWiePTDNPp8uph9GpImeZCKXiQH6pUvxoQ+LXmyTQ2mrtxD6+QUvlmyQ0PSJE9R0YvkUFREGH1vrM73/VpSuVRBEr9YykMfLGTXQQ1Jk7xBRS/iJdXLFGbc4815/qbazNu0n7YDU/lo3lYNSRO/U9GLeFF4mPFQy8pMTYynfoVi/P3bFdw5ah6bNSRN/EhFL+IDFUoU4KOHm/D6bfVYvfsw7QelMjxlI2fPnfd3NAlBKnoRHzEzbm9cgelJCSTUiOHVyWu4ddgcVu3SkDTJXSp6ER8rUyQfI+5rxNC7G7L70Ak6vz2LN6eu1ZA0yTUqepFcYGZ0qleWaYkJdK5/JW/9uIFOQ2axaKuGpInvqehFclHxglEk316fDx5szInT5+g2fA4vfbeSY6c0JE18J8dFb2bhZrbEzCZ6rlc2s/lmtt7MvjCzqJzHFAku19UszZTEeO5rWon3Z2+h3aBUZq5P93csCVLe2KPvD6y+4PprwEDnXHXgAPCwF9YhEnQKRUfwjy51GftYM6LCw7jvvQU8PW4ph45rSJp4V46K3szKA52Adz3XDbgBGOdZZAxwS07WIRLsmlQuwaT+reh1XVW+WryT1gNT+GHFr/6OJUEkp3v0g4Cngd9fHFwSOOic+/2A4w6gXA7XIRL08kWG80z7Wox/ogUxhaJ5/ONF9P5kEXuPnPR3NAkC2S56M7sJ2OucW3ThzRdZ9KLv/zaznmaWZmZp6ek6NikCULdcUcb3acFT7WoyffVe2iSn8tUiDUmTnMnJHn0LoLOZbQE+J+OQzSCgmJlFeJYpD+y62J2dcyOdc3HOubiYmJgcxBAJLpHhYTxxfTUm9WtFtdKFePLLpfR4fyE7Dhz3dzQJUNkueufcc8658s65WOBO4Efn3D3AT0A3z2I9gPE5TikSgqqVLsSXjzXjpc51SNuyn3YDU/lw7hYNSZPL5ovX0T8DJJnZBjKO2b/ng3WIhISwMKNH81imDIinYaXiPD9+JXeMnMvG9KP+jiYBxPLCsb+4uDiXlpbm7xgieZpzjq8W7+Tlias4ceYc/W+sTs/4KkSG632PocrMFjnn4jJbTt8hIgHCzOjWqDzTkuJpfVVp3piylluGzmbFzkP+jiZ5nIpeJMCULpyPYfc0Yvi9Ddlz+BRdhs7m9R/WcPKMhqTJxanoRQJU+7plmZGUQNcG5Rj280Y6DplJ2pb9/o4leZCKXiSAFS0QyRvdr+HDh5pw6sx5uo+YywvjV3BUQ9LkAip6kSAQXyOGqYnx9GgWy4fzttJuYCop6/RGRMmgohcJEgWjI3ixcx2+fKwZ+SLD6DF6AUljf+Hg8dP+jiZ+pqIXCTJxsSX4vl8r+lxfjQm/7KJ1cgqTlu/2dyzxIxW9SBDKFxnOX9rVZHyfFlxRNB+9P1nM4x8tYu9hDUkLRSp6kSBW58qifNu7Bc+0r8WPa/fSOjmFsWnbNSQtxKjoRYJcRHgYva6ryg/9W1HriiI8PW4Z949ewPb9GpIWKlT0IiGiSkwhPu/ZlJe71GHx1gO0G5TK+7M3c05D0oKeil4khISFGfc1i2VqUgKNY0vw0ner6D58Dhv2HvF3NPEhFb1ICCpXLD8fPNiYgXdcw6Z9x+g4eBZv/7ieM+fOZ35nCTgqepEQZWbc2qA805MSaFOnDP+euo6b35rF8h0akhZsVPQiIa5UoWiG3t2QEfc1Yv+x09wybDavTtaQtGCiohcRANrVuYJpSQl0a1ie4Skb6TB4JvM3/ebvWOIFKnoR+Y+i+SN5rVs9PnnkWs6eP88dI+fx929XcOTkGX9HkxxQ0YvI/2hRrRRTBsTzUIvKfDw/Y0jaT2v3+juWZJOKXkQuqkBUBM/fXJuvejWnYHQED76/kMQvfmH/MQ1JCzQqehG5pIYVizOxX0v63Vid75buok1yChOX7dIYhQCioheRTEVHhJPUpgbf9W1JueL56fPpEnp+tIg9GpIWEFT0IpJlV5Utwte9mvPXjrVIXZdO6+QUvli4TXv3eZyKXkQuS0R4GD3jqzJlQDy1yxbhma+Wc8+789n2m4ak5VUqehHJlthSBfns0ab869a6LNtxiHaDUnlvloak5UUqehHJtrAw455rKzEtKZ5mVUvy8sRV3PbOHNbt0ZC0vERFLyI5VrZoft7rEcfgO+uz9bdjdBoykyEz1nP6rIak5QUqehHxCjOjS/1yTE9KoEPdsiRPW0fnt2exdPtBf0cLeSp6EfGqkoWiGXJXA969P46Dx89w67DZ/L9JqzlxWkPS/EVFLyI+0bp2GaYmxXNH44qMTN1Eh8GpzN2oIWn+oKIXEZ8pki+SV7pezaePXosD7ho1j79+s5zDGpKWq1T0IuJzzauW4of+8TzaqjKfL9hG2+RUflyzx9+xQka2i97MKpjZT2a22sxWmll/z+0lzGyama33fCzuvbgiEqjyR4Xzf51q83XvFhTNH8lDH6TR//Ml/Hb0lL+jBb2c7NGfBZ50zl0FNAWeMLPawLPADOdcdWCG57qICAD1KxTju74tSWxdg0nLd9NmYCoTlmpImi9lu+idc7udc4s9l48Aq4FyQBdgjGexMcAtOQ0pIsElKiKM/q2rM7FvKyqUKEC/z5bw6Idp7D50wt/RgpJXjtGbWSzQAJgPlHHO7YaMXwZA6T+5T08zSzOztPT0dG/EEJEAU/OKwnzdqzl/63QVszbso21yKp/O38Z5jVHwqhwXvZkVAr4CBjjnDmf1fs65kc65OOdcXExMTE5jiEiACg8zHmlVhSkD4qlbrih//WY5d787jy37jvk7WtDIUdGbWSQZJf+Jc+5rz817zKys5/NlAZ1/TEQyValkQT599Fpe7Xo1K3cepv3gVEalbtKQNC/IyatuDHgPWO2cS77gUxOAHp7LPYDx2Y8nIqHEzLizSUWmJSXQqnoM/5q0mq7DZrP2Vw1JywnL7jPdZtYSmAksB36fXPRXMo7TjwUqAtuA7s65/Zd6rLi4OJeWlpatHCISnJxzfL98Ny+MX8nhk2fofV01el9fleiIcH9HyzPMbJFzLi7T5fLCS5pU9CLyZ/YfO83LE1fxzZKd1ChTiNduq0eDinp7DmS96PXOWBHJ00oUjGLgHfUZ/UAcR06epes7c3h54iqOnz7r72gBQ0UvIgHhhlplmJoYzz3XVuS9WZtpP2gmczbs83esgKCiF5GAUThfJP+85Wo+79mU8DDj7nfn8+xXyzh0QkPSLkVFLyIBp2mVkkzu34rHEqowNm07bQemMG2VhqT9GRW9iASkfJHhPNfhKsY/0ZLiBaJ49MM0+ny6mH0akvY/VPQiEtCuLl+U7/q25Mk2NZi6cg+tk1P4ZskODUm7gIpeRAJeZHgYfW+szvf9WlKlVEESv1jKQx8sZNdBDUkDFb2IBJHqZQrz5ePNeeHm2szbtJ+2A1P5aN7WkB+SpqIXkaASHmY82KIyUxPjqV+hGH//dgV3jprHpvSj/o7mNyp6EQlKFUoU4KOHm/B6t3qs2X2YDoNnMjxlI2fPnc/8zkFGRS8iQcvMuD2uAtOTEriuZgyvTl7DLcNms2pXlieqBwUVvYgEvdJF8jH83kYMu6chvx46See3Z/Hm1LWcOnvO39FyhYpeREKCmdHx6rJMT0qgS/1yvPXjBjoNmcWirQf8Hc3nVPQiElKKFYjizduvYcxDTThx+hzdhs/hpe9WcuxU8A5JU9GLSEhKqBHDlMR47mtaifdnb6HdoFRmrg/O81er6EUkZBWKjuAfXery5ePNiIoI4773FvDUl0s5dDy4hqSp6EUk5DWOLcGkfq3ofV1Vvl6yk9YDU/hhxa/+juU1KnoRETKGpD3dvhbjn2hBTKFoHv94Eb0/WcTeIyf9HS3HVPQiIheoW64o4/u04Kl2NZm+ei9tklMZtyiwh6Sp6EVE/iAyPIwnrq/GpH6tqF66EH/5cik93l/IjgPH/R0tW1T0IiJ/olrpQox9rBkvda5D2paMIWlj5mwJuCFpKnoRkUsICzN6NI9lamI8cbEleGHCSm4fMZeNATQkTUUvIpIF5YsXYMyDjXmz+zWs33uUDoNnMvSnDZwJgCFpKnoRkSwyM25rVJ7pSQm0vqo0b0xZS5e3Z7Ni5yF/R7skFb2IyGWKKRzNsHsaMfzehqQfPUWXobN57Yc1nDyTN4ekqehFRLKpfd2yTE9MoGuDcrzz80Y6Dp7Jwi37/R3rf6joRURyoGiBSN7ofg0fPdyE0+fO0334XJ4fv4KjeWhImopeRMQLWlWPYcqAeB5sEctH87bSbmAqKevyxpA0Fb2IiJcUjI7ghZvrMO7x5uSPCqfH6AUkjf2FA8dO+zWXil5ExMsaVSrO9/1a0veGakz4ZRdtBqYwafluv41R8EnRm1l7M1trZhvM7FlfrENEJC+LjgjnybY1mdCnJWWL5qf3J4t5/ONF7D2c+0PSvF70ZhYODAU6ALWBu8ystrfXIyISCGpfWYRvejfn2Q61+HltOq2TUxibtj1X9+59sUffBNjgnNvknDsNfA508cF6REQCQkR4GI8nVGVy/1bUKluEp8ct4773FrB9f+4MSfNF0ZcDtl9wfYfnNhGRkFYlphCfP9qUf95Sl1+2H6TtwFS+W7rL5+v1RdHbRW77n79RzKynmaWZWVp6et54CZKIiK+FhRn3Nq3E1MR4WlQrReVSBX2/Th885g6gwgXXywP/8yvLOTfSORfnnIuLiYnxQQwRkbzrymL5ebdHHHXLFfX5unxR9AuB6mZW2cyigDuBCT5Yj4iIZEGEtx/QOXfWzPoAU4BwYLRzbqW31yMiIlnj9aIHcM5NAib54rFFROTy6J2xIiJBTkUvIhLkVPQiIkFORS8iEuRU9CIiQc78NTbzv0KYpQNbs3n3UsA+L8YJBNrm0KBtDg052eZKzrlM33GaJ4o+J8wszTkX5+8cuUnbHBq0zaEhN7ZZh25ERIKcil5EJMgFQ9GP9HcAP9A2hwZtc2jw+TYH/DF6ERG5tGDYoxcRkUsImKLP7ITjZhZtZl94Pj/fzGJzP6V3ZWGbk8xslZktM7MZZlbJHzm9KasnljdEICFiAAADR0lEQVSzbmbmzCzgX6GRlW02s9s9X+uVZvZpbmf0tix8b1c0s5/MbInn+7ujP3J6i5mNNrO9ZrbiTz5vZjbE8/+xzMwaejWAcy7P/yNj3PFGoAoQBSwFav9hmd7AcM/lO4Ev/J07F7b5eqCA53KvUNhmz3KFgVRgHhDn79y58HWuDiwBinuul/Z37lzY5pFAL8/l2sAWf+fO4TbHAw2BFX/y+Y7AZDLO0NcUmO/N9QfKHn1WTjjeBRjjuTwOuNHMLnZaw0CR6TY7535yzv1+duF5ZJzNK5Bl9cTyLwOvAydzM5yPZGWbHwWGOucOADjn9uZyRm/LyjY7oIjnclEucpa6QOKcSwX2X2KRLsCHLsM8oJiZlfXW+gOl6LNywvH/LOOcOwscAkrmSjrfuNyTrD9Mxh5BIMt0m82sAVDBOTcxN4P5UFa+zjWAGmY228zmmVn7XEvnG1nZ5heBe81sBxnntuibO9H85nJ/3i+LT0484gNZOeF4lk5KHkCyvD1mdi8QByT4NJHvXXKbzSwMGAg8kFuBckFWvs4RZBy+uY6Mv9pmmlld59xBH2fzlaxs813AB865N82sGfCRZ5vP+z6eX/i0vwJljz4rJxz/zzJmFkHGn3uX+lMpr8vSSdbNrDXwf0Bn59ypXMrmK5ltc2GgLvCzmW0h41jmhAB/Qjar39vjnXNnnHObgbVkFH+gyso2PwyMBXDOzQXykTETJlhl6ec9uwKl6LNywvEJQA/P5W7Aj87zLEeAynSbPYcxRpBR8oF+3BYy2Wbn3CHnXCnnXKxzLpaM5yU6O+fS/BPXK7Lyvf0tGU+8Y2alyDiUsylXU3pXVrZ5G3AjgJldRUbRp+dqytw1Abjf8+qbpsAh59xubz14QBy6cX9ywnEz+weQ5pybALxHxp93G8jYk7/Tf4lzLovb/AZQCPjS87zzNudcZ7+FzqEsbnNQyeI2TwHamtkq4BzwlHPuN/+lzpksbvOTwCgzSyTjEMYDgbzjZmafkXHorZTneYcXgEgA59xwMp6H6AhsAI4DD3p1/QH8fyciIlkQKIduREQkm1T0IiJBTkUvIhLkVPQiIkFORS8iEuRU9CIiQU5FLyIS5FT0IiJB7v8DxQt9TWal6IwAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7ff9e0d58240>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "5RHkYEr-QSwJ"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}