{
  "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_ODE.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wVdYiDLFHEEL"
      },
      "source": [
        "# 常微分方程式の数値解法"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "MO6-C742mcJ-"
      },
      "outputs": [],
      "source": [
        "#使用するモジュールのimport\n",
        "import numpy as np\n",
        "from matplotlib import pyplot as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JxV3JSWOd-W8"
      },
      "source": [
        "## 常微分方程式"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AGBSpCwiHIgY"
      },
      "source": [
        "多くの自然現象や社会現象は微分方程式として定式化される。\n",
        "\n",
        "身近な例としてはニュートンの運動方程式にはじまり、\n",
        "金融や群集動態、感染症のモデルなどにも当てはまる。\n",
        "\n",
        "このノートでは、微分方程式のうちとくに一変数関数の場合である  \n",
        "常微分方程式(ordinary differential equation, ODE)について  \n",
        "いくつかの代表的な数値解法を解説するとともに、  \n",
        "アルゴリズムを自分で実装してみよう。\n",
        "\n",
        "例えば1階の常微分方程式は、一般に以下の形で書くことが出来る:  \n",
        "$dx(t)/dt = f(x,t)$ (Eq. 1)\n",
        "\n",
        "これは、ある質点の時刻$t$での座標を$x(t)$と呼ぶことにすると、  \n",
        "質点の位置の変化が$f(x,t)$の関数として書けるような状況に相当する。\n",
        "\n",
        "注) 別に$t$は時間である必要はないし、  \n",
        "$x(t)$は座標以外の量(金融商品の価格など)何でも構わない。\n",
        "\n",
        "\n",
        "一般に常微分方程式が与えられた時、まずはその解が変数分離などを用いて  \n",
        "閉じた形で書き下せないかと考えるわけだが、  \n",
        "実際の多くの問題ではそうは上手く行かず、数値計算が必要となる。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "duBNxeB_cUMU"
      },
      "source": [
        "以下では、主として$dx/dt = f(x,t)$を考えるが、  \n",
        "高階の微分を含む場合も、\n",
        "\n",
        "$$\n",
        "\\left[\n",
        "\\begin{matrix} \n",
        "x_1 \\\\ \n",
        "x_2 \\\\\n",
        "\\vdots \\\\\n",
        "x_r\n",
        "\\end{matrix}\n",
        "\\right]\n",
        "= \n",
        "\\left[\n",
        "\\begin{matrix} \n",
        "x \\\\ \n",
        "dx/dt \\\\\n",
        "\\vdots \\\\\n",
        "d^{r-1}x/dt^{r-1}\n",
        "\\end{matrix} \\right]\n",
        "$$\n",
        "\n",
        "と置き換えれば、連立1階常微分方程式:\n",
        "\n",
        "$$\n",
        "\\left[\n",
        "\\begin{matrix} \n",
        "x_2 \\\\ \n",
        "x_3 \\\\\n",
        "\\vdots \\\\\n",
        "x_{r}\n",
        "\\end{matrix}\n",
        "\\right]\n",
        "= \n",
        "\\frac{d}{dx}\n",
        "\\left[\n",
        "\\begin{matrix} \n",
        "x_1 \\\\ \n",
        "x_2 \\\\\n",
        "\\vdots \\\\\n",
        "x_{r-1} \n",
        "\\end{matrix} \\right]\n",
        "$$\n",
        "とみなすことが出来る。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "brvU1jS2W0NK"
      },
      "source": [
        "\n",
        "常微分方程式の解を求めるためには、 \n",
        "適切な条件(初期条件や境界条件)を設定しなければならない。\n",
        "\n",
        "そのうち、変数$t$の最初の値(初期値,初期条件)$t_0$と  \n",
        "$x(t=t_0)$の値が与えられた問題を初期値問題と呼ぶ。  \n",
        "※$t_0$は問題によって0だったり、$-\\infty$だったりする。\n",
        "\n",
        "\n",
        "解の存在性や下記のリプシッツ条件など適切な条件の元での  \n",
        "解の一意性の証明等については教科書に譲る。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o2JwiqJAd2LL"
      },
      "source": [
        "###リプシッツ条件"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ipc62VyEeFXo"
      },
      "source": [
        "関数$f$が$U$(たとえば$R^n$の部分集合) 上でリプシッツ(Lipshitz)連続であるとは、  \n",
        "正実数$K$が存在し、任意の$x,y \\in U$に対し、  \n",
        "$||f(x)-f(y)|| \\leq K ||x −y||$が成立することである。  \n",
        "以下では$||\\cdot||$はユークリッドノルムとし、$K$をリプシッツ定数と呼ぶ。\n",
        "\n",
        "この$K$が存在しているとき、  \n",
        "前述の常微分方程式の初期値問題は一意の解を持つ。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3BMukvRtgdZI"
      },
      "source": [
        "## 代表定期な数値解法"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-pwbBdUBgjMr"
      },
      "source": [
        "### オイラー(Euler)法"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JrC4_IHYgu8-"
      },
      "source": [
        "1階の常微分方程式$dx(t)/dt=f(x,y)$を数値的に解く上で  \n",
        "最も基本的な考え方は時刻(便宜上そう呼ぶ)を  \n",
        "細かく分割した上で、$x$の任意時刻$t$での近似値を  \n",
        "初期値から初めて逐次的に求めていくことである。  \n",
        "\n",
        "以下では、その最も単純な例であるオイラー法を説明する。\n",
        "\n",
        "微分の定義:\n",
        "$dx(t)/dt = \\lim_{h\\to 0} \\frac{x(t+h) - x(t)}{h}$から  \n",
        "微小変分$h$を十分小さく取れば、  \n",
        "$x_{i+1} = x_{i} + h f(x,t_i)$という近似式が得られる。\n",
        "\n",
        "つまり、$i=0$の初期値$(x_0,t_0=0)$から順番に上の漸化式を用いて、  $i$番目の時刻$t_i$での座標$x_i$を推定していく。\n",
        "\n",
        "※一般に微小変分$h$は各ステップ$i$に依存しても良いが、簡単のため共通とした。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VgupCOzCjPwO"
      },
      "source": [
        "例: $dx/dt = -2x^2t$, 初期条件$x(t=0)=1.0$のもとで  \n",
        "$0 < t \\leq 1.0$での$x$の値を予想してみよう:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 225
        },
        "id": "lI-2SaBkHCXa",
        "outputId": "be62334d-2414-425b-a2a3-29a2efee59ff"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAADQCAYAAABP/LayAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXDc533f8fezi2OJG8QNLMD7AC+QBK3DsiRK1EFd1EHRsTxxbdeJ2mmcpkmbjj3N1K4zHtvJpLHTpHEVRbWTaeKE1EHqpG7LliVZBEmQBCge4oVd3Adxn7vf/rEQTVEkBEpc/ADs5zWDAX67j374Yp4B9cHzPL/ncWaGiIiIiEwtn9cFiIiIiCQihTARERERDyiEiYiIiHhAIUxERETEAwphIiIiIh5QCBMRERHxQJLXBVyu/Px8mz9/vtdliIiIiHysmpqadjMruNh7My6EzZ8/nz179nhdhoiIiMjHcs6dvtR7mo4UERER8YBCmIiIiIgH4hbCnHOPOedanXOHLvG+c879lXPuuHPugHNufbxqEREREZlu4jkS9hNg8wTv3wEsGf94GPjbONYiIiIiMq3EbWG+mb3hnJs/QZN7gX+w2AnibzvncpxzJWbWFK+aJuNQuJuf/OoUfufw+Rx+H+e+9jmH/9xn8Pt8BJJ9pCb5CST7CCT5CST7SU3yEUj2n3tvToqfrEASGYEk5iT7cc55+SOKiIjINODl05FlQMN516Hx1z4SwpxzDxMbLaOioiKuRXX0j/Cr4+1EzIhEIWpGJGpEzYhGjYgZ0Sjj79tl39/vc2SkJpGRmkRm4IOPZDJSk8iek0xuegp56SnkpqcwNy2F3PRk5qankJuWQiDZH4efWERERLwwI7aoMLNHgEcANmzYcPnJ5zLcuLSAX31z06TaRqPGSCTK8GiUobEIQ6MRhkajDI/FPseuIwyORugdGqNveIzeoVH6hsboHR6LvTY0RmvvEO+3jdEzOMrZwVHsEj9heoqf3PQU8jNSKcpKpTAzEPucFaAoK/Z1UWaAnLRkjbaJiIhMc16GsDBQft51cPy1GcPncwR8sSnIbJKvyD3HIlG6B0fpGhihs3+Uzv5hOvs/uI59tPcNc7K9n7dPdNI9OPqRe6T4fRRkplKSHSCYO4ey3DkEc9Moy5lDMHcOpTlzNKomIiLiMS9D2C7g6865nwFXA91erwebDpL8PvIyUsnLSJ1U+6HRCK09w7T2DtHSM0xLzxAtvUO0dA/R1D3EntNdPH2g6SNTpwWZqbGAljOH8rlpzM9LY0F+BvPz0yjISNVImoiISJzFLYQ55/4Z2AjkO+dCwLcgNlxkZj8GngPuBI4DA8BX41XLbBZI9lORl0ZFXtol24xForT0DhPqHCB8dpBQ1yDhrkHCZwc5FO5md10zo5HfhLT0FD/z89OZn5/Ogrzxz/mxkDY3PWUqfiwREZFZL55PRz70Me8b8Hvx+v7yG0l+H2U5sVGvixmLRGk8O8TJjn5Otfdzsr2fUx391IW7eeFQ84dG0eamp7C4MIOlRRksLcoc/zqTvPQUjZ6JiIhchhmxMF/iK8nvOzeaduPSD58xOhqJEuoa5FR7P++39fF+Wx9HW/rYub+R3qGxc+1y05JZUpTJkvFQtrQokxUlWWSnXZm1ciIiIrONQphMKNnvY0F+Ogvy07lpeeG5182M1t5hjrb0cqylj2OtvRxt6WNX7YfDWVnOHFaUZrGiJOvc52DuHI2aiYhIwlMIk0/EOTe+LUaA65f8ZvTsg3D2XnMv9Y091Df1UN/YzcuHW85tvZEZSGJFSRaV48FsdVk2SwozSPLrKFMREUkcCmFyRZ0fzs6f2hwYGeNIcy/1TT0cbuqhvrGHf3m3gcHRCACBZB+rSrNZE8yhqjybqmAO8/LSNGImIiKzlrNL7Qw6TW3YsMH27NnjdRlyBUSixqmOfg6GuqkNneVAqJtD4W6Gx6IAZM9JZk0we/wjh7XlORRlBTyuWkREZPKcczVmtuFi72kkTDzj9zkWFWSwqCCD+9aVAbEnNY+29I2HsrPUNnTz45+fOPeEZkl2gPXzctkwL5fqeblUlmSRrGlMERGZgTQSJtPe0GiEusYeahvOsq/hLDWnOmnsHgJgTrKfqvJsqsdD2fqKXHLStJeZiIhMDxoJkxktkOw/F7I+0Hh2kL1nuthzqou9Z7o+NFq2uDDj3EjZNQvz9DSmiIhMSxoJk1lhYGSM2obu8WDWSc3pLnrGt8oozQ5w9cI8rlk4l6sX5GnBv4iITBmNhMmsl5aSxLWL8rh2UR4A0ahxrLWPd0528M6JTt442saT+2LnwxdlpXLNwjyuXhALZgvy0xXKRERkymkkTBKCmfF+Wx9vnejknRMdvHOyk7beYSB2mPnVC+Zy7aI8rl9cMOE5nCIiIpdjopEwhTBJSGbGifZ+3jnRyTsnO3j7RActPbFQVj53Dp9bnM91i/P57KJ8HVouIiKfmEKYyMf4IJS9ebydXx5r5633O+gdjq0pW1maxecW5/O5Jfl8Zv5cAsl+j6sVEZGZQiFM5DKNRaIcDHfzy2Pt/PJ4O3vPdDEaMVKSfGyYl8t1i/O5fkk+q0qz8fm0nkxERC5OIUzkUxoYGePXJztjI2XHOzjc1ANAXnoKNywt4MalBdywtEBTlyIi8iF6OlLkU0pLSWLjskI2LisEoL1vmF8ea+f1I638fPzJS+dgTTCHjUsL2LisgDXBHPwaJRMRkUvQSJjIpxSNGgfD3bx+pI3Xj7ayv+EsZpCblswN44Hs+iUF5Gekel2qiIhMMU1Hikyhrv4R3jjWxs+PtvHG0Tba+0ZwDlaXZbNxWSG3VBZqLZmISIJQCBPxSDRq1DX28PqRVl4/2sa+M11EDQozU9lUWcQtlYVctzhfT1yKiMxSCmEi00Rn/wivH2nl5cMtvHG0nb7hMQLJPj63OJ9NlUVsWl5IYVbA6zJFROQK8SyEOec2Az8C/MCjZvb9C96fBzwGFACdwG+bWWiieyqEyWwxMhblnZMdvFzfwsuHWwmfHQSgKpgdC2SVhawoydKRSiIiM5gnIcw55weOArcCIeBd4CEzqz+vzXbgGTP7qXPuZuCrZvalie6rECazkZlxpKWXVw638lJ9C7Wh2OL+0uwAmyqLuG1lEdcszCPZ7/O6VBERuQxehbBrgW+b2e3j198EMLPvndemDthsZg0u9ud+t5llTXRfhTBJBG29w7z2Xmza8hfH2hkcjZAVSOKWyiJuW1nMjUsLmJOidWQiItOdV/uElQEN512HgKsvaFMLPEBsyvJ+INM5l2dmHXGsS2TaK8hM5fOfKefznylnaDTCG0fb2F3XwsuHW3hiX5hAso8blxZw+8piNi0vIjst2euSRUTkMnm9Wet/Af7aOfcV4A0gDEQubOScexh4GKCiomIq6xPxXCDZz20ri7ltZTGjkSi/PtnJ7rpmXqxrYXddC0k+x7WL8rh9ZTG3rSjSwn4RkRnC0+nIC9pnAO+ZWXCi+2o6UiQmGjVqQ2fZXdfC7rpmTrb34xysK89h86pibl9ZzLy8dK/LFBFJaF6tCUsitjB/E7ERrneBL5pZ3Xlt8oFOM4s6574LRMzsv090X4UwkY8yM4619rH7UDMv1DVT1xg727KyJIu7Vhdz5+oSFhZkeFyliEji8XKLijuBHxLbouIxM/uuc+47wB4z2+WcexD4HmDEpiN/z8yGJ7qnQpjIx2voHGB3XTPPH2qm5nQXoEAmIuIFbdYqksCaugd57mAzzx1sOhfIlhdncveaEgUyEZE4UwgTEUCBTERkqimEichHKJCJiMSfQpiITOhSgeyeqlK2VJVSPjfN4wpFRGYmhTARmbSm7kGeP9jMMwca2XvmLABry3PYUlXK3WtKtA+ZiMhlUAgTkU+koXOAZw40sau2kcNNPfgcXLMwjy1VpdyxqkQ79YuIfAyFMBH51I639rJrfyO7ahs51TFAst9x49IC7qkq5ZbKItJTvT6AQ0Rk+lEIE5Erxsw4FO5hV22Yp2ubaO4ZYk6yn02VhWypKuXGZQWkJulwcRERUAgTkTiJRo13T3Wyq7aR5w420TUwSlYgic2ritlSVca1i/Lw+5zXZYqIeEYhTETibjQS5ZfH23m6tpEX61roGx4jPyOVe6pKuH9dGavLsnFOgUxEEotCmIhMqaHRCK+918rO/Y28+l4rI5EoCwvSuX9tGfetK9OWFyKSMBTCRMQz3YOjPH+wiSf3hXnnZCcAG+blct+6Mu5eU0JOWorHFYqIxI9CmIhMC+Gzg+zcH+bJvWGOtfaR7HfctKyQ+9eVcdPyQgLJWtAvIrOLQpiITCtmRn1TD0/tC7NzfyOtvcNkBpK4a3UJ960r46r5c/FpQb+IzAIKYSIybUWixlvvd/DEvhC7DzXTPxKhNDvAvevKuH9dGUuLMr0uUUTkE1MIE5EZYWBkjJfqW3hqX5g3jrUTiRorSrJ4YH0Z91SVUqQjk0RkhlEIE5EZp71vmGdqG3lyfyO1DWfxObhucT4PrC/j9pXFpKVoh34Rmf4UwkRkRjvR1sdT+8I8uT9MQ+cg6Sl+7lxdwgPrg1y9QOvHRGT6UggTkVnhgx36n9gb5tmDTfQNjxHMncMD68p4YH2Q+fnpXpcoIvIhCmEiMusMjkR4sb6ZHTUh3jzeTtSgel4uW9cHuWtNCdlzkr0uUUREIUxEZrfm7iGe2h/m8ZoQx1r7SEnycduKIrauD3L9knyS/D6vSxSRBKUQJiIJwcw4GO7m8ZoQu2ob6RoYpSAzlfvWlvLA+iCVJVlelygiCcazEOac2wz8CPADj5rZ9y94vwL4KZAz3uYbZvbcRPdUCBORyRgZi/LakVYerwnx2pFWRiOx7S62Vge5d20p+RmpXpcoIgnAkxDmnPMDR4FbgRDwLvCQmdWf1+YRYJ+Z/a1zbgXwnJnNn+i+CmEicrk6+0d4uraRx/eGOBDqxu9zbFxawNbqIJsqC0lN0nFJIhIfE4WweG60cxVw3MxOjBfxM+BeoP68NgZ8MD+QDTTGsR4RSVBz01P48mfn8+XPzudoSy+P7w3x1L4wr7zXSvacZO5eU8LW6iDrynNwTttdiMjUiOdI2IPAZjP7nfHrLwFXm9nXz2tTArwI5ALpwC1mVnORez0MPAxQUVFRffr06bjULCKJIxI13jzezuN7Q+yua2ZoNMrCgnQerA6ydX1Qu/OLyBXh1XTkZELYH43X8BfOuWuBvwdWmVn0UvfVdKSIXGm9Q6M8d7CJHTUh3j3Vhc/BDUsLeLA6yC2VRQSSNV0pIp+MV9ORYaD8vOvg+Gvn+xqwGcDM3nLOBYB8oDWOdYmIfEhmIJnf+kwFv/WZCk6297OjpoEn9ob5+j/tI3tOMveuLeXB6iCry7I1XSkiV0w8R8KSiC3M30QsfL0LfNHM6s5r8zzwL2b2E+dcJfAKUGYTFKWRMBGZCh9MV26viU1XjoxFWVaUybYNQe5bV6anK0VkUrzcouJO4IfEtp94zMy+65z7DrDHzHaNPxH5d0AGsUX6/9XMXpzongphIjLVugdHebq2ke01IWobzpLkc2xcVsi2DUFuXl5IsjaDFZFL0GatIiJXyLGWXnbUhHhiX5i23mHy0lO4d20Z2zZoM1gR+SiFMBGRK2wsEuWNY21s3xPi5cMtjEaMVWVZbKsuZ0tVKbnpKV6XKCLTgEKYiEgcdfaPsHN/mB01Ieoae0jx+7hlRSHbqst1dqVIglMIExGZIvWNPWyvaWDn/kY6+0cozEzl/vVlbKsuZ3FhhtflicgUUwgTEZliI2NRXn2vlR01Dbx2pI1I1FhXkcOD1UHuqSolK5DsdYkiMgUUwkREPNTWO8xT+8Jsr2ngaEsfqUk+Nq8q5sHqIJ9dlI/fp73HRGYrhTARkWnAzDgY7mb7nhA794fpGRqjNDvA1vGjkubnp3tdoohcYQphIiLTzNBohJcPt7B9T4hfHGsjanDV/Lk8uCHIXatLSE+N54EmIjJVFMJERKax5u4hHt8b4vGaECfa+0lL8XPHqhK2bQhy1fy5+DRdKTJjfaoQ5pxbYWb1F7y20cxev3IlTp5CmIjMVmbG3jNdbN8T4pkDTfQNj1ExN42t64NsrS4jmJvmdYkicpk+bQg7BPwj8GdAYPzzBjO79koXOhkKYSKSCAZHIrxQ18SOmhBvHu8A4LOL8ti2IcjmlSXMSfF7XKGITManDWHpwA+AaiAT+H/AD8wseqULnQyFMBFJNKGuAR6vCbNjbwMNnYNkpCZx95rYdOX6ilyc03SlyHQ1UQibzMrPUWAQmENsJOykVwFMRCQRBXPT+INblvD7Ny/m16c6x5+ubORn7zawMD/93NOVxdkBr0sVkcswmZGwWmAn8KdAPvBjYMTMtsW/vI/SSJiICPQNj/Hcwdh05a9PduJz8LklBWyrDnLriiICyZquFJkOPu105AYz23PBa18ys3+8gjVOmkKYiMiHne7o5/GaEI/vDRM+O0hWIIkta0vZVl3OmmC2pitFPKQtKkREEkA0avzq/Q521DTw/KFmhseiLC3K4MHqIPetK6MwU9OVIlNNIUxEJMH0DI3y7IEmtu9pYO+Zs/h9jo1LC9i2IcjNy4tISfJ5XaJIQlAIExFJYMdb+3h8b4gn9oZo6RkmNy2Ze9eWsW1DkJWl2V6XJzKrKYSJiAhjkSi/ON7OjpoQL9W1MBKJUlmSxbbx6cq56Slelygy6yiEiYjIh5wdGOHp2ka214Q4EOom2e+4eXkh26rLuXFZAcl+TVeKXAkKYSIicklHmnvZUdPAk/vCtPeNkJ+Ryv3rStm2oZylRZlelycyoymEiYjIxxqNRPn5kTa21zTwyuFWxqJGVTCbB6uDbKkqIzst2esSRWYcz0KYc24z8CPADzxqZt+/4P2/BG4av0wDCs0sZ6J7KoSJiMRfR98wO/fHpisPN/WQ4vdx68oitlUHuX5JAX6f9h4TmQxPQphzzg8cBW4FQsC7wENmVn+J9r8PrDOzfzvRfRXCRESm1qFwNztqQuzcH6ZrYJSirFQeWB/kweogiwoyvC5PZFrzKoRdC3zbzG4fv/4mgJl97xLtfwV8y8xemui+CmEiIt4YHovw2nutbN8T4vWjbUSixvqKHLZtKOfuNSVkBjRdKXIhr0LYg8BmM/ud8esvAVeb2dcv0nYe8DYQNLPIRd5/GHgYoKKiovr06dNxqVlERCantXeIp/aF2b4nxLHWPgLJPjavLGbbhnKuXZiHT9OVIsDEISxpqou5hC8AOy4WwADM7BHgEYiNhE1lYSIi8lGFmQEevmERv3v9Qg6Eutle08Cu/Y08tb+Rspw5bF1fxoPV5VTkpXldqsi0Fc8QFgbKz7sOjr92MV8Afi+OtYiISBw456gqz6GqPIc/uWsFL9W3sL0mxP967Th/9epxrlowl23VQe5cXUJ66nT5u19keojndGQSsYX5m4iFr3eBL5pZ3QXtlgMvAAtsEsVoTZiIyPTX1D3IE3vD7KgJcbK9n7QUP3esKmHr+jKu0XSlJBAvt6i4E/ghsS0qHjOz7zrnvgPsMbNd422+DQTM7BuTuadCmIjIzGFm7D3TxfY9IZ490ETv8Bil2QHuX1/GA+v1dKXMftqsVUREPDc0GuGl+hae2BvijWPtRKJGVXkOD64v4+41peTq7EqZhRTCRERkWmntHWLX/kYe3xvmcFPPubMrH1gf5KZlhaQk6exKmR0UwkREZNqqb+zhyX0hntzXSHvfMLlpyWypKuWB9UHWBLNxTuvHZOZSCBMRkWlvLBLlF8fbeWJvmBfrmhkei7K4MIMH1pdx39oySnPmeF2iyGVTCBMRkRmlZ2iU5w408cTeML8+1Ylz8NlFeWxdH+T2lcXa7kJmDIUwERGZsc50DPDEvhBP7A1zpnOAtBQ/m1cVs3V9kGsW5ukwcZnWFMJERGTGMzNqTnfx+N4QzxxoondojJLsAPevi213sbhQ213I9KMQJiIis8rQaISXD7fwxN4wPx8/TLwqmM3W6iD3aLsLmUYUwkREZNZq6x1m5/4wT+wNUz++3cXGZYU8sK6Mm5YXEkj2e12iJDCFMBERSQiHm3p4Ym+Infsbae0dJjOQxF2rS7h3bRlXL5ir45JkyimEiYhIQolEjbfe7+DJfWFeONRE/0iE0uwAW9aWcf+6MpYVZ3pdoiQIhTAREUlYgyMRXjrcwlP7wrxxtI2xqLG8OJP715WxZW0pJdnaf0ziRyFMREQE6Ogb5tmDTTy5L8y+M2dxDq5ZkMf968rYvLqYrECy1yXKLKMQJiIicoFT7f08tT/Mzv2NnGzvJyXJx62VRdy3rowblxbo/Eq5IhTCRERELsHMqA1189S+ME/XNtLRP0JOWjJ3rS7h/nVlVM/L1fmV8okphImIiEzCaCTKL4+189T+MLvrmhkajVI+dw73rS3j3rVl2hBWLptCmIiIyGXqGx7jxbpmntwX5s3j7UQNVpdlc9+6Mu6pKqEwM+B1iTIDKISJiIh8Cq09Q+yqbWTn/kYOhrvxOfjsony2VJVy+6pisudoQb9cnEKYiIjIFXK8tZed+xvZVdvI6Y4BUvw+blxWwJaqUm6pLGJOinbol99QCBMREbnCzIwDoW521TbydG1sh/60FD+3riji3rWlfG6xnrAUhTAREZG4ikSNd0528HRtI88dbKZ7cJSctGTuWFXClqpSrlowF7+OTEpInoUw59xm4EeAH3jUzL5/kTafB74NGFBrZl+c6J4KYSIiMp2NjEX5xbE2dtU28lJ9CwMjEYqyUrl7TSlbqkpZE8zWlhcJxJMQ5pzzA0eBW4EQ8C7wkJnVn9dmCfCvwM1m1uWcKzSz1onuqxAmIiIzxcDIGK8cbmVXbSM/P9LGSCTK/Lw07qmKBbIlRTrDcrbzKoRdC3zbzG4fv/4mgJl977w2fwYcNbNHJ3tfhTAREZmJugdG2V3XzM7aMG+930HUYHlxJveujW15EcxN87pEiYOJQlhSHL9vGdBw3nUIuPqCNksBnHNvEpuy/LaZvRDHmkRERDyRnZbM5z9Tzuc/U05r7xDPHmhiV20jP3jhPX7wwntUz8tlS1Upd64uoSAz1etyZQrEM4RN9vsvATYCQeAN59xqMzt7fiPn3MPAwwAVFRVTXaOIiMgVVZgZ4KvXLeCr1y2goXPg3BOW39pVx/94uo6rF+Rx15oS7lhVTF6GAtls5fV05I+Bd8zs/45fvwJ8w8zevdR9NR0pIiKz1dGWXp450MQzBxo50daP3+e4dmEskG1eWUxueorXJcpl8mpNWBKxhfmbgDCxhflfNLO689psJrZY/8vOuXxgH7DWzDoudV+FMBERme3MjPeae3l2PJCd6hjA73Nctzifu1eXcPvKYrLTtEv/TODlFhV3Aj8ktt7rMTP7rnPuO8AeM9vlYs/o/gWwGYgA3zWzn010T4UwERFJJGZGXWMPzx5s4tkDTZzpHCDZ7/jc4nzuWlPKrSuKdGzSNKbNWkVERGYBM+NguHt8hKyJ8NlBUvw+bliaz11rSrilsojMgALZdKIQJiIiMsuYGfsbzvLsgSaePdhEU/cQKUk+Ni4t4K41JWyqLCIj1evn70QhTEREZBaLRo19DWd55kAjzx1soqVnmNQkHzcvL+SuNSXcvLyQtBQFMi8ohImIiCSIaNSoOdN1boSsrXeYQLKPjUsLuWN1MTcvL9SU5RRSCBMREUlAkajx65OdPH+oiRcONdPaO0yK38f1S/K5Y3UJt1YW6SnLOFMIExERSXCxKcsunjvYzAuHmgmfHSTJ57h2UR53ri7hthVF2hg2DhTCRERE5Bwz40Com+cPNfP8oSZOdwzgc3D1gjzuWF3M7SuLKcoKeF3mrKAQJiIiIhdlZhxu6uX5Q008f6iZ4619OAfVFblsXlXMHatLKMuZ43WZM5ZCmIiIiEzKsZbe8RGyZg439QBQFczmjtWxsyzn5aV7XOHMohAmIiIil+1Ue/+5KcsDoW4AKkuyuHNVMZtXFbO4MIPY4TdyKQphIiIi8qk0dA6wuy42QlZzuguABfnp3LayiNtWFLOuPAefT4HsQgphIiIicsW09AzxYn0LL9Y189b7HYxFjYLMVG5dUcTtK4u5dmEeKUk+r8ucFhTCREREJC66B0d5/Ugru+uaef1IGwMjETJTk7hpeSG3rSxi47LChD4+SSFMRERE4m5oNMKbx9t5sa6Flw+30NE/Qorfx3WL87h9ZTGbKosoyEysvcgmCmGJG01FRETkigok+9lUWcSmyiIiUaPmdBe765rZXdfMa0cO4txBNszL5bYVsb3IKvLSvC7ZUxoJExERkbj6YC+yF+ub2V3Xcm7ri+XFmdy2spjbVhSxsjRrVj5pqelIERERmTY+eNLyxfoW9pzqJGpQljOHWyoLuWVFEVcvmD0L+xXCREREZFrq6Bvm5cMtvFTfyi+PtzE0GiUjNYkblxawqbKQm5YVkpue4nWZn5hCmIiIiEx7gyMRfvV+Oy8fbuGVw6209g7jc7Bh/tzYKFllEQsLMrwu87IohImIiMiMEo0aB8PdvHy4hZcPt55bR7YwP51bVhSxaXkh1fNySfJP72lLhTARERGZ0UJdA7xyuJWXD7fw9okORiNGTloyNy8rZFNlETcszSczkOx1mR+hECYiIiKzRu/QKL841s7L9S28eqSVswOjJPsd1yzM45bKIjZVFhLMnR7bX3gWwpxzm4EfAX7gUTP7/gXvfwX4cyA8/tJfm9mjE91TIUxEREQ+MBaJsvfM2fFpyxZOtPUDse0vbl5eyKbKQtaW5+L36FxLT0KYc84PHAVuBULAu8BDZlZ/XpuvABvM7OuTva9CmIiIiFzKibY+XjncykuHW6g53UUkauSmJbNxWSE3LS/kxiUFZKdN3bSlVzvmXwUcN7MT40X8DLgXqJ/wvxIRERH5hBYWZLCwIIPfvWEh3QOjvHGsjVffa+X1I608uS+M3+eonpfLzcsL2byymPn56Z7VGs8QVgY0nHcdAq6+SLutzrkbiI2a/aGZNVzYwDn3MPAwQEVFRRxKFRERkdkmOy2Ze6pKuaeqlEjU2N/QxSwpltsAAAXKSURBVKvvtfLqe218//n3APj3Ny7yrD6vz458GvhnMxt2zv074KfAzRc2MrNHgEcgNh05tSWKiIjITBcbAZtL9by5/PHty2k8O0gg2e9pTfHcXCMMlJ93HeQ3C/ABMLMOMxsev3wUqI5jPSIiIiIAlObMYa7HO/HHM4S9Cyxxzi1wzqUAXwB2nd/AOVdy3uUW4HAc6xERERGZNuI2HWlmY865rwO7iW1R8ZiZ1TnnvgPsMbNdwH90zm0BxoBO4CvxqkdERERkOtFmrSIiIiJxMtEWFdP7wCURERGRWUohTERERMQDM2460jnXBpyO87fJB9rj/D3k8qlfph/1yfSkfpl+1CfT01T0yzwzK7jYGzMuhE0F59yeS83finfUL9OP+mR6Ur9MP+qT6cnrftF0pIiIiIgHFMJEREREPKAQdnGPeF2AXJT6ZfpRn0xP6pfpR30yPXnaL1oTJiIiIuIBjYSJiIiIeCChQ5hzbrNz7ohz7rhz7hsXeT/VOfcv4++/45ybP/VVJp5J9MsfOefqnXMHnHOvOOfmeVFnIvm4Pjmv3VbnnDnn9BTYFJhMvzjnPj/++1LnnPunqa4x0Uzi368K59xrzrl94/+G3elFnYnEOfeYc67VOXfoEu8759xfjffZAefc+qmqLWFDmHPOD/wNcAewAnjIObfigmZfA7rMbDHwl8APprbKxDPJftkHbDCzNcAO4M+mtsrEMsk+wTmXCfwB8M7UVpiYJtMvzrklwDeB68xsJfCfprzQBDLJ35U/Af7VzNYBXwD+99RWmZB+Amye4P07gCXjHw8DfzsFNQEJHMKAq4DjZnbCzEaAnwH3XtDmXuCn41/vADY559wU1piIPrZfzOw1MxsYv3wbCE5xjYlmMr8rAH9K7A+VoaksLoFNpl9+F/gbM+sCMLPWKa4x0UymTwzIGv86G2icwvoSkpm9AXRO0ORe4B8s5m0gxzlXMhW1JXIIKwMazrsOjb920TZmNgZ0A3lTUl3imky/nO9rwPNxrUg+tk/Gh+/LzezZqSwswU3md2UpsNQ596Zz7m3n3ESjAfLpTaZPvg38tnMuBDwH/P7UlCYTuNz/71wxSVPxTUTiwTn328AG4Eava0lkzjkf8D+Br3hcinxUErEplo3ERozfcM6tNrOznlaV2B4CfmJmf+Gcuxb4R+fcKjOLel2YTL1EHgkLA+XnXQfHX7toG+dcErGh444pqS5xTaZfcM7dAvw3YIuZDU9RbYnq4/okE1gFvO6cOwVcA+zS4vy4m8zvSgjYZWajZnYSOEoslEl8TKZPvgb8K4CZvQUEiJ1fKN6Z1P934iGRQ9i7wBLn3ALnXAqxBZK7LmizC/jy+NcPAq+aNlaLt4/tF+fcOuD/EAtgWuMSfxP2iZl1m1m+mc03s/nE1ultMbM93pSbMCbzb9hTxEbBcM7lE5uePDGVRSaYyfTJGWATgHOuklgIa5vSKuVCu4B/M/6U5DVAt5k1TcU3TtjpSDMbc859HdgN+IHHzKzOOfcdYI+Z7QL+nthQ8XFii/q+4F3FiWGS/fLnQAawffw5iTNmtsWzome5SfaJTLFJ9stu4DbnXD0QAf7YzDSaHyeT7JP/DPydc+4PiS3S/4r+uI8v59w/E/tjJH98Ld63gGQAM/sxsbV5dwLHgQHgq1NWm/peREREZOol8nSkiIiIiGcUwkREREQ8oBAmIiIi4gGFMBEREREPKISJiIiIeEAhTEQSmnMuxzn3H7yuQ0QSj0KYiCS6HEAhTESmnEKYiCS67wOLnHP7nXN/7nUxIpI4tFmriCQ059x84BkzW+VxKSKSYDQSJiIiIuIBhTARERERDyiEiUii6wUyvS5CRBKPQpiIJDQz6wDedM4d0sJ8EZlKWpgvIiIi4gGNhImIiIh4QCFMRERExAMKYSIiIiIeUAgTERER8YBCmIiIiIgHFMJEREREPKAQJiIiIuIBhTARERERD/x/8OZ6MjjlM9gAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 720x216 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "# f(x,t)のトイモデルとして dx/dt = -2tx^2 というものを考える\n",
        "def f(x,t): \n",
        "    return -  2.0*  x * x * t\n",
        "#時刻tを間隔hで細かく分割する\n",
        "h = 1.e-4\n",
        "tr = np.arange(0.0,1.0,h) \n",
        "#初期条件\n",
        "x0 = 1.0 \n",
        "\n",
        "#求解\n",
        "x = x0\n",
        "xs = [x0]\n",
        "for t in tr:\n",
        "    x += h * f(x,t)\n",
        "    xs += [x]\n",
        "\n",
        "#描画\n",
        "fig = plt.figure(figsize=(10,3))\n",
        "plt.xlabel(\"t\");plt.ylabel(\"x\")\n",
        "plt.plot(tr,xs[:-1])\n",
        "plt.show()\n",
        "plt.close()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_xzmhI0Rq7A_"
      },
      "source": [
        "実は今の初期条件のもとでこの微分方程式は  \n",
        "閉じた形$x(t)=x_0/(x_0t^2+1)$で解が与えられるので、  \n",
        "真の解と数値解法による近似解を比較できる。\n",
        "\n",
        "差分のlogを取ってみると...以下のようになる:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 225
        },
        "id": "dMy_G8hjq6i-",
        "outputId": "acdfa5e0-6bf7-4e04-922a-8d9885e98d63"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAADQCAYAAAC3HE1FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZCcd33n8fe3u6fnvmc0ujWyRpYs2WCbsWNMHBscjDEIJ+EIBnJh8JKFEDYhtbBLLSRAmWwqpEIgh3ZxHJLFJjhUyrdNuBRjgy3ZxpZsWZZlHaNz7nump7u/+8fTI41kaaYlzdPd0/15VXU93U8/ep6v9Xikj37XY+6OiIiIiBSeSL4LEBEREZHTU1ATERERKVAKaiIiIiIFSkFNREREpEApqImIiIgUKAU1ERERkQIVy3cBYWhpafH29vZ8lyEiIiIyp23btvW4e+vpvivKoNbe3s7WrVvzXYaIiIjInMxs35m+U9eniIiISIEqqqBmZpvMbPPg4GC+SxERERE5b0UV1Nz9Pne/rb6+Pt+liIiIiJy3ogpqIiIiIsWkKCcTiJS6ZCrNRDLNxFRqxitNIpVmKpkmmXamUmmmUk4yFexPppxkOk0is2/6+6lUmlTaSbuTdki74w4+y+e0AzjpdPAZIGJGJAJmRsTAyGzNiJhhBhEj8376u8yvM6MsasSiEcqikeB9JEIsGuwvi0aIRSInjokYZbEIsUjmu6gRj0aoKItSHjuxjUX1b1URKWwKaiJ5lE47o4kko5MpRiaTjEwmGZ1MMjwRbEcTJ97P/H4skWJyKs1EMsV4IsVEMghiE5n3Uymf1zrNIDojUM0MUGZgQCRimX1AJoRNfw/gM0JcEOocz/wenPQdM485sQ1DNGJUxCKUl0WPb8tP+XzyNkJ1PEZVPEZVPEpVeZTqeIzK+IztKfsqyiLY9G+CiMhZUlATmQeptNM/lqB/NMHA+BQDY1MMjk8xMJbIbKcy+4PP0/uGJqbwLEJINGLUlMeOvyriUSrLIjRVx6moj1IZj1JRFqE8lnkfi1IZD1qOKmJRKuJB0CiLRYhHMy1NsQhlx1ulIidapjKtT7HoiRapaCS/QSOddpLpoMVvKulMpYMWwKnUzNbBE62C0y2BM49JJNNMZgLtzO104A22aSanUse3Q+NTHJtKkZhunUymGUskmZhKZ127GVSVRakqj1Edj1IZD7Y1FTFqK8qorYgFr/ITn2tmvK+rKMscG6NMLYAiJUdBTeQMEsk0x4Yn6B1J0DMySe9Igu7Mtmdkkt7RSXqGE/SOTtI3mjhjq0/EoK6yjIbKMuqr4jRWxVndUk19ZRn1lcFfxtUzQlhN+YzPmb+0y2Ol3SoTiRjxiBEnAvF8VxME8/GpFGOTSUYTKcYSQSvn6GSS8UTqpH0njgn2jU4G277RBPt6xxiemGJoIkkiOXf4qyiLUFNeRl0muNVUxDL/H8VpqAr+H2uoKjt5X1UZDZVxteyJLFAKalKSRiaTHBmcCF5DExwZHM9sJ45ve0YSp/211fEozTXltNTEWdlcxeWrGmmpidNSU05jdfz4X5YNlXHqq8qoLY8RyXOLlMyvmS2c82UymWJkIujqDl5TDE/OeJ/ZjkwmGZpxzJHBieMttMlZ+ojj0Qj1J4W5+PFQN72vsTpOU3Wc5upymqrjNFaVaRyfSJ4VVVAzs03Apo6OjnyXInk2OpnkQP8YXX3jHOgf40Bm29U/Tlf/GMMTydf8moaqMhbXVbC4voKLl9azuL6CtroKWmvKac4EsZaacirj0Tz8F0mxK49FKa8J/hFwLtydsUTqRBd7prv9RNd7Zl+mW/7gwDgvHBpkYHyKsUTqjOdtqCqjqTpOU1UmxNUE26bqcpqrp9+feFWU6edDZD6ZZzNAZoHp7Ox0PUKq+A2OT/Fqzyh7ukeCbc8oXX1jHOgfp2/05NawyrIoK5oqWdFYxfLGSpY0VB4PZdNb/QUjpSqRTDMwnqB/dOp4V37faILekcTx98f3jSboH0uQOkPr3XSLc2ttOa015bTUxmmtqQg+z3i11MQpj+lnTgTAzLa5e+fpviuqFjUpPqm0s79vjJePDmdC2WgmlI2c1DUZMVjeWMWq5iretrSeFU2VLG+sYkVjJSuaqmiujmt8jsgZxGMRFtVWsKi2Aqid8/h02hmamApCWya8zQx3vaOT9IxMsqdnhJ+/Okn/2NRpz1NfWXY8tLXWBq3XJwW6mnIW1ZXTVBXX8AEpWQpqUhDcne6RSV46MnzidXSYXUeHT5ph11ITDMS/fn0bq1uruaClmgtaq1nZVE08prE0IrkQiRgNVXEaquLQOvfxiWSa3tFJuoeDANc9POOV+fx81wDdw5OMnqYbtixqLKo90QLeVlfB4vryYJtpEW+rU6u4FCcFNck5d2df7xjPHxwMXl2DvHR0+KTuypaactYtruEDV65i/eJa1rbVcEFrDfWVZXmsXETORTwWYUl9JUvqK+c8dnQyeVKYOzo0wZGhzHZwghcPD/Gjl46ddlzd9DjT6QDXlgl2S+qDMLe0oZK6ipha12VBUVCTULk7Xf3jPNc1yHMHB9ieCWZDmcH88WiE9UtqeetFbaxbXMv6xbWsW1x7zgOqRWRhq84sT7OqufqMx7g7w5NJjs6YpT09YzsIdhPsODRE7+jka9YprCmPsayhkqUNQXBb2lDJ8sbK4+/bass101UKioKazKupVJoXDg2xdV8/W/f2sXVfP93Dk0DQfbF+cR3vfP1SLllWzyXL6rmwrVZdliJyVsyMuooy6irKWNt25jF1U6k0x4YnOTI4weHBcQ4PTHBwYJyDA+McGhjn2QMDrxk/F40Yi+sqjge5ZZkAt2xGqKuex2VZROai/9vkvExMpdi2r58nXull674+fnFgkPGpoEtieWMlb1rTzBvam7h0eQMXLq7RLC8RyZmyaIRlmZAFjac9ZnQyyeHBcQ4OTHBoYJyD/UGIOzgwztP7+3ngucOvWZ+uuTrOiqYqVjRVsbKpkpVNVaxoDD4vqa9Qi5zMKwU1OSuptPP8wUF+uruHn+7uYeu+fhLJNNGIsWFJHb95xQquaG+is72RtrqKfJcrIjKr6vIYHYtq6Vh0+pa5VNrpHp483hLX1T/Ggb5gbcZfHBjgwecPn7RUSSxiLG3IhLemykyYC4LcyqYqGqrKNEZOzkrBBzUzuw74IrADuNvdf5zXgkrQsaEJfvTSMX648xiPv9J7fLHY9Ytr+a2rVvHLHS1csbppXldpFxEpBNGIBbNN6yt4w6rXtsolU2kOD04E4a1/jP19Y+zvG+dA3xiP7jhK7ylrOtaWx1jZXEV7SzBrvb25mvaWala3VNOoECenEerfrGZ2B/BO4Ji7Xzxj/43AXwNR4P+6+1dmOY0DI0AF0BViuZLh7uw4NMR/vHiUH+48xnNdgwAsra/gna9bwtVrWnjjmmZaNOBfREpcLBo53g16OtNPSdnfGyzGvb93lL29Y2w/OMjD24+c1BpXVxFjdWsNqzNBbvWMIKcZ76Ur1CcTmNmvEISsb00HNTOLAruAtxIEr6eAWwhC2+2nnOLDQI+7p82sDfiqu39wruvqyQRnL512nt7fz/3PHebh7Uc4MjSBGVy2ooHrL2rj+osWsa6tVv/aExGZJ4lkmq7+Mfb2Bot57+0dZW/PGK/2jHJocPykGavN1XHaM8FtzaJqOlpr6FhUw8qmKo2JKwJ5ezKBu28xs/ZTdl8J7Hb3PZni7gZudvfbCVrfzqQfUBPOPHJ3njkwwAPPHebB5w9zeHCCeCzCdRe28umN67huXatazUREQhKPRbigNVgj8i3rT/5uYirFgb4x9vSMsrcnCHGv9ozy2O5u/u3pE51L8WiE9pYqOhbV0NFaw5pFQYBb01qjBYCLRD4GFS0DDsz43AX80pkONrPfAN4GNABfn+W424DbAFauXDkvhRargwPj/Nu2Lu7Z1sX+vjHi0Qi/cmEr//3G9Vx/0SJqK9TELiKSTxVlUda21Z52+ZHhiSle6R5l97GR468XDw/z8PYjTPekmgUz76db3qZfa9tqqdOf8QtKwY/+dvfvAd/L4rjNwGYIuj7DrmuhmZhK8egLR/nu1gM8trsHd7h6TTOfvH4tN2xs0w+uiMgCUVtRxqUrGrh0RcNJ+yeTKfb2jJ0IcN3B9vFXeplMnngU37KGyuOLiwcLjddxQWs1ZepCLUj5CGoHgRUzPi/P7DtvZrYJ2NTR0TEfpysKBwfG+Zef7eOuJ/czMDbFsoZKPvmWtbznDcvPOPhVREQWnvJY9Hj4mimVdg72j/PysWF2znie8k92dR9fI64saqxprZkR3oIAt6S+QmOT8yzUyQQAmTFq98+YTBAjmExwPUFAewr4gLvvmK9rlvpkAnfnqb393Pn4qzyy4yjuzg0bFvOhq1Zx9ZpmIhH90ImIlLpEMs2enhF2Hp4OcEO8dGSYQ4MTx4+prYhx0eI6NiytY+PSOi5eVk/Hohq1vs2zvE0mMLO7gOuAFjPrAj7v7t80s08AjxDM9LxjvkJaqbeouTs/ePEYX//Rbp49MEB9ZRkfuWY1v3XVKpY3qvVMREROiMcirF9cx/rFdSftHxyfYtfRYXYeHmLnkWFePDzEd546cPypM8Gvq2Xj0vrj4W394lpNXghJ6C1q+VBqLWrptPPwjiP8zQ938+LhIZY3VvJfrl3Dey5fTmVcPzgiInJ+Umnn1Z5RdhwaZMehIbYfDLaD48GzUqMRo6O1ho3L6ti4tJ7XLa/n4qX1+jsoS7O1qCmoLWDuzo93dfPnD+1k55FhLmit5uPXdfCuS5eqWVpERELl7nT1j78mvB0bngSC8LaurZbXr2jgshUNvH5FAx2Laohq+M1rlExQm9H1+dGXX3453+WE6rmuAW5/cCdP7OllVXMVf3zDOt5xyRL9AIiISF4dG5rg+YODPHtg4Phr+tGD1fEor1sehLbpmauL6/Vc6JIJatOKuUWtd2SS2x/ayT3bumiujvPJ69dyy5UricfUgiYiIoUnnXZe7R3l2f0D/KIrCG4vHh5iKhXkj8V1wXNUO9sbuaK9ifWLa0vuaQslE9SKuUUtnXa+s/UAX3loJ6OTST5yzQV8/M1rtDitiIgsOBNTKV44PMSz+4Pgtm1fPwcHxgGoKY9x2coGrmhvorO9kUtXNFAVL/hlX89LyQS1acXWora3Z5RPf/cXbN3Xzy+tbuJLv3bxaVerFhERWagODoyzdW8fW/f289TePl46Oow7xCLGxqV1dLY3cUV7E1dd0ERDVTzf5c4rBbUFyt359pP7+fIDLxKLGP9r00beffkyLT4oIiJFb3B8iqf397N1bx9P7e3n2QMDJJJpzODipfVcvaaZN65p5or2JqrLF3aLm4LaAtQ/muDT3/0FP9h5jF/uaOEv3vs6ltRX5rssERGRvJhMpniua5DHd/fy+Cs9PLN/gEQqTSxiXLqiIRPcWrhsZcOCW9OtZIJasYxR235wkI/9yzaODU3y2ZvW8ztvbNfTBERERGYYT6TYtq+fn77Sw+Ov9PJ81wBph/JYhCtXN3HdukVct66VC1qqC74nqmSC2rSF3KL23a0H+J//vp2W6jh/+6E3vOahuyIiIvJaQxNTPLmnj5++0sNPdnWzp3sUgBVNlVx7YSvXXbiIN65pLshuUgW1BcDd+avv7+JrP9zNmzqa+dr7L6O5pjzfZYmIiCxIB/rG+PGubn7y0jEef6WXsUSKeDTCFasbefO6Rbx1QxurmqvzXSagoFbwplJpPvu957lnWxfv61zOl3/9Ej1ZQEREZJ5MJlNs3dvPT3Z18+OXjrHr6AgA69pquWFjG2/buJiNS+vy1kVaMkFtIY5RSyTTfPzbT/P9F47yqV9dyx9ev7bg+9JFREQWsgN9Yzz6wlEe3XGEp/b2kXZYWl/BDRsXc8OGNq5Y3ZTTBpOSCWrTFkqL2lQqzcf/39M8+sJR/uzmjfz2G9vzXZKIiEhJ6RtN8B8vHuXRHUf5z5e7mUymaagq48aNi7npkiW8cU1z6KFNQa0ApdLOH9z1NA8+f4QvbNrA775pdb5LEhERKWljiSRbdnXz0PYj/McLRxlNpPjbD17OTZcsCfW6swW1wpv6UCK+eP8LPPj8ET73josU0kRERApAVTzGjRcv4caLlzAxleInu7q5Zm1LXmtSUMuDbz72Knc+vpdbf3k1H7nmgnyXIyIiIqeoKIvyto2L810GmlqYY//5cjdfeuAF3raxjf9x00X5LkdEREQKWFEFNTPbZGabBwcH813KaR0ZnOBTdz/L2kU1/NVvXkpUTxsQERGRWRRVUHP3+9z9tvr6+nyX8hrJVJpP3vUM41PBwMSquHqdRUREZHZKCznyzcde5cm9fXz1fa+nY1FtvssRERGRBaCoWtQK1SvdI/zl93dxw4Y2fv2yZfkuR0RERBYIBbWQuTuf/bfnqSyL8qVfu1hPHRAREZGsKaiF7P7nDvPk3j4++/b1LKqryHc5IiIisoAoqIVoYirFVx7ayYYldby3c0W+yxEREZEFpqiCWqEtz3Hn43s5ODDO595xkZbiEBERkbNWVEGtkJbnGEsk2bxlD9de2MrVHfl9/ISIiIgsTEUV1ArJ3U8eoG80wR+8pSPfpYiIiMgCpaAWgkQyzeYte7hydROd7U35LkdEREQWKAW1EDyy4whHhib4/WvX5LsUERERWcAU1ELw7Z/vZ3ljJdde2JrvUkRERGQBU1CbZ3u6R3hiTy+3XLmSiGZ6ioiIyHlQUJtn393WRTRivLdzeb5LERERkQXunIKamd0/34UUA3fngecO86aOFhbV6ikEIiIicn7OGNTM7M8z2/ee5uuPhlbRArb94BD7+8Z4xyWL812KiIiIFIHZWtRusuAJ4p899Qt3PxxeSSczs4iZfdnM/sbMfidX1z0X9z9/iFjEuGGDgpqIiIicv9mC2sNAP/A6Mxua8Ro2s6FsTm5md5jZMTPbfsr+G83sJTPbbWafmeM0NwPLgSmgK5vr5ssPXjzGVRc001gdz3cpIiIiUgRmC2qfc/cG4AF3r5vxqnX3uizPfydw48wdZhYFvgG8HdgA3GJmG8zsEjO7/5TXImAd8Li7/xHw+2f7H5grhwbG2X1sREtyiIiIyLyJzfLdE8DlQFatZ6fj7lvMrP2U3VcCu919D4CZ3Q3c7O63A+889Rxm1gUkMh9T51pL2B57uQeAay7Ucz1FRERkfswW1OJm9gHgajP7jVO/dPfvneM1lwEHZnzuAn5pluO/B/yNmV0DbDnTQWZ2G3AbwMqVK8+xtHO35eVuWmvLWddWm/Nri4iISHGaLah9DPgg0ABsOuU7JwhQoXP3MeDWLI7bDGwG6Ozs9LDrOuXa/GxPH2/qaCaYfyEiIiJy/s4Y1Nz9MeAxM9vq7t+cx2seBFbM+Lw8s++8mdkmYFNHR8d8nC5rBwfG6RmZ5A2rGnN6XRERESlus62j9pbM234z+41TX+dxzaeAtWa22sziwPuBe8/jfMe5+33uflt9ff18nC5rz+wfAOCyFQpqIiIiMn9m6/q8Fvghr+32hCy7Ps3sLuA6oCUzKeDz7v5NM/sE8AgQBe5w9x1nW/gZrpeXFrVnDwxQHouwfonGp4mIiMj8MfecDufKic7OTt+6dWvOrvfuv3scA+75/atzdk0REREpDma2zd07T/fdGVvUzOyPZjupu3/1fAubb/loUXN3dh4e4r2dK+Y+WEREROQszLbgbW3m1Umw0OyyzOtjBOurFZx8jFE7ODDOaCLF2raanF1TRERESsNssz7/FMDMtgCXu/tw5vMXgAdyUt0C8PLREQAu1PppIiIiMs9ma1Gb1saJJwOQed8WTjnnx8w2mdnmwcHBnF3zpaPDAFy4SEFNRERE5lc2Qe1bwJNm9oVMa9rPCZ7hWXDy0fW56+gwi2rLqa8qy9k1RUREpDTMtjwHAO7+ZTN7CLgms+v33P2ZcMtaOPb1jnFBa3W+yxAREZEiNGdQA3D3p4GnQ65lQerqH+NX1rbmuwwREREpQtl0fS4YuR6jNplMcXRokuWNVTm5noiIiJSWogpquR6jdmhgAoDljZU5uZ6IiIiUlqIKarnW1T8GKKiJiIhIOM4pqJnZ5vkuZCHq6h8HYHmTuj5FRERk/p1ri9o/zGsV8yTXY9QOD05gBm215Tm5noiIiJSWswpqZhYxszp33xZWQecj12PUekcmaayKE4uqB1lERETm35wJw8y+bWZ1ZlYNbAdeMLM/Cb+0wtczMklLTTzfZYiIiEiRyqYpaIO7DwG/BjwErAZ+K9SqFoiekQTN1er2FBERkXBkE9TKzKyMIKjd6+5TgIdb1sLQOzJJi8aniYiISEiyCWr/AOwFqoEtZrYKGAqzqIWiZyShrk8REREJzZxBzd2/5u7L3P0mD+wD3pyD2s5aLmd9TkylGJlM0lKjFjUREREJRzaTCZrN7Gtm9rSZbTOzvwZyM63yLOVy1mfvaAKA5mq1qImIiEg4sun6vBvoBt4NvCfz/jthFrUQDI1PAVBfWZbnSkRERKRYxbI4Zom7f3HG5y+Z2W+GVdBCMTyRBKCmIpvfQhEREZGzl02L2qNm9v7MYrcRM3sf8EjYhRW64YmgRa22Qi1qIiIiEo4zNgeZ2TDBMhwGfAr4l8xXEWAE+HTo1RWwkcmgRa1WLWoiIiISkjOmDHevzWUhC81QpuuztlxBTURERMKRVcows0ZgLVAxvc/dt4RV1Lkys03Apo6OjtCvpa5PERERCVs2y3N8BNhCMC7tTzPbL4Rb1rnJ5fIcIxNJYhGjokwPZBcREZFwZJMy/hC4Atjn7m8GLgMGQq1qARieSFJTEcPM8l2KiIiIFKlsgtqEu08AmFm5u+8E1oVbVuEbnUxSo/FpIiIiEqJskkaXmTUA/w5838z6gX3hllX4JpIpKsqi+S5DREREiticQc3dfz3z9gtm9iOCx0c9HGpVC8DEVFrj00RERCRUZ9V35+4/CauQhWZiKkVFTC1qIiIiEh41CZ2jyWSacrWoiYiISIiUNM6RWtREREQkbApq50gtaiIiIhK2gl9fwsyuAT5IUOsGd786zyUBalETERGR8IXaJGRmd5jZMTPbfsr+G83sJTPbbWafme0c7v6f7v4x4H7gn8Ks92yoRU1ERETCFnaL2p3A14FvTe8wsyjwDeCtQBfwlJndC0SB20/59R9292OZ9x8Abg253qxNTKUoV4uaiIiIhCjUoObuW8ys/ZTdVwK73X0PgJndDdzs7rcD7zzdecxsJTDo7sMhlntWJqfSWvBWREREQpWPvrtlwIEZn7sy+2ZzK/CPsx1gZreZ2VYz29rd3X2eJc4unXYSqTTlMXV9ioiISHgWRNJw98+7++NzHLPZ3TvdvbO1tTXUehKpNIDGqImIiEio8pE0DgIrZnxentl33sxsk5ltHhwcnI/TnVEy7QCURRTUREREJDz5SBpPAWvNbLWZxYH3A/fOx4nd/T53v62+vn4+TndGqVQQ1KIRC/U6IiIiUtrCXp7jLuAJYJ2ZdZnZre6eBD4BPAK8CPyru++Yp+vlqEUt6PpUUBMREZEwhT3r85Yz7H8QeDCE690H3NfZ2fnR+T73TKm0WtREREQkfEU1yCpXLWopD4JaTEFNREREQlRUQS1XY9SSGqMmIiIiOVBUQS1X1PUpIiIiuVBUQS3XXZ8KaiIiIhKmogpqOVueIz09Rq2ofvtERESkwChpnIMTY9TyXIiIiIgUNUWNc3BijJp++0RERCQ8RZU0tDyHiIiIFJOiCmq5G6OmJxOIiIhI+IoqqOWK1lETERGRXFBQOwdankNERERyoaiCWs7GqKU1Rk1ERETCV1RBLWePkNKTCURERCQHiiqo5UpKY9REREQkBxTUzoHGqImIiEguKKidAz1CSkRERHKhqJJGriYTnBijFuplREREpMQVVdTI1WSCE9T1KSIiIuEpqqAmIiIiUkwU1M6BZyYTiIiIiIRJQe08mHo+RUREJEQKaiIiIiIFSkHtPKhBTURERMKkoCYiIiJSoIoqqOVqHTXNJRAREZFcKKqglut11EyzCURERCRERRXUcsVRk5qIiIiET0HtPKg9TURERMKkoCYiIiJSoBTUzoEmE4iIiEguKKidB80lEBERkTApqImIiIgUKAW1c6CuTxEREckFBbXzYJr3KSIiIiGK5buAuZjZSuBrQB+wy92/kueStIqaiIiI5ESoLWpmdoeZHTOz7afsv9HMXjKz3Wb2mTlOcwlwj7t/GLgstGLPgSYTiIiISJjCblG7E/g68K3pHWYWBb4BvBXoAp4ys3uBKHD7Kb/+w8DPgHvM7MPAP4dcr4iIiEjBCDWoufsWM2s/ZfeVwG533wNgZncDN7v77cA7Tz2HmX0a+HzmXPcA/3i6a5nZbcBtACtXrpy3/4bTcc0mEBERkRzIx2SCZcCBGZ+7MvvO5GHgk2b298DeMx3k7pvdvdPdO1tbW+elUBEREZF8KvjJBO6+HXhPNsea2SZgU0dHR7g1hXp2ERERkUA+WtQOAitmfF6e2Xfe3P0+d7+tvr5+Pk43J00mEBERkTDlI6g9Baw1s9VmFgfeD9ybhzpEREREClrYy3PcBTwBrDOzLjO71d2TwCeAR4AXgX919x3zdL1NZrZ5cHBwPk53Zur7FBERkRwIe9bnLWfY/yDwYAjXuw+4r7Oz86Pzfe7TMfV9ioiISIiK6hFSOWtRExEREckBK8Y1wcysG9gX8mVagJ6QryFnT/el8OieFCbdl8Kje1KYcnFfVrn7adcWK8qglgtmttXdO/Ndh5xM96Xw6J4UJt2XwqN7UpjyfV+KqutTREREpJgoqImIiIgUKAW1c7c53wXIaem+FB7dk8Kk+1J4dE8KU17vi8aoiYiIiBQotaiJiIiIFCgFtTmY2Y1m9pKZ7Tazz5zm+3Iz+07m+5+bWXvuqywtWdyTPzKzF8zsOTP7gZmtykedpWau+zLjuHebmZuZZreFLJt7Ymbvy/y87DCzb+e6xlKUxZ9hK83sR2b2TObPsZvyUWcpMbM7zOyYmW0/w/dmZl/L3LPnzOzyXNWmoDgr7EQAAAPXSURBVDYLM4sC3wDeDmwAbjGzDaccdivQ7+4dwF8Bf57bKktLlvfkGaDT3V8H3AP879xWWXqyvC+YWS3wh8DPc1th6cnmnpjZWuCzwJvcfSPwqZwXWmKy/Fn5HMHjFS8jeB723+a2ypJ0J3DjLN+/HVibed0G/F0OagIU1OZyJbDb3fe4ewK4G7j5lGNuBv4p8/4e4HrTs6XCNOc9cfcfuftY5uPPgOU5rrEUZfOzAvBFgn/MTOSyuBKVzT35KPANd+8HcPdjOa6xFGVzXxyoy7yvBw7lsL6S5O5bgL5ZDrkZ+JYHfgY0mNmSXNSmoDa7ZcCBGZ+7MvtOe0zmgfODQHNOqitN2dyTmW4FHgq1IoEs7kumq2CFuz+Qy8JKWDY/KxcCF5rZT83sZ2Y2W4uCzI9s7ssXgA+ZWRfBc7H/IDelySzO9u+eeRPqQ9lF8snMPgR0Atfmu5ZSZ2YR4KvA7+a5FDlZjKAr5zqCluctZnaJuw/ktSq5BbjT3f/SzN4I/LOZXezu6XwXJrmnFrXZHQRWzPi8PLPvtMeYWYygmbo3J9WVpmzuCWb2q8D/BN7l7pM5qq2UzXVfaoGLgR+b2V7gKuBeTSgIVTY/K13Ave4+5e6vArsIgpuEJ5v7civwrwDu/gRQQfC8ScmfrP7uCYOC2uyeAtaa2WozixMM6rz3lGPuBX4n8/49wA9di9OFac57YmaXAf9AENI05iY3Zr0v7j7o7i3u3u7u7QRjB9/l7lvzU25JyObPr38naE3DzFoIukL35LLIEpTNfdkPXA9gZhcRBLXunFYpp7oX+O3M7M+rgEF3P5yLC6vrcxbunjSzTwCPAFHgDnffYWZ/Bmx193uBbxI0S+8mGIj4/vxVXPyyvCd/AdQA383M69jv7u/KW9ElIMv7IjmU5T15BLjBzF4AUsCfuLt6BEKU5X35Y+D/mNl/I5hY8LtqAAiXmd1F8I+WlszYwM8DZQDu/vcEYwVvAnYDY8Dv5aw23XsRERGRwqSuTxEREZECpaAmIiIiUqAU1EREREQKlIKaiIiISIFSUBMREREpUApqIiJzMLMGM/uv+a5DREqPgpqIyNwaAAU1Eck5BTURkbl9BVhjZs+a2V/kuxgRKR1a8FZEZA5m1g7c7+4X57kUESkxalETERERKVAKaiIiIiIFSkFNRGRuw0BtvosQkdKjoCYiMgd37wV+ambbNZlARHJJkwlERERECpRa1EREREQKlIKaiIiISIFSUBMREREpUApqIiIiIgVKQU1ERESkQCmoiYiIiBQoBTURERGRAqWgJiIiIlKg/j/KjitIl7twiwAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 720x216 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "def fexact(t):\n",
        "    return x0/(x0* t*t + 1.0)\n",
        "xe = np.array([ fexact(t) for t in tr])\n",
        "fig = plt.figure(figsize=(10,3))\n",
        "plt.xlabel(\"t\"); plt.ylabel(\"abs. diff.\")\n",
        "plt.yscale(\"log\")\n",
        "plt.plot(tr,abs(xe-np.array(xs[:-1])))\n",
        "plt.show()\n",
        "plt.close()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q7jGOD0zssEL"
      },
      "source": [
        "Taylor展開するだけで簡単に調べられるようにEuler法は局所打ち切り誤差が$O(h^2)$、  \n",
        "大局打切り誤差が$\\sim O(h^2)\\times O(h^{-1})\\sim O(h)$の手法となっている。  \n",
        "つまり1桁精度を上げたいなら$h$を1/10、つまり計算量を10倍にしなくてはならない。\n",
        "\n",
        "上の$h$を変えてチェックしてみよう。  \n",
        "\n",
        "大局誤差が$h$の$p$乗に比例する手法を一般に$p$次公式と呼び  \n",
        "上のEuler法は1次公式となる。\n",
        "\n",
        "Euler法に少し工夫を加えた修正Euler法は、2次公式であり  \n",
        "$h$を1/10すると、大局誤差は1/100にできる。  \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6ru34y_hvFV6"
      },
      "source": [
        "### 4次のルンゲクッタ法 (RK4)\n",
        "\n",
        "次に4次公式の代表例である4次のルンゲクッタ法を紹介しよう。  \n",
        "古典的ルンゲ=クッタ法などとも呼ばれるこの方法は、  \n",
        "幾つかの$(x,t)$に対して$f(x,t)$の値を計算し、  \n",
        "打切り誤差が互いにキャンセルしあうように  \n",
        "うまく重み付けをして次ステップの$x$の値を求めていく方法であり  \n",
        "様々なところで用いられている。\n",
        "\n",
        "基本的には2変数関数のTaylor展開をやれば更新式を導出できるが\n",
        "煩雑なので導出については省略することにして結果だけ示すと\n",
        "\n",
        "$i$番目のステップでの$x,t$の値$x_i,t_i$が所与のとき、\n",
        "以下の4点での$f(x,t)$の(近似)値を用意して  \n",
        "$ f_1 = f(x_i,t_i),\n",
        "f_2 = f(x_i+\\frac{h}{2} f_1,t_i+\\frac{h}{2}),$\n",
        "$f_3 = f(x_i+\\frac{h}{2} f_2,t_i+\\frac{h}{2}),\n",
        "f_4 = f(x_i+h f_3,t_i+h)\n",
        "$\n",
        "\n",
        "$x_{i+1} = x_{i} + \\frac{h}{6} (f_1+2f_2+2f_3 + f_4 )$\n",
        "とすることで、4次公式が得られる。\n",
        "\n",
        "また、前述のような\"重み付け\"の選択には自由度があり、  \n",
        "$ f_1 = f(x_i,t_i),\n",
        "f_2 = f(x_i+\\frac{h}{3} f_1,t_i+\\frac{h}{3}),$\n",
        "$\n",
        "f_3 = f(x_i-\\frac{h}{3}f_1 + h f_2,t_i+\\frac{2h}{3}),\n",
        "f_4 = f(x_i+hf_1-hf_2+h f_3,t_i+h)\n",
        "$\n",
        "\n",
        "$x_{i+1} = x_{i} + \\frac{h}{8} (f_1+3f_2+3f_3 + f_4 )$\n",
        "\n",
        "という公式も知られている。\n",
        "\n",
        "\n",
        "前者を実装して、先程の例でEuler法と比べてみると...\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 225
        },
        "id": "bpzo_HH03zsy",
        "outputId": "2256e450-8361-4754-e1a8-332e05ff5372"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAADQCAYAAAC+9+0/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXicZbnH8e89M9mbrukCLSUFClgKFCj7YkWwwGEVRFZFqijIEQFRELUFUVBcjkAVUDgsHpayqAUru4Bs0iJby1oKbVMK3Zfsycxz/nhmMpPJJJkks2Sa3+e6pjPv82538jYz9zzba845RERERKT/C+Q7ABERERFJjxI3ERERkQKhxE1ERESkQChxExERESkQStxERERECoQSNxEREZECEcp3ALlQVVXlqqur8x2GiIiISLdeeeWVNc65kanWDYjErbq6mgULFuQ7DBEREZFumdnSztapqVRERESkQChxExERESkQStxERERECoQSNxEREZECMSAGJ4j0R5GIoyUSoTXsaA3HX7eEI7RGHK3hCC1hR2skQkv0dTgSfTiHc45wBMIRR8S5tudItDwS3a79+qTy6OvE8ohzOOdjdNHX0cXo6/j6tm06We+iZbEjOBffhrZ1nhkEzDAD8M+WWA6YXxldZwQs+jq6nmi53yf+2qI7GfFzGBAIGAEzggGiz/7R9tqMQCBpfawsuuyPQcK2yfvT6X6h6PZFgQChYPR1MEAwYJn7jyYiWxQlbjKghCOOhpYwDc3RR4t/NLWEaQ5HaGqJ0NQaoTkcpqklklAWprnVr4s/2pc1t4ajzz7Rap+M+QQsMTELR1z3AedQwHyCYWZtSQ+QkEBFl6P/xFILS0iE2paj+0VLktbHE6rY8WIJnsMnkbFkL5Y0Ohcrj6Z8yeUJ25KwHIlu7/rXr7pbZhAKGKFAwD8HjVAw/roo4JO71GXx5K8oaAQDAYqi5cFAgKJg9LhBi57DKA4FKApGH6EAJcEARSF/nOJoWXEwkLCdUZK4T3RdcXRdKKjGHJFsKcjEzcymAT8FFgH3OOeezmtAknHOOZpaI9Q2tVLb2EptUyubo8+1TS3UNoWpbWylvrmVhuYw9QnJWH1LmMbmMPUtrdQ3x177dU2tkV7HZAbFwQAloQDFoSAloQAlRf7DqqTILw8qCVFcHmj78Cxq+2CNf2D6Dzb/ujjUcb3fN2HbgP/gLAok1AjFanA6qdFJpzy5BmkgiNUgxhK6iHNEInSsgWxXK0mHmsuOtZy02z++vmONaOL+sQTe16ZGa1hjr6O1rr721bXVwrZGXLsa2XAkcZsIDS2xY0YSjh9pO09rJEI44UtFaxa+QASMtqSvLdkLWTSxi5fFk0JrS/xKQkFKivzfWUnC31nK191tGwooiZQtTs4TNzO7FTgaWOWcm5xQfgTwOyAI/Mk5d00Xh3FALVAK1GQxXOmDxpYwmxpa2NDQwob6FjbUN7OxoYWNseWGZjY2tLK5sYW6domZT9bS+UAxg/KiIGXFIcqKA5QXhSgtDlJeFGRUZSllxUHKioKUR5/LihNfh6LP8Tf64ugHQHEo0PbGHysrClpbrZMUpljtYABdxxgXTSBbwr62uDmaEMZqjtvKWn15S9jXMPvm+3gNc3PCPu2PFaGl1dEcjh8nsby+oaVdWayGu6k1QmNLmL7mlcGAJfwdd57gJSaBZUVBSouDlIaCbe8hZUV+fex9pKwoSGn0EV8OUBoKDpgvQpIf+ahxuw24AbgjVmBmQWA2cDg+EZtvZnPxSdzVSfufDfzLOfeMmY0GfgOcnoO4BzTnHLVNraytbWZNbRNraptZW9fUthx7TkzKGls6r90KGAwtL2ZIWRGVpSEGlYQYP7ycQaUhKktCVJSE2l4PKg0xqKSIQSUhKkuj66KP0qKAkimRPjDzNbxFwQDlxfmOpqPWcPvuCbHuDLHkLjHRi61vDieXR2hqCXd6nNqm1rbtG1siNLb2rYa+JBSgLCHxKy0KUlYUiD77pDCW6JUllJWG/BfL8pIQFcVByotDlBcHqSiJvy4vDlEcUi3iQJbzxM0596yZVScV7wMsds4tATCze4DjnHNX42vnOrMeKMlGnANJQ3OYTzY18snGRj7Z1MDKjY18urHRP29qZPXmJtbUNdPcyZvYkLIiRgwqZkRFMeOHl7PbuKK2pGxIWRFDy4sYWlbM0HK/PKS8iEHFIX0rFZFuhYK+ubMiD+/0kYjvshHrC9sY7XLRGF1uaA7T2BqhMaG/bGxdvMzXHMb23djQEl2OxPvbtoR7FFdR0BISuSAVJfGkrrw4SEVxiPKSYFtZWxJYEl2XsBz7ElxeHNSX4ALRX/q4jQWWJyzXAPt2trGZfRGYDgzF196l2uYc4ByA8ePHZyzQQtTYEqZmfQPL19WzLOGxfF09n2xqZEN9S4d9BpeG2GpIGaOHlLLDqEqqBhVTNajEJ2iDShhR4ZeHVxTr25+IbJECAfPNoMXBrJ4n1qc3lsTVN7dS1xSmLtqHt645TH1TK3XNYRqaW9uW65vD1Df77eqbwqza3Eh9U3y5rrk17abmgEFFSWIrR4hBpUV+OaGssjTVclF8n5KQRkVnWX9J3HrEOfcg8GA329wM3AwwderUAhtT1nPOOVZtbuL9T2t5f9VmFq+q5f1VtSxb65OzRGVFQcYPL2eb4WXsXT2cMUNKGTO4lK2GlPrXQ0opLy7I/xoiIgXHzNr6yw3L4HFjCWF9c5i6plYaWvxzbDmW9HUYABZ93tjQwor19W1ldc3p1QyWFwfbkrvKkhCDy4r8o7SorSVmcFko/jqhvLI0pAEl3egvn84rgG0SlsdFyySFlnCEdz/ZzMIVG3lzxUbeXrmJ91fVsrmxtW2bwaUhJo6u5KCJVYwfXh5N1Pxz1aBiVYmLiGzhEhPC4RV978AYjrguEr2WDonf5ug2mxpaWLG+gU2NfnBaS7jrupRBJaG2JC6e6MWTuyFlRQyrKGZ4eTHDKooYXlHMsPJiSouyWzPaX/SXxG0+MNHMJuATtlOA0/IbUv+xanMj8z9cz8sfruW15Rt4+5PNbf3NKktCfGbrwRw/ZSwTRw9ih5GD2GH0IEYOKlFyJiIiGRMMGINLfQ1Zbznn59Lc1NDaNsvApujzxoaWtuQutn5TQwvL1tW3ra/votavvDjIsGgyN6y8uC2hG15RHE/0yqNJX/RRVIC1e/mYDuRuYBpQZWY1wEzn3C1mdj7wKH4k6a3OuUW5jq2/qG1q5bn31/DMe6v495J1LFlTB/j/lLuNG8JZB1QzeewQdh07hG2Hl6uTv4iIFASz2MCKEGOGlPZ4/5ZwJDp7QTPr6lpYV9fM+vpm/1zXzLr6ZjbU+/Jl6+pZV9fcrjUq2fCKYkYOKqGq0j+PrCyhKvoce1QNKmF4eXG/+aw1V2hTivfC1KlT3YIFC/IdRpfW1jbx9zdX8vhbn/LvJetoDkeoLA2xT/Vw9t1uOPtMGMEuWw8uyG8HIiIi+dISjrA+IaFbX9fM2rrY1FZNrN7sH2tqm1m1uTHlVFbBgLUNyjts0mguOnzHrMZsZq8456amWtdfmkoHpJZwhMff+pQHXqnhmfdW0xpxbDeygrMOrObQnUex17bDlKiJiIj0QVEwwKjKUkZVdl/D55yjrjkcTeQSk7r466I817wpccuDTY0t3PPyMm57/iM+3tjI6MElzDhoAifsOZadxwzOd3giIiIDkpm1TWsyoaoi3+GkpMQth5paw/z5pWVc/9T7bKhvYb/thnPFcZM5dOdRmvdGREREuqXELUdeXbaei+97nSWr6zhwhxH84Iid2W3c0HyHJSIiIgVEiVuWOee44anF/PaJ99hqSBn/e9beTNtppKbqEBERkR5T4pZFLeEIF895nbmvf8xxU7bmquMnU9mH+W9ERERkYFPiliXhiOPCe1/j4TdWcsn0nThv2vaqZRMREZE+UeKWJb99/D0efmMllx25M9/87Pb5DkdERES2AJokLAte/GAts59ezMlTxylpExERkYxR4pZhreEIs+YuYtywMmYdu0u+wxEREZEtiBK3DHvojY9599PNXHbkZygvVku0iIiIZI4Stwy7/YWlbDeygiMnj8l3KCIiIrKFUeKWQYtX1fLa8g2cse+2GkEqIiIiGafELYOeeudTAKartk1ERESyQIlbBj3z3mp2Gl3J2KFl+Q5FREREtkBK3DLEOcebNRvZq3pYvkMRERGRLZQStwypWd/ApsZWdtl6cL5DERERkS1UQc5XYWYHA6fj45/knDsgzyGxeFUtADuPqcxzJCIiIrKlynmNm5ndamarzGxhUvkRZvaumS02s0u7OoZz7l/OuW8BDwO3ZzPedH28sQGArdW/TURERLIkHzVutwE3AHfECswsCMwGDgdqgPlmNhcIAlcn7X+2c25V9PVpwIxsB5yOlRsaCQaMUZWl+Q5FREREtlA5T9ycc8+aWXVS8T7AYufcEgAzuwc4zjl3NXB0quOY2Xhgo3NucyfrzwHOARg/fnxmgu/Cyo2NjKosIRjQ/G0iIiKSHf1lcMJYYHnCck20rCszgP/tbKVz7mbn3FTn3NSRI0dmIMSubWxoZlh5cdbPIyIiIgNXQQ5OAHDOzcx3DInqmsJUlATzHYaIiIhswfpLjdsKYJuE5XHRsoJR19xKRUnB5sEiIiJSAPpL4jYfmGhmE8ysGDgFmJvnmHqkrqmVimIlbiIiIpI9+ZgO5G7gRWAnM6sxsxnOuVbgfOBR4G1gjnNuUa5j6ws1lYqIiEi25WNU6amdlM8D5uU4nIypa1JTqYiIiGRXf2kqLXhN4QjFIf06RUREJHuUaWSQoTncREREJHuUuImIiIgUCCVuIiIiIgVCiZuIiIhIgVDiJiIiIlIglLhlist3ACIiIrKlU+KWQaZBpSIiIpJFStxERERECoQSNxEREZECocRNREREpEAocRMREREpEErcMsRpWKmIiIhkmRK3DNKgUhEREckmJW4iIiIiBUKJm4iIiEiBKMjEzcwmmdkcM/uDmZ2U73hEREREcqFXiZuZPdzbE5rZrWa2yswWJpUfYWbvmtliM7u0m8McCVzvnDsX+EpvYxEREREpJJ0mbmb2i+jzl1Ks/kYfznkbcETSuYLAbHxCNgk4NVqrtquZPZz0GAXcCZxiZtcCI/oQS8Y4DSoVERGRLAt1se6oaM3XZcB9iSuccyt7e0Ln3LNmVp1UvA+w2Dm3BMDM7gGOc85dDRzdyaG+HU34Hky10szOAc4BGD9+fG/D7RHdq1RERESyqavE7RFgPTDIzDYllBvgnHODMxjHWGB5wnINsG9nG0cTvx8CFcC1qbZxzt0M3AwwdepU1YeJiIj0cy0tLdTU1NDY2JjvUHKitLSUcePGUVRUlPY+XSVuP3LOXWJmf3POHdf38DLHOfcR0do0ERER2TLU1NRQWVlJdXU1toU3YznnWLt2LTU1NUyYMCHt/boanPBi9HlTF9tkygpgm4TlcdEyERERGSAaGxsZMWLEFp+0AZgZI0aM6HHtYlc1bsVmdhpwgJl9MXmlcy5l37Jemg9MNLMJ+ITtFOC0DB5fRERECsBASNpievOzdlXj9i3gYGAocEzSo7MBA90ys7vxtXk7mVmNmc1wzrUC5wOPAm8Dc5xzi3p7jnxQJzoREZHCFwwGmTJlStvjmmuu6XL72267jfPPPz9H0XVR4+acew54zswWOOduydQJnXOndlI+D5iXqfPkg+lupSIiIgWtrKyM1157LWvHb21tJRTqqsGza13N43Zo9OV6M/ti8qPXZxQREREpMNXV1axZswaABQsWMG3atA7brF69mhNPPJG9996bvffem+effx6AWbNmceaZZ3LggQdy5pln9imOrlK+zwJP4ZtGkzk6mT9NREREpK+ueGgRb32c2fGRk7YezMxjdulym4aGBqZMmdK2fNlll/HlL385reNfcMEFXHjhhRx00EEsW7aM6dOn8/bbbwPw1ltv8dxzz1FWVtb7H4Cum0pnRp+/1qcziIiIiBSIvjSVPvHEE7z11ltty5s2baK2thaAY489ts9JG3SRuJnZRV3t6Jz7TZ/PLiIiIpJCdzVjuRYKhYhEIgCdTuERiUR46aWXKC0t7bCuoqIiI3F0Naq0MvqYCpyLv7vBWPxo0z0zcnYRERGRAlBdXc0rr7wCwAMPPJBymy984Qtcf/31bcvZGOTQaeLmnLvCOXcFfjLcPZ1zFzvnLgb2AnJz888C4nSXeRERkYIX6+MWe1x66aUAzJw5kwsuuICpU6cSDAZT7nvdddexYMECdtttNyZNmsSNN96Y8fjSGY86GmhOWG6OlkmSATRnoIiIyBYpHA6nLD/44IN57733OpSfddZZnHXWWQBUVVVx7733dthm1qxZGYsvncTtDuBlM/tLdPl44LaMRSAiIiIiaek2cXPO/czM/oG/iwLA15xzr2Y3LBERERFJltbUvc65/wD/yXIsIiIiItKFrkaVioiIiEg/osQtQzSmVERERLJNiVsGaVCpiIiIZFOvEjczuznTgYiIiIjkWzAYZMqUKUyePJljjjmGDRs2APDRRx8xefLktu3++Mc/stdee7F+/fq2sl//+teYWdvN6LOhtzVuN2U0ChEREZF+IHav0oULFzJ8+HBmz57dYZs777yT66+/nkcffZRhw4YBsHz5ch577DHGj8/uPQp6lLiZWcDMBjvnXslWQCIiIiL9wf7778+KFSvalc2ZM4drrrmGxx57jKqqqrbyCy+8kF/+8pdYlmfj73Y6EDO7C39/0jAwHxhsZr9zzl2b1cji598OuBwY4pw7qbMyERER2YL841L45M3MHnPMrnDkNWltGg6HefLJJ5kxY0Zb2dKlSzn//PN59dVXGTNmTFv53/72N8aOHcvuu++e2XhTSKfGbZJzbhP+jgn/ACYAZ6ZzcDO71cxWmdnCpPIjzOxdM1tsZpd2dQzn3BLn3IzuyvJNtyoVEREpfLF7lY4ZM4ZPP/2Uww8/vG3dyJEjGT9+PHPmzGkrq6+v5+c//zlXXnllTuJLZwLeIjMrwiduNzjnWsws3TTlNuAG/G2zADCzIDAbOByoAeab2VwgCFydtP/ZzrlVaZ4r/3SzUhERkcxIs2Ys02J93Orr65k+fTqzZ8/mO9/5DgDl5eXMmzePgw8+mFGjRnH66afzwQcf8OGHH7bVttXU1LDnnnvy8ssvt6uVy5R0ErebgI+A14FnzWxbYFM6B3fOPWtm1UnF+wCLnXNLAMzsHuA459zVwNHphS0iIiKSPeXl5Vx33XUcf/zxnHfeeW3lo0aN4pFHHmHatGlUVVUxffp0Vq2K1zFVV1ezYMGCdv3fMqnbplLn3HXOubHOuaOctxT4XB/OORZYnrBcEy1LycxGmNmNwB5mdllnZSn2O8fMFpjZgtWrV/chXBERERmI9thjD3bbbTfuvvvuduUTJkxg7ty5nH322bz88ss5jSmdwQkjgJnAQfgbBDwHXAmszW5onnNuLX5wRJdlKfa7GbgZYOrUqeqBJiIiIt2qra1tt/zQQw+1vV64MN5lf/fdd+8w4hT8fG/ZlM7ghHuA1cCJwEnR1/f24ZwrgG0SlsdFy0RERESkC+kkbls5537qnPsw+rgKGN2Hc84HJprZBDMrBk4B5vbheCIiIiIDQjqJ22Nmdkp08t2AmZ0MPJrOwc3sbuBFYCczqzGzGc65VuD86DHeBuY45xb19gfoTzSmVERERLKp0z5uZrYZ36fNgO8Cf46uCgC1wPe6O7hz7tROyucB83oarIiIiGzZnHNZv/tAf+F6MQlsp4mbc66yT9GIiIiI9EBpaSlr165lxIgRW3zy5pxj7dq1lJaW9mi/dOZxw8yGAROBtqM7557t0ZlEREREujBu3DhqamoYKNN4lZaWMm7cuB7tk850IF8HLsCP/nwN2A/fb+3QXsQoIiIiklJRURETJkzIdxj9WjqDEy4A9gaWOuc+B+wBbMhqVAWmN23UIiIiIj2VTuLW6JxrBDCzEufcO8BO2Q2rMG3hzfEiIiKSZ+n0casxs6HAX4HHzWw9sDS7YYmIiIhIsm4TN+fcCdGXs8zsn8AQ4JGsRiUiIiIiHaQ1qjTGOfdMtgIRERERka6l08dNRERERPoBJW4iIiIiBUKJWwZoNhARERHJBSVuGWS6zbyIiIhkkRI3ERERkQKhxE1ERESkQChxExERESkQStxERERECkS/T9zMbDszu8XM7k8o+4yZ3Whm95vZufmMD0CDSkVERCQXspq4mdmtZrbKzBYmlR9hZu+a2WIzu7SrYzjnljjnZiSVve2c+xZwMnBg5iPvHd1kXkRERLIp2zVutwFHJBaYWRCYDRwJTAJONbNJZrarmT2c9BjV2YHN7Fjg78C87IUvIiIi0n/06F6lPeWce9bMqpOK9wEWO+eWAJjZPcBxzrmrgaN7cOy5wFwz+ztwV2YiFhEREem/8tHHbSywPGG5JlqWkpmNMLMbgT3M7LJo2TQzu87MbqKTGjczO8fMFpjZgtWrV2cwfBEREZH8yGqNWyY459YC30oqexp4upv9bgZuBpg6darGD4iIiEjBy0eN2wpgm4TlcdGyguV0s1IRERHJgXwkbvOBiWY2wcyKgVOAuXmII+M0qFRERESyKdvTgdwNvAjsZGY1ZjbDOdcKnA88CrwNzHHOLcpmHCIiIiJbgmyPKj21k/J5aBoPERERkR7p93dOEBERERFPiZuIiIhIgVDilgEaUyoiIiK5oMQtg3SvUhEREckmJW4iIiIiBUKJm4iIiEiBUOImIiIiUiCUuImIiIgUCCVuGaBblYqIiEguKHHLINOwUhEREckiJW4iIiIiBUKJm4iIiEiBUOImIiIiUiCUuImkq34dNKyHSDhe1toMK9/wjwW3wsYav37V21C7On+xiogMVHVrOpY1rIfWJnj9Xpg1BH6zC7w1t/Nj1K+D5fPbHysS8cfJs1C+AxDJq79/D+b/sWP5TkfBu/Ng8omw8IHU+47bG2rm9+x84/aBrz/u3wBwEAjG1zkHLtK+bKCIROCdh2DsXvD01fDqnxNWGhz9W3j4u1BUARe/A6WD8xaqSFoiEYi0QqjYLy95GrbeE95/DB6Y4cuKyqGlPr7P6ffDxMNzHmrBWv8R/G733u27qQbmnBlfnnAIHHABjNsLflHd9b6Dx8JFb/XuvBmgxC0DnG4z70UiULcKBo2G1kYIlfbuBq6L/gpv/RUW/QVOfwC2ngKPXApv3ufXn/EgjN8fmmuhdAiEStrvv/o9eO43sP2h8OA3evezvDvPP3eWtEHPkzaAmpf9t73ubLU7TDoealdB1Q6w25d9TV7Z0O73jUSiv/+SjklguAX+fRM8djmc/RgM3w6CRVBSmZ2EMdwKy1/yP0ftp/DCDTBiO/iv38K6JXDXl9I4iPNJG0BLHVyzjX897TI45BKf7FoQAmpAkH7gldvhoe+kt21i0gbwfyfFX39lLozd0/9t5sqmj+E3n4FAEZz/Mny6CO49w6/b/3yY/rPcxRKJQP1a/54XLOq4/saD4JM3M3e+D5/1j3TsckLmztsL5vr5JGRmth1wOTDEOXdStGwa8FNgEXCPc+7pro4xdepUt2DBgqzF2NQaZqcfPcIl03fi25/bIWvnyYtwCzx0ga8ufv9RX7b3N+C/fuVfN9fBwgdh7vmp9z/2ev/G091/9OZ6WP0OLLglqbYlh3Y8ElYtgg3L2pd//iewx5n+5ygqgzfv99+YS4fCUdfCZ45pn6S2NkMgBO89Am/9DTavhJE7w0fP+eP31defhHFTO5b/+jPQsM4nbclmbfTNAonfMJPtejKcmKL2MWbtB1BcAeVVPrHeZh8YOj6+fsNyX8MwfALc9zVY9GD6P1OysmFwwk3+Wsz7Xnr7nH6/j2/8/r37wiDSV6/dDX/9Vs/2KSqHI67xX6g6e3848AKYOgOeugrenANnzYPqA/seb7LuvlQOGQ8XvN7zL0mv3QWDt4Y7joM9vwI7/RfsOB2W/xu23sN/yVz2kv/5PvqXf79c/U58/52P9u81lVvB5k/gpdntj3/eS1C1o0+EN62EEdv7L6If/ssnYwf8t0/+isra7xcJ+/M8dVX8yzpAxSi45H3/urXJJ7Qlg/37S1Fpz372XjKzV5xzKd7os5y4mdmtwNHAKufc5ITyI4DfAUHgT865a9I41v0JidtngUuBT4GrnHOLu9pXiVsv1a+DX07I7DFDpfC9931TVzo1T9n2jad881yurHrbv1E3rIebP9tx/fE3pvfG/+M18MYcWPm6T+KevRbWvNf3+CZ8Fr6aot/Hilfgj4emd4zyKqhP0cekK197BLbdP71tV74BNx3c+foJh8BXH+rZ+UXS4RxcEa31Lq6E1gY4+GL43A/9B/xVo/y6ETvAEb+A8fvB386DERPh8z/2NdDBLhq6PnoOnrzSJzTdOel/YfRkePIK2O1kmHRc6u3CrfDTEfDZSyHc7FtFqg+B3b/sY372Wv/oiRNv8bWB//y5/2I7dLx/L3rlNhh/gP8iXzwIJhzsW06yZc+vwrHXZe/4eZTPxO0QoBa4I5a4mVkQeA84HKgB5gOn4pO4q5MOcbZzblV0v8TELeCci5jZaOA3zrnTu4pjQCRuzsHbc33frFTVyj2xfD7ccljv9h02AdZ/2LfzA/xknf/G1LQZmmph0Kh4U96/b4bFT/g3y8Tkp3QING70yc+O033tz6BR7Y/buMknTl29eeZD02b/RpdcU9S4ydferf8I/pBmYgOw77n+Zy+p7FhjdeEiGDLOn7O1CSzQPkE/5BI49EdQswD+9Ple/0htLl3mf+eRsP9ZAkF49Iew+6mw1W69P25LA9x2tL/uHzzZcf05T/tv87nUuNH//Uzs5d+P9G9XDPNN8935yfrMNN0vfQH+98j0tp10HJxwc/saoeUvwy097DN38XtQOdonXCWVsEP0//K878PLN/XsWL015QyfjG6zD/zlm77lItmB34XDr8hNPHmQt8QtevJq4OGExG1/YJZzbnp0+TIA51xy0pZ8nLbELaGsGLgruTzZFp24hVvh9/vC2oRKx+qDYenzHZaYTrQAABD+SURBVN9gfrSqY3+wVJJrwr7/IZQP969bm+Gqke3Xn3qvb6Kbcppfjv2fuuNYX4Pz1E/T+1mqdvLnOf2+3PbrKCQN67vuODtzQ8fkb/ET8OcT/esfr+k8sY81AXdl33Ph9bvh5Nt9s0dXZm3sen2mPfYjeOH6juUn3gK7dvkW0b3Nn/j+NrcdDUf+wg9aiX2RaK6Hpk2+KeWX2/laGIh/ADbXwa3TffPQdtN8LUwmm3JX/Af++Dnfp/PMLNZu9MXaD+D6Pf3r5P8XDevhps/CcbN9Dc2G5b5JLBiC9Ut9V4Mh4/wXm9pVvnYnF81VdWvgjuN9M+WuJ/m/jwe/3v1+31sMg0Z2v11PxH5/F7wOZcPj/TxTOfq3fvsXb+jZOb50u+9bO7yLVpZ//AD+fWPPjgv+mkfC8b+ZSASWveC/OD7/O//3cMaDnfezdc4/Ni6HO4+Hb/wzvf6+Bay/JW4nAUc4574eXT4T2Nc5l7KTlJmNAH6Gr6H7k3PuajP7IjAdGAr8IVUfNzM7BzgHYPz48XstXbo0wz9ZXF4St3ArvHEP/O3b6e8zdi84bQ5UVPnlmlfg0zdhr7Pi2zz+E/+HFPNfv4G9u/kw74mHLvDV6TGX1cDT18ChP85Z34GC9/CFfuqRAy+Az8/0b2hmmRlc0FXz9XcXwtAuPjAg/q38wrdgyNi+x9Mb4Rb4aVX7ssOvhN1P8zWX2+zds+M118PPt+pY/oOl8MkbcPsxne/7zWfh5s+BS5hCZtuD4Gt/T+/czsGjl/vf5Yu/h+N/7/vs7HICLHvRJxPLX4pv//Un/f+P2KCdrabAxC/A+H3TO19f1K/zg4L2+7aveQyEYN0H/vxXJHzIfvUhuPtUP7hojzM679N68bvw651Sr9v5aPjyn7PblzHxb2Gfb7avbUpMPiNhWPykrxk652kYtm32YoqpWws3TIVjfuf72F7RTRJz+Sd+AE9shGtis+7nf+Kbe9MVifgmV/C19X+/CF690zfbfus5//voby0aBaigE7dMyHaNW2NLmJ1/nOXEraXBv8FNneGr4P96Hrz2fz0/zvDt4Duv+jfWa6Idy3+y3n+w3HoErEj4PWWjxiSxxujQH8MhaXY8l9yoXQ2/Svg//L33OzY3F4KupgnoaR+4bPTFPPux9JKpB74eH03dF7mo/fzrt+G1HA8s2uFwWPw4/Hht5pOFzq57V7XW+VK3Fv5zu+/vluxbz8OYyR3Lm+v9AKI9zsh+fNJjXSVu+Rg/vwJI/No+LlpW8LI6kO1nY3xfpVunwyOXpZe0zXjCJ2UVCR+865bADfvEkzaAK4f54ycmbZd8kLnYE5UN8x8iMzcoaeuPBo301yb2KMSkDWBYtY9/+xQDKj581n8or3o7XlazwE970JI0IveJWT0/98wNqct/lDAh861f6Hz/ujXw55N8LJlI2sA3p2bKytf9wJhEznWftI3pZX/GrRIS8JE7t1+3+HH//I9Luj5G3RqYf0u8G0df9LekDaBiBBx8ka/hBV/zOnODf69NlbQBFJcraStQ+ajPnA9MNLMJ+ITtFOC0PMTR/zSs96P3Whp9glNc4Sf6S2xCrHnZPxIlzq+z8AG4/2zfRBJrErpwUft+aWve7XjuSGv75bJhff95uqLpGvqvLeXamMHJd8LVnTTZ/n4/P9Ju6z3g7i/7shev94MzWpt9X77nfhvfPlZr9cgPO05HELP1Hv68n7sc/vmzjvsmeu0u2PGIeP9R8P3hrt3ev76+B6Odv7e4fU1pspr5fhRgsg3Lff+xnlzze8+EDUuhamJ88MeHz3S/38m3+y+d7z0SLzvoQjhsVvvt3n8C/u/E+PKp9/ipJGI2roDfTmq/z6DR3ce87AXfTWSHw/x7bMWIzrdf/1H89dZ7wsfRxPf7GRh4lU1b7Z77vqWSc9keVXo3MA2owk/dMdM5d4uZHQX8D34k6a3OuazO6perptLvH7ET503rQ1NpZ1XzEw7pfGJAC8LMdX07frLDroCDvpvetiKFoCfNnbM2wi+3bz+lyWUroGSQf73mfd+/KOaUu2Hnozoep6XB12TvfiqcEO3Q/cvt/CCH5PN1F2dsm+T1sfJIGK6MJoD7ngtHXgP3ndV+KobkD/RPF8EfDvBziO13burzpowlIYbY6MlY2Rf/1HkH/p4kFLHjHXu9n/eruzi6O37itlU7+S+v6W6vREjyoKum0qzWuDnnTu2kfB4wL9W6QrbNuhdgVvRb8un3+1mwdzgMvnCVn3tq9y/7qvpnf+UHAHzjKd80Ne/7flLFznQ1m/P3l6Qf4PafTz1tQrJ9v5n+MUUKwUXv+Mk5nfOJU1dNla/f03EeuljSBr6m6St/i4+qTZW0gZ/s87tvwqAx8bLvvgk/37r9djWv+GM+/z+pj3N6wt07LlkCGz7yAzBG7xIvDwThnGfg/cfjHc1PuNnfXWL2Pn65fl372r2Hol/OHrnUj5ItG+5j2O9cX9sf88Z9Phm75ANfYxUqi4+c/fWO/meKad7sa6U+XejPN2i0b0XYLsWchV25rAaWvgg7dDEdzfkL/L2B7zzeLy/6K+xyfMftXk6aVDrW4hAb1BOz5Bnf8pDYvH50J9dEJI809CODjnnjv+MLsVuXLH7CP8Anbp+8Af+8yi/P3hsmnwQL7+/dCaf9sGdDopNHBO79jdT36QxpdKdsYQYnjgrtplb8L0lfXA5Icfui6kNgzK7dj8ZLvLME+PnsRu/qR3PH/OlQP5Hof25PfYzE0bkVIzpv4tt6in/EhIr9bPIxj/4wXvMH7btczPkq7PVVP3VPw/r2tzaK1aA9+ytfk5eobnX7yVtjTb8TDomXpTuxcqKSStixi+QafLI7IuFa3vdV2CVF7Vhnd95Y9iJse0B8+Y5j/fP5r8TLOpvUViSPlLhlSBVpVKdvWO5vnZEoudkkXdsfCtN+0LN9Rkz0z7GmFGifuF26XDfvFkn2hRTzEAYCfuqDnjKDc5/r2MzXWdIGfsqF3jLzNWkN63yfvUnH+5rA5NuiLXvB3xMXfOLWFted8dfvPOxn3Y/VtsX869fx15Uppk3JpuS+eQ98HaZfHZ9HbeXrne+7+En/u1j2kp/KImbzyvjrxBpKkX5CiVuGLChNo4/I/6QY3ZN8X8yubDUFVr7mX8dms+6J2K2dxu+Xen3yfdxEtlRTTu/ddDr5UNbH5KEhoQ9sbBBGKv+5wz9v/sQ/v3ADPHZ5fP3G5f7RlXwPbHnzPv+I9Uu7KaHmLxBqPwjrX7/yj2S3H53dGEX6SIlbvtWtbr88Zld/K4/kGeynng1H/crPv+YiUN7FiKjObLt/x1m9L10GgSL/Dbw/DnMXyYbjZvupEEqH+hHUzbW+5uj1u9tv98OPU+/fV5d84KfouOtLHdd97R++b9iHz/g7j/R1Fv7h2/lpgNI1rNo/L3uxZ+eZ0uWdB/Pv8k/h4Qs6n/A32Qk5ur2TSA8pccu35NtSjdwZdvlix8Rt+Pa+A3Jfq+6TPwRKo002xeV9O65IITFr378J/LQYyRI76WdSRVXnfbhKBsOI7f0jE3r6M7REm0Jbm3q23+Ctu98mn4Ih36c43cStNAsTL4tkgBK3fGuubb/sXOqbE2fidkYi0rlAP6lxznRzY0/7yL1+l3/0+Dz96D3qp6MgnCLxzHdTrkgG5OPOCVue5ObOvtjpyPjrxFmtt9knc+cQkY6Sp6yY0MMpLPqrXU/O/DEPuqhjWU/vAZsppSlG1qdK2qD9KNTuJI7IFelHlLhlQlNd349x0Tu+38uu0WlEfrQKjrne3ybn+x/GBxaISHaM3w9KEprHDpuV/XOelmr+xgzXCu3/bT//W7KvPeJvPv6997s/xqikOxUcNhN+8FH7sjGd3Bs2277/IVz4VvfbgG8Ov/zT9I6bqaZqkQxTU2l/MThpGH2oxD8HiiGkIekiOVExApqiIxJz0ayWiwFBZqnnfhu9ix9Jns5o8pIU0wRl+7Z46QoEur+nbuLvuUjzVEphU42biEhKeeoPlat+WBnvS9ef+4/159hEekaJm4hITBbv3dwvztdb/TopExlYlLiJiKSSt2QlV+ftyXnS2bYfJ3dKPGULosRNREREpEAocRMRSUl93Hq0bb+u1erPsYn0jLlC6WPRB2a2Gliag1NVAWtycB5Jn65J/6Tr0v/omvQ/uib9Uy6uy7bOuZT3uxsQiVuumNkC59zUfMchcbom/ZOuS/+ja9L/6Jr0T/m+LmoqFRERESkQStxERERECoQSt8y6Od8BSAe6Jv2Trkv/o2vS/+ia9E95vS7q4yYiIiJSIFTjJiIiIlIglLj1kJkdYWbvmtliM7s0xfoSM7s3uv7fZlad+ygHnjSuy0Vm9paZvWFmT5rZtvmIcyDp7pokbHeimTkz0+i5HEjnupjZydG/l0VmdleuYxxo0nj/Gm9m/zSzV6PvYUflI86BxMxuNbNVZrawk/VmZtdFr9kbZrZnrmJT4tYDZhYEZgNHApOAU81sUtJmM4D1zrkdgN8Cv8htlANPmtflVWCqc2434H7gl7mNcmBJ85pgZpXABcC/cxvhwJTOdTGzicBlwIHOuV2A7+Y80AEkzb+VHwFznHN7AKcAv89tlAPSbcARXaw/EpgYfZwD/CEHMQFK3HpqH2Cxc26Jc64ZuAc4Lmmb44Dbo6/vBz5v1q+nFN8SdHtdnHP/dM7VRxdfAsblOMaBJp2/FYCf4r/cNOYyuAEsnevyDWC2c249gHNuVY5jHGjSuSYOGBx9PQT4OIfxDUjOuWeBdV1schxwh/NeAoaa2Va5iE2JW8+MBZYnLNdEy1Ju45xrBTYCI3IS3cCVznVJNAP4R1Yjkm6vSbRpYRvn3N9zGdgAl87fyo7Ajmb2vJm9ZGZd1TpI36VzTWYBZ5hZDTAP+O/chCZd6OnnTsaEcnESkf7CzM4ApgKfzXcsA5mZBYDfAGflORTpKIRv/pmGr5l+1sx2dc5tyGtUA9upwG3OuV+b2f7AnWY22TkXyXdgknuqceuZFcA2CcvjomUptzGzEL5ae21Oohu40rkumNlhwOXAsc65phzFNlB1d00qgcnA02b2EbAfMFcDFLIunb+VGmCuc67FOfch8B4+kZPsSOeazADmADjnXgRK8ffLlPxJ63MnG5S49cx8YKKZTTCzYnwn0blJ28wFvhp9fRLwlNNkednW7XUxsz2Am/BJm/rsZF+X18Q5t9E5V+Wcq3bOVeP7HR7rnFuQn3AHjHTew/6Kr23DzKrwTadLchnkAJPONVkGfB7AzD6DT9xW5zRKSTYX+Ep0dOl+wEbn3MpcnFhNpT3gnGs1s/OBR4EgcKtzbpGZXQkscM7NBW7BV2MvxndsPCV/EQ8MaV6Xa4FBwH3RsSLLnHPH5i3oLVya10RyLM3r8ijwBTN7CwgDlzjn1GqQJWlek4uBP5rZhfiBCmepQiC7zOxu/BeYqmjfwplAEYBz7kZ8X8OjgMVAPfC1nMWmay8iIiJSGNRUKiIiIlIglLiJiIiIFAglbiIiIiIFQombiIiISIFQ4iYiIiJSIJS4iYj0gpkNNbPz8h2HiAwsStxERHpnKKDETURySombiEjvXANsb2avmdm1+Q5GRAYGTcArItILZlYNPOycm5znUERkAFGNm4iIiEiBUOImIiIiUiCUuImI9M5moDLfQYjIwKLETUSkF5xza4HnzWyhBieISK5ocIKIiIhIgVCNm4iIiEiBUOImIiIiUiCUuImIiIgUCCVuIiIiIgVCiZuIiIhIgVDiJiIiIlIglLiJiIiIFAglbiIiIiIF4v8BgOMboFFSbPsAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 720x216 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "#求解(RK4)\n",
        "xsRK4 = [x0]\n",
        "x = x0\n",
        "for t in tr:\n",
        "    f1 = f(x,t)\n",
        "    f2 = f(x+0.5*h*f1, t+0.5*h)\n",
        "    f3 = f(x+0.5*h*f2, t+0.5*h)\n",
        "    f4 = f(x+h*f3, t+h)\n",
        "    x += h*(f1 + 2*f2 + 2*f3 + f4)/6.0\n",
        "    xsRK4 += [x]\n",
        "#描画\n",
        "fig = plt.figure(figsize=(10,3))\n",
        "plt.xlabel(\"t\");plt.ylabel(\"abs. diff.\")\n",
        "plt.yscale(\"log\")\n",
        "plt.plot(tr,abs(xe-np.array(xs[:-1])),label=\"Euler\")\n",
        "plt.plot(tr,abs(xe-np.array(xsRK4[:-1])),label=\"RK4\")\n",
        "plt.legend()\n",
        "plt.show()\n",
        "plt.close()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WqKW52Ncbg7c"
      },
      "source": [
        "となり、高精度で解が得られていることが分かる。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uIqPPSsJb-KD"
      },
      "source": [
        "## 積分を用いた解法"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SUuqmECsuB4l"
      },
      "source": [
        "常微分方程式 $\\frac{dx}{dt} = f(x,t)$は形式的には  \n",
        "$x(t') = x(t_0) + \\int^{t'}_{t_0} f(x(t),t,) dt$ と書くことが出来る。\n",
        "\n",
        "右辺の積分を数値的に計算することで、  \n",
        "元の微分方程式の解を求めることを考えよう。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GSy9sdWUdyXI"
      },
      "source": [
        "### Adams-Bashforth法"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0mWaZoxBd1m5"
      },
      "source": [
        "\n",
        "\n",
        "以下では特に、Adams-Bashforth法(AB法)と呼ばれる手法を考える。\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "svbcn9pIM-n1"
      },
      "source": [
        "まずは、2次のAB法を考えよう。\n",
        "\n",
        "$t_i$と$t_{i+1}$の間の$f(x(t),t)$を、現在居る$i$番目のステップと  \n",
        "1個前の$i-1$番目のステップを用いて、1次のLagrange補間で近似してみる。\n",
        "\n",
        "$\n",
        "\\begin{align} \n",
        "x_{i+1} &= x_i + \\int^{t_{i+1}}_{t_i} P_1(t) dt  \\\\\n",
        "P_1(t) &= f_{i-1}L_{i-1}(t) + f_i L_i(t) \\\\\n",
        "L_{i-1} &= \\frac{t-t_i}{t_{i-1}-t_i}, L_{i} = \\frac{t-t_{i-1}}{t_{i}-t_{i-1}}\n",
        "\\end{align}$\n",
        "\n",
        "各ステップの刻み幅$t_{i}-t_{i-1}$ for $\\forall i $を一定$h$としよう。  \n",
        "このもとで上の積分を評価すると、  \n",
        "$x_{i+1} = x_i + \\frac{h}{2} (-f_{i-1}+3f_i)$\n",
        "という更新式を得る。\n",
        "\n",
        "この手法の局所打切り誤差は(またTaylor展開して)$O(h^3)$であり、  \n",
        "２次の公式となる。一般的な分類に倣えば、この手法は過去の2ステップの情報を利用する、  \n",
        "陽(explicit)解法かつ２段法(2-step method)となっている。\n",
        "\n",
        "\n",
        "同様にして、一般に$m$ステップのAB法($n-1$次のLegendre補間を使用)は$m$次公式を与える。  \n",
        "たとえば4次公式は$x_{i+1}=x_i + \\frac{h}{24}(-9f_{i-3} + 37f_{i-2} -59f_{i-1} +55f_{i})$\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "40VvY7-xRbXC"
      },
      "source": [
        "**Adams-Bashforth法とRunge-Kutta法の違い**\n",
        "\n",
        "AB法では、過去に評価した$f$の値を用いるため、計算量の増加を防ぐことができる。  \n",
        "つまり、各$f$の評価は１回だけで良い。  \n",
        "一方で、RK4では各ステップで$f$を４回計算する必要がある。\n",
        "\n",
        "これは$f$の評価が複雑な場合  \n",
        "(例えば解きたい微分方程式が行列(連立)になっていて、$f$に行列演算を含む場合など)  \n",
        "Runge-Kutta法よりも、高速に計算できるかもしれない。\n",
        "\n",
        "一方で、$m$次のAB法では過去の情報を利用するため、  \n",
        "はじめの$m-1$回については、RK法などで予め求めておく必要がある。  \n",
        "また、Runge-Kutta法と比較すると、誤差が大きかったり、  \n",
        "数値的に不安定になりやすいことも知られている。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 225
        },
        "id": "v8lSgfj-jvu5",
        "outputId": "45d3597a-12ff-43bf-eaa5-b24b1f4b6d0c"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAADQCAYAAAC+9+0/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5xcdb3/8ddnZrZveiGQEDb0EiBA6FUUErw0BRVEREC5lvxAUJqNYKGoWEC8gMKlKEEIXgiIdBBBWqKUhBpasrT0stk2O/P9/fGdvrOzs7tTdrLv5+MxzDnf0z6bw85+5tuOOecQERERkcEvUO4ARERERCQ/StxEREREKoQSNxEREZEKocRNREREpEIocRMRERGpEErcRERERCpEqNwBlMLYsWNdU1NTucMQERER6dWCBQtWOOfGZds2JBK3pqYm5s+fX+4wRERERHplZu/1tE1NpSIiIiIVQombiIiISIVQ4iYiIiJSIZS4iYiIiFSIITE4QWQwikYd4WiUroijK5JcDkeidEUdXZEo4YijKxolHFuORGMv53DOEYlCJOqIOpd4j8bKo7H90rdnlMeWU8ujzuGcj9HFlmOrseXk9sQ+PWx3sbL4GZxL7kNim2cGATPMAPy7pZYD5jfGthkBiy3HthMr98ckly12kJG8hgGBgBEwIxgg9u5fiWUzAoGM7fGy2Lo/Byn7Zh5Pj8eFYvtXBQKEgrHlYIBgwAr3P5qIbFSUuMmQEok62sIR2jpjr7B/dYQjdEaidISjdHRF6YxE6AhHU8oidHb5bclXellnVyT27hOt9GTMJ2CpiVkk6noPuIQC5hMMM0skPUBKAhVbj/0nnlpYSiKUWI8dFyvJ2J5MqOLniyd4Dp9ExpO9eNLoXLw8lvJllqfsS8p6NLa/G1z/1L0yg1DACAUC/j1ohILJ5aqAT+6ylyWTv6qgEQwEqIqVBwMBqoKx8wYtdg2jOhSgKhh7hQLUBANUhfx5qmNl1cFAyn5GTeoxsW3VsW2hoBpzRIqlIhM3MzsE+AmwCLjNOfd4WQOSgnPO0dEVpaWji5b2Llo6ulgfe2/pCNPSEaGlvYvWzi7aOiO0piRjreEI7Z0RWsNdtHbGl/22jq5ov2Myg+pggJpQgOpQkJpQgJoq/8eqpsqvN9aEqK4PJP54ViX+sCb/YPo/bH65OtR9uz82Zd+A/8NZFUipEYrX4PRQo5NPeWYN0lAQr0GMJ3RR54hG6V4DmVYrSbeay+61nKQdn9zevUY09fh4Au9rU2M1rPHlWK2rr311iVrYrqhLq5GNRFP3idIWjp8zmnL+aOI6XdEokZQvFV1F+AIRMBJJXyLZC1kssUuWJZNCSyR+NaEgNVX+96wm5fcs63Jv+4YCSiJlo1PyxM3MbgCOBJY556amlM8EfgsEgT865y7LcRoHtAC1QHMRw5UBaA9HWNcWZk1bmDWtYda0drK2Lcza+HpbJ2vbuljfHmZDWmLmk7V8/qCYQX1VkLrqEHXVAeqrQtRWB6mvCjJ+WC111UHqqoLUx97rqlOXQ7H35Ad9dewPQHUokPjgj5dVBS1R6ySVKV47GED3Mc7FEshwxNcWd8YSwnjNcaKsy5eHI76G2TffJ2uYO1OOST9XlHCXozOSPE9qeWtbOK0sXsPd0RWlPRxhoHllMGApv8c9J3ipSWBdVZDa6iC1oWDiM6Suym+Pf47UVQWpjb2S6wFqQ8Eh80VIyqMcNW43Ar8Dbo4XmFkQuBo4DJ+IPW9m8/BJ3KUZx58G/NM59w8z2wT4FXBSCeIe0pxztHR0sbKlkxUtHaxo6WTlho7Eevw9NSlrD/dcuxUwGFlfzYi6KobVhmisCTF5dD2NtSGG1YRoqAkllhtrQzTWVNFYE2JYbWxb7FVbFVAyJTIAZr6GtyoYoL663NF01xVJ754Q784QT+5SE7349s5IZnmUjnCkx/O0dHQl9m8PR2nvGlgNfU0oQF1K4ldbFaSuKhB790lhPNGrSymrDfkvlvU1IRqqg9RXh6ivDtJQk1yurw5RHVIt4lBW8sTNOfeEmTVlFO8FLHbOvQ1gZrcBxzjnLsXXzvVkNVBTjDiHkrbOCB+ta+ejte18tK6ND9e28/Hadv++rp3l6ztYsaGTzh4+xEbUVTGmsZoxDdVMHl3PLpOqEknZiLoqRtZXMbKumpH1fn1EfRWN1SF9KxWRXoWCvrmzoQyf9NGo77IR7wvbHuty0R5bb+uM0N4VpT2lv2x8W7LM1xzGj13bFo6tR5P9bcORPsVVFbSURC5IQ00yqauvDtJQHaK+JpgoSySBNbFtKevxL8H11UF9Ca4Qg6WP20Rgacp6M7B3Tzub2WeBGcBIfO1dtn3OAM4AmDx5csECrUTt4QjNq9tYuqqVJSmvpata+WhdO2taw92OGV4bYtMRdWwyopatxw9jbGM1YxtrfILWWMOYBr8+uqFa3/5EZKMUCJhvBq0OFvU68T698SSutbOLDR0RNsT68G7ojNDa0cWGzghtnV2J9dbOCK2dfr/WjgjL1rfT2pFc39DZlXdTc8CgoSa1lSNEY22VX08pG1abbb0qeUxNSKOii2ywJG594pz7K/DXXva5DrgOYPr06RU2pqzvnHMsW9/Bmx+38Oay9Sxe1sKby1pYstInZ6nqqoJMHl3P5qPr2LNpNBNG1DJheC2bjqj1yyNqqa+uyP81REQqjpkl+suNKuB54wlha2eEDR1dtIX9e3w9nvR1GwAWe1/bFub91a2Jsg2d+dUM1lcHE8ndsJoQw+uq/Ku2KtESM7wulFxOKR9WG9KAkl4Mlr/O7wObp6xPipVJFuFIlNc/Ws/C99fy8vtrefXDdby5rIX17V2JfYbXhthmk2EcsM1YJo+ujyVq/n1sY7WqxEVENnKpCeHohoF3YIxEXY5EL9wt8Vsf22ddW5j3V7exrt0PTgtHctelNNaEEklcMtFLJncj6qoY1VDN6PpqRjVUMbqhmlH11dRWFbdmdLAYLInb88A2ZjYFn7CdAHyxvCENHsvWt/P8O6t57p2VvLB0Da9+tD7R32xYTYgdNhvOsdMmss0mjWw9rpGtN2lkXGONkjMRESmYYMAYXutryPrLOT+X5rq2rsQsA+ti72vbwonkLr59XVuYJataE9tbc9T61VcHGRVL5kbVVycSutEN1clErz6W9MVeVRVYu1eO6UDmAIcAY82sGbjIOXe9mc0CHsCPJL3BObeo1LENFi0dXTz55gr+8cYynn17FW+v2AD4/yl3mTSCr+zXxNSJI9h54gi2GF2vTv4iIlIRzOIDK0JMGFHb5+PDkWhs9oJOVm0Is2pDJ6tbO/37hk5WtXayptWXL1nVyqoNnWmtUZlGN1QzrrGGscP8+7hhNYyNvcdfYxtrGF1fPWj+1pqrtCnF+2H69Olu/vz55Q4jp5UtHfzt5Q956JWPefbtVXRGogyrDbFX02j23nI0e00Zw06bDa/IbwciIiLlEo5EWZ2S0K3e0MnKDfGprTpYvt6/VrR0smx9e9aprIIBSwzK+9SOm3DOYdsWNWYzW+Ccm55t22BpKh2SwpEoD73yMXcuaOYfbyynK+rYclwDX9m/iUO3H88eW4xSoiYiIjIAVcEA44fVMn5Y7zV8zjk2dEZiiVxqUpdcripzzZsStzJY1x7mtueWcONT7/LB2nY2GV7D6QdM4TO7T2T7CcPLHZ6IiMiQZGaJaU2mjG0odzhZKXEroY6uCH96ZglXPfoma1rD7LPlaC4+ZiqHbj9e896IiIhIr5S4lch/lqzmO3e8yNvLN7D/1mM4f+b27DJpZLnDEhERkQqixK3InHP87tHF/PrhN9h0RB3/+5U9OWS7cZqqQ0RERPpMiVsRhSNRvnP7i8x78QOOmbYZPz12KsMGMP+NiIiIDG1K3IokEnWc/ZcXuPelDzl3xnZ885CtVMsmIiIiA6LErUh+/dAb3PvSh1x4xPb898FblTscERER2QhokrAiePqtlVz9+GI+P32SkjYREREpGCVuBdYViTJ73iImjapj9tE7lTscERER2YgocSuwe176gNc/Xs+FR+xAfbVaokVERKRwlLgV2E3/eo8txzVwxNQJ5Q5FRERENjJK3Apo8bIWXli6hi/tvYVGkIqIiEjBKXEroEdf+xiAGaptExERkSJQ4lZA/3hjOdttMoyJI+vKHYqIiIhshJS4FYhzjpeb17JH06hyhyIiIiIbKSVuBdK8uo117V3stNnwcociIiIiG6mKnK/CzA4ETsLHv6Nzbr8yh8TiZS0AbD9hWJkjERERkY1VyWvczOwGM1tmZgszymea2etmttjMLsh1DufcP51zXwfuBW4qZrz5+mBtGwCbqX+biIiIFEk5atxuBH4H3BwvMLMgcDVwGNAMPG9m84AgcGnG8ac555bFlr8InF7sgPPx4Zp2ggFj/LDacociIiIiG6mSJ27OuSfMrCmjeC9gsXPubQAzuw04xjl3KXBktvOY2WRgrXNufQ/bzwDOAJg8eXJhgs/hw7XtjB9WQzCg+dtERESkOAbL4ISJwNKU9eZYWS6nA//b00bn3HXOuenOuenjxo0rQIi5rW3rZFR9ddGvIyIiIkNXRQ5OAHDOXVTuGFJt6IjQUBMsdxgiIiKyERssNW7vA5unrE+KlVWMDZ1dNNRUbB4sIiIiFWCwJG7PA9uY2RQzqwZOAOaVOaY+2dDRRUO1EjcREREpnnJMBzIHeBrYzsyazex051wXMAt4AHgVuN05t6jUsQ2EmkpFRESk2MoxqvTEHsrvA+4rcTgFs6FDTaUiIiJSXIOlqbTidUSiVIf0zykiIiLFo0yjQELj7mJJ+zPlDkNEREQ2YkrcCiQ48l88ue5X5Q5DRERENmJK3EREREQqhBI3ERERkQqhYZAF9sbqNzD880oT72aJ5eRbch8z61YW3zf1HD0dF5fteqn79OVcAQsQsEBiW8ACBAiAQYDYNtNzWUVEREpJiVuBHTfvuHKHUFKGdUvwsiZ7seV4UpnrmMT21H3Nsh6fbb9sywELELRgYrmv6z1tCwaCGJZ9PeD3M4xgIJg4zixl/35cM3U9GAj69yzLIQslE3Al2SIiGwUlbgV2xcFX4HAAiXdcctm59G0OlyiLy7Zvrn0yr9dt3xz7JGLMcj2HI+qiRF00UZa2TLTHsmzboy6adt6ezhl1UXAkltOOSY2L2H6py7FjuqJdacdEXCTxs6Qup65HXCRt31zrlSZkoUTy2ONyIJR3MphtOddxqe+hQCiReMavmXr9XLGELOS3BfyxifXYuRIv676PklcR2RgocSuww5sOL3cIUgJ5JXrRCI7YejSWIBJNLscSzrT1XpLLbMlmJBrx71mWu1xXIpYu10UkGkkktvFzxJd7O09npJO2aFvinPkcl1lWTqmJaigQoipQlZY8piZ+aeupyeJAj88z2Qxa0J8/I1lNu2ZsW3y/gKnLsshQoMRNpB/iTZCSv3jtasRFsiaOvSWD8W1d0a70ZdeVLIumlyWOi0YIR8OJ5DXzXKnlmcd3Rjtp62rLenw4Gk5LilOPL7XURLAqUJWWEFYFq5LLsW1p+/S2HjtHVaAqkUhm7p9ZVhWoyhlHKBBbT4lJtaIivVPiJiIlEe/bFyRIdbC63OEUVbzmNS05zJIspiabiUQwNTmNRgi7cLdkMzVhTD1/fFvqe3yfcCSctn9ntJPWrtbsx8SOC0eS54i6aNH/3eIJXV4JZTzhC4aosqoe98+ZUGYkkL3tn7YcrOqWBOvLnJSCEjcRkQIzM58IEKImWFPucAoiXiOameSlJXsZyWNqAhh24ezbs5SFo+Fu56iE5DPexJ2Z0GUmeT0uB3vfpy/7Jmo0sySZ8X2DgWBR/02k8JS4iUhFC4fDNDc3097eXu5QSqa2tpZJkyZRVVVVsmsGLEB1sLqia0szayl7Sj57Wk9NKLMtx5PE3vaLJ6HtXe05901NVovFsIImgr3t21PtZa/7puwTsqHdrK7ETUQqWnNzM8OGDaOpqWlIfJg751i5ciXNzc1MmTKl3OFUlGDAN9VXWi2ocy49QcyWMMZrNrPUiKYt55Nc9nC+cDRMa6SVcGd6AtrTNYup19rIfJPGHElrT83jExomsNOYnYr68+X82ct2ZRGRAmhvbx8ySRv4ZtgxY8awfPnycociJWJmPmkIlq6GdaDi/TwzaxnzqYnMq4azhwQ0W1lHVwcboht6PV84Gk6bIqsnM5tm8ouDf1GCf8XslLiJSMUbKklb3FD7eaXyJPp5BkLUUVfucPKWGIGeoyayIdRQ1hgrMnEzsx2B2cBK4BHn3NzyRgQuUs/2jQeWOwwRERHpp2AgOOgHbPRr7LKZ3dvfC5rZDWa2zMwWZpTPNLPXzWyxmV3Qy2mOAK5yzn0D+HJ/YxERKYRgMMi0adMSr8suuyzn/jfeeCOzZs0qUXQisjHpscbNzC53zp1vZp9zzt2RsflrA7jmjcDvgJtTrhUErgYOA5qB581sHhAELs04/jTgFuAiMzsaGDOAWAqn92ZxEdlI1dXV8cILLxTt/F1dXYRCFdlAIiIFlqvG7dPmO1JcmLnBOfdhfy/onHsCWJVRvBew2Dn3tnOuE7gNOMY597Jz7siM17LY61vABcCKbNcxszPMbL6ZzS9ZJ151OxGRFE1NTaxY4T+i5s+fzyGHHNJtn+XLl3Pcccex5557sueee/LUU08BMHv2bE4++WT2339/Tj755FKGLSKDWK6vcPcDq4FGM1uXUm6Ac84NL2AcE4GlKevNwN497WxmTcD3gAYg69AO59x1wHUA06dPV32YyBBw8T2LeOWDdb3v2Ac7bjaci47KPfS/ra2NadOmJdYvvPBCvvCFL+R1/rPOOouzzz6bAw44gCVLljBjxgxeffVVAF555RWefPJJ6uoqp3O3iBRXrsTtB865c83sbufcMSWLKA/OuXeBM8odh4gIDKyp9OGHH+aVV15JrK9bt46WlhYAjj76aCVtIpImV+L2NLA7UNivr9m9D2yesj4pViYikrfeasZKLRQKEY36xyz19GSHaDTKM888Q21tbbdtDQ3lnXZARAafXH3cqs3si8B+ZvbZzFeB43ge2MbMpphZNXACMK/A1xARKammpiYWLFgAwJ133pl1n8MPP5yrrroqsV7MQQ4iUvlyJW5fBw4ERgJHZbyO7O8FzWwOvjZvOzNrNrPTnXNdwCzgAeBV4Hbn3KL+XkNEpJTifdzirwsu8DMaXXTRRZx11llMnz6dYDD73FBXXnkl8+fPZ5dddmHHHXfkmmuuKWXoIlJhzLnc/fZjidX1JYqnKKZPn+7mz59f1GtMvWFvth92EHM/V77HYIgMRa+++io77LBDucMouaH6c4sMBWa2wDk3Pdu2XPO4HeqcexRYna1p1Dn31wLGuFEwzQciIiIiRZRrcMLBwKP4ptFMDlDiJiIiIlJCPSZuzrmLYu+nli4cEREREelJrqbSc3Id6Jz7VeHDEREREZGe5GoqHRZ73w7Yk+T0HEcBzxUzKBERERHpLldT6cUAZvYEsLtzbn1sfTbwt5JEJyIiIiIJueZxi9sE6ExZ74yViYgIEAwGmTZtGlOnTuWoo45izZo1ALz77rtMnTo1sd8f/vAH9thjD1avXp0ou+KKKzCzxMPoRURyySdxuxl4zsxmx2rbngVuLGZQIiKVJP6s0oULFzJ69GiuvvrqbvvccsstXHXVVTzwwAOMGjUKgKVLl/Lggw8yefLkUocsIhWq18TNOfcz4FRgdex1qnPu0mIHJiJSifbdd1/efz/9Ucu33347l112GQ8++CBjx45NlJ999tn8/Oc/x0xzQIpIfnINTkhwzv0b+HeRYxERGZi/XwAfvVzYc07YGY64LK9dI5EIjzzyCKeffnqi7L333mPWrFn85z//YcKECYnyu+++m4kTJ7LrrrsWNl4R2ajl01Qqecn96DAR2XjFn1U6YcIEPv74Yw477LDEtnHjxjF58mRuv/32RFlrayuXXHIJP/7xj8sRrohUsLxq3CQ/auwQKbM8a8YKLd7HrbW1lRkzZnD11Vdz5plnAlBfX899993HgQceyPjx4znppJN46623eOeddxK1bc3Nzey+++4899xzabVyIiKZlLiJiBRIfX09V155Jcceeyzf/OY3E+Xjx4/n/vvv55BDDmHs2LHMmDGDZcuWJbY3NTUxf/78tP5vIiLZ9Kup1MyuK3QgIiIbg912241ddtmFOXPmpJVPmTKFefPmcdppp/Hcc5rDXET6p781btcWNAoRkQrW0tKStn7PPfcklhcuXJhY3nXXXbuNOAU/35uISD76lLiZWQBodM4tKFI8IiIiItKDXptKzexWMxtuZg3AQuAVMzu3+KElrr+lmV1vZnNzlYmIiIhs7PLp47ajc24dcCzwd2AKcHI+JzezG8xsmZktzCifaWavm9liM7sg1zmcc287507vrUxERERkY5dP4lZlZlX4xG2ecy5M/pOW3QjMTC0wsyBwNXAEsCNwopntaGY7m9m9Ga/xef8kg4FmPxcREZEiyqeP27XAu8CLwBNmtgWwLp+TO+eeMLOmjOK9gMXOubcBzOw24JjYY7SOzC9sERERkaEnn2eVXumcm+ic+7Tz3gM+MYBrTgSWpqw3x8qyMrMxZnYNsJuZXdhTWZbjzjCz+WY2f/ny5QMIV0RERGRwyGdwwhgzu9LM/m1mC8zst8CIEsQGgHNupXPu6865reIPt89WluW465xz051z08eNG1eqcEVkiLrrrrswM1577TXAT/FRV1fHtGnT2HXXXdlvv/14/fXX045ZsmQJjY2N/PKXvyxHyCJSgfLp43YbsBw4Djg+tvyXAVzzfWDzlPVJsbIKp2eVigxlc+bM4YADDkibeHerrbbihRde4MUXX+SUU07hkksuSTvmnHPO4Ygjjih1qCJSwfJJ3DZ1zv3EOfdO7PVTYJMBXPN5YBszm2Jm1cAJwLwBnE9EpKxaWlp48sknuf7667ntttuy7rNu3TpGjRqVWL/rrruYMmUKO+20U6nCFJGNQD6DEx40sxOA22PrxwMP5HNyM5sDHAKMNbNm4CLn3PVmNit2jiBwg3NuUZ8jH4Q0plSkvC5/7nJeW/VaQc+5/ejtOX+v83Puc/fddzNz5ky23XZbxowZw4IFCxgzZgxvvfUW06ZNY/369bS2tvLss88CPtG7/PLLeeihh9RMKiJ90mONm5mtN7N1wNeAW4HO2Os24Ix8Tu6cO9E5t6lzrso5N8k5d32s/D7n3LaxPmo/G/iPISJSPnPmzOGEE04A4IQTTkg0l8abSt966y1+85vfcMYZ/qNz9uzZnH322TQ2NpYtZhGpTD3WuDnnhpUyEBGRgeqtZqwYVq1axaOPPsrLL7+MmRGJRDAzvvWtb6Xtd/TRR3PqqacC8OyzzzJ37lzOO+881qxZQyAQoLa2llmzZpU8fhGpLHk9q9TMRgHbALXxMufcE8UKSkSkUsydO5eTTz6Za6+9NlF28MEHs3Tp0rT9nnzySbbaaisA/vnPfybKZ8+eTWNjo5I2EclLr4mbmX0VOAs/+vMFYB/gaeDQ4oYmIjL4zZkzh/PPT6/pO+6447j00ksTfdycc1RXV/PHP/6xTFGKyMYinxq3s4A9gWecc58ws+2BS3o5ZkhxTlOBiAxVjz32WLeyM888kzPPPDOv42fPnl3giERkY5bPdCDtzrl2ADOrcc69BmxX3LBEREREJFM+NW7NZjYSuAt4yMxWA+8VNywRERERydRr4uac+0xscbaZPYZ/3NX9RY1KRKQPnHOYDZ2ZFNU9Q2ToymtUaZxz7h/FCkREpD9qa2tZuXIlY8aMGRLJm3OOlStXUltb2/vOIrLR6VPiJrmZnp0gUnKTJk2iubmZ5cuXlzuUkqmtrWXSpEnlDkNEykCJm4hUtKqqKqZMmVLuMERESiKfUaXSC3U3ERERkVJQ4iYiIiJSIZS4iYiIiFQIJW4iIiIiFUKJm4iIiEiFUOImIiIiUiEGfeJmZlua2fVmNjelbAczu8bM5prZN8oZH4AGlYqIiEgpFDVxM7MbzGyZmS3MKJ9pZq+b2WIzuyDXOZxzbzvnTs8oe9U593Xg88D+hY+8f4bApO0iIiJSRsWucbsRmJlaYGZB4GrgCGBH4EQz29HMdjazezNe43s6sZkdDfwNuK944YuIiIgMHkV9coJz7gkza8oo3gtY7Jx7G8DMbgOOcc5dChzZh3PPA+aZ2d+AWwsT8UCowVRERESKqxx93CYCS1PWm2NlWZnZGDO7BtjNzC6MlR1iZlea2bX0UONmZmeY2Xwzm1+6ZxiqrVRERESKZ9A/q9Q5txL4ekbZ48DjvRx3HXAdwPTp01UdJiIiIhWvHDVu7wObp6xPipVVLKeHlYqIiEgJlCNxex7YxsymmFk1cAIwrwxxiIiIiFSUYk8HMgd4GtjOzJrN7HTnXBcwC3gAeBW43Tm3qJhxiIiIiGwMij2q9MQeyu9D03iIiIiI9Mmgf3KCiIiIiHhK3EREREQqhBK3AtCYUhERESkFJW4FpOl3RUREpJiUuImIiIhUCCVuIiIiIhVCiVtBqbFUREREikeJm4iIiEiFUOJWAHpUqYiIiJSCErcCMrWUioiISBEpcRMRERGpEErcRERERCqEEjcRERGRCqHETSRfraugbTVEI8myrk748CX/mn8DrG3225e9Ci3LyxeriMhQtWFF97K21dDVAS/+BWaPgF/tBK/M6/kcratg6fPp54pG/XnKLFTuAETK6m/fhef/0L18u0/D6/fB1ONg4Z3Zj520JzQ/37frTdoLvvqQ/wDAQSCY3OYcuGh62VARjcJr98DEPeDxS+E/f0rZaHDkr+Heb0NVA3znNagdXrZQRfISjUK0C0LVfv3tx2Gz3eHNB+HO031ZVT2EW5PHnDQXtjms5KFWrNXvwm937d+x65rh9pOT61MOgv3Ogkl7wOVNuY8dPhHOeaV/1y0AJW4F4PSYeS8ahQ3LoHET6GqHUG3/htouugteuQsW/R+cdCdsNg3uvwBevsNv/9JfYfK+0NkCtSMgVJN+/PI34MlfwVaHwl+/1r+f5fX7/MskGXIAABZZSURBVHtPSRv0PWkDaH7Of9vrzaa7wo7HQssyGLs17PIFX5NXN7L3Y6PR2L9/TfckMBKGZ6+FB78Ppz0Io7eEYBXUDCtOwhjpgqXP+J+j5WP41+9gzJbwX7+GVW/DrZ/L4yTOJ20A4Q1w2eZ++ZAL4aBzfbJrQQioAUEGgQU3wT1n5rdvatIG8Ofjk8tfngcTd/e/m6Wy7gP41Q4QqIJZz8HHi+AvX/Lb9p0FM35WuliiUWhd6T/zglXdt19zAHz0cuGu984T/pWPnT5TuOv2g7lBPgmZmW0JfB8Y4Zw7PlZ2CPATYBFwm3Pu8VznmD59ups/f37RYuzoirDHzXux68gZ/PmzlxTtOmURCcM9Z/nq4jcf8GV7fg3+65d+uXMDLPwrzJuV/fijr/IfPL39j97ZCstfg/nXZ9S2lNC2R8CyRbBmSXr5J38Eu53sf46qOnh5rv/GXDsSPv0L2OGo9CS1qxMCIXjjfnjlblj/IYzbHt590p9/oL76CEya3r38ih2gbZVP2jLNXuubBVK/YWba+fNwXJbax7iVb0F1A9SP9Yn15nvByMnJ7WuW+hqG0VPgjlNh0V/z/5ky1Y2Cz1zr78V9383vmJPm+vgm76u5eaQ8XpgDd329b8dU1cPMy/wXqp4+H/Y/C6afDo/+FF6+Hb5yHzTtP/B4M/X2pXLEZDjrxb5/SXrhVhi+Gdx8DOz+Zdjuv2DbGbD0WdhsN/8lc8kz/ud795/+83L5a8njtz/Sf9YM2xTWfwTPXJ1+/m8+A2O39Ynwug9hzFb+i+g7//TJ2H7/zyd/VXXpx0Uj/jqP/jT5ZR2gYTyc+6Zf7urwCW3NcP/5UlXbt5+9n8xsgXMuywd9kRM3M7sBOBJY5pybmlI+E/gtEAT+6Jy7LI9zzU1J3A4GLgA+Bn7qnFuc69jSJW4z+fNnS/iNpNhaV8HPpxT2nKFa+O6bvqkrn5qnYvvao755rlSWveo/qNtWw3UHd99+7DX5ffD/cAW8dDt8+KJP4p74Bax4Y+DxTTkYTsnS7+P9BfCHQ/M7R/1YaM3SxySXU++HLfbNb98PX4JrD+x5+5SD4JR7+nZ9kXw4BxfHar2rh0FXGxz4HfjE9/wf+J+O99vGbA0zL4fJ+8Dd34Qx28Anf+hroIM5GrrefRIe+bFPaHpz/P/CJlPhkYthl8/Djsdk3y/SBT8ZAwdfAJFO3yrSdBDs+gUf8xO/8K++OO56Xxv42CX+i+3Iyf6zaMGNMHk//0W+uhGmHOhbTopl91Pg6CuLd/4yKmfidhDQAtwcT9zMLAi8ARwGNAPPAyfik7hLM05xmnNuWey41MQt4JyLmtkmwK+ccyflimNIJG7OwavzfN+sbNXKfbH0ebj+U/07dtQUWP3OwK4P8KNV/htTx3roaIHG8cmmvGevg8UP+w/L1OSndgS0r/XJz7YzfO1P4/j087av84lTrg/PcuhY7z/oMmuK2tf52rvV78L/5JnYAOz9Df+z1wzrXmN19iIYMclfs6sDLJCeoB90Lhz6A2ieD3/8ZL9/pIQLlvh/82jE/yyBIDzwPdj1RNh0l/6fN9wGNx7p7/tbj3Tffsbj/tt8KbWv9b8/2/Tz90cGt4tH+ab53vxodWGa7t/7F/zvEfntu+Mx8Jnr0muElj4H1/exz9x33oBhm/iEq2YYbB37f/m+8+C5a/t2rv6a9iWfjG6+F/zff/uWi0z7fxsOu7g08ZRB2RK32MWbgHtTErd9gdnOuRmx9QsBnHOZSVvmeRKJW0pZNXBrZnmmjTpxi3TB7/eGlSmVjk0HwntPdf+A+cGy7v3BssmsCTvvHagf7Ze7OuGn49K3n/gX30Q37Yt+Pf7/1M1H+xqcR3+S388ydjt/nZPuKG2/jkrStjp3x9mL1nRP/hY/DH86zi//cEXPiX28CTiXvb8BL86Bz9/kmz1ymb029/ZCe/AH8K+rupcfdz3snPMjonfrP/L9bW48Eo643A9aiX+R6GyFjnW+KeXnW/paGEj+AezcADfM8M1DWx7ia2EK2ZT7/r/hD5/wfTpPLmLtxkCsfAuu2t0vZ/5/0bYarj0Yjrna19CsWeqbxIIhWP2e72owYpL/YtOyzNfulKK5asMKuPlY30y58/H+9+OvX+39uO8uhsZxve/XF/F/v7NehLrRyX6e2Rz5a7//07/r2zU+d5PvWzs6RyvL38+HZ6/p23nB3/NoJPk7E43Ckn/5L45P/db/Pnzprz33s3XOv9YuhVuOha89ll9/3wo22BK344GZzrmvxtZPBvZ2zmXtJGVmY4Cf4Wvo/uicu9TMPgvMAEYC/5Otj5uZnQGcATB58uQ93nvvvQL/ZEllSdwiXfDSbXD3t/I/ZuIe8MXboWGsX29eAB+/DHt8JbnPQz/yv0hx//Ur2LOXP+Z9cc9Zvjo97sJmePwyOPSHJes7UPHuPdtPPbL/WfDJi/wHmllhBhfkar7+9kIYmeMPBiS/lZ/9CoyYOPB4+iMShp+MTS877Mew6xd9zeXme/btfJ2tcMmm3cvPfw8+egluOqrnY//7CbjuE+BSppDZ4gA49W/5Xds5eOD7/t/y6d/Dsb/3fXZ2+gwsedonE0ufSe7/1Uf8/x/xQTubToNtDofJe+d3vYFoXeUHBe3zLV/zGAjBqrf89S9O+SN7yj0w50Q/uGi3L/Xcp/U7r8MV22Xftv2R8IU/FbcvY+rvwl7/nV7blJp8RiOw+BFfM3TG4zBqi+LFFLdhJfxuOhz1W9/H9uJekpjvf+QH8MRHuKY2637yR765N1/RqG9yBV9b/7dz4D+3+Gbbrz/p/z0GW4tGBaroxK0Qil3j1h6OMP2WIidu4Tb/ATf9dF8Ff9c34YU/9/08o7eEM//jP1gvi3Us/9Fq/4flhpnwfsq/UzFqTFJrjA79IRyUZ8dzKY2W5fDLrZPr332ze3NzJcg1TUBf+8AVoy/maQ/ml0zd+dXkaOqBKEXt513fghdKPLBo68Ng8UPww5WFTxZ6uu+5aq3LZcNK+PdNvr9bpq8/BROmdi/vbPUDiHb7UvHjkz7LlbiVY/z8+0Dq1/ZJsbKKV9RxbD+b4Psq3TAD7r8wv6Tt9Id9UtaQ8od31dvwu72SSRvAj0f586cmbee+VbjYU9WN8n9ELlqjpG0wahzn7038VYlJG8CoJh//VlkGVLzzhP+jvOzVZFnzfD/tQThjRO7Ds/t+7YvWZC//QcqEzDcc3vPxG1bAn473sRQiaQPfnFooH77oB8akcq73pG1CP/szbpqSgI/bPn3b4of8+9/PzX2ODSvg+euT3TgGYrAlbQANY+DAc3wNL/ia14vW+M/abEkbQHW9krYKVY76zOeBbcxsCj5hOwH4YhniGHzaVvvRe+F2n+BUN/iJ/lKbEJuf869UqfPrLLwT5p7mm0jiTUJnL0rvl7bi9e7Xjnalr9eNGvjPk4umaxi8NpZ7YwafvwUu7aHJ9vf7+JF2m+0Gc77gy56+yg/O6Or0ffme/HVy/3it1f3f6z4dQdxmu/nrfuL78NjPuh+b6oVbYduZyf6j4PvD/WIrv3xVH0Y7f3dxek1ppubn/SjATGuW+v5jfbnnfzkZ1rwHY7dJDv545x+9H/f5m/yXzjfuT5YdcDZ8anb6fm8+DH8+Lrl+4m1+Kom4te/Dr3dMP6Zxk95jXvIv301k60/5z9iGMT3vv/rd5PJmu8MHscT3vAIMvCqmTXctfd9SKblijyqdAxwCjMVP3XGRc+56M/s08Bv8SNIbnHNF7RhWqqbSaSNn8qeBNJX2VDU/5aCeJwa0IFy0amDnz/Spi+GAb+e3r0gl6Etz5+y18POt0qc0ufB9qGn0yyve9P2L4k6YA9t/uvt5wm2+JnvXE+EzsQ7dP9/SD3LIvF5vccb3ydweL49G4MexBHDvb8ARl8EdX0mfiiHzD/rHi+B/9vNziO3zjezXzRpLSgzx0ZPxss/+secO/H1JKOLnO/oqP+9Xb3H0dv7Ufcdu57+85ru/EiEpg1xNpUWtcXPOndhD+X3Afdm2VbJh7R8kf+FPmutnwd76U3D4T/3cU7t+wVfVP/FLPwDga4/6pqn7zvOTKvYk12zO572df4BbfTL7tAmZ9v7v/M8pUgnOec1PzumcT5xyNVW+eFv3eejiSRv4mqYv350cVZstaQM/2ee3X4bGCcmyb78Ml2yWvl/zAn/Op36T/TwnpTy949y3Yc27fgDGJjslywNBOOMf8OZDyY7mn7nOP13i6r38euuq9Nq9e2Jfzu6/wI+SrRvtY9jnG762P+6lO3wydu5bvsYqVJccOXvFtv5niutc72ulPl7or9e4iW9F2DLLnIW5XNgM7z0NW+eYjmbWfP9s4FuO9euL7oKdju2+33MZk0rHWxzig3ri3v6Hb3lIbV4/sod7IlJGGvpRQFuueDS5En90yeKH/Qt84vbRS/DYT/361XvC1ONh4dz+XfCQ7/VtSHTmiMA9v5b9OZ0hje6Ujczw1FGhOZoUwY8OTLVflscXNR0EE3bufTRe6pMlwM9nt8nOfjR33B8P9ROJ/vum7OdIHZ3bMKbnJr7NpvlXXKjazyYf98D3kjV/kN7l4vZTYI9T/NQ9bavTH20Ur0F74pe+Ji/VhuXpk7fGm36nHJQsy3di5VQ1w2DbHMk1+GR3TMq9vOMU2ClL7VhPT95Y8jRssV9y/eaj/fusBcmynia1FSkjJW4FklcPkTVL/aMzUmU2m+Rrq0PhkPP7dsyYbfx7vCkF0hO3C5bq4d0imQ7PMg9hIOCnPugrM/jGk92b+XpK2sBPudBfZr4mrW2V77O347G+JjDzsWhL/uWfiQs+cUvEdUty+bV7/az78dq2uH9ekVwelmXalGLK7Jt351dhxqXJedQ+fLHnYxc/4v8tljzjp7KIW/9hcjm1hlJkkFDiViB1luX5kJl+k2V0T+ZzMXPZdBp8+IJfjs9m3RfxRztN3if79sznuIlsrKad1L/pdMqhboDJQ1tKH9j4IIxs/n2zf1//kX//1+/gwe8nt69d6l+5lHtgy8t3+Fe8X9q1KTV/gVD6IKx//tK/Mt10ZHFjFBkgJW4F1K+PrA3L09cn7Owf5ZE5g/300+DTv/Tzr7ko1OcYEdWTLfbtPqv3BUsgUOW/gQ/GYe4ixXDM1X4qhNqRfgR1Z4uvOXpxTvp+3/sg+/EDde5bfoqOWz/Xfdupf/d9w975h3/yyEBn4R+9pZ8GKF+jmvz7kqf7dp1pOZ88WH7f/xjuPavnCX8zfaZEj3cS6SMlbuWW+ViqcdvDTp/tnriN3sp3QB5o1X3mH4HaWJNNdf3AzitSSczS+zeBnxYjU2on/UJqGNtzH66a4TBmK/8qhL7+DOFYU2hXR9+OG75Z7/uUUzDk+xTnm7jVFmHiZZECUOJWIG2BAO39aSbobElfdy77w4kL8TgjEelZYJDUOBe6ubGvfeRevNW/+nydQfQZ9ZPxEMmSeJa7KVekAMrx5ISNzqpVbwDwl+EFeDD6dkckl1Nntd58r4GfW0R6ljllxZQ+TmExWO38+cKf84Bzupf19RmwhVKbZWR9tqQN0keh9iZ1RK7IIKLErQBaW/OcADeXc17z/V52jk0j8oNlcNRV/jE5572THFggIsUxeR+oSWke+9Ts4l/zi9nmbyxwrdC+3/Lzv2U69X7/8PHvvtn7OcZnPKngUxfB+e+ml03o4dmwxXbeO3D2K73vA745/Psf53feQjVVixSYmkoHi+EZw+hDNf49UA0hDUkXKYmGMdARG5FYima1UgwIMss+99smO/mR5PmMJq/JMk1QsR+Ll69AoPdn6qb+O1dpnkqpbKpxExHJqkz9oUrVD6vgfekGc/+xwRybSN8ocRMRiSvis5sHxfX6a1AnZSJDixI3EZFsypaslOq6fblOPvsO4uROiadsRJS4iYiIiFQIJW4iIlmpj1uf9h3UtVqDOTaRvjFXKX0sBsDMlgPvleBSY4EVJbiO5E/3ZHDSfRl8dE8GH92TwakU92UL51zW590NicStVMxsvnNuernjkCTdk8FJ92Xw0T0ZfHRPBqdy3xc1lYqIiIhUCCVuIiIiIhVCiVthXVfuAKQb3ZPBSfdl8NE9GXx0Twanst4X9XETERERqRCqcRMRERGpEErc+sjMZprZ62a22MwuyLK9xsz+Etv+rJk1lT7KoSeP+3KOmb1iZi+Z2SNmtkU54hxKersnKfsdZ2bOzDR6rgTyuS9m9vnY78siM7u11DEONXl8fk02s8fM7D+xz7BPlyPOocTMbjCzZWa2sIftZmZXxu7ZS2a2e6liU+LWB2YWBK4GjgB2BE40sx0zdjsdWO2c2xr4NXB5aaMcevK8L/8BpjvndgHmAj8vbZRDS573BDMbBpwFPFvaCIemfO6LmW0DXAjs75zbCfh2yQMdQvL8XfkBcLtzbjfgBOD3pY1ySLoRmJlj+xHANrHXGcD/lCAmQIlbX+0FLHbOve2c6wRuA47J2OcY4KbY8lzgk2aDekrxjUGv98U595hzrjW2+gwwqcQxDjX5/K4A/AT/5aa9lMENYfncl68BVzvnVgM455aVOMahJp974oDhseURwAcljG9Ics49AazKscsxwM3OewYYaWabliI2JW59MxFYmrLeHCvLuo9zrgtYC4wpSXRDVz73JdXpwN+LGpH0ek9iTQubO+f+VsrAhrh8fle2BbY1s6fM7Bkzy1XrIAOXzz2ZDXzJzJqB+4D/V5rQJIe+/t0pmFApLiIyWJjZl4DpwMHljmUoM7MA8CvgK2UORboL4Zt/DsHXTD9hZjs759aUNaqh7UTgRufcFWa2L3CLmU11zkXLHZiUnmrc+uZ9YPOU9Umxsqz7mFkIX629siTRDV353BfM7FPA94GjnXMdJYptqOrtngwDpgKPm9m7wD7APA1QKLp8fleagXnOubBz7h3gDXwiJ8WRzz05HbgdwDn3NFCLf16mlE9ef3eKQYlb3zwPbGNmU8ysGt9JdF7GPvOAU2LLxwOPOk2WV2y93hcz2w24Fp+0qc9O8eW8J865tc65sc65JudcE77f4dHOufnlCXfIyOcz7C58bRtmNhbfdPp2KYMcYvK5J0uATwKY2Q74xG15SaOUTPOAL8dGl+4DrHXOfViKC6uptA+cc11mNgt4AAgCNzjnFpnZj4H5zrl5wPX4auzF+I6NJ5Qv4qEhz/vyC6ARuCM2VmSJc+7osgW9kcvznkiJ5XlfHgAON7NXgAhwrnNOrQZFkuc9+Q7wBzM7Gz9Q4SuqECguM5uD/wIzNta38CKgCsA5dw2+r+GngcVAK3BqyWLTvRcRERGpDGoqFREREakQStxEREREKoQSNxEREZEKocRNREREpEIocRMRERGpEErcRET6wcxGmtk3yx2HiAwtStxERPpnJKDETURKSombiEj/XAZsZWYvmNkvyh2MiAwNmoBXRKQfzKwJuNc5N7XMoYjIEKIaNxEREZEKocRNREREpEIocRMR6Z/1wLByByEiQ4sSNxGRfnDOrQSeMrOFGpwgIqWiwQkiIiIiFUI1biIiIiIVQombiIiISIVQ4iYiIiJSIZS4iYiIiFQIJW4iIiIiFUKJm4iIiEiFUOImIiIiUiGUuImIiIhUiP8P0bmUYk1dJw0AAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 720x216 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "def AB4(x0,m=4):\n",
        "    xsAB = [x0]\n",
        "    x = x0\n",
        "    fs = [ 0.0 for i in range(m)]\n",
        "    for i,t in enumerate(tr):\n",
        "        if len(xsAB) < m : #RK4\n",
        "            f1 = f(x,t)\n",
        "            f2 = f(x+0.5*h*f1, t+0.5*h)\n",
        "            f3 = f(x+0.5*h*f2, t+0.5*h)\n",
        "            f4 = f(x+h*f3, t+h)\n",
        "            x += h*(f1 + 2*f2+ 2*f3 + f4)/6.0\n",
        "            fs[0] = f1; fs[1] = f2; fs[2] = f3; fs[3] = f4\n",
        "        else:\n",
        "            f1,f2,f3 = fs[1:]\n",
        "            f4 = f(x,t)\n",
        "            x += h *(-9*f1 +37*f2 -59*f3 +55*f4 )/24.0 \n",
        "            fs[0] = f1; fs[1] = f2; fs[2] = f3; fs[3] = f4\n",
        "        xsAB += [x]\n",
        "    return xsAB\n",
        "\n",
        "xsAB4 = AB4(x0)\n",
        "#描画\n",
        "fig = plt.figure(figsize=(10,3))\n",
        "plt.xlabel(\"t\");plt.ylabel(\"abs. diff.\")\n",
        "plt.yscale(\"log\")\n",
        "plt.plot(tr,abs(xe-np.array(xs[:-1])),label=\"Euler\")\n",
        "plt.plot(tr,abs(xe-np.array(xsRK4[:-1])),label=\"RK4\")\n",
        "plt.plot(tr,abs(xe-np.array(xsAB4[:-1])),label=\"AB4\")\n",
        "plt.legend()\n",
        "plt.show()\n",
        "plt.close()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dYfoXPujX-vw"
      },
      "source": [
        "### Adams-Moulton法"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2VU4odEJYEE7"
      },
      "source": [
        "AB法は、現在$i$と過去のステップの情報のみを使用したが、  \n",
        "$i+1$ステップの情報を使う(陰解法)ことで、より高精度の結果を得ることができる。\n",
        "\n",
        "3次のAB法($i-2,i-1,i$の情報)に加えて$i+1$番目の情報を用いて  \n",
        "3次のLagrange補間多項式$P_3(t)=\\sum^2_{j=-1}f_{i-j}L_{i-j}$のもと  \n",
        "$x_{i+1}=x_i + \\sum^2_{j=-2} f_{i+j}  \\int^{t_{i+1}}_{t_i} L_{i+j}(t) dt $\n",
        "という  \n",
        "更新式を考える。愚直に上と同様に各積分を評価すると  \n",
        "$x_{i+1}=x_i + \\frac{h}{24}(f_{i-2}-5f_{i-1}+19f_{i} +9f_{i+1}) $\n",
        "という式が得られる。\n",
        "\n",
        "これは、４次精度の公式になっており、3ステップ(4次)のAdams-Moulton法と呼ばれる。  \n",
        "右辺にある$f_{i+1}$の評価に関する計算コストの分だけ、3次のAB法に比べて増えることになる。\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yI-0E3olcM4m"
      },
      "source": [
        "### Adams-Bashforth-Moulton法"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eLY3gCIokyI1"
      },
      "source": [
        "AM法とAB法を組み合わせることで、使いやすい高精度な陽解法を作ろう  \n",
        "というのが予測子修正子法の発想である。\n",
        "\n",
        "4ステップのAB法を予測子, 4ステップのAM法を修正子として採用すると  \n",
        "以下の更新式が得られる:\n",
        "\n",
        "$\\begin{align}\n",
        "x^{P}_{i+1} &= x_i + \\frac{h}{24}(-9f_{i-3}+37f_{i-2}-59f_{i-1}+55f_i) \\\\\n",
        "x_{i+1} &= x_i + \\frac{h}{24}(f_{i-2}-5f_{i-1}+19f_i + 9f(t_{i+1},x^P_{i+1}))\n",
        "\\end{align}\n",
        "$\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 225
        },
        "id": "xu2CrHFXFh_p",
        "outputId": "7c978fc2-0411-46a9-c1cb-cb48ceaee61c"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAADQCAYAAAC+9+0/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZwcVb3+8c+3l9myLyQBQkjYAwECBlREiQsSvMhy8SKIiIjiAj8RFAW9XuIGuOCGeAUBAe8liOiFgMgiiBhZgxIhhCWBkAwCCdln7+X8/ji9VPd09/QkvUxnnvfr1emqU6eqzkxnup8+VXXKnHOIiIiIyNAXqncDRERERKQ8Cm4iIiIiDULBTURERKRBKLiJiIiINAgFNxEREZEGoeAmIiIi0iAi9W5ALUycONFNnz693s0QERERGdCTTz75pnNuh0LLhkVwmz59OosXL653M0REREQGZGavFFumQ6UiIiIiDULBTURERKRBKLiJiIiINAgFNxEREZEGMSwuThAZipJJRyyZJJ5wxBPZ6VgiSTzpiCeSxBKOeDJJLDWdSKYezuGcI5GERNKRdC7znEyVJ1P1cpfnlaemg+VJ53DOt9GlplOzqens8kydIstdqiy9Beeydcgs88wgZIYZgH+2YDlgfmFqmRGy1HRqOalyv0522lIrGdl9GBAKGSEzwiFSz/6RmTYjFMpbni5LzfttEKibvz5F14uk6kdDISLh1HQ4RDhklfuPJiLbFQU3GVYSSUd3LEF3X+oR84/eWIK+RJLeWJLeeJK+RILeWDJQlqAv7pdlH7llffFE6tkHrdww5gNYMJglkm7gBtdQyHzAMLNM6AECASo1n/onHS0sEIQy86n1UiV5y7OBKr29dMBz+BCZDnvp0OhcujwV+fLLA3UJzCdT9d3Q+lUPyAwiISMSCvnnsBEJZ6ejIR/uCpdlw180bIRDIaKp8nAoRDSc2m7YUvswmiIhouHUIxKiORwiGvHbaUqVNYVDgXpGc3Cd1LKm1LJIWAdzRKqlIYObmc0FvgUsBW52zj1Y1wZJxTnn6I0n6eiN09ETp6M3zpbUc0dvjI7eBB09cbr64nT3JegKhLGuWIKevgRdsThdfelpv6w3ntzqNplBUzhEcyREUyRMcyREc9R/WDVH/fzI5ghNbaHMh2c088Ga/cD0H2x+uinSf7lfN1A35D84o6FAj1C6B6dIj0455fk9SMNBugcxHeiSzpFM0r8HMqdXkn49l/17OclZP7u8f49ocP10gPe9qake1vR0qtfV9766TC9sPOlyemQTyWCdJN2x9DaTge0nM/uJJ5MkAl8q4lX4AhEyMqEvE/Yilgp22bJsKLRM8GuOhGmO+r+z5sDfWcHpgepGQgqRst2peXAzs+uAY4A1zrlZgfJ5wE+AMHCNc+6yEptxQAfQArRXsbmyDXpiCTZ3x9jYHWNjV4yNXX1s6o6xKT3f3cem7jhbemJ05gQzH9bK+UAxg7ZomNamCK1NIdqiEVqawrRFw0wa1UJrU5jWaJi21HNrU3A6knrOvtE3pT4AmiKhzBt/uiwatkyvkzSmdO9gCL2OaS4VIGMJ31vclwqE6Z7jTFncl8cSvofZH77P9jD3BdbJ3VaSWNzRl8huJ1je1R3LKUv3cPfGk/TEEmxrrgyHLPB3XDzgBUNgazRMS1OYlkg48x7SGvXL0+8jrdEwLalHdj5ESyQ8bL4ISX3Uo8fteuBnwI3pAjMLA1cCR+KD2BNmthAf4i7NW/8TwF+dc38xs8nAD4FTa9DuYc05R0dvnHUdfbzZ0cubHX2s6+zNzKefg6GsJ1a8dytkMLatiTGtUUa1RBjZHGHa+DZGtkQY1RxhRHMkMz2yJcLI5igjmyOMakktSz1aoiGFKZFtYOZ7eKPhEG1N9W5Nf/FE7ukJ6dMZ0uEuGPTSy/sS+eVJemOJotvp6I1n6vfEkvTEt62HvjkSojUQ/FqiYVqjodSzD4XpoNcaKGuJ+C+Wbc0RRjSFaWuK0NYUZkRzdrqtKUJTRL2Iw1nNg5tz7iEzm55XfCiw3Dn3EoCZ3Qwc55y7FN87V8wGoLka7RxOuvsSvL65h9c39fD65m5e29TDG5t6/PPmHtZu6eXNzj76iryJjWmNMmFkExNGNDFtfBsHTI1mQtmY1ihj26KMbW1ibJufH9MWZWRTRN9KRWRAkbA/3DmiDu/0yaQ/ZSN9LmxP6pSLntR8d1+CnniSnsD5sull2TLfc5hed1N3LDWfzJ5vG0sMql3RsAWCXJgRzdlQ19YUZkRThLbmcKYsEwKbU8sC8+kvwW1NYX0JbhBD5Ry3nYHVgfl24K3FKpvZvwNHAWPxvXeF6pwFnAUwbdq0ijW0EfXEErRv6Gb1+i5WBR6r13fx+uYeNnbF+q0zuiXCjmNamTymhT0mjWLiyCYmjmz2AW1kMxNG+PnxI5r07U9EtkuhkPnDoE3hqu4nfU5vOsR19cXp7E3QmTqHt7MvQVdvnM6+BN198cx8V1+Crj5fr6s3wZotPXT1Zuc7++JlH2oOGYxoDh7liDCyJernA2WjWgrNR7PrNEd0VXSVDZXgNijOud8Dvx+gztXA1QBz5sxpsGvKBs85x5otvbz4RgcvrtnC8jUdvLimg1XrfDgLao2GmTa+jV3Gt3LI9PFMGdPClNEt7DimxU+PaaGtqSH/a4iINBwzy5wvN66C200Hwq6+BJ29cbpj/jk9nw59/S4ASz1v6o7x6oauTFlnX3k9g21N4Uy4G9UcYXRr1D9aopkjMaNbI9npQPmologuKBnAUPl0fhXYJTA/NVUmBcQSSZ5/fQvPvLqJp1/dxLLXNvPimg629MQzdUa3RNhz8igO33Mi08a3pYKaf544skld4iIi27lgIBw/YttPYEwkXYmgF+sX/Lak6mzujvHqhm429/iL02KJ0n0pI5sjmRCXDXrZcDemNcq4EU2Mb2ti3Igo40c0Ma6tiZZodXtGh4qhEtyeAPY0sxn4wHYy8JH6NmnoWLOlhyde3sDjL6/jqdUbWfb6lsz5ZqOaI8zcaTTHz96ZPSePZI8dRrLH5JHsMLJZ4UxERComHDJGt/gesq3lnB9Lc3N3PDPKwObU86buWCbcpZdv7o6xan1XZnlXiV6/tqYw41JhblxbUybQjR/RlA16banQl3pEG7B3rx7DgSwA5gITzawduNg5d62ZnQPcg7+S9Drn3NJat22o6OiNs+jFN/nLC2t47KX1vPRmJ+D/Ux4wdQwfP2w6s3Yew/47j2HX8W06yV9ERBqCWfrCighTxrQMev1YIpkavaCP9Z0x1nf2saGrzz939rG+q4+NXb581fou1nf25RyNyjd+RBM7jGxm4ij/vMOoZiamntOPiSObGd/WNGQ+a8012pDiW2HOnDlu8eLF9W5GSes6evnD069x37Nv8NhL6+lLJBnVEuHQ6eN5627jOXTGBPbbaXRDfjsQERGpl1giyYZAoNvQ2ce6zvTQVr2s3eIfb3b0sWZLT8GhrMIhy1yU9759J3P+kXtVtc1m9qRzbk6hZUPlUOmwFEskue/ZN/jdk+385YW1xJOO3XYYwcffMZ337DOJt+w6TkFNRERkG0TDISaNamHSqIF7+JxzdPYlUkEuGOqy09E697wpuNXB5p4YNz++iuv/tpJ/beph8uhmzjx8BiccvDP7TBld7+aJiIgMS2aWGdZkxsQR9W5OQQpuNdQbT/A/j67iigdeZGNXjLftNp5vHDeL9+wzSePeiIiIyIAU3GrkH6s28MXfLuGltZ28Y48JfGXePhwwdWy9myUiIiINRMGtypxz/OyB5fzoTy+w45hWfvXxQ5i79w4aqkNEREQGTcGtimKJJF+8ZQkLl/yL42bvxLePn8WobRj/RkRERIY3BbcqSSQd5/3mKe7852tccNTefG7u7uplExERkW2i4FYlP7rvBe7852tcdPQ+fPqI3evdHBEREdkOaJCwKnhkxTqufHA5J82ZqtAmIiIiFaPgVmHxRJL5C5cydVwr84/dr97NERERke2IgluF3fHPf/H8G1u46OiZtDXpSLSIiIhUjoJbhd3w8CvstsMIjp41pd5NERERke2MglsFLV/TwVOrN/LRt+6qK0hFRESk4hTcKuiB594A4Cj1tomIiEgVKLhV0F9eWMvek0ex89jWejdFREREtkMKbhXinOPp9k28Zfq4ejdFREREtlMKbhXSvqGbzT1x9ttpdL2bIiIiItuphhyvwszeCZyKb/++zrnD6twklq/pAGCfKaPq3BIRERHZXtW8x83MrjOzNWb2TF75PDN73syWm9mFpbbhnPurc+4zwJ3ADdVsb7n+takbgJ10fpuIiIhUST163K4HfgbcmC4wszBwJXAk0A48YWYLgTBwad76n3DOrUlNfwQ4s9oNLsdrG3sIh4xJo1rq3RQRERHZTtU8uDnnHjKz6XnFhwLLnXMvAZjZzcBxzrlLgWMKbcfMpgGbnHNbiiw/CzgLYNq0aZVpfAmvbeph0qhmwiGN3yYiIiLVMVQuTtgZWB2Yb0+VlXIm8KtiC51zVzvn5jjn5uywww4VaGJpm7r7GNfWVPX9iIiIyPDVkBcnADjnLq53G4I6exOMaA7XuxkiIiKyHRsqPW6vArsE5qemyhpGZ1+cEc0Nm4NFRESkAQyV4PYEsKeZzTCzJuBkYGGd2zQonb1xRjQpuImIiEj11GM4kAXAI8DeZtZuZmc65+LAOcA9wDLgFufc0lq3bVvoUKmIiIhUWz2uKj2lSPldwF01bk7FdPbqUKmIiIhU11A5VNrwehNJmiL6dYqIiEj1KGlUSGSH21jV82i9myEiIiLbMQW3CgmPfZhFm39Y72aIiIjIdkzBTURERKRBKLiJiIiINAhdBllhL2x4AcPfrzTzbJaZzj5l65hZv7J03eA2iq2XVmh/wTqD2VbIQoQslFkWshAhQmAQIrXMdF9WERGRWlJwq7ATF55Y7ybUlGH9Al7BsJeaTofKUutklgfrmhVcv1C9QtMhCxG2cGZ6sPPFloVDYQwrPB/y9QwjHApn1jML1N+KfQbnw6Gwfy4wHbFINoArZIuIbBcU3Crs8iMux+EAMs+47LRzucscLlOWVqhuqTr5++tXt0SdTBsL7M/hSLokSZfMlOVMkyxaVmh50iVztltsm0mXBEdmOmedYLtI1QtOp9aJJ+M56yRcIvOzBKeD8wmXyKlbar7RRCySCY9Fp0ORssNgoelS6wWfI6FIJnim9xncf6m2RCzil4X8upn51LYyD+tfR+FVRLYHCm4V9v7p7693E6QGygp6yQSO1HwyFRBJZqdTgTNnfoBwWShsJpIJ/1xgOu7imbbEXZxEMpEJtultpKcH2k5foo/uZHdmm+Wsl19WT8GgGglFiIaiOeExGPxy5oNhcVvXLzNshi3st58XVnP2mVqWrhcynbIsMhwouIlshfQhSClfunc14RIFg+NAYTC9LJ6M5067eLYsmVuWWS+ZIJaMZcJr/raC5fnr9yX76I53F1w/lozlhOLg+rUWDILRUDQnEEbD0ex0allOnYHmU9uIhqKZIJlfP78sGoqWbEcklJoPtEm9oiIDU3ATkZpIn9sXJkxTuKnezamqdM9rTjgsEBaDYTMTBIPhNJkg5mL9wmYwMAa3n14WfE7XiSViOfX7kn10xbsKr5NaL5bIbiPpklX/vaUDXVmBMh34whGiFi1av2SgzAuQA9XPmQ5H+4VgfZmTWlBwExGpMDPzQYAIzeHmejenItI9ovkhLyfs5YXHYACMuVjh5QXKYslYv200QvhMH+LOD3T5Ia/odHjgOoOpm+nRLBAy03XDoXBVfydSeQpuFbLba471o+rdChGR6ghZiKZwU0P3lub3UhYLn8Xmg4Gy0HQ6JA5ULx1Ce+I9JesGw2q1GFbRIDhQ3WK9lwPWDdSJ2PA+rK7gViGXXZ+gJwqcXe+WiIhIIeGQP1TfaL2gzrncgFgoMKZ7Ngv0iOZMlxMui2wvlozRlegi1pcbQIvts5oG7I0sNzSWCK3FDo9PGTGF/SbsV9Wfr+TPXrc9b4daYrD2p1f4mfS3AbPMoLvpMjPLLseK1jfLW4bl1ssUFaiX2VbuwL+YDVA/NRkKgYV8WchS8wYWwkKWmSY1nakfSm0/s37utiyUXhZcL7utgvtJlWW2GwoBRbaVXze9rXAIwn4MNcJhX1dEpAGYmQ8N4Wi9m1K29Hme+b2M5fREltXDWSSAFirrjffSmewccHuxZCxniKxi5k2fx/eP+H4NfouFKbhV2Js//3m9myDlCochFPIhLh3mCswTDmGhsH+2wLJwOBU2w9k6g1zfh8/g+untBdoTDmXrWGp74TCEwtlAGo6kyiNYJLV+6kE4kq0XiWTXDdQrtCxTHg5DsWUKwCJSQOY8z1CEVlrr3ZyyZa5AL9ETOSIyoq5tbMjgZmb7AvOBdcD9zrlb69si2NxqLNt/Gh+/8W4gMAiuc/6Rns57zmT7QnXyyrLj6g6wzbwBeEtvMztIMOTVS/rBazPTST/ILS7ZbxnO4ZIusyxn3jm/btL5fSSTqWV+uUvVL7mfnG2l9klgu/n7yWw7vc2EL08mcImkf04mIZFalkj6+ulliWR2vUSx9R0kEv3XTyRwfX249M+SSO0rGdh2wo/5llmWSOCcb0+wfv/1q39l36CYZUKdhULZgJcTHNPhLy9gButFUgEzZxuBepGwD7SRMBaJ+pAZSQfKSP/5qC+j2LJ0GI1EfRBN100ts0gktTxdN5pdL91zKyLbnXAoPOQv2Niq4GZmdzrnjtnKda8DjgHWOOdmBcrnAT8BwsA1zrnLSmzmaOAK59xfzWwhUPfglq//4cgi9WrQFtl+uHRQDga+ZBIXj/vpRBIScVwi4cuSSVw8kSoLLkv4EBlP4BLpeqltFFyW2kY8HVQL1Iun2xSoF0/gEtn9Z+tll7lEHNcd8/Xi8dxtJFL7SgT2Efc/A7HUOvUIs4FQlxPyIpFsuCwQAC0aDJN5daORVM9nYD4QNvstS8+nQ2qmDVEsmlcejebUtVQZwbLw0P6wEhGvaHAzs+86575iZv/hnPtt3uJPbcM+rwd+BtwY2FcYuBI4EmgHnkgFsjBwad76nwB+DVxsZscCE7ahLSINJX2OHuGwQn+KSyZ94Es/EglcLJYKl76MTHkcFw8uS2TnY/FUOIwXWZaaz4THVN2Cy/LmU21IdnVnQ3WwvQXqZgJqokZ3nAiFsoEymgp/kWDQS5VHojnhz9cPBMZg/WBoTJVl1ikWMiMRLNqUrZ8JmZFs+Azsh2CZDt3LMFCqx+0DZnYhcBGQE9ycc69t7Q6dcw+Z2fS84kOB5c65lwDM7GbgOOfcpfjeuULOTgW+3xdaaGZnAWcBTJs2bWubKyJDnIVC0NSENTXuMBWlZHpWg6E0nu1VzcynHzFf5sNgflk8Wx6LpYKsL8tsJxbcVmA7Odvw6yQ7eiBYP7hOLDgfr00ADYRPi0ZTQbNIyIzkBcBoXshML88PoMFQ2xSYDu4v2pQNpNHcOpkQGo2m6il0yuCUCm53AxuAkWa2OVBugHPOja5gO3YGVgfm24G3FqucCn5fBUYABS/tcM5dDVwNMGfOnIEvExERGYIyF6lEG+eKwkIyPaM54S8dCPvywmcs2zsaXKdU0BwgZPYPq34+2dPdL2Rm1snbTlXDZ/oCoUJBLxoIg0XqZHpECwXMQJ3+4TGvfrBHM1M3tawpb3s6vF4XpYLbfzrnLjCz251zx9WsRWVwzq0k1ZsmIiJD3/bQM5rt5fRhMxseY4FwGcsGxFjM14vH/TmZwTp9efXjscw0RetkA2qyp9vXjRXZf6A3tmrMcgIiTQMEx/R0U4ngGKhTtCezQG9nvzAaLbDN9OH4Br+4qFRwewQ4GNhcok6lvArsEpifmioTEREZEjJXSjc34w/4DH05PZ1Fwp0PiAWCY7pOX6F1U3VzQmqh7fdlAmays7N4GwJlxKo7eG8mDJYKl4GASF64bD1oNuM/8pGqtrGUUsGtycw+AhxmZv+ev9A5V/Dcsq30BLCnmc3AB7aTgfr9VkRERLYDjdjT6ZzLO0ReKETG+vdi9guRRcJo0aAZqJ+qk+zqxsW35NQPjxlT199PqeD2GeBUYCzwwbxljiIXBQzEzBYAc4GJZtYOXOycu9bMzgHuwV9Jep1zbunWbF9EREQal5llesSkv6LBzTm3CFhkZoudc9dWaofOuVOKlN8F3FWp/dSDa+zD5iIiIjLElRrH7T3OuQeADTU4VLpdMI2sJSIiIlVU6lDpEcAD9D9MCttwqFREREREtk6pQ6UXp57PqF1zRERERKSYUodKzy+1onPuh5VvjoiIiIgUU+pQ6ajU897AIcDC1PwHgcer2SgRERER6a/UodJvAJjZQ8DBzrktqfn5wB9q0joRERERySjnrraTgb7AfF+qTERERERqqNSh0rQbgcfN7P9S88cD11etRSIiIiJS0IDBzTn3HTP7I/DOVNEZzrl/VLdZIiIiIpKvnB43nHN/B/5e5baIiIiISAnlnOMmZTDn6t0EERER2c4puImIiIg0CAU3ERERkQaxVcHNzK6udENEREREpLSt7XG7qqKtEBEREZEBlXVVaZqZhYCRzrknq9QeERERESliwB43M7vJzEab2QjgGeBZM7ug+k3L7H83M7vWzG4tVSYiIiKyvSvnUOm+zrnN+Dsm/BGYAZxWzsbN7DozW2Nmz+SVzzOz581suZldWGobzrmXnHNnDlQmIiIisr0rJ7hFzSyKD24LnXMxoNxBy64H5gULzCwMXAkcDewLnGJm+5rZ/mZ2Z95jUtk/yVBgVu8WiIiIyHasnHPcrgJWAkuAh8xsV2BzORt3zj1kZtPzig8FljvnXgIws5uB45xzlwLHlNdsERERkeFnwB4359xPnXM7O+c+4LxXgHdvwz53BlYH5ttTZQWZ2QQz+wVwkJldVKyswHpnmdliM1u8du3abWiuiIiIyNAwYI+bmU0ALgYOxx8iXQR8E1hX3aZ5zrl1wGcGKiuw3tXA1QBz5szR/ahERESk4ZVzjtvNwFrgROBDqenfbMM+XwV2CcxPTZU1NJ3dJiIiItVWTnDb0Tn3Lefcy6nHt4HJ27DPJ4A9zWyGmTUBJwMLt2F7Q4a69URERKSayrk44V4zOxm4JTX/IeCecjZuZguAucBEM2sHLnbOXWtm56S2EQauc84tHXTLt1EsFqO9vZ2enp6KbC/6w58xKxph2bJlFdletbW0tDB16lSi0Wi9myIiIiJlMucK9xOZ2RZ8J5IBI4BkalEI6HDOja5JCytgzpw5bvHixTllL7/8MqNGjWLChAlYBYbx6Fz6DN1tzUycsec2b6vanHOsW7eOLVu2MGPGjHo3R0RERALM7Enn3JxCy4oeKnXOjXLOjU49h5xzkdQj1EihrZienp6KhbZGY2ZMmDChYr2NIiIiUhtl3avUzMYBewIt6TLn3EPValStDMfQljacf3YREZFGVc69Sj8JPIQ/J+0bqef51W3W8BAOh5k9e3bmcdlll5Wsf/3113POOefUqHUiIiIy1JTT43YucAjwqHPu3Wa2D3BJdZvVWIqdJziQ1tZWnnrqqQq3JisejxOJlNWpKiIiIg2gnOFAepxzPQBm1uycew7Yu7rNGt6mT5/Om2++CcDixYuZO3duvzpr167lxBNP5JBDDuGQQw7hb3/7GwDz58/ntNNO4x3veAennXZaLZstIiIiVVZOd0y7mY0FbgPuM7MNwCvVbVZtfeOOpTz7r7Juv1pUsquTRChEtMUHrn13Gs3FH9yv5Drd3d3Mnj07M3/RRRfx4Q9/uKz9nXvuuZx33nkcfvjhrFq1iqOOOiozFMmzzz7LokWLaG1t3cqfRkRERIaiAYObc+6E1OR8M/szMAa4u6qtGia25VDpn/70J5599tnM/ObNm+no6ADg2GOPVWgTERHZDg3qBCjn3F+q1ZB6GqhnbCDOObqeXVqxcdwikQjJpB82r9iQHclkkkcffZSWlpZ+y0aMGLHNbRAREZGhp5xz3KTGpk+fzpNPPgnA7373u4J13v/+93PFFVdk5qt5kYOIiIgMDQpudZQ+xy39uPDCCwG4+OKLOffcc5kzZw7hcLjguj/96U9ZvHgxBxxwAPvuuy+/+MUvatl0ERERqYOit7zanhS65dWyZcuYOXNmRbZf6UOltVLJ34GIiIhUxlbd8kpEREREhhYFNxEREZEGoeAmIiIi0iAU3EREREQahIKbiIiISIMY8sHNzHYzs2vN7NZA2Uwz+4WZ3Wpmn61n+0RERERqparBzcyuM7M1ZvZMXvk8M3vezJab2YWltuGce8k5d2Ze2TLn3GeAk4B3VL7ltREOh5k9ezazZs3igx/8IBs3bgRg5cqVzJo1K1Pvl7/8JW95y1vYsGFDpuzyyy/HzDI3oxcREZHtX7V73K4H5gULzCwMXAkcDewLnGJm+5rZ/mZ2Z95jUrENm9mxwB+Au6rX/OpK36v0mWeeYfz48Vx55ZX96vz617/miiuu4J577mHcuHEArF69mnvvvZdp06bVuskiIiJSR1UNbs65h4D1ecWHAstTPWl9wM3Acc65p51zx+Q91pTY9kLn3NHAqdX7CWrn7W9/O6+++mpO2S233MJll13Gvffey8SJEzPl5513Ht/73vcws1o3U0REROpoUDeZr5CdgdWB+XbgrcUqm9kE4DvAQWZ2kXPuUjObC/w70EyRHjczOws4Cxi4Z+qPF8LrT5f/E/TjaOnqoilk0NLmi6bsD0dfVtbaiUSC+++/nzPPzB4RfuWVVzjnnHP4xz/+wZQpUzLlt99+OzvvvDMHHnjgNrRXREREGlE9gtugOOfWAZ/JK3sQeHCA9a4GrgZ/y6sqNW+bpO9V+uqrrzJz5kyOPPLIzLIddtiB8ePHc8stt3DeeecB0NXVxSWXXMK9995bryaLiIhIHdUjuL0K7BKYn5oqq58ye8aKco6eZ5fSPaKZidPLv1dp+hy3rq4ujjrqKK688ko+//nPA9DW1sZdd93FO9/5TiZNmsSpp57KihUrePnllzO9be3t7Rx88ME8/rQ3ygEAABzUSURBVPjjOb1yIiIisn2qR3B7AtjTzGbgA9vJwEfq0I4ho62tjZ/+9Kccf/zxfO5zn8uUT5o0ibvvvpu5c+cyceJEjjrqKNasyZ72N336dBYvXpxz/puIiIhsv6o9HMgC4BFgbzNrN7MznXNx4BzgHmAZcItzbmk129EIDjroIA444AAWLFiQUz5jxgwWLlzIJz7xCR5//PE6tU5ERESGgqr2uDnnTilSfhcNPIxHpXR0dOTM33HHHZnpZ57JDn134IEH9rviFPx4byIiIjJ8DPk7J4iIiIiIp+AmIiIi0iAU3CpgSI41IiIiItsdBTcRERGRBqHgJiIiItIgFNxEREREGoSCW53ddtttmBnPPfcc4If4aG1tZfbs2Rx44IEcdthhPP/88znrrFq1ipEjR/KDH/ygHk0WERGROlFwq7MFCxZw+OGH5wy8u/vuu/PUU0+xZMkSTj/9dC655JKcdc4//3yOPvroWjdVRERE6kzBrY46OjpYtGgR1157LTfffHPBOps3b2bcuHGZ+dtuu40ZM2aw33771aqZIiIiMkTU416lQ853H/8uz61/bpu2kejsJBkOEX2uFYB9xu/DVw79Ssl1br/9dubNm8dee+3FhAkTePLJJ5kwYQIrVqxg9uzZbNmyha6uLh577DHAB73vfve73HfffTpMKiIiMgypx62OFixYwMknnwzAySefnDlcmj5UumLFCn784x9z1llnATB//nzOO+88Ro4cWbc2i4iISP2oxw0G7BkbSNI5up9dSveIZiZO37OsddavX88DDzzA008/jZmRSCQwM84+++ycesceeyxnnHEGAI899hi33norX/7yl9m4cSOhUIiWlhbOOeecbWq/iIiINAYFtzq59dZbOe2007jqqqsyZUcccQSrV6/Oqbdo0SJ23313AP76179myufPn8/IkSMV2kRERIYRBbc6WbBgAV/5Sm5P34knnsill16aOcfNOUdTUxPXXHNNnVopObrWgxk0j4ZQ2JfF+2Bt6vzIVxfDnu+HUTvCmy9A20QYuUP92isiMhx1vgkjJuaWdW+AaBssvQ3+7ywYPRXmXQr7Hlt4G13rYd0KGD8ju61kEno3Qeu4wuvUiDm3/d9pc86cOW7x4sU5ZcuWLWPmzJkV2f7WHCodCir5O2hYf/gSPPHL/uV7fwCevwtmnQjP/K7wulMPgfYnBre/qYfCJ+/zbwC4bAAEcA5cMrdsuEgm4bk7YOe3wIOXwj/+J7DQ4JgfwZ1fgOgI+OJz0DK6bk0VKUsyCck4RJr8/EsPwk4Hw4v3wu/O9GXRNoh1Zdc59VbY88iaN7VhbVgJPzmwMtua8S447FyY+hb47vTSdUfvDOc/W5n9FmFmTzrn5hRaph43qZxkEjrXwMjJEO+BSIvvoRqspbfBs7fB0v+DU38HO82Guy+Ep3/rl3/09zDt7dDXAS1jINKcu/7aF2DRD2H398DvP7V1P8vzd/nnYqENBh/aANofh/ljBq6344Gw7/HQsQYm7gEHfBiSCWgdO/C6yWTq99/cPwQmYvDYVXDv1+AT98L43SAcheZR1QmMiTisftT/HB1vwMM/gwm7wb/9CNa/BDf9RxkbcT60AcQ64bJd/PTci+BdF/iwa2EI6VorGQKevAHu+Hx5dYOhDeB/P5Sd/thC2Plg/7dZK5v/BT+cCaEonPM4vLEUfvNRv+zt58BR36ldW5JJ6Frn3/PC0f7Lf3E4vP505fb38kP+UY79TqjcfrfCkO9xM7PdgK8BY5xzH0qVzQW+BSwFbnbOPVhqG+pxK6ys30EiBnec67ueX7zHlx3yKfi31HAkfZ3wzO9hYZFz7Y69wr/xDPQfva/LH3JcfG1eb0sN7XU0rFkKG1fllr/3v+Cg0/zPEW2Fp2/135hbxsIHvg8zP5gbUuN9EIrAC3fDs7fDltdgh31g5SK//W31yfthaoEvYpfPhO71PrTlm78Jnl0It5xWfLv7nwQnFuh9TFu3AppG+EPAz94GuxwKY6dll29c7XsYxs+A354BS39f/s+Ur3UcnHCVfy3u+lJ565x6q2/ftLdv3RcGkW311AK47TODWyfaBvMu81+oir0/vONcmHMmPPBtePoW+PhdMP0d297efAN9qRwzDc5dMvgvSU/dBKN3ghuPg4M/Bnv/G+x1FKx+DHY6yH/JXPWo//lW/tW/X64NDNG1zzH+vWbUjrDldXj0ytztf+5RmLiXD8KbX4MJu/svoi//1Yexw/6fD3/R1tz1kgm/nwe+nf2yDjBiElzwop+O9/pA2zzav79EWwb3s2+lUj1uVQ1uZnYdcAywxjk3K1A+D/gJEAaucc5dVsa2bg0EtyOAC4E3gG8755aXWlfBrbABfwdd6+F7Myq700gLfOlFf6irnJ6navvUA/7wXK2sWebfqLs3wNVH9F9+/C/Ke+P/+pvwz1vgtSU+xD30fX9e3baacQScvrB/+atPwi/fU9422iZC15uD2+8Zd8Ouby+v7mv/hKveWXz5jHfB6XcMbv8i5XAOvpHq9W4aBfFueOcX4d1f9R/w357kl03YA+Z9F6a9DW7/HEzYE977dd8DHS5xoGvlIrj/mz7QDORDv4LJs+D+b8ABJ8G+xxWul4jDtybAERdCos8fFZn+Ljjww77ND33fPwbjxGt9b+CfL/FfbMdO8+9FT14P0w7zX+SbRsKMd/ojJ9Vy8Olw7E+rt/06qmdwexfQAdyYDm5mFgZeAI4E2oEngFPwIe7SvE18wjm3JrVeMLiFnHNJM5sM/NA5d2qpdgyL4OYc9GzygcjK+zZU9Hew+gm49n1b145xM2DDy1u3btB/rfffmHq3QG8HjJyUPZT32NWw/E/+zTIYflrG+N/B8b/w3+aScb9eUM9mH5xKvXnWQ+8W/0aX31PUs9n33m1YCf9dZrABeOtn/c/ePKp/j9V5S2HMVL/PeK///xIM6O+6AN7zn9C+GK5571b/SBkXrvK/82TC/yyhMNzzVTjwFNjxgK3fbqwbrj/Gv+4r7u+//KwH/bf5WurZ5P9+9tzKvx8Z2r4xzh+aH8h/bajMoftXHoZflXl7w32PgxOuzu0RWv04XDvIc+a++AKMmuwDV/Mo2CP1f/muL8PjV5Vet1Jmf9SH0V0Ohf/7tD9yke8dX4Ajv1Gb9tRB3YJbaufTgTsDwe3twHzn3FGp+YsAnHP5oS1/O5ngFihrAm7KL8+3XQc353wvTqI3W9Y00p//lW/HA3NCXdHfQX5P2JdfhrbxfjreB9/Ou1LylN/4Q3SzP5JtE8CNx/oenAe+Vd7PMnFvv59Tf1vb8zoaSfeG0ifOXryxf/hb/if4nxP99NffLHy+CGQPAZfy1s/CkgVw0g3+sEcp8zeVXl5p9/4nPHxF//ITr4X9S75FDGzL6/58m+uPgaO/6y9aSX+R6OuC3s3+UMr3dvO9MJD9AOzrhOuO8oeHdpvre2EqeSj31b/DL9/tz+k8rYq9G9ti3Qq44mA/nf//onsDXHUEHHel76HZuNofEgtHYMMr/lSDMVP9+1rHGt+7U4vDVZ1vwo3H+8OU+3/I/338/pMDr/el5ZW/mjz9+zt3CbSOz57nWcgxP/L1H/nZ4PbxHzf4z4jxJY6y/PEr8NgvBrdd8K95MpH9m0kmYdXD/ovj337i/x4++vvi59k65x+bVsOvj4dP/bm8830b2FALbh8C5jnnPpmaPw14q3Ou4ElSZjYB+A6+h+4a59ylZvbvwFHAWOC/C53jZmZnAWcBTJs27S2vvPJKzvKGD27O+bCUfz5WKdG27MnowLKnlzCz7x/wlo9n69z3X/4PKe3ffgiHDPBhPhh3nOu709MuaocHL4P3fL1m5w40vDvPg8XX+Q+U917s/y+YVebiglKHr7/wDIwt8YEB2W/l5z0LY3be9vZsjUQMvpU3FMCR34QDP+J7Lnc5ZHDb6+uCS3bsX/6VV+D1f8INHyy+7qcfgqvfDS6RLdv1cDjjD+Xt2zm452v+d/nIz+H4n/tzdvY7AVY94sPE6kez9T95v///kb5oZ8fZfoiaaW8tb3/bomu9vyjobWf7nsdQBNav8Pv/RuBD9vQ7YMEp/svlQR8tfk7rF5+Hy/cuvGyfY+DD/1PdcxmDfwuHfjq3tykYPpMJWH6/7xk660EYt2v12pTWuQ5+Ngc++BN/ju03BggxX3vdX8CTvsI1eFj3vf/lD/eWK5n0h1zBdwT84Xz4x6/9YdvPLPK/j6F2RKMBNXRwq4Sq97glHd3Lti643XbbbZxwwgksW7aMffbZh5UrVzJz5kz23ntvnHOMGDGCX/3qV+y95548eM9C3v2BE/jlL3/JJ088ErrX89Qzz3PQUafw/a9/gS995mOZ7V7+i1/zpW/9iLVP38/E8YExZ8LNMHlfSMZZ9tj9zLznJN+t7xJw3Tw/FllaNXpMgj1G7/k6vKvME8+lNjrWwg/2yM5/6cX+h5sbQalhAgZ7Dlw1zsX8xL3lhanffTJ7NfW2qEXv521nw1M1vrBojyNh+X3w9XWVDwvFXvdSvdb10rkO/n6DP98t32f+BlNm9S/v6/IXEB300eq3TwatVHCrx/XzrwLBr+1TU2XD0oIFCzj88MMz9ymF7L1KlyxZwumnn84ll1wCry+BzrXMmrkXt9x0o+9tAxbcfjcH7rtXzjZXv/o69z70CNOm7QKT9/fffNMSvf7QavAy6m+Og+9MyQ1tF6yoys9L6zj/IXLxRoW2oWjkDv61ST8aMbQBjJvu2797gQsqXn7IfyivWZYta1/shz2I5V2R+6f5g9/3xRsLl//n2uz0de8vvn7nm/A/H/JtqURoA384tVJeW+IvjAlybuDQNmUrz2fcMRDAd9gnd9ny+/zzHy8ovY3ON+GJa7OncWyLoRbaAEZMgHee73t4wfe8XrzRv9cWCm0ATW0KbQ2qHsHtCWBPM5uROkftZKDAZWzbv46ODhYtWsS1117LzTff7E+m793ie796O6Cvi80bNzBubLYbfNedJtHT1cEba9fhnOPuPz/M0Ucd6QcE3OkgGDed8+Zfzvcu/gpmIX8oYfJ+uTsuNFxEMp47X+2RoTVcw9Blln00MjM46dfFl//8bfDnS+H5u/1FGMvugEdS58jF+/x4XIt+lK0/f5N/vO3swtsD/zdoBu/+Wm75/E3Zw1RpT93kDy8G9XXC93f3geSKQVzt/KWSF9YXH3Nw4+rBh5nfnObHR/zXP7JlL/9l4PVOugH2mpdbdvh52d9r+nFq3tiJp9ycXXb2Y/4wfL6Rkwdu8x/Oh2UL/e+4c13p+htWZqd3Ojg7/eUKXHhVTTse6H9PU+c0/t+vFFXVA9FmtgCYC0w0s3bgYufctWZ2DnAP/krS65xzFRjcauu9fskl9C57buCKRTgg2dVJMhSis8WPE9M8cx+mfPWrJde7/fbbmTdvHnvttRcTJkzgyXt/w4RxY1nx0svMPvgtbOnsoqu7h8fuuTVnvQ/92/v47Z33cdCsfTh4/31oHp09n+f2ex9i5z3248B3Hw+kBi21kP9ACb7RlvK+bwzP0ftl+9M8MvWhX+Sw11/yRiJ64Nv+qtofzswd0uSiwEGBOWfkjiN18gLY5wO52zniy37sqO9M8VfPprVN8Bc5ANz2Wf8cPIx5yU7Z6c3t2el0nfyfI12eDJxD99bPwtGXwW8/nh2K4Y9fhrd+OnfdN5bCfx/mxxB722cp28bU+cJXz81ePZm+UOXfryl+Av/43eAjvxl4+8Erco+9wo//FTRm5/6v6YOXwtwLi29z1cP++ZaP+Yug3ny+9OHj4GH2s/48cJtFaqiqwc05d0qR8ruAuwota2wuezVnvM8HpeZRvjcs1u2vmHTOjyDf8QYLbvpfzv3852HDSk4++nAW3HYP55zxYXbfdSpP3XczAL+5/R7OOu8i7v7f7AfFSR88kg9/9kKeW76SU874LA8/6sf86erq4pJLLuHee+8t3LzmUb5HbyD5b/Aije785/zgnM754FTqUOWSm/uPQ9c8Mjs9cU/42O3ZsJIf2tKirfCFp2HklGzZF57ODWcA7U/6bf7tx4W3E+yBuuAl2LjSX4AR7EkPheGsv8CL92VPND/han93iSsP9fNd67NXhwPckfpid/eF/irZ1vG+DW/7rB9oNO2fv/Vh7IIV/p6NkdbslbOX7+V/prS+Lb5X6o1n/P5GTvZjAO5WYMzCUi5qh1cegT1KDEdzzmLY1O6vMgR/x5X9ju9f7/G8QaXffN4/py/qSXvpL/7IQ/Dw+jFFXhOROtKlHzBgz9hA0hcnxFocY1vzDjn2bsmOAN023ge4La+xfsMmHnjgzzz9zyUYjkQiiZlx9sdPyln92PcfwRnn555wOmXSRKKRCPctWsxPrrkpE9xWrFjByy+/zIEH+m+L7e3tHHzwwTz++ONMmTIFwnmHaooNuRHR1Z2ynRkdvCp0j6LVAH91YNBhBW5fNP1dMGX/ga/GC95ZAvyV3ZP3hzcCYeea9/iBRP9+Q+FtBK/OHTHBPwrZabZ/pEWa/Gjyafd8FU4IDOXQ/nh2+pbT4S2n+6F7ujfk3too3YP20A98T15Q59rcwVv3muff52a8K1tW7sDKQc2jYK8S4Rp82J0QeC1/ezrsV6AXrdidN1Y9Arselp2/MXWz8XOezJYVG9RWpI4U3Gop3gfJGAC3/uFPnHbiB7jqx9/131KBI078JKv/9XrOKoueeIrdp0/tt6lvfvWLrOltIhzOHtLcf//9WbNmTWZ++vTpLF68mIkTU4dS04FsxA5+XKRNy8hx4WrdvFsk3/sLjEMYCvmhDwbLDD67qP8hz2KhDcoeULvo/lrH+4uZlizw97+NtvY/z3XVw/6euOCDW6ZdgXMEn7vTj7qf7m1L++vl2elRBYZNqab887h+90k46tLsOGqvLSm+7vL7/e9i1aN+KIu0La9lp4M9lCJDhIJbLQXuQ7fgtnv4ytmnZ8fDAU78wHu59Ge/YsUr7cw+8mScczQ1Rbnm+1/3FQI3Uz/snXMHf8VftM0/Bw+D5CxvLVwusr2ZfSo89b/1bkV5WrcxPHQHLoBY8OHi9f5+o3/ekvry+PDP4N7ARRabVvtHKfU+If7p3/pH+vy1qwI9f6FI7kVYf/2Bf+S74ZjqtlFkGym41cmfb73aTwTeSD5/5il8/tOf8IEsfQJwWttE5h53GnM/+BHA5QzxMX/+/IL7WLlyZW5B80j/zTJ4OfuFqyAU9d/Ah+Jl7iLVcNyVfiiElrH+Cuq+Dt9ztGRBbr2v/qs6+79ghR+i46b/6L/sjD/6c8Ne/ou/88i2jsI/fjdY/1L59cdN98+rHhncfmaXvPNg/X3tDbjz3OID/uY7oUa3dxIZJAW3usu7FD/a4j9I8oNbpNl/m93WQSbzw1lL6pBNU9u2bVekkZjlnt8E/vSBfMV6p7fViInFz+FqHg0TdvePShjszxBLHQqN95auly//6s+hJhyBWR8qP7i1VGHgZZEKUHCrkNDWjuuYf8NiR/0PN4gMR6Eh0uNc6b//wZ4jt+Qm/xj0fobQEELfmpR7/+Y0vbfKdqAeA/Bud+K9nQA0d1fgTSF4cUDwxNhqffMXES9/yIoZgxzCYqja/6SB6wzW4ef3LxvsPWArpaXAfToLhTbIvQp1IBP3GriOSB0M6+BWqfu0JhOxbd/I5P38+WfpsLbjgTBmmn+evH/Fg1u171Er0nCmvQ2aA4fH3je/+vv8yC0FCivcK/T2s/34b/nOuNvffPxLLw68jUn75s6/72L4ysrcsilF7g1bbV9+ufDdFPLrgD8c/rU3yttupQ5Vi1TYsA1uLS0trFu3bugEmHBT7vln6dtVWajiN092zrFu3TpaWjRem0iO4BhptTisVosLgswKj/02eT9/JXk5V6c3FxgmqNq3xStXKDTwzxD8PUf1vieNbdie4zZ16lTa29tZu3btwJUHEOvtgnX+xtLRLfEBaheRP6ZalbW0tDB1aoGTsUUkpU7nQ9XqPKyKn0s3lM8fG8ptExmcYRvcotEoM2bMqMi2Xnz4NuJnXwTAzJO3cviAUvfNE5HaqHUP/FDp8R/IkA5lIsPLsD1UKiJSUt3CSq32O5j9lFN3CIc7BU/Zjii4iYiIiDQIBTcRkYJ0jtug6g7pXq2h3DaRwbEhc1VlFZnZWuCVAStuu4nAmzXYj5RPr8nQpNdl6NFrMvToNRmaavG67OqcK3i/u2ER3GrFzBY75+bUux2SpddkaNLrMvToNRl69JoMTfV+XXSoVERERKRBKLiJiIiINAgFt8q6ut4NkH70mgxNel2GHr0mQ49ek6Gprq+LznETERERaRDqcRMRERFpEApug2Rm88zseTNbbmYXFljebGa/SS1/zMym176Vw08Zr8v5Zvasmf3TzO43s13r0c7hZKDXJFDvRDNzZqar52qgnNfFzE5K/b0sNbObat3G4aaM969pZvZnM/tH6j3sA/Vo53BiZteZ2Roze6bIcjOzn6Zes3+a2cG1apuC2yCYWRi4Ejga2Bc4xcz2zat2JrDBObcH8CPgu7Vt5fBT5uvyD2COc+4A4Fbge7Vt5fBS5muCmY0CzgUeq20Lh6dyXhcz2xO4CHiHc24/4As1b+gwUubfyn8CtzjnDgJOBn5e21YOS9cD80osPxrYM/U4C/jvGrQJUHAbrEOB5c65l5xzfcDNwHF5dY4DbkhN3wq812xIDym+PRjwdXHO/dk515WafRSYWuM2Djfl/K0AfAv/5aanlo0bxsp5XT4FXOmc2wDgnFtT4zYON+W8Jg4YnZoeA/yrhu0blpxzDwHrS1Q5DrjReY8CY81sx1q0TcFtcHYGVgfm21NlBes45+LAJmBCTVo3fJXzugSdCfyxqi2SAV+T1KGFXZxzf6hlw4a5cv5W9gL2MrO/mdmjZlaq10G2XTmvyXzgo2bWDtwF/L/aNE1KGOznTsVEarETkaHCzD4KzAGOqHdbhjMzCwE/BD5e56ZIfxH84Z+5+J7ph8xsf+fcxrq2ang7BbjeOXe5mb0d+LWZzXLOJevdMKk99bgNzqvALoH5qamygnXMLILv1l5Xk9YNX+W8LpjZ+4CvAcc653pr1LbhaqDXZBQwC3jQzFYCbwMW6gKFqivnb6UdWOiciznnXgZewAc5qY5yXpMzgVsAnHOPAC34+2VK/ZT1uVMNCm6D8wSwp5nNMLMm/EmiC/PqLAROT01/CHjAabC8ahvwdTGzg4Cr8KFN5+xUX8nXxDm3yTk30Tk33Tk3HX/e4bHOucX1ae6wUc572G343jbMbCL+0OlLtWzkMFPOa7IKeC+Amc3EB7e1NW2l5FsIfCx1denbgE3OuddqsWMdKh0E51zczM4B7gHCwHXOuaVm9k1gsXNuIXAtvht7Of7ExpPr1+LhoczX5fvASOC3qWtFVjnnjq1bo7dzZb4mUmNlvi73AO83s2eBBHCBc05HDaqkzNfki8Avzew8/IUKH1eHQHWZ2QL8F5iJqXMLLwaiAM65X+DPNfwAsBzoAs6oWdv02ouIiIg0Bh0qFREREWkQCm4iIiIiDULBTURERKRBKLiJiIiINAgFNxEREZEGoeAmIrIVzGysmX2u3u0QkeFFwU1EZOuMBRTcRKSmFNxERLbOZcDuZvaUmX2/3o0RkeFBA/CKiGwFM5sO3Omcm1XnpojIMKIeNxEREZEGoeAmIiIi0iAU3EREts4WYFS9GyEiw4uCm4jIVnDOrQP+ZmbP6OIEEakVXZwgIiIi0iDU4yYiIiLSIBTcRERERBqEgpuIiIhIg1BwExEREWkQCm4iIiIiDULBTURERKRBKLiJiIiINAgFNxEREZEG8f8BHVbNEIGorL8AAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 720x216 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "def ABM(x0,m=4):\n",
        "    xsABM = [x0]\n",
        "    x = x0\n",
        "    fs = [ 0.0 for i in range(m)]\n",
        "    for i,t in enumerate(tr):\n",
        "        if len(xsABM) < m : #RK4\n",
        "            f1 = f(x,t)\n",
        "            f2 = f(x+0.5*h*f1, t+0.5*h)\n",
        "            f3 = f(x+0.5*h*f2, t+0.5*h)\n",
        "            f4 = f(x+h*f3, t+h)\n",
        "            x += h*(f1 + 2*f2+ 2*f3 + f4)/6.0\n",
        "            fs[0] = f1; fs[1] = f2; fs[2] = f3; fs[3] = f4\n",
        "        else:\n",
        "            f1,f2,f3 = fs[1:]; f4 = f(x,t)\n",
        "            ## 予測子\n",
        "            xp = x + h *(-9*f1+37*f2-59*f3+55*f4) / 24.0            \n",
        "            f5 = f(xp,t+h)\n",
        "            ## 修正子\n",
        "            x += h *(f2 -5*f3 +19*f4 +9*f5 )/24.0 \n",
        "            fs[0] = f1; fs[1] = f2; fs[2] = f3; fs[3] = f4\n",
        "        xsABM += [x]\n",
        "    return xsABM\n",
        "\n",
        "xsABM = ABM(x0)\n",
        "fig = plt.figure(figsize=(10,3))\n",
        "plt.xlabel(\"t\");plt.ylabel(\"abs. diff.\");plt.yscale(\"log\")\n",
        "plt.plot(tr,abs(xe-np.array(xs[:-1])),label=\"Euler\")\n",
        "plt.plot(tr,abs(xe-np.array(xsRK4[:-1])),label=\"RK4\")\n",
        "plt.plot(tr,abs(xe-np.array(xsAB4[:-1])),label=\"AB4\")\n",
        "plt.plot(tr,abs(xe-xsABM[:-1]),label=\"ABM4\")\n",
        "plt.legend();plt.show();plt.close()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6DUB18TOpj-9"
      },
      "source": [
        "1桁程度、AB法よりも精度が改善していることがわかる。"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "include_colab_link": true,
      "name": "Python_misc_ODE.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}