{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.5"
    },
    "colab": {
      "name": "Интегрирование уравнений движения.ipynb",
      "provenance": [],
      "collapsed_sections": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ggaOOzeSKJCf"
      },
      "source": [
        "# Интегрирование уравнений движения"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R85WF0AKKJCk"
      },
      "source": [
        "На предыдущих лекциях было показано, что уравнения, определяющие движение тела, могут быть представлены в виде (смотри первую лекцию):\n",
        "\n",
        "$$\n",
        "\\begin{cases}\n",
        "\\frac{d\\vec{r}}{dt} = \\vec{v}(t)\\\\\n",
        "\\frac{d\\vec{v}}{dt} = m^{-1}\\vec{R}(t, \\vec{r})\n",
        "\\end{cases}\n",
        "\\label{equations}\\tag{1}\n",
        "$$\n",
        "\n",
        "Начальные условия при этом:\n",
        "\n",
        "$$\n",
        "\\begin{cases}\n",
        "\\vec{r}(t_0) = \\vec{r}_0\\\\\n",
        "\\vec{v}(t_0) = \\vec{v}_0\n",
        "\\end{cases}\n",
        "\\label{init_cond}\\tag{2}\n",
        "$$\n",
        "\n",
        "Данная запись эквивалентна следующей:\n",
        "\n",
        "$$\n",
        "\\begin{cases}\n",
        "\\frac{d\\vec{y}}{dt} = \\vec{f}(t, \\vec{y})\\\\\n",
        "\\vec{y}(t_0) = \\vec{y}_0\n",
        "\\end{cases}\n",
        "\\label{normal_task}\\tag{3}\n",
        "$$\n",
        "\n",
        "где\n",
        "\n",
        "$$\n",
        "\\vec{y}(t) = \\begin{pmatrix}\n",
        "\\vec{r}(t)\\\\\n",
        "\\vec{v}(t)\\\\\n",
        "\\end{pmatrix} = \n",
        "\\begin{pmatrix}\n",
        "x(t)\\\\\n",
        "y(t)\\\\\n",
        "z(t)\\\\\n",
        "v_x(t)\\\\\n",
        "v_y(t)\\\\\n",
        "v_z(t)\\\\\n",
        "\\end{pmatrix}\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\vec{f}(t, \\vec{y}) = \\begin{pmatrix}\n",
        "\\vec{v}(t)\\\\\n",
        "m^{-1}\\vec{R}(t, \\vec{x})\\\\\n",
        "\\end{pmatrix} = \\begin{pmatrix}\n",
        "v_x(t)\\\\\n",
        "v_y(t)\\\\\n",
        "v_z(t)\\\\\n",
        "m^{-1}R_x(t, \\vec{x})\\\\\n",
        "m^{-1}R_y(t, \\vec{x})\\\\\n",
        "m^{-1}R_z(t, \\vec{x})\\\\\n",
        "\\end{pmatrix}\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\vec{y}_0 = \\begin{pmatrix}\n",
        "\\vec{r}_0\\\\\n",
        "\\vec{v}_0\\\\\n",
        "\\end{pmatrix} = \\begin{pmatrix}\n",
        "\\vec{x}_0\\\\\n",
        "\\vec{y}_0\\\\\n",
        "\\vec{z}_0\\\\\n",
        "\\vec{v_x}_0\\\\\n",
        "\\vec{v_y}_0\\\\\n",
        "\\vec{v_z}_0\\\\\n",
        "\\end{pmatrix}\n",
        "$$\n",
        "\n",
        "\n",
        "Для решения данной задачи был предложен \"Метод Эйлера\", позволяющий рекурсивно выразить решение в последующий момент вемени $t_{i + 1}$ через предыдущий $t_{i}$ с помощью данных формул: \n",
        "\n",
        "$$\n",
        "\\begin{cases}\n",
        "    \\vec{y}_{i + 1} = \\vec{y}_i + h * \\vec{f}(t_i, \\vec{y}_i)\\\\\n",
        "    \\vec{y}_0 = \\vec{y}_0\n",
        "\\end{cases}\n",
        "$$\n",
        "\n",
        "где  h - шаг по времени, $t_0, t_1, ..., t_{N - 1}$ - разбиение времени, $t_i = t_0 + i * h$. Здесь для расчета последующего шага используется информация лишь о состоянии объекта на предыдущем шаге. А что, если учитывать больше информации?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VnCd2ku2KJCl"
      },
      "source": [
        "## Метод Симпсона"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EBKtcJt4KJCm"
      },
      "source": [
        "Как и в предыдущих случаях, будем полагать, что дана система уравнений и начальное условие (3), а также некоторое разбиение времени $t_i = t_0 + h*i$, где h - шаг разбиеня.\n",
        "Метод сипсона выглядит следующим образом:\n",
        "\n",
        "\n",
        "$$\n",
        "\\vec{y}_{i + 1} = \\vec{y}_i + h * (\\frac{1}{6}\\vec{k_1} + \\frac{2}{3}\\vec{k_2} + \\frac{1}{6}\\vec{k_3})\\\\\n",
        "\\vec{k_1} = \\vec{f}(t_i, \\vec{y_i})\\\\\n",
        "\\vec{k_2} = \\vec{f}(t_i + \\frac{h}{2}, \\vec{y_i} + \\frac{1}{2}h\\vec{k_1})\\\\\n",
        "\\vec{k_3} = \\vec{f}(t_i + h, \\vec{y_i} + h * (-\\vec{k_1} + 2\\vec{k_2}))\\\\\n",
        "\\vec{y}_0 = \\vec{y}(0)\n",
        "$$\n",
        "\n",
        "Данный метод гораздо точнее, нежели разобранный ранее метод Эйлера. Реализуем его для уравнения $y' = 4 t^3$ при условии $y(0) = 0$. Решение данной задачи $y = t^4$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "47uYP4daKJCn",
        "outputId": "e035c9f9-b9dd-49f0-95bf-b5626cadbf8d"
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "h = 0.1\n",
        "\n",
        "def f(t, y):  # Функция правой части уравнения\n",
        "    return 4 * t**3  # Здесь может быть другая функция\n",
        "\n",
        "\"\"\"\n",
        "Метод Эйлера\n",
        "\"\"\"\n",
        "t_eul = [0]\n",
        "y_eul = [0]\n",
        "for i in range(11):\n",
        "    y_eul.append(y_eul[-1] + h * f(t_eul[-1], y_eul[-1]))\n",
        "    t_eul.append(i * h)\n",
        "    \n",
        "\"\"\"\n",
        "Метод Симпсона\n",
        "\"\"\"\n",
        "t_sim = [0]\n",
        "y_sim = [0]\n",
        "for i in range(11):\n",
        "    k_1 = f(t_sim[-1], y_sim[-1])\n",
        "    k_2 = f(t_sim[-1] + h / 2, y_sim[-1] + h / 2 * k_1)\n",
        "    k_3 = f(t_sim[-1] + h, y_sim[-1] + h * (2 * k_2 - k_1))\n",
        "    y_sim.append(y_sim[-1] + h * (k_1 / 6 + 2 / 3 * k_2 + 1 / 6 * k_3))\n",
        "    t_sim.append(i * h)\n",
        "    \n",
        "plt.plot(np.linspace(0, 1, 100),np.linspace(0, 1, 100)**4)\n",
        "plt.scatter(t_eul, y_eul, c='r')\n",
        "plt.scatter(t_sim, y_sim, c='b')\n",
        "plt.show()\n",
        "    "
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAj7klEQVR4nO3de3xU9Z3/8ddnZnIhEK4JyC0JVxW8ARHvgl1bqe1Kra1Vg60WTaXbdn/bq7vptl0tvbjb7q5b1GbXS5V47VpLV6xtFaqoIChyBwkhCSFcEi4JyeQ+398fEzTEhEySmUxm8n4+HvPInDNn5nxOEt58c873fL/mnENERGKfJ9oFiIhIeCjQRUTihAJdRCROKNBFROKEAl1EJE74orXjtLQ0l5WVFa3di4jEpHfeeafSOZfe0WtRC/SsrCw2bNgQrd2LiMQkMyvp7DWdchERiRMKdBGROKFAFxGJEwp0EZE4oUAXEYkTXQa6mT1iZofNbGsnr5uZ3W9mhWa22cxmh79MEZHYV1AAWVng8QS/FhSE9/NDaaE/Biw4zeufBKa1PnKBB3tflohIfCkogNxcKCkB54Jfc3PDG+pdBrpz7jXg6Gk2WQg87oLWAsPNbGy4ChQRiQd5eeD3O4ZfuZOEtGoA/P7g+nAJxzn08cC+Nstlres+wsxyzWyDmW2oqKgIw65FRGJDaSkkTTjKsEv2kDim+pT14dKnF0Wdc/nOuWznXHZ6eod3roqIxKWMDBhy/j4CDT78u8aesj5cwhHo+4GJbZYntK4TEZFW37+nicFnHaB22zhcsxeAlBRYujR8+whHoK8Avtja2+VioMo5dyAMnysiEjeSpu3HfAFSKzIwg8xMyM+HnJzw7aPLwbnM7ClgPpBmZmXAD4EEAOfcQ8BK4FqgEPADt4evPBGR+PD0+n3MHDeUFzcOi9g+ugx059zNXbzugL8LW0UiInFm6/4qtpVXc+/CmRHdj+4UFRGJsKfXl5Lk83DdBR12AAwbBbqISATVNbbw+43lfOrcsQwblBDRfSnQRUQi6A+byznR0MwXLpzY9ca9pEAXEYmgJ9eVMnX0EOZOGhnxfSnQRUQiZFt5Fe/tO87NczMws4jvT4EuIhIhT64rJdHn4YbZkb0YepICXUQkAmobmvn9e+V8+ryxDE9J7JN9KtBFRCJgxaZyahqaybkojIO1dEGBLiISAU+uK+XMManMzhjRZ/tUoIuIhNnmsuNs2V/FLRf1zcXQkxToIiJh9sRbJQxK8HJ9H10MPUmBLiISRsf9jazYVM5nZo1naHJk7wxtT4EuIhJGz20oo6E5wK0XZ/b5vhXoIiJhEgg4lq8rITtzBDPGDe3z/SvQRUTC5PXCSkqO+Ln1kr5vnYMCXUQkbJ54q5i0IYksOOeMqOxfgS4iEgZlx/y8svMwN12YQZLPG5UaFOgiImHwxNoSDLi5D+8MbU+BLiLSS3WNLTyzfh/XzDyD8cMHRa0OBbqISC/9/r39HPc3cdulWVGtQ4EuItILzjkee7OYs8cO7ZNJLE5HgS4i0gtri46y8+AJbr80q0/HbemIAl1EpBd+82YxI1ISuO6CcdEuRYEuItJTZcf8/Gn7QW6am0FyQnS6KralQBcR6aHH3yrBzFgUhXFbOqJAFxHpgdqGZp56u5RPnhPdroptKdBFRHrguQ37OFHfzOLLJ0W7lA8o0EVEuqkl4HjkjWJmZwxnVh9OMdcVBbqISDf9ZcchSo/6WXz55GiXcgoFuohINz28Zi/jhw/impljuvfGggLIygKPJ/i1oCCsdSnQRUS6YUtZFW/vPcrtl2Xh83YjQgsKIDcXSkrAueDX3NywhnpI1ZjZAjPbZWaFZnZ3B69nmNkqM9toZpvN7NqwVSgi0o/89+tFDEnyceOFE7v3xrw88PtPXef3B9eHSZeBbmZeYBnwSWAGcLOZzWi32feBZ51zs4CbgAfCVqGISD9RdszPi1sOcPPcid2fALq0tHvreyCUFvpcoNA5V+ScawSeBha228YBJyfQGwaUh61CEZF+4uE1ezHg9st60FUxo5Nx0jtb3wOhBPp4YF+b5bLWdW39CFhkZmXASuDrHX2QmeWa2QYz21BRUdGDckVEoqPK38Qz6/dx3fnjGNeTG4mWLoWUlFPXpaQE14dJuC6K3gw85pybAFwLPGFmH/ls51y+cy7bOZednp4epl2LiETe8nUl+BtbuPPKHnZVzMmB/HzIzASz4Nf8/OD6MPGFsM1+oO3Z/wmt69paDCwAcM69ZWbJQBpwOBxFiohEU31TC4++UcyV09M5e+zQrt/QmZycsAZ4e6G00NcD08xskpklErzouaLdNqXA3wCY2dlAMqBzKiISF17YuJ/Kmga+0tPWeR/pMtCdc83A14CXgR0Ee7NsM7N7zOy61s2+BdxpZpuAp4DbnHMuUkWLiPSVloDj168Vcc74oVw6ZVS0yzmtUE654JxbSfBiZ9t1P2jzfDtwWXhLExGJvj9uPcjeyloeyJkd9RmJuqI7RUVEOuGc44HVhUxOG8w1M8+IdjldUqCLiHTi9d2VbCuv5ivzJuP19O/WOSjQRUQ69cDqQs4Ymsz1syZEu5SQKNBFRDrwbukx1hYd5Y4rJpHoi42ojI0qRUT62AOrChmeksDNc8N3a36kKdBFRNrZVl7FX3Yc5suXTWJwUkidAfsFBbqISDu/erWQ1CQfX7o0K9qldIsCXUSkjfcPneClrQe57bIshg3q5hC5UaZAFxFpY9mqQlISvT0bIjfKFOgiIq32Vtbyh03l3HpxJiMHJ0a7nG5ToIuItFq2qpBEn4c7rujfg3B1RoEuIgIUV9byu437uWVuJumpSdEup0cU6CIiwP2v7ibBa9w1PzZb56BAFxGhqKKGFzbuZ9FFmYxOTY52OT2mQBeRAe+/Xg2eO//KvCnRLqVXFOgiMqDtqajh9+/t54uXZMXsufOTFOgiMqDd/8puknxecvv59HKhUKCLyIC182A1KzaV86VLs0gbEtutc1Cgi8gA9os/vc+QRB93zYv91jko0EVkgCkogKwsSB53nD9vP8RFwyczPCX27grtiAJdRAaMggLIzYWSEhh25S5a/Ik8/S+TKCiIdmXhoUAXkQEjLw/8fkiaeIRBWZVUvTWF2iofeXnRriw8FOgiMmCUlgI4hs/bSfOJZE5szGyzPvYp0EVkwMjIgEFTD5E8/jhVb0yDFu8H6+OBAl1EBox7fhxg5FW7aDoymJrNEwBISYGlS6NcWJgo0EVkwEg6swzfyBoStp+F4SEzE/LzIScn2pWFR+zMfioi0gt1jS38+593MytjOM+vGYNZtCsKP7XQRWRAeOzNYg5W13P3grOweExzFOgiMgAcrW3kgdWFfOys0Vw0eVS0y4kYBbqIxL37X9lNbUMzd3/yrGiXElEhBbqZLTCzXWZWaGZ3d7LNjWa23cy2mdmT4S1TRKRniipqWL62hJvmZjB9TGq0y4moLi+KmpkXWAZ8HCgD1pvZCufc9jbbTAP+EbjMOXfMzEZHqmARke742Us7SfJ5+Ierp0e7lIgLpYU+Fyh0zhU55xqBp4GF7ba5E1jmnDsG4Jw7HN4yRUS6b13REf60/RBL5k+J+ckrQhFKoI8H9rVZLmtd19Z0YLqZvWFma81sQUcfZGa5ZrbBzDZUVFT0rGIRkRAEAo6lK3cwdlgyiy+Pj+FxuxKui6I+YBowH7gZ+G8zG95+I+dcvnMu2zmXnZ6eHqZdi4h81PMb97O5rIrvXHMmgxK90S6nT4QS6PuBiW2WJ7Sua6sMWOGca3LO7QXeJxjwIiJ9rqahmfv+uJMLJg7nMxe0P6EQv0IJ9PXANDObZGaJwE3AinbbvECwdY6ZpRE8BVMUvjJFREL3wKpCDp9o4Id/OwOPJz5vIupIl4HunGsGvga8DOwAnnXObTOze8zsutbNXgaOmNl2YBXwHefckUgVLSLSmdIjfv7n9b18dvZ4ZmWMiHY5fSqksVyccyuBle3W/aDNcwd8s/UhIhI1S1dux+c1vrcgvm8i6ojuFBWRuPH67gpe3naIv7tqKmOGJke7nD6nQBeRuNDYHOCHK7aRNSqFO66YFO1yokLD54pIXHh4zV6KKmp59PYLSfINjG6K7amFLiIx70BVHf/16m4+PmMMV505cEceUaCLSMz78Ys7aAk4fvDpGdEuJaoU6CIS09bsruTFzQdYMn8KE0emRLucqFKgi0jMqm9q4Z9/v5WsUSncNW9KtMuJOgW6iMSsB1bvYW9lLT/+zLkkJ4R4IbSgALKywOMJfi0oiGSJfUq9XEQkJu2pqOGh1XtYeME4Lp+WFtqbCgogNxf8/uBySUlwGSAnJzKF9iG10EUk5jjn+OcXtpKU4CHvU2eH/sa8vA/D/CS/P7g+DijQRSTmPP/uft7cc4TvLjiL0anduCO0tLR762OMAl1EYkplTQP3vridOZkjyJmb0b03Z3SyfWfrY4wCXURiyr/8YTv+hhZ+fsO53R8ad+lSSGnXtTElJbg+DijQRSRmvLLjEH/YVM7XPjaVqaNTu/8BOTmQnw+ZmWAW/JqfHxcXREG9XEQkRtQ0NPP9F7YyfcyQ3vU5z8mJmwBvTy10EYkJP1m5g0PV9fzshvNI9Cm6OqLvioj0e2t2V/LkulLuuGIyswfYLETdoUAXkX6tpqGZ7/3vZianDeabH58e7XL6NZ1DF5F+7acrd1BeVcdv77ok9Nv7Byi10EWk31qzu5KCdaXccfkk5mSOjHY5/Z4CXUT6paq6Jr7z201MTh/Mtz5xZrTLiQk65SIi/dKPVmzj8IkGnl9yqU61hEgtdBHpd1ZuOcDvNu7na1dN5fyJw6NdTsxQoItIv3K4up68323hvAnD+NrHpka7nJiiQBeRfiMQcHz7t5vxN7bwyxsvIMGriOoOfbdEpN949M1iXnu/gu9/egZTRw+JdjkxR4EuIv3C9vJqfv7STq4+ewyLLoqP4Wz7mgJdRKKurrGFbzy9keEpCdz3ufMw6+awuAKo26KI9AP3vridwsM1PLF4LiMHJ0a7nJilFrqIREVBAWRlwZAZ5Ty5rpQr0iZzxbT0aJcV0xToItLnCgogNxf2V9Uy8pot1O8fzvP/ciYFBdGuLLaFFOhmtsDMdplZoZndfZrtbjAzZ2bZ4StRROJNXh74G1pIW/guLmBUrpiFv8ZDXl60K4ttXQa6mXmBZcAngRnAzWY2o4PtUoG/B9aFu0gRiS+lpTDiqh0knVHNkZXn0VKd8sF66blQWuhzgULnXJFzrhF4GljYwXb3Aj8H6sNYn4jEoYzL9zN0TgnV6ydRV3jGh+vVW7FXQgn08cC+Nstlres+YGazgYnOuRdP90FmlmtmG8xsQ0VFRbeLFZHY9/6hEyRctoWm8hEcW33WB+tTUmDp0igWFgd6fVHUzDzAL4FvdbWtcy7fOZftnMtOT9fVbJGBpqahmbuWv8OwwV6+/7HZZE70YAaZmZCfH7dzN/eZUPqh7wcmtlme0LrupFTgHGB1680AZwArzOw659yGcBUqIrHNOcd3f7uJ4spaCu64mEumJLPktmhXFV9CaaGvB6aZ2SQzSwRuAlacfNE5V+WcS3POZTnnsoC1gMJcRE7x4F/3sHLLQb634CwumTIq2uXEpS4D3TnXDHwNeBnYATzrnNtmZveY2XWRLlBEYt+qXYf515d38bfnjyP3ysnRLiduhXTrv3NuJbCy3bofdLLt/N6XJSLxoriylr9/aiNnnTGU+27QOC2RpDtFRSRiquubuPPxDXg8Rv6tcxiUqKnkIkmBLiIR0RJwfOOpjeytrOWBnNlMHJkS7ZLinkZbFJGI+MnKHazeVcFPrj+XS6ekRbucAUEtdBEJu6ffLuXhNXu5/bIsbtFkFX1GgS4iYfX67gq+/8JWrpyeTt61Z0e7nAFFgS4iYbPzYDVLlr/L1NFDWHbLLHya5LlP6bstImFxqLqe2x9dz+AkL4/efiGpyQnRLmnAUaCLSK+dqG/i9kfXU13XxCO3XcjYYYOiXdKApEAXkV5paG7hK0+8w/uHTrAsZzYzxw0L7Y0n56DzeIJfNV1Rr6nbooj0WCDg+Oazm3hzzxF+eeP5zD9zdGhvPDkHnd8fXC4pCS6DhlzsBbXQRaRHnHPc83/beXHzAf7p2rP47OwJob85L+/DMD/J70dz0PWOAl1EeuTf/7Kbx94s5o7LJ5F75ZTuvbmzueY0B12vKNBFpNv+5/Ui7n9lNzdmTyDvUz3oa97ZXHOag65XFOgi0i3PrC/lxy/u4FPnjuWnn+3h6IlLlwbnnGtLc9D1mgJdREL2u41l3P38FuZNT+ffv3ABXk8Ph8LNyQnOOZeZieagCx/1chGRkKzYVM63nt3ExZNG8dCiOST6etkezMlRgIeZWugi0qWVWw7wD8+8R3bmSB6+LVvjmvdTCnQROa2VWw7wjac2csHE4Txy+4WkJOoP+/5KgS4inVqxqZyvP7WR8ycO57HbL2RIksK8P9NPR0Q69Py7ZXz7uU1kZ43kkdsU5rFALXQR+Ygn15Xyrec2cfHkUWqZxxD9lETkFA/9dQ8/e2knV52ZzoOL5pCcoAugsUItdJEB7NQBDx1f+sVOfvbSTj593lh+fWu2wjzGqIUuMkCdMuChBThx1jb+WlFK9oiJ/OdN5/b8piGJGrXQRQaokwMemq+F9M+8S+oFpVS9OYUNDyrMY5Va6CIDVGkpeJIbSb9hA0njj3H0zzM58W4WVcrymKVAFxmgMmbU0nTJenzD6qj8/Wz8u8YG12vAw5ilQBcZgDYUHyVl4QZqauDQ0xfRsH8koAEPY53OoYsMMC9s3M8t/72OMSMS+fb5l3GGb6QGPIwTaqGLDBAtAce/vryLh/66h7mTRvLrRXMYMTiRb3w52pVJuITUQjezBWa2y8wKzezuDl7/ppltN7PNZvaKmWWGv1QR6akT9U3c+fgGHvrrHm65KIPliy9ixODEaJclYdZlC93MvMAy4ONAGbDezFY457a32WwjkO2c85vZEuA+4AuRKFhEuqfw8Alyn3iHkiN+7lk4k1svzuzZLEPS74XSQp8LFDrnipxzjcDTwMK2GzjnVjnnTk7hvRboxvTfIhIpf9x6gIW/eoMqfxPLF1/EFy/JUpjHsVDOoY8H9rVZLgMuOs32i4GXOnrBzHKBXIAM9Y0SiZimlgD/9vIufv1aERdMHM6Di2YzdtigaJclERbWi6JmtgjIBuZ19LpzLh/IB8jOznbh3LeIBJUfr+PrT23knZJjLLo4g3/+9AySfBqTZSAI5ZTLfmBim+UJretOYWZXA3nAdc65hvCUJyLd8erOQ3zq/tfZeaCa+2+exY8/c+7pw/zU0bmCyxKzQmmhrwemmdkkgkF+E3BL2w3MbBbwa2CBc+5w2KsUkdOqb2rhpyt38Ju3Sjh77FAeyJnNpLTBp3/TKaNzASUlwWVQZ/QYZc51febDzK4F/gPwAo8455aa2T3ABufcCjP7C3AucKD1LaXOuetO95nZ2dluw4YNvSpeRGDXwRN846mN7Dp0gi9fNonvLjgztGFvs7KCId5eZiYUF4e7TAkTM3vHOZfd4WuhBHokKNBFeqcl4Hh4TRH/9vL7DB3k498+fz7zzxwd+gd4PNDRv38zCATCV6iE1ekCXXeKisSg0iN+vv3cJt4uPso1M8ew9PpzSRuS1L0PycjouIWuHmgxS4EuEkNaAo7fvFnMv768C5/H+MXnz+ezs8f3rG/50qWnnkMHjc4V4xToIjGi8HAN3/vfzbxTcoyrzkxn6fXnMm54L/qWn7zwmZcXHBw9IyMY5rogGrMU6CL9XH1TCw+u3sODq/cwKNHLL288n+tn9bBV3l5OjgI8jijQRfqxNwsryXthK3sra1l4wTjyPnU2o1OTo12W9FMKdJF+oKDg1DMf3/lhHduTdvDi5gNkjkrhicVzuWJaerTLlH5OgS4SZafc3+Nt4djYvdy3tZCERMc/XD2dr8ybHFq/chnwFOgiUZaXB36/I+XsA4yYtxPfsDr8u8Zg78/g73+aEu3yJIYo0EWi7FDgKGMW7SB5/HEaDw3l0MrzqC9NQ6PcSncp0EWiZOfBau774y7G3HKY5hNJVK48j9qtE8AFk1z390h3KdBF+ljh4Rr+85Xd/N/mclKTfFxzxlksfzALf/WH58l1f4/0hAJdpI/sqahh2auFvPDefpITvNw1bwpfuXIyw1MSuTJd9/dI7ynQRSJsx4Fqlq0q5MUtB0jyebjjisnkXjn5lLFXdH+PhIMCXSQCnHOsLTpK/mt7WLWrgiFJPpbMm8KXL5/U8SBa7Tuiq4kuPaBAFwmjxuYAL209wMNr9rK5rIq0IYl86+PT+eIlWQxLSej4TZpoQsJE46GLhEFlTQNPv13K42+VcPhEA5PTBrP4ikncMHtC1zcFaaIJ6QaNhy4SAc453t57lOXrSvnj1gM0tTiumJbGz284j3nT0/F4QuxIXlravfUinVCgi3RTZU0Dz79bxjPr97GnopbUZB+LLs4k56JMpo4e0v0P1EQTEiYKdJEQNDS3sGrnYf733f2s2nmY5oBjTuYI7rthCn97/jgGJfZirBVNNCFhokAXadW+o8m9P3ZMueQIf9h0gJVbDlBV10R6ahJfvnwSN2ZPYOro1PDsWBNNSJjooqgIbTqa1DmSJhwlZfpBBp99AO/gBlISvXxixhiunz2By6aMwuf1RLtcGcB0UVTkNOqbWsj7VSXJVxxi5LRDeAc3EmjyUF+UTkrleN55ZXTvTqmI9BEFugxI5cfrWL2rglW7DrNmdyXMa2Fwg5e6otH4d42lrigd1+TDDAYlRrtakdAo0GVA8Dc2s67oKK/vrmRNYQXvH6oBYPzwQXxuzgSe+NkY9m0YCS2ntsTV0URiiQJd4lJtQzMbS4+ztugIa4uOsKnsOE0tjkSfh7lZI/n8nInMPzOdqaOHYGac9ac15K6bhZ/BH3xGCrUsvXYjcHn0DkSkGxToEhfKj9exsfQ475YeY0PxUbaWV9MScHg9xrnjh7H48slcNnUUF2aN7PDOzZyVi4BLyeMnlJJBBqUs5Z/IWfkmUNzXhyPSI+rlIv1OV+NUHa1tZOv+Krbsr2LTvuNsLqviYHU9AEk+DxdMHM6FWSPJzhrBnMwRpCZ3MoZKWx4PdPRvwQwCgTAdmUjvqZeLxIxTxqmyAPtP+Pn6fdX8+eAJfGnVbCuv5kBV/QfbT0obzNxJI5mVMZzZGSM4e+xQEn096FaouzUlDijQJepqGpoprqxlT0UNeU/VknJ1DcNG1ZAwshbzBVvHqw8ZUy0Y3jPHDWXmuGGcM25Y5yMYdpfu1pQ4oECXToVriO6mlgAHq+opP17HvmN17DvqZ98xP6VH/BQf8VNZ0/DBtm4mJFYNoqkylbqidJoqU2msSKX5yGCKmyP466q7NSUO6By6dKigAHK/3Iy/8cMQTUlsJv8R3wcZ55zjuL+JypoGKk40cPhEA4eq6zl8ooGDVfUcqKrjQFU9h6rrCbT5NTODM4YmM3FkCpNGDSYzLYWsUYOZnD6YBVnHKWn46GmOTG8Zxc0TIn3YIv1er8+hm9kC4D8BL/A/zrmftXs9CXgcmAMcAb7gnCvuTdEd+erM1eRvv5wWvHhpIXfGGh7YNj/cu+lQwVfXkJefRWnLODK85SzNLSbngch3Z+uL/QYCjtrGZk7UBx/V9U3k/aIam+YjdVATnuRGvIOa8AxqJO+PdSyvDHCkpoGjtY00Bz7aIBiU4GXssGTGDk/mkimjGD98UPAx4sOvSb6O77xc2rCEXPI/2n2w5XtAQViPWyTedBnoZuYFlgEfB8qA9Wa2wjm3vc1mi4FjzrmpZnYT8HPgC+Es9KszV/Pg9nlAcIzpFnzB5ZmrIx7qBV9dQ+6DH/ZRLmmZQO6DI4A1YQnXQMDR4hwtAUdzwNHS4mgKBHjue+v4bsFM6lOT8XlrKfcM5au/m0KJdzXz/985NDYHgo+WAA1NARqaAzQ0t1DfFKC+qYX65hbqG1uoa2rB33jy0Yy/sYXahmZqG4JfaxqbP9rB4xOQ1vrUOQjUJxCoS6ShLoHxw5M5b/wwRg5JJG1IEumpSaQNSWTM0GRGpyYxJMmHWYhjgbeTk/kGlNz50e6DmW/25lssMiB0ecrFzC4BfuScu6Z1+R8BnHM/bbPNy63bvGVmPuAgkO5O8+HdPeXis2Za8DH43H0MnVv0YX04Zsz4aHiEciKpbXmnbO9Ofbr3/SaaTv7f12ZXCdbMhEneD8Iw4BzOgcMRcMHPD7jg+kAg+FqLc63LfBDikZLo85CS6GVQgpfkBC8piV4GJ/oYlOhlSJKPwUleBif5SE1OIDXJx5BkH0OTExg6yMet8yrYX5cZDPKGBHDBA8+kmGKXFbGaPzIdGwQvTubn63y2CL0/5TIe2NdmuQy4qLNtnHPNZlYFjAIq2xWSC+QCZHSzO1gLwT/RA3UJNFWeOonAtDEdv8cIoZVoHT49pYW59TUXfLU1ex2AMxpxzL3KwIL78ljw/LBheDwfLnvNMLMPnns9wWWvJ7js8Rg+j+H1ePB6wOfx4PMaX1sCgYAXAoYLeHDNHggYBOCVP/tI9HlIan0En3tJ8nlITgh+DXnGnA4sbbyX3OM/xc+HA5mkUMvSUb8E7u/x53ZJFydFes45d9oH8DmC581PLt8K/KrdNluBCW2W9wBpp/vcOXPmuO7w0uSCf/yf+vDS1K3P6YlM774O953p3ReX+3XOObd8uVuecJvLZK8zWlwme93yhNucW7488vsWkU4BG1wnuRrKHRj7gYltlie0rutwm9ZTLsMIXhwNm9wZa/joiRTXuj6yluYWk0LtKetSqGVpbnFc7heAnBxyHr2a4sz5BMxHceZ8ch69Wi1lkf6ss6Q/+SB4WqYImAQkApuAme22+TvgodbnNwHPdvW53W2hO+fckhmrWlvqAeelyS2Zsaqn/8l12/Ilr7tM775ga9W7zy1f8npc71dE+idO00IPqR+6mV0L/AfBbouPOOeWmtk9rR+8wsySgSeAWcBR4CbnXFGnH4j6oYuI9ESv+6E751YCK9ut+0Gb5/XA53tTpIiI9I4mRxQRiRMKdBGROKFAFxGJEwp0EZE4EbXRFs2sAuhgRoGQpNHuLtQBQMc8MOiYB4beHHOmcy69oxeiFui9YWYbOuu2E690zAODjnlgiNQx65SLiEicUKCLiMSJWA30/GgXEAU65oFBxzwwROSYY/IcuoiIfFSsttBFRKQdBbqISJzo14FuZgvMbJeZFZrZ3R28nmRmz7S+vs7MsqJQZliFcMzfNLPtZrbZzF4xs8xo1BlOXR1zm+1uMDNnZjHfxS2UYzazG1t/1tvM7Mm+rjHcQvjdzjCzVWa2sfX3+9po1BkuZvaImR02s62dvG5mdn/r92Ozmc3u9U47G1c32g+CQ/XuASbz4TjsM9pt81VOHYf9mWjX3QfHfBWQ0vp8yUA45tbtUoHXgLVAdrTr7oOf8zRgIzCidXl0tOvug2POB5a0Pp8BFEe77l4e85XAbGBrJ69fC7xEcPbLi4F1vd1nf26hzwUKnXNFzrlG4GlgYbttFgK/aX3+W+BvrKfTzfcPXR6zc26Vc+7kDMprCc4gFctC+TkD3Av8HKjvy+IiJJRjvhNY5pw7BuCcO9zHNYZbKMfsgKGtz4cB5X1YX9g5514jOD9EZxYCj7ugtcBwMxvbm33250DvaHLq8Z1t45xrBk5OTh2rQjnmthYT/B8+lnV5zK1/ik50zr3Yl4VFUCg/5+nAdDN7w8zWmtmCPqsuMkI55h8Bi8ysjOD8C1/vm9Kiprv/3rsU0gQX0v+Y2SIgG5gX7Voiycw8wC+B26JcSl/zETztMp/gX2Gvmdm5zrnj0Swqwm4GHnPO/cLMLgGeMLNznHOBaBcWK/pzC71fTE7dx0I5ZszsaiAPuM4519BHtUVKV8ecCpwDrDazYoLnGlfE+IXRUH7OZcAK51yTc24v8D7BgI9VoRzzYuBZAOfcW0AywUGs4lVI/967oz8H+npgmplNMrNEghc9V7TbZgXwpdbnnwNeda1XG2JUl8dsZrOAXxMM81g/rwpdHLNzrso5l+acy3LOZRG8bnCdcy6WJ6QN5Xf7BYKtc8wsjeApmNPO09vPhXLMpcDfAJjZ2QQDvaJPq+xbK4AvtvZ2uRiocs4d6NUnRvtKcBdXia8l2DLZA+S1rruH4D9oCP7AnwMKgbeBydGuuQ+O+S/AIeC91seKaNcc6WNut+1qYryXS4g/ZyN4qmk7sIXgxOtRrzvCxzwDeINgD5j3gE9Eu+ZeHu9TwAGgieBfXIuBu4C72vyMl7V+P7aE4/dat/6LiMSJ/nzKRUREukGBLiISJxToIiJxQoEuIhInFOgiInFCgS4iEicU6CIiceL/A3zb2ia4ojE5AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LauoKDPZKJCp"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}