{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "view-in-github"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/SotaYoshida/Lecture_DataScience/blob/main/notebooks/Python_misc_NewtonsMethod.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aLf3PqMgkFNC"
      },
      "source": [
        "# ニュートン法によるN次元多項式の求根\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0ZIFvSR9kOo9"
      },
      "source": [
        "$n$次元多項式$f(x) = c_0 + c_1 x + c_2 x^2 + \\ldots + c_nx^n$が与えられた時、  \n",
        "$f(x)=0$となるような解(根)をどのように求めたら良いだろうか？\n",
        "\n",
        "$n=2$までの場合については、中学までで習うようによく知られた公式がある。  \n",
        "$n=3,4$の場合についてもそれぞれ、[カルダノの公式]や[フェラリの解法]として知られる  \n",
        "代数的な求根法が存在する。  \n",
        "なお、5次以上の代数方程式には代数的な解法※が存在しないことが知られている(アーベル-ルフィニの定理)  \n",
        "(※方程式の係数の有限回の四則演算および冪根操作で解を表示すること)\n",
        "\n",
        "\n",
        "方程式の解や積分値などが代数的/解析的に解けない場合は、もちろん数値計算の出番となる。  \n",
        "\n",
        "このノートブックでは、ニュートン法(あるいはニュートン・ラフソン法)として知られる求根アルゴリズムを紹介する。  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tPW10YYsnOHs"
      },
      "source": [
        "## ニュートン法"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZnosTDdTnPFL"
      },
      "source": [
        "ニュートン法による関数$f(x)=0$の求根アルゴリズム自体は非常にシンプルで\n",
        "\n",
        "1. 初期値$x \\in \\mathbb{R} $ (または$x \\in \\mathbb{C})$を決める\n",
        "2.  上の初期値を$x_0$とでも呼ぶことにして、  \n",
        "    $ x_{t+1} = x_t - f(x_t) / f'(x_t) $  と$x$を更新し  \n",
        "    $x$の値が収束するまで更新を繰り返す\n",
        "\n",
        "というものである。  \n",
        "\n",
        "$\\clubsuit$  進んだ注  \n",
        "高次元のベクトル$\\boldsymbol{x}$に拡張することもできる。   \n",
        "$ \\boldsymbol{x}_{t+1} = \\boldsymbol{x}_{t} - \\partial f(\\boldsymbol{x}_t)^{-1} f(\\boldsymbol{x}_t)$\n",
        "$\\partial f(\\boldsymbol{x}_t)$はヤコビ行列になるが、  \n",
        "数値計算上の困難が有り、実際に解く場合には様々な工夫が必要となる。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lP5pa0eSpqfz"
      },
      "source": [
        "### ニュートン法のアルゴリズムの実装例\n",
        "\n",
        "$n=5$の多項式の例を用いて、ニュートン法のアルゴリズムを実装してみよう。  \n",
        "実数係数の多項式の場合、非ゼロの係数を持つ最大の$n$が奇数次ならば、  \n",
        "関数$f(x)$は必ず実軸を横切る(つまり、$f(x)=0$の解が存在する。  \n",
        "\n",
        "まず係数を適当に...\n",
        "$c_0 = 1.0, c_1 = 1.0, c_2 = 2.0, c_3 = 3.0, c_4 = 4.0, c_5 = 5.0$\n",
        "とでもして、プロットしてみよう。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 276
        },
        "id": "mDpCwZF0kEnj",
        "outputId": "c5eb065f-87bf-42d0-fc38-1c57f88b8848"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAABIAAAAEDCAYAAABNtH+LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZhcVZ3/8c+prfd9TTrdSSedPSELTQIhIYQQQIwgCigiI4oyzrjxqOOM+pPx0dFZXEGdAQZRQGQRhQz7biCQfU93lk7Snd73tXqr7fz+qCJ0MGFLJ5Wuer+ep6i659669Q33VFfXp+89x1hrBQAAAAAAgNjliHYBAAAAAAAAOLUIgAAAAAAAAGIcARAAAAAAAECMIwACAAAAAACIcQRAAAAAAAAAMY4ACAAAAAAAIMZFLQAyxtxjjGk1xux5D9teYIzZZowJGGOuftu6zxhjqiK3z5y6igEAAAAAAMamaJ4B9HtJl73HbWsl3SjpjyMbjTHZkv5V0mJJiyT9qzEma/RKBAAAAAAAGPuiFgBZa1+V1DmyzRgzxRjzrDFmqzHmNWPMjMi2NdbaXZJCb9vNpZJesNZ2Wmu7JL2g9x4qAQAAAAAAxAVXtAt4m7skfdFaW2WMWSzpvyVd9A7bF0mqG7FcH2kDAAAAAABAxBkTABljUiUtkfQnY8ybzQnRqwgAAAAAACA2nDEBkMKXo3Vba+e/j+c0SLpwxPIESX8dxZoAAAAAAADGvDNmGnhrba+kamPMNZJkwua9y9Oek3SJMSYrMvjzJZE2AAAAAAAARERzGvgHJa2XNN0YU2+MuUnS9ZJuMsbslFQh6crItucYY+olXSPpTmNMhSRZazsl/VDS5sjtB5E2AAAAAAAARBhrbbRrAAAAAAAAwCl0xlwCBgAAAAAAgFMjKoNA5+bm2kmTJkXjpQEAAAAAAGLS1q1b2621ecdbF5UAaNKkSdqyZUs0XhoAAAAAACAmGWOOnGgdl4ABAAAAAADEOAIgAAAAAACAGEcABAAAAAAAEOMIgAAAAAAAAGIcARAAAAAAAECMIwACAAAAAACIcQRAAAAAAAAAMY4ACAAAAAAAxKW9Tb267cUqWWujXcopRwAEAAAAAADizrbaLn3izvV6cFOtOvt90S7nlCMAAgAAAAAAceWNg+369N0blZXi0Z++eJ5yUhOiXdIp54p2AQAAAAAAAKfLC5Ut+tIft6k0J0X337RI+emJ0S7ptCAAAgAAAAAAcWHNjgZ9/ZGdmlOUoXs/e44ykz3RLum0IQACAAAAAAAx7/4NR3Trmj1aXJqtuz9zjlIT4isSia9/LQAAAAAAiCvWWv3ixSrd/lKVVs7I12+uX6hEtzPaZZ12BEAAAAAAACAmBYIhfW9NhR7cVKtryyfox1fNlcsZn/NhEQABAAAAAICYM+QP6qsPbtfzlS360oop+uYl02WMiXZZUUMABAAAAAAAYkrPgF+fv2+zthzp0vc/Mks3nl8a7ZKijgAIAAAAAADEjPquAX32d5tV09GvX123QKvPGh/tks4IBEAAAAAAACAm7Krv1k33btGQP6h7P7tIS8pyo13SGYMACAAAAAAAjHkvVrboKw9uV3aKRw98frGmFaRFu6QzCgEQAAAAAAAY0+5bX6Pv/1+FZo/P0G9vLFd+WmK0SzrjjFoAZIxxStoiqcFau3q09gsAAAAAAHA8wZDVj5/eq9+uq9bFM/N1+3ULlOzhXJfjGc3/K1+TtFdS+ijuEwAAAAAA4G/0Dfn11Qe365X9bbpxySR9b/UsOR3xO837u3GMxk6MMRMkfVjS3aOxPwAAAAAAgBOp7RjQx/77Db1a1a5/++gcff+K2YQ/72K0zgD6paRvSTrhCEvGmJsl3SxJJSUlo/SyAAAAAAAgnmw83KEv/mGrQla6/3PM9PVenfQZQMaY1ZJarbVb32k7a+1d1tpya215Xl7eyb4sAAAAAACIMw9vrtWnf7tRWSkePf6l8wl/3ofROAPofElXGGMul5QoKd0Y8wdr7adHYd8AAAAAACDO+QIh/fDJSt2/4YiWTc3Vrz+1UBlJ7miXNaacdABkrf22pG9LkjHmQknfJPwBAAAAAACjobVvSP/4h23acqRLf3/BZP3TpdPlco7KkMZxhbnRAAAAAADAGWlbbZf+4Q9b1TPo1+3XLdAV88ZHu6Qxa1QDIGvtXyX9dTT3CQAAAAAA4s9Dm2p165oKFWQk6C//cL5mjU+PdkljGmcAAQAAAACAM8agL6hb1+zRn7bWa9nUXN3+yQXKSvFEu6wxjwAIAAAAAACcEarb+/UPf9iqfc19+upFZfraxdPkdJholxUTCIAAAAAAAEDUPb27Sd96dJdcTqPfffYcrZieH+2SYgoBEAAAAAAAiBpfIKT/eGaf7nm9WvOLM/Wb6xeqKDMp2mXFHAIgAAAAAAAQFTXt/frqQ9u1q75HNy6ZpO9cPlMeF1O8nwoEQAAAAAAA4LRbs6NB331sj5wOozs+fbYum1MY7ZJiGgEQAAAAAAA4bQZ8Af3rmgr9aWu9yidm6bbrFnDJ12lAAAQAAAAAAE6LPQ09+upD21Xd3q+vXFSmr62cKpeTS75OBwIgAAAAAABwSgVDVne+eki/eOGAslM8euDzi7VkSm60y4orBEAAAAAAAOCUqesc0Dce2alNNZ368Nxx+tFVc5SZ7Il2WXGHAAgAAAAAAIw6a60e296gW9dUSJJ+fu08XbWgSMaYKFcWnwiAAAAAAADAqGrrG9b/e3y3nqto0aJJ2frZtfNUnJ0c7bLiGgEQAAAAAAAYNU/uatT3Ht+jfl9Q//KhGfrCsslyOjjrJ9oIgAAAAAAAwEnr8A7r1jUVemp3k+ZNyNDPrp2nsvy0aJeFCAIgAAAAAADwgVlr9dTuJv3rmgr1DQX0rcum6+Zlk5ne/QxDAAQAAAAAAD6Qpp5Bfe/xPXpxb6vmFmXop9fM0/RCzvo5ExEAAQAAAACA9yUUsnpgU63+85l9CoRC+u7lM/XZ8ydx1s8ZjAAIAAAAAAC8Zwdb+/Ttv+zW5pouLS3L1Y+vmquSHGb4OtMRAAEAAAAAgHc16Avq169U6a5XDyvZ49JPrj5LV589QcYww9dYQAAEAAAAAADe0cv7WnTrmgrVdw3qYwuL9J3LZyo3NSHaZeF9IAACAAAAAADH1dg9qB88UalnK5pVlp+qh24+V+dOzol2WfgACIAAAAAAAMAxhvxB3f3aYf3mlUOysvqnS6frC8smy+NikOexigAIAAAAAABIkqy1eqGyRT98qlJ1nYP60JxCfefymSrOZpDnsY4ACAAAAAAA6GCrVz94slKvHmjT1PxUPfD5xTq/LDfaZWGUEAABAAAAABDHOrzDuu2lKj2wsVbJHqduXT1LN5w3UW4nl3vFEgIgAAAAAADi0JA/qHvfqNGvXz6oAX9Qn1pUolsunqocZveKSQRAAAAAAADEkVDI6sndTfqvZ/epvmtQF83I13cun6Gy/LRol4ZTiAAIAAAAAIA48VpVm/7z2X3a09CrGYVp+sNNi7V0KuP8xAMCIAAAAAAAYtzOum7913P79PrBDhVlJunn187TlfOL5HSYaJeG04QACAAAAACAGLW/uU+3vXRAT+9uVnaKR7eunqXrzy1RgssZ7dJwmhEAAQAAAAAQY/Y19+r2l6r09O5mpSa49NWLyvSFCyYrLdEd7dIQJQRAAAAAAADEiL1N4eDnmT3h4OcrF5XppqWlykz2RLs0RNlJB0DGmGJJ90kqkGQl3WWtve1k9wsAAAAAAN6bysZw8PNsxVtn/HyO4AcjjMYZQAFJ37DWbjPGpEnaaox5wVpbOQr7BgAAAAAAJzAy+Ekj+ME7OOkAyFrbJKkp8rjPGLNXUpEkAiAAAAAAAE6BbbVdunPtIT1X0RIOflZO1U3nlyojmTF+cHyjOgaQMWaSpAWSNh5n3c2SbpakkpKS0XxZAAAAAABiXjBk9eLeFv3vq4e15UiX0hNd+trKqfocwQ/eg1ELgIwxqZL+LOkWa23v29dba++SdJcklZeX29F6XQAAAAAAYtmQP6hHt9brt+uqVd3erwlZSfrXj8zSteXFSklgbie8N6PSU4wxboXDnwestX8ZjX0CAAAAABDPOrzDum/9Ed2/4Yg6+306a0KGfv2pBbpsdqFcTke0y8MYMxqzgBlJv5W011r785MvCQAAAACA+HW4zau711Xrz1vrNRwI6eKZ+frCsslaVJqt8Fdw4P0bjTOAzpd0g6TdxpgdkbbvWGufHoV9AwAAAAAQ84Ihq5f3teoPG47o1ao2uZ0OfXxhkW5aWqqy/LRol4cYMBqzgK2TRAQJAAAAAMD71No3pIc31enBTbVq7BlSflqCvnLRVN1w7kTlpSVEuzzEEEaLAgAAAADgNLLWav3hDj2woVbPVTQrELJaWparWz8ySytnFsjN+D44BQiAAAAAAAA4DXoG/frz1no9sPGIDrX1KyPJrc8smaTrF5docl5qtMtDjCMAAgAAAADgFAmFrDYc7tCj2+r19O4mDflDmlecqZ9cfZY+Mm+8Et3OaJeIOEEABAAAAADAKDvc5tVftjXose0NaugeVFqCS1ctKNL1iydqTlFGtMtDHCIAAgAAAABgFPQM+vXkrkb9eWu9ttV2y2GkpVPz9K3LpuvS2YWc7YOoIgACAAAAAOADCgRDeq2qXY9uq9cLlS3yBUKaVpCqb39ohj66oEgF6YnRLhGQRAAEAAAAAMD7EgiGtLG6U0/tbtKze5rV2e9TVrJbn1pUoo8vnKA5RekyxkS7TOAYBEAAAAAAALyLQDCkTdWdenJ3k57b06yOfp+SPU6tnFmg1WeN04rp+fK4mL4dZy4CIAAAAAAAjiMYstpY3aGndoXP9Ono9ynJ7dTKmflafdY4LZ+WryQP4/pgbCAAAgAAAAAgYsgf1BuH2vXi3lY9X9Gsdm849LloZr5Wzx2nC6cT+mBsIgACAAAAAMS11r4hvby3VS/ubdXrB9s16A8q2ePUiun5+nDk8i5CH4x1BEAAAAAAgLhirVVlU69e2tuql/a2aGd9jyRpfEairj57glbOzNe5k3OYth0xhQAIAAAAABDzeof8Wn+oQ69Vtenlva1q7BmSJM0rztQ3Vk3TypkFmjkujdm7ELMIgAAAAAAAMccfDGlHXbdeq2rXuqo27azvUTBklexx6vyyXH3t4qlaMSNf+WmJ0S4VOC0IgAAAAAAAY561Vofa+rWuqk3rDrZrw+FOeYcDchhp7oRM/cPyKVo6NVcLS7KYrh1xiQAIAAAAADDmvBn4bKru1OaaTm083HH0sq6S7GRdMX+8lpXlasmUXGUku6NcLRB9BEAAAAAAgDNeMGS1t6lXm6o7j4Y+Hf0+SVJuaoIWl2brS2U5WlaWp5Kc5ChXC5x5CIAAAAAAAGec/uGAdjf0aOuRLm2u6dTWmi71DQckSROykrR8ep4Wl2ZrUWmOJuUkM3gz8C4IgAAAAAAAURUMWVW19mlHbbd21IVvB1r6FLLh9VPzU/WR+eO1uDRb50zK1vjMpOgWDIxBBEAAAAAAgNPGWqumniHtqu+JhD1d2l3fo35fUJKUnujSvOJMXTKrQPNLMjW/OEvZKZ4oVw2MfQRAAAAAAIBTIhiyqm73qqKxVxWNvaps7FVFY4+6BvySJJfDaNb4dH387AmaX5yp+cWZKs1N4XIu4BQgAAIAAAAAnLT+4YCqWr2qbOxVZVOPKhp7ta+pT4P+8Jk9HqdD0wpTdcmsQs0uStfs8RmaPT5diW5nlCsH4gMBEAAAAADgPRv0BXWw1av9LX2qaunTgZY+HWjxqqF78Og2aQkuzRyfrk8uKj4a9JTlp8rtdESxciC+EQABAAAAAI5hrVVnv0/V7f063N6vw239OtgaDnrqugZkI4Mze5wOTc5L0dkTs3TdomJNLUjTzMJ0FWcncRkXcIYhAAIAAACAONU/HFB1e/8xt8Pt/apu86p3KHB0O7fTqDQ3RXMnZOjjCydoWkGqphakaVJOslyc1QOMCQRAAAAAABCjAsGQmnqGVNc5oNrOAdV1Dai2c1B1nQOq7xpQu9d3zPZFmUkqzU3RlfOLVJqbotK8FE3OTVFRZhJBDzDGEQABAAAAwBjlHQ6oqXtQjT1DR+8buwfV2D2ouq4BNXYPKRiyR7d3OoyKMpNUnJ2kVbMKNCErWZMjQc/E7BQleRiQGYhVBEAAAAAAcIbxB0Nq9w6rtXdYbX3Dau0L3zf3DqmpZ1BN3UNq7BlU34jLtCTJGCk/LUHjM5O0oDhLV85LVnF2koqzk1WclaxxGYmcyQPEKQIgAAAAADgNrLXqHQq8Fex4h9XaO6Q277Daji6H7zv7fcfdR3aKR+MzE1WSk6xzJ2drXGaSxmUkanzkviA9kZm2ABwXARAAAAAAvE+hkFXvkF9dA3519vvUPeBTZ79PXQM+dQ341fXm435/pC3cPvJyrDd5nA7lpSUoLy1BJTnJKp+Upby0BOWnJUbuw+tyUxPkcRHuAPhgCIAAAAAAxJ1gyMo7FFDvkD98Gwyob8iv3qHI/dFlv/oi2/UNBdQ7GN6me8Cn42Q5kiSXwygrxaOsZLeykj2akpeqrBSPslPCy7mpCccEOxlJbqZMB3DKjUoAZIy5TNJtkpyS7rbW/sdo7BcAAABA/AqFrIYCQQ35QxryByO3UKQtfOsfDmrAF3jr3hdU//CxywPDkfsR2w34gu/6+klup9KTXEpLdCs90aWsZI8m5qQoLdGl7GSPMpPdyk7xRMIeT7gtxa20BBeBDoAzzkkHQMYYp6TfSFolqV7SZmPM/1lrK0923wAAAABOLWutQlYKhEIKhqwCIatQyMoftPIHQ0dvvsBby75gKLw+8Lblo9u+bTkYkj9gjwY3w/6QhkcGO8eEPCEN+4MaDoSf9355nA4lJziV4nEp2eNUcoJLKR6nMpM9SklwKtkTXk5JcCk9ya20RJfSIwFPWqL7aOCTluhiLB0AMWU0zgBaJOmgtfawJBljHpJ0paS4CIB+//vfa/78+Zo/f76CwaDuv/9+LVy4UGeddZb8fr8eeOABlZeXa86cORoaGtJDDz2kxYsXa+bMmRoYGNAjjzyi8847T9OnT5fX69Wjjz6qpUuXqqysTD09PXrsscd0wQUXaPLkyerq6tKaNWt04YUXatKkSWpvb9eTTz6plStXqri4WK2trXr66ae1atUqFRUVqbm5Wc8++6wuu+wyFRYWqqGhQS+88IIuv/xy5efnq66uTi+99JJWr16t3Nxc1dTU6K9//auuvPJKZWVl6fDhw3r11Vd11VVXKSMjQwcPHtS6det09dVXKzU1Vfv379f69et17bXXKjk5WXv37tXGjRv1yU9+UomJidqzZ4+2bNmi66+/Xm63W7t27dK2bdt0ww03yOl0aseOHdqxY4duvPFGSdLWrVtVUVGhv/u7v5Mkbd68WQcOHND1118vSdqwYYOqq6t13XXXSZLeeOMN1dfX69prr5UkrVu3Ts3Nzbr66qslSWvXrlVHR4c+9rGPSZJeeeUV9fb26sorr5QkvfjiixocHNRHPvIRSdLzzz8vv9+vD3/4w5KkZ599VpJ02WWXSZKeeuopud1uXXLJJZKkJ554QklJSbr44oslSWvWrFF6erpWrFghSfrLX/6inJwcLV++XJL06KOPqrCwUEuXLpUkPfLII5owYYKWLFkiSXrwwQdVWlqqc889V5L0wAMPaNq0aTrnnHMkSffdd59mz56ts88+m75H34uLvldaWqrOzi793//9n5ZdsFwTSkrU3t6u5555WudfcKHGjZ+gttZWvfLSczrvgouUX1Co1pYWvbH2JS2+YIWycwvU2tyoza+v1TnLViojO0ctTY3auXGdFi69SKkZ2WptrFPltg2at/RiJadmqK3xiKp2btH8ZauUmJKm1vojOrxnq+Ytu0SepBS11lWrpnK7zrrgMnkSk9Ry5JBq9+/S/OWXy+lOUPORA2qoqtBZyz8sh8ut5ur9ajq0V3NXrJbD4VTz4X1qqd6nuSvCfaHpUKXa6w5p9vKPyFqrpkMV6mo8ohlLL5dk1Vy1Wz2tDZq6JNwXmg7slLezRVMWr5K1UtOBHRrs6VDpOSslSc37tmnI26NJ5eG+0LR3i3yD/Zq4MNwXGis2KeAbVsmCZZKkhj0bFAoGVTzvfElS/a71kqQJZ50nSarb+bocTqeK5oT7Ru321+T0JGj8rHDfqN22Vu6kVI2bGe4bR7a+ooTUTBVOXyBJqtnyspLSc1QwbZ4kqXrTi0rJzld+2VmSpMMbn1da7njlTZkjSTq0/lllFJYot3SWJOng608rq2iycibNkCRVrXtS2SXTlFMyTTYU0sE3nlbOxOnKLp6qUCCgQxueVW7pTGUVTVHQ79Phjc8rb/JsZY4vVWB4SNWbX1R+2VxlFE6Uf2hANVteVsHUeUovKJZv0KsjW/+qwmkLlJpXpOH+XtXteFWF089Wau44DXm7Vb9zncbNPEcp2QUa7O1Uw+43NH7WYiVn5Wmwp0MNe9araM55SsrI0UBXmxorN6po7hIlpWerv7NFTXs3a8K8pUpMzZS3vUnN+7eqeP4F8qSkqa+tQa0Hdqh44XJ5klLV11qv1qqdKjl7hdyJyeptrlXbod2aWL5SroRE9TTWqL26QpMWrZLT7VF3w2F11OxV6eJL5XC51FV3UJ21+zX5vA/JOBzqrD2grroqTTk//HOmo2afuhurNfm8D8laqaOmUn0tdZq4+NLw+sN75G1vVMk5qyRJ7Yd2a7C7VRPOXilZqf3QTg33dmr8ggvD66t2aLi/R+Pnhfta24FtCgz1q2BuuK+179+ioH9Y+XPCP3fa9m2WDQaVN+tcWUntezdKknJnLApvv3eD5HAqZ3q4r7VVvCGnO0FZU8+WtVbtFa/LlZiijCnzw9vveU2u5AxllIb7VvvuV+VOy1L6pLmy1qpj91p5MnKVVjI7/O/b9Yo8mYVKLZ4ZXt75khJyi5RSNF3WSp07XlRi/kQljZ8qSerc/rwSCycroWCKQsGgene/LHfBZLnzShUK+DWwd62cBWVy5pQo6B9WoOp1mbwy2cwJsoEh2cPrFcqbqlD6OFnfkFy1m+TPnaZAar7kG1BCw1YN5UyTPzlPxtevlOYd8mZPky8pR2bIq4z2XerJmq7hhCw5h/uU1Vmh9ozpGnClyzXco4Le/apPmSavI1WJvh4VDx5UladMXkeKUvzdmuKv1i7HFPUpSRnBbs20ddoQnKxem6B82605zka95p+sfutRkaNHc11NWuubokG5Vezo1mxXs17xTdGw3Jro6NJMV4te8pXJL5dKnZ2a7mzVC76pCsqpyc4OTXO26TnfNFk5VOZsV5mzXS+HZsnjdGias00T1KG9KfOU6HKqKNCgcf4OdRQuUqLbobTeGiUMtstRtlSJbqdM635Zb4cK569Qotupnprd8vV1as75lyrB7VBd5VYN9XXr4suvUIrHqS3r16nf2/e2z9zed/7M9fKZy+97/L4X732vrKxM8WA0AqAiSXUjluslLX77RsaYmyXdLEklJSWj8LIAgLez1ioYshr2h9TWNyxfMKT+4fBsI3saeuQPhtTZPyx/e7+G97VqOBBSS++Qeuu6VeOolS8YUmPngJoPtOmNnir5gyENdPTr8M5GPVbnUTBk5W7zau/GWt2z1yoYDCqjpU/bXjus3q1DCgUCGtfeq3UvVqnz1R7ZgF+T+3r0/JOVan2mVY6gT7OGuvVYw041qU6u4LAW2i7dX71VDaEqJYSGtMTVqbuqNqoxuFepZlhL3R26o2qjWkKVSjdDWuLu1J0HN6k1VKlMM6hz3V26855NarcpyjYDWuTu0h2/26JOm6xc069yd7fuuHeLum2S8h1eLXR16477t6nXJqrA0acFrh7d8Ydt8toEjXP0al5k+a0vQj367/u3jfgi1Kvf3L91xBehXv36vi0jvgj16vb7toz4ItSr2+7dMuKLUJ9+ee8WSdI0Z5tKnb365X3h5RnOVhU7e/TLP2yVJM1ytmics1e/eGCbJGmOq1l5Dq9+cWi7JGmuq0nZZkC/OBhenudqUoYZ0s+rwssLXM1KMT79bH94+WxXixJMQD/dF14+x9UqpwlpQ2V4eZG7VZK0qSK8fK67TUHr0Obd4eUl7nYNW5e27tohSVrq7lS/9Wr7zvDyBe4u9dhB7dwevuxhubtLnXZYu7eFB8lY4elWW21Ae7aE/6K/0tOjplqrys3hKYxXeXpVV1+vfRvDM99c6ulTdWOdDmwYkiRd5vHqYFOtDr4xIKOQLvV49XRTrQ6/3i+nglrl8eqppiOqDvbJrYBWerx6orFGR0I9SpBfKzxerWmsVl2oS0nya7nHq8caD6sh1KEU49Myt1d/bjykplBbpO959afGg2oJtUT6nlePNFapNdQU6XtePdRwQO22IdL3vPrjc/tH9D2v/vDs/hF9z6v7nt43ou959fun9o7oe1797om9I/qeV79dUzmi73n1v4/vGdH3vLrzL7tH9D2v/ufPu0b0Pa9+86edI/qeV79+eOeIvtenXz+yc0Tf69Wv/rTzmL73q0d3HdP3bq/eFel7zco2A7rt8O5I32tRhhnS7Yd2R/peq1KMT+sO7o70vTYlmIDeqNoT6Xvt4b53YE+k73WE+97+ikjf6wz3vX0Vkb7XFe57eysifa9b/XZA2ysrI32vRz12WDsrKiN9r0ed1q/deyojfa9XbY0h7dntiPS9PjU1GlWGy9Mqj1d1zS3at9NG+l6/qlubdWBnMNL3BnSwrVkHdwQifW9AB9pbdDgYkMsEtcozpAOdLaqzQXlMQMtdQ6rqblejsUowfp1nfDrk7VSLwyhJPp0d8unQULc6HA4la1hzAwHVtPSpy+lWsh3UTF9AtZ0D6nMlKDk0oAR/SK19Pg0MDSk5OKy0YEgDvoB8CsklyWGkZI9TiYkeJfkTlBhwqDQ3RUrKUMJQSEntbp1blCNHUrpcAy45mlv10dLxciWly/S1yN/QrS9ML5U7KVWBriYN1PXqG3OnKSEpWcMd9equ8erSRWcpKTlZvc3VatP+WVYAABpGSURBVKoa1KcuPEcpyYlqOXJQNft8uuXDy5ScmKCaqkrtr9ilH39ilZI8blXu2aXdu3bqjhs/JOmtL+G3/92Fkt78Eh7U9deHv4Ru2OBXdbVP110RDuveeKNT9fXDuvb8UknSOl+tmpsHtWJGviRpbX2iOvxuFWUmSZKcDi67AoATMdaeYOSy97oDY66WdJm19vOR5RskLbbWfvlEzykvL7dbtmw5qdcFgLHGWqvhQEje4YD6hwPyDofHH/AOBzToC58SP+g/dpyDQV/4tPhBX2S8A99bp8m/uS7c9tbySf5YP4bLYeR2OuR2hu+dDiOXw8jpNHI5Riwfcx9pd56g/c1l59+2OxxGTmPkMJIxRo7IY4fDyBi9tfw360a2v+25ZuRz33rOyPV6816SkfTmsA1GZsTjN//zt+1vjvMQ2dXRx28uHdtuTvgaOkF7+Pnv/Bqn0+ke1iIaX+dO99gd0fk3jnx9c9z2E24/YsGcaJsRffZ42+o9vP6x+37vr/k3+x+FfRqjY36WmKNtBA4AgDOHMWartbb8eOtG4wygBknFI5YnRNoAIGb4AqHIDCF+9URm/3jrcXimEO+wX/3DkYEnfQF533w8Iuw53tSvJ+IwUqLbqSS3U4lupxLdjsh9uC0r2X3M8sj1CS6HElwOeVwOuZ3H3ntGLjsd8riMPE6n3C4jj9Mhd2SbNwMfAAAAAGPfaARAmyVNNcaUKhz8fFLSp0ZhvwAw6vzBkLoGfOrs96nT61Nn5HGH16fuAd8Jw51B/zvPFOJxOpSa6FJKZNDJlASXMpLcKspMPLqcmuBScoJTqQmuY9pSEpxK8oRDnCS3UwmRe7fT8JdlAAAAAKPipAMga23AGPNlSc8pPA38PdbaipOuDADeA2utugb8au0bUkvvsNr6htXV71NHv+/ofWf/sLoG/OrwDqt3KHDCfaUnupSZ7FF6Ung2kLL81PCsIEnhMCc9ya30RHfkcaQtMdye6Haexn81AAAAALw/o3EGkKy1T0t6ejT2BQBSONjpHvCrpW9Irb3DaukdUmvfsFp7w0HPm+1vDnT8dm6nUXaKR9kpCcpOcasoK1nZye7wcqpH2ckeZad4lJPqUVayR1nJbrmY6hUAAABAjBqVAAgA3q++Ib8augfV2D2ohq5BNXQPHV1u7hk6YbCTnuhSQXqiCtITtbg0RfnpiSpIT1BBeqLy0xKUl5ag7BSPUhNcXD4FAAAAABEEQABOiUFfULWdAzrS0a/azgHVdg5Egp7wre9tl2K5nUbjMpI0PjNRi0qzlZ+eoIK0xPB9euLRx1xqBQAAAADvHwEQgA+sZ8CvQ+1e1XYM6EjHQCTo6deRjgG19g0fs21aoksTspI1IStJi0uzNT4zSeMzk1SUlaSizCTlpSbIwYxTAAAAAHBKEAABeEfBkFVD16AOtXlH3Pp1uM2rdq/vmG0L0xNVkpOs5dPyNDEnWSU5KSrJTtbE7GRlJru5JAsAAAAAooQACICkcNBT09Gv/c192t/cp6rWPh1q7Vd1R798gbfG4slO8WhybopWzijQlPwUleamalJOsoqzk7k8CwAAAADOUARAQJyx1qqtb1j7IkHPvuY+7W/pVVWLV8ORoMcYaWJ2ssryU7V8ep6m5KVoSl6qJuelKjvFE+V/AQAAAADg/SIAAmKYtVZHOga0u6FHexp6tLuhR3ubetU14D+6TV5agmYUpumGcydqemGaZhSmqyw/VUkezuYBAAAAgFhBAATEiFDkEq6RYU9FY+/R2bY8ToemF6bpklmFmjEuTdML0zS9IE05qQlRrhwAAAAAcKoRAAFjVId3WNtru7Wttkvbaru0p6FX3uFI2ONyaGZhmq6YN15zijI0tyhD0wrS5HE5olw1AAAAACAaCICAMSAYstrf3KettV3afiQc+NR0DEiSXA6jWePTddWCIs0pStecSNjjdhL2AAAAAADCCICAM9CgL6itR7q0sbpDW490aWddt/p9QUlSbqpHC0qy9IlzSnT2xCzNLcpgvB4AAAAAwDsiAALOAIO+oLbVdmn9oQ5tONyhnfXd8getnA6jGYVp+vjZE7SwJEsLS7JUnJ0kY0y0SwYAAAAAjCEEQEAUvBn4bDgcDnx21L0V+MwpytDnlpbq3Mk5Kp+YpbREd7TLBQAAAACMcQRAwGlgrdW+5j69eqBNaw+0aUtNl3zBEIEPAAAAAOC0IAACTpHuAZ9eq2rXqwfa9GpVm1p6hyVJMwrTdOP5k3TeFAIfAAAAAMDpQQAEjJJgyGpnfbfW7g+f5bOrvlshK2UkubV0aq6WT8vTBVPzVJiRGO1SAQAAAABxhgAIOAlD/qDWVbXrhcoWvbi3RR39PhkjzS/O1Fcumqrl0/M0b0KmnA4GbQYAAAAARA8BEPA+dQ/49PK+Vj1f0aK1B9o06A8qLcGlFTPydfGsAi0ry1VWiifaZQIAAAAAcBQBEPAeNHQP6oWKZj1f2aKN1Z0KhqwK0hN09dkTdMnsAi0uzZHH5Yh2mQAAAAAAHBcBEHACTT2DenJnk57Y1ahd9T2SpKn5qfri8sm6ZFah5hZlyMGlXQAAAACAMYAACBihwzusp/c064mdjdpc0ylrpbMmZOjbH5qhS2YXqjQ3JdolAgAAAADwvhEAIe71Dfn1XEWLntjZqHUH2xUMWU3NT9XXL56m1fPGE/oAAAAAAMY8AiDEJV8gpJf2tmjNjka9vL9VvkBIE7KSdPMFk3XFvPGaUZgmY7i8CwAAAAAQGwiAEFf2Nffqkc31enxHgzr7fcpLS9CnFpXoivnjtaA4k9AHAAAAABCTCIAQ83oG/XpiZ6Me2VKnXfU9cjuNLplVqGvKJ2jZ1Dw5GcgZAAAAABDjCIAQk0Ihqw3VHXpkc52e2dOs4UBIMwrTdOvqWfrogiJlp3iiXSIAAAAAAKcNARBiSmvfkB7eVKdHttaprnNQaYkuXVM+QZ8oL9GconQu8QIAAAAAxCUCIIx51lptr+vWvW/U6OndTfIHrZZMydE3L5muS2cXKtHtjHaJAAAAAABEFQEQxqwhf1BP7mrSfetrtKu+R2kJLn363Im64dyJmpyXGu3yAAAAAAA4YxAAYcxp7B7UAxuP6MFNders96ksP1U/vHK2rlo4QakJdGkAAAAAAN6Ob8sYE6y12ljdqXvfqNHzlS2y1mrlzALduGSSlkzJYWwfAAAAAADeAQEQzmjBkNXzFc36n7WHtKu+R5nJbn1+Wak+vXiiirOTo10eAAAAAABjwkkFQMaYn0j6iCSfpEOSPmut7R6NwhDfhgNBPb69QXeuPazD7f2alJOsH181V1ctKFKSh0GdAQAAAAB4P072DKAXJH3bWhswxvynpG9L+ueTLwvxyjsc0B83HtFv11WrpXdYc4rS9ZtPLdRlcwrldHCZFwAAAAAAH8RJBUDW2udHLG6QdPXJlYN41e4d1u9fr9F962vUOxTQkik5+uk187S0LJfxfQAAAAAAOEmjOQbQ5yQ9fKKVxpibJd0sSSUlJaP4shjLGroHdefaQ3p4c518wZAum12oLy6fonnFmdEuDQAAAACAmPGuAZAx5kVJhcdZ9V1r7ZrINt+VFJD0wIn2Y629S9JdklReXm4/ULWIGS29Q/rNKwf14KZaSdLHFkzQzcsna0peapQrAwAAAAAg9rxrAGStvfid1htjbpS0WtJKay3BDt5Rh3dYd6w9pPvWH1EwZHVNebG+clGZxmcmRbs0AAAAAABi1snOAnaZpG9JWm6tHRidkhCLegb8uuu1Q/rd6zUa8gd11YIJ+trKqSrJYSp3AAAAAABOtZMdA+jXkhIkvRAZqHeDtfaLJ10VYkbfkF+/e71G//vaYfUNBbT6rHG65eJpKsvnUi8AAAAAAE6Xk50FrGy0CkFsGfQFdd/6Gt2x9pC6BvxaNatAX181TTPHpUe7NAAAAAAA4s5ozgIGKBiyenRrnX7+wgG19A5r+bQ8fX3VNGb1AgAAAAAgigiAMGrWHmjTvz+9V/ua+zS/OFO/um6hFpVmR7ssAAAAAADiHgEQTlplY6/+/Zm9eq2qXSXZyfrNpxbq8rmFiowLBQAAAAAAoowACB9Yh3dYP3vhgB7aVKu0RLe+t3qWPn1uiRJczmiXBgAAAAAARiAAwvvmD4Z03/oj+uWLBzTgC+ozSybplpXTlJHsjnZpAAAAAADgOAiA8L6sPdCmHzxRoUNt/Vo2NVe3rp6lqQVp0S4LAAAAAAC8AwIgvCcN3YP64ROVeraiWZNykvXbz5Trohn5jPMDAAAAAMAYQACEd+QLhHT3usP61UsHZWX1T5dO1+eXlTLODwAAAAAAYwgBEE7ojUPt+t7je3SorV+rZhXo1tWzVJydHO2yAAAAAADA+0QAhL/R2e/Tj57aqz9vq1dxdpLuubFcF80oiHZZAAAAAADgAyIAwlHWWj2+o0E/fHKvegf9+tKKKfrKRVOV6OZyLwAAAAAAxjICIEiSajsG9N3Hd+u1qnbNL87Uf3x8rmYUpke7LAAAAAAAMAoIgOJcMGT1+zdq9JPn9snlcOgHV87W9Ysnyulgdi8AAAAAAGIFAVAcO9Tm1bce3aWtR7p00Yx8/eiqORqXkRTtsgAAAAAAwCgjAIpDwZDVPeuq9dPn9yvR7dTPr52nqxYUyRjO+gEAAAAAIBYRAMWZw21efeNPO7W9tlurZhXoRx+do/z0xGiXBQAAAAAATiECoDhhrdX9G47ox0/vVYLLqds+OV9XzBvPWT8AAAAAAMQBAqA40NwzpH96dKdeq2rX8ml5+q+rz1IBZ/0AAAAAABA3CIBi3BM7G/X/Ht8jXyCkf/voHF2/uISzfgAAAAAAiDMEQDHKOxzQrY/v0V+2N2h+caZ+8Yn5Ks1NiXZZAAAAAAAgCgiAYtCOum597aHtqusc0C0XT9WXV5TJ5XREuywAAAAAABAlBEAxJBSyuuPVQ/r58wdUkJ6oh//+PJ0zKTvaZQEAAAAAgCgjAIoRrb1DuuXhHXrjUIc+PHecfnzVXGUku6NdFgAAAAAAOAMQAMWA1w+262sPbVf/cFD/+fG5ura8mIGeAQAAAADAUQRAY1gwZHX7S1W6/eUqleWl6sEvLNTUgrRolwUAAAAAAM4wBEBjVFvfsG55eLteP9ihjy0s0r99dI6SPRxOAAAAAADwt0gMxqCNhzv05Qe3q3fQr/+6+ixdW14c7ZIAAAAAAMAZjABoDLHW6rfrqvXvz+zTxOxk3X/TIs0oTI92WQAAAAAA4AxHADRGDPgC+uc/79YTOxt16ewC/fSaeUpLZJYvAAAAAADw7giAxoCa9n79/f1bVdXap29dNl3/sHwKs3wBAAAAAID3jADoDPfKvlZ99aHtcjqM7v3cIi2bmhftkgAAAAAAwBjjGI2dGGO+YYyxxpjc0dgfwuP93LH2kD5372aVZCfriS8vJfwBAAAAAAAfyEmfAWSMKZZ0iaTaky8HkjTkD+pf/rxLj+9o1OqzxuknV89TkscZ7bIAAAAAAMAYNRqXgP1C0rckrRmFfcW9lt4h3XzfFu2s79E3L5mmL60oY7wfAAAAAABwUk4qADLGXCmpwVq7891CCmPMzZJulqSSkpKTedmYtbOuW1+4b4u8wwHddcPZumR2YbRLAgAAAAAAMeBdAyBjzIuSjpdEfFfSdxS+/OtdWWvvknSXJJWXl9v3UWNceGZ3k255eIfy0hL0l5uWaEZherRLAgAAAAAAMeJdAyBr7cXHazfGzJVUKunNs38mSNpmjFlkrW0e1SpjmLVWd756WP/xzD4tLMnU//5duXJSE6JdFgAAAAAAiCEf+BIwa+1uSflvLhtjaiSVW2vbR6GuuOAPhvT/Htujh7fUafVZ4/TTa+Yp0c1gzwAAAAAAYHSNxiDQ+AB6Bv36xwe26vWDHfryijJ9fdU0ORwM9gwAAAAAAEbfqAVA1tpJo7WvWNfYPajP3LNJNR39+snVZ+ma8uJolwQAAAAAAGIYZwCdZvuae3XjPZvVPxzQvZ9dpCVludEuCQAAAAAAxDgCoNNo/aEO3Xz/FiV7nHrki+dp5jhm+gIAAAAAAKceAdBp8uSuRn394Z0qyUnWvZ9bpKLMpGiXBAAAAAAA4gQB0Glwz7pq/fCpSp1dkqW7P1OuzGRPtEsCAAAAAABxhADoFLLW6ucvHNCvXj6oS2cX6LZPLmCadwAAAAAAcNoRAJ0ioZDV95+o0H3rj+gT5cX60VVz5HI6ol0WAAAAAACIQwRAp4A/GNI3/7RTa3Y06uYLJuvbH5ohY0y0ywIAAAAAAHGKAGiUDfqC+tIft+nlfa36p0un6x8vnEL4AwAAAAAAoooAaBT1Dfl1071btLmmU//20Tn69LkTo10SAAAAAAAAAdBo6Rn06zP3bNLuhh798hPzdeX8omiXBAAAAAAAIIkAaFR09fv06d9u1IGWPv339Qt16ezCaJcEAAAAAABwFAHQSWrrG9an796o6o5+3XVDuVbMyI92SQAAAAAAAMcgADoJzT1D+tTdG9TYPajf3XiOzi/LjXZJAAAAAAAAf4MA6ANq6hnUJ+/aoPa+Yd33ucVaVJod7ZIAAAAAAACOiwDoA8pIcmtqfqp+8Yn5WliSFe1yAAAAAAAATogA6ANK9rh092fOiXYZAAAAAAAA78oR7QIAAAAAAABwahEAAQAAAAAAxDgCIAAAAAAAgBhHAAQAAAAAABDjCIAAAAAAAABiHAEQAAAAAABAjCMAAgAAAAAAiHEEQAAAAAAAADHOWGtP/4sa0ybpyGl/4VMjV1J7tItAVHDs4xvHP35x7OMXxz6+cfzjF8c+vnH849dYPfYTrbV5x1sRlQAolhhjtlhry6NdB04/jn184/jHL459/OLYxzeOf/zi2Mc3jn/8isVjzyVgAAAAAAAAMY4ACAAAAAAAIMYRAJ28u6JdAKKGYx/fOP7xi2Mfvzj28Y3jH7849vGN4x+/Yu7YMwYQAAAAAABAjOMMIAAAAAAAgBhHAAQAAAAAABDjCIDeI2PMNcaYCmNMyBhT/rZ13zbGHDTG7DfGXDqi/bJI20FjzL+c/qpxKhhjHjbG7IjcaowxOyLtk4wxgyPW3RHtWjG6jDHfN8Y0jDjGl49Yd9yfA4gdxpifGGP2GWN2GWMeM8ZkRtp578cBPtPjhzGm2BjzijGmMvK739ci7Sf8DEBsifx+tztynLdE2rKNMS8YY6oi91nRrhOjyxgzfcT7e4cxptcYcwvv/dhljLnHGNNqjNkzou2473UTdnvk94BdxpiF0av8g2MMoPfIGDNTUkjSnZK+aa1988NglqQHJS2SNF7Si5KmRZ52QNIqSfWSNku6zlpbeZpLxylkjPmZpB5r7Q+MMZMkPWmtnRPdqnCqGGO+L8lrrf3p29qP+3PAWhs87UXilDHGXCLpZWttwBjzn5Jkrf1n3vuxzxjjFJ/pccMYM07SOGvtNmNMmqStkj4q6Vod5zMAsccYUyOp3FrbPqLtvyR1Wmv/IxICZ1lr/zlaNeLUivzcb5C0WNJnxXs/JhljLpDklXTfm7/Hnei9Hgn+viLpcoX7xW3W2sXRqv2D4gyg98hau9dau/84q66U9JC1dthaWy3poMJfAhdJOmitPWyt9Ul6KLItYoQxxij8y+CD0a4FUXeinwOIIdba5621gcjiBkkTolkPTis+0+OItbbJWrst8rhP0l5JRdGtCmeAKyXdG3l8r8KhIGLXSkmHrLVHol0ITh1r7auSOt/WfKL3+pUKB0XWWrtBUmbkDwZjCgHQySuSVDdiuT7SdqJ2xI5lklqstVUj2kqNMduNMWuNMcuiVRhOqS9HTvu8Z8Tp37zf48/nJD0zYpn3fmzjPR6nImf4LZC0MdJ0vM8AxB4r6XljzFZjzM2RtgJrbVPkcbOkguiUhtPkkzr2j7y89+PHid7rMfG7AAHQCMaYF40xe45z4698ceY99oXrdOwHQ5OkEmvtAklfl/RHY0z66awbJ+9djv3/SJoiab7Cx/tnUS0Wo+69vPeNMd+VFJD0QKSJ9z4Qg4wxqZL+LOkWa22v+AyIJ0uttQslfUjSlyKXiRxlw2NoMI5GjDLGeCRdIelPkSbe+3EqFt/rrmgXcCax1l78AZ7WIKl4xPKESJv+f3v382JjFMdx/P3JyMKWZEnZKxulpPwuKVZj4UdZUKwtbJSVLGzYiR3KgkyI/APyK6XBAlEk/gAli6/F80xd08yYYow59/3a3HvPfe7t1H2+5zz323POd4Z2/ed+dy4kGQH2AusGPvMd+N4/f5rkLd1+UE/msKv6y2Y7DiS5CNzuX840DmgBmUXsHwJ2AZv7iwJjfzgY40MmyWK65M+VqroBUFVfBt4fnAPUmKr61D9+TXKTbhnolyQrq+pzv+zj67x2UnNpJ/BsIuaN/aEzXaw3cS3gHUB/bgwYTbIkySpgDfCIboPINUlW9Vnk0f5YtWEL8LqqPk40JFnebxhHktV058K7eeqf5sCkdb57gImKAdONA2pIkh3ACWB3VX0baDf22+ecPkT6Pf4uAa+q6txA+3RzgBqSZGm/+TdJlgLb6H7rMeBgf9hB4Nb89FD/wC93+Rv7Q2e6WB8DDvTVwNbTFQL6PNUX/M+8A2iWkuwBzgPLgTtJnlfV9qoaT3IdeEm3JODYROWfJMeB+8Ai4HJVjc9T9/X3TV4XDLAROJ3kB13FuKNVNXlTMS1sZ5OspbsV9D1wBGCmcUBNuQAsAR50/w95WFVHMfab11d+c04fHhuA/cCLJM/7tpPAvqnmADVnBXCzH+dHgKtVdS/JY+B6ksPAB7pCIGpMn/Tbyq/xPeX1nxa+JNeATcCyJB+BU8AZpo71u3QVwN4A3+iqwy04loGXJEmSJElqnEvAJEmSJEmSGmcCSJIkSZIkqXEmgCRJkiRJkhpnAkiSJEmSJKlxJoAkSZIkSZIaZwJIkiRJkiSpcSaAJEmSJEmSGvcTIpMUKhx0SIAAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 1440x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light",
            "tags": []
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "import numpy as np \n",
        "from matplotlib import pyplot as plt\n",
        "\n",
        "\n",
        "## 多項式の係数を定義\n",
        "cs = np.array( [1.0, 1.0, 2.0, 3.0, 4.0, 5.0] )\n",
        "\n",
        "## 愚直に多項式とその微分を定義\n",
        "def xpoly(x):\n",
        "    return np.array([1.0,x,x**2,x**3,x**4,x**5])\n",
        "def xpoly_prime(x):\n",
        "    return np.array([0.0,1.0,2*x,3*x**2,4*x**3,5*x**4])\n",
        "\n",
        "# 適当な区間で点を作成して描画してみる\n",
        "xr = np.arange(-100,100,1.0) \n",
        "yr = np.array([ np.dot(cs,xpoly(x)) for x in xr])\n",
        "fig = plt.figure(figsize = (20,4))\n",
        "ax = fig.add_subplot(111)\n",
        "ax.plot(xr,yr)\n",
        "ax.plot(xr,0.0*yr,linestyle=\"dotted\",color=\"gray\") # y=0\n",
        "plt.show()\n",
        "plt.close() \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1anohi6wsU2l"
      },
      "source": [
        "ゼロ点($f(x)=0$となる$x$)があることはわかるが、このスケールだと具体的な値はよくわからない。\n",
        "\n",
        "次にニュートン法のアルゴリズムに対応する関数を作って解を求めてみよう"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "dqnujbsFspn8",
        "outputId": "ceb388d5-e79c-4190-bce1-b6cf812fabe3"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "it 1 x -16.03067071498258 y -15383219.0 -5041023.966714532\n",
            "it 2 x -12.854892481957197 y -5041023.966714532 -1651968.106196959\n",
            "it 3 x -10.313885990777683 y -1651968.106196959 -541379.1850391383\n",
            "it 4 x -8.28061578625216 y -541379.1850391383 -177429.79752234442\n",
            "it 5 x -6.653441486353502 y -177429.79752234442 -58155.11457333112\n",
            "it 6 x -5.3510406877452255 y -58155.11457333112 -19063.39609222064\n",
            "it 7 x -4.308351359326742 y -19063.39609222064 -6250.001804089965\n",
            "it 8 x -3.4733344532483637 y -6250.001804089965 -2049.456302435744\n",
            "it 9 x -2.8044032112313464 y -2049.456302435744 -672.1354918462077\n",
            "it 10 x -2.268410685784332 y -672.1354918462077 -220.3982199780121\n",
            "it 11 x -1.8391438506400222 y -220.3982199780121 -72.1807419252783\n",
            "it 12 x -1.4963673552041776 y -72.1807419252783 -23.526420007144324\n",
            "it 13 x -1.2256340399442445 y -23.526420007144324 -7.546944476833993\n",
            "it 14 x -1.0192792149298004 y -7.546944476833993 -2.3017219833384583\n",
            "it 15 x -0.878192974042394 y -2.3017219833384583 -0.6001347172980255\n",
            "it 16 x -0.8072662373727413 y -0.6001347172980255 -0.09757282508715664\n",
            "it 17 x -0.7905468880551323 y -0.09757282508715664 -0.00434806118170894\n",
            "it 18 x -0.789729888238273 y -0.00434806118170894 -9.879715826832669e-06\n",
            "it 19 x -0.7897280233719703 y -9.879715826832669e-06 -5.135225578101199e-11\n",
            "it 20 x -0.7897280233622771 y -5.135225578101199e-11 5.551115123125783e-16\n",
            "x -20.0 => -0.7897280233622771  after  20 iterations\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "-0.7897280233622771"
            ]
          },
          "execution_count": 36,
          "metadata": {
            "tags": []
          },
          "output_type": "execute_result"
        }
      ],
      "source": [
        "def Newton(cs,xini,tolxrel=1.e-6,toly=1.e-6,maxit=1000):\n",
        "    xp = xpoly(xini)\n",
        "    xpp = xpoly_prime(xini)\n",
        "    x = xini\n",
        "    itnum = 0\n",
        "    while True:\n",
        "        y = np.dot(cs,xp)\n",
        "        yp = np.dot(cs,xpp)        \n",
        "        delta = y / yp\n",
        "        x += - delta\n",
        "        xp = xpoly(x)\n",
        "        xpp = xpoly_prime(x)\n",
        "        ynew = np.dot(cs,xp)\n",
        "        itnum += 1\n",
        "        print(\"it\", itnum, \"x\",x, \"y\", y,ynew)\n",
        "        if abs(delta/x) < tolxrel and abs(ynew) <toly:\n",
        "            break\n",
        "        if  itnum == maxit:\n",
        "            print(\"x not converged\")\n",
        "            return None\n",
        "    print(\"x\", xini, \"=>\", x, \" after \",itnum, \"iterations\")\n",
        "    return x \n",
        "\n",
        "x_initial = -20.0\n",
        "Newton(cs,x_initial)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lg6vMncbwyfa"
      },
      "source": [
        "$x=-20.0$から始めると、20回のiteration(反復)で、  \n",
        "yの値が$5.e-16 \\simeq 0$の点が求められている事がわかる。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_yFWkaZxxBqB"
      },
      "source": [
        "今のようにうまくいく例もある一方で、関数や初期値によっては解に収束しない場合があるので注意が必要"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dPK_KIGcyuod"
      },
      "source": [
        "# LICENSE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q943wB7Z4DYK"
      },
      "source": [
        "\n",
        "Copyright (C) 2021 Sota Yoshida\n",
        "\n",
        "[ライセンス:クリエイティブ・コモンズ 4.0 表示 (CC-BY 4.0)](https://creativecommons.org/licenses/by/4.0/deed.ja)"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "authorship_tag": "ABX9TyMWxHIsatxKmpa/oePk47i/",
      "collapsed_sections": [],
      "include_colab_link": true,
      "name": "Python_misc_NewtonsMethod.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3.9.13 64-bit",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.9.13"
    },
    "vscode": {
      "interpreter": {
        "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}