{"nbformat":4,"nbformat_minor":0,"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.8.1"},"colab":{"name":"06-ode_04_odes_in_various_areas.ipynb","provenance":[],"collapsed_sections":[],"toc_visible":true}},"cells":[{"cell_type":"markdown","metadata":{"id":"_hNhUJmaJJqJ"},"source":["# 経済や生物で使う微分方程式"]},{"cell_type":"markdown","metadata":{"id":"A2uQh0WuJJqK"},"source":["## マルサスの『人口論』と微分方程式\n","あなたが高校生以上なら社会の授業で次のような言葉を聞いたことがあるかもしれません。\n","\n","> 人口は制限されなければ幾何級数的に増加するが、\n","> 生活資源は算術級数的にしか増加しないので生活資源は必ず不足する。\n","\n","これは有名なマルサスの『人口論』での議論です。\n","幾何級数はいわゆる等比数列の和で、\n","算術級数は等差数列の和です。\n","マルサスは経済学者なので**経済の話**です。\n","経済の話でこういう数学ネタが出てくるわけです。\n","\n","ちなみに大学の経済学だと経済の話で微分積分が本当に出てきます。\n","ときどき経済学部の入試で数学が受験科目に入っていることがあります。\n","実際に大学に入ってから使うからです。\n","それも高校の理系のレベルを越えた数学です。\n","統計データをいろいろいじらないといけないので**統計学**が必要で、\n","そちらでも割といろいろ数学が必要です。\n","\n","統計学については「[プログラミングで数学を 中高数学虎の穴](https://phasetr.com/lp/mpgh1/)」で入門的な解説をしています。\n","興味があればぜひ受講してみてください。"]},{"cell_type":"markdown","metadata":{"id":"obnkmPOm1i_g"},"source":["## 微分方程式を立てる\n","話を元に戻しましょう。\n","マルサスの人口論に合わせて人口の増え方を考えます。\n","ここではもっと一般的に生物の集団だと思いましょう。\n","まずは生物の種類が 1 種類だとして考えます。\n","生物っぽく人の数というよりも一般に個体数と呼ぶことにして、\n","時刻 $t$ での個体数を $x (t)$ とします。\n","\n","個体数の増加率を「単位時間あたりの個体増加量」と「個体数の比」として定義しましょう：数学的には $x'(t) / x (t)$ です。\n","増加率が一定値 $\\alpha$ だとすると「単位時間あたりの個体増加量」が $x'(t)$、\n","「個体数」が $x(t)$ であり、\n","これが $\\alpha$ に等しいことは $x'(t) / x (t) = \\alpha$ と書けます。\n","この分母を払うと次の微分方程式が出てきます。\n","\n","\\begin{align}\n"," \\frac{d x (t)}{dt}\n"," =\n"," \\alpha x (t)。\n","\\end{align}\n","\n","この微分方程式を**マルサスモデル**と言います。\n","この方程式の解は $x (t) = x (0) e^{\\alpha t}$ です：あなたが指数関数の導関数をご存知なら実際に解になっていることはすぐわかるでしょう。\n","\n","何はともあれ、これで指数関数が出てきました。\n","等比数列は一般項が $a_n = 2^{n}$ のように指数関数で書けている数列だから、\n","これを足し上げていけばまさに幾何級数です。\n","放射性物質の崩壊とは指数の肩の符号が変わっているだけで、\n","それ以外は同じ式です。\n","放射性物質の崩壊と生物の個体数の変化を追う方程式が同じ形をしているわけです。"]},{"cell_type":"markdown","metadata":{"id":"BDipZs2aJJqL"},"source":["### 注意\n","先程こう書きました。\n","\n","> 時刻 $t$ での個体数を $x (t)$ とします。\n","\n","本当はここはかなり微妙です。\n","本来 $x (t)$ は整数です。\n","しかし実はここではいったん実数だと思って処理しています。\n","このあたりの処理は物理や化学でもよく出てきます。\n","ここではとてもやりきれませんが、この近似の合理化もきちんとやっておかないといけないことでしょう。\n","本当は時間についてもいろいろ議論があります。"]},{"cell_type":"markdown","metadata":{"id":"djCiB3CpJJqL"},"source":["## シミュレーション\n","シミュレーションの結果は放射性物質の崩壊のときと基本的に同じです。\n","いちおうきちんとやっておきましょうか。\n","まず微分方程式を次のように近似します。\n","\n","\\begin{align}\n"," \\frac{x_{n+1} - x_{n}}{h}\n"," =\n"," \\alpha x_n、 \\quad\n"," x_{n+1}\n"," =\n"," \\alpha x_n + h x_n. \\tag{1}\n","\\end{align}\n","\n","この漸化式を解いてみましょう。"]},{"cell_type":"code","metadata":{"scrolled":true,"id":"sU7K9ZzYJJqM","outputId":"5d3e571a-d643-4ea9-dc93-08746c1df228"},"source":["import numpy as np\n","import matplotlib.pyplot as plt\n","\n","def malthus_euler(nt, init = 10):\n","    dt = 2 / (nt - 1)\n","    # 初期条件設定\n","    x = np.zeros(nt)\n","    x[0] = init\n","\n","    for i in range(1, nt):\n","        x[i] = x[i-1] + c * dt * x[i-1] # ここで方程式を解いている\n","\n","    return x\n","\n","# パラメータ設定\n","c = 1\n","init = 1\n","nt = 101\n","\n","# 近次解\n","x_approx = malthus_euler(nt, init)\n","plt.plot(np.linspace(0, 2, nt), x_approx)\n","\n","# 厳密解\n","t = np.linspace(0, 2, nt)\n","x_exact = init * np.exp(c * t)\n","plt.plot(t, x_exact)\n","\n","# 凡例\n","plt.legend(['approximation', 'exact'])\n","# 描画\n","plt.show()"],"execution_count":null,"outputs":[{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deVyU5f7/8dfFJgiIoLgiIikuILtr7tsxU1NL09xNTa1OttppOWnLqX5p27GTWWmnNLNUzErNzH0XDPddUXEHlX2duX5/QH49hjIoM/cAn+fj4UNg7pl5M9y+vbnmuq9baa0RQghhvxyMDiCEEOL2pKiFEMLOSVELIYSdk6IWQgg7J0UthBB2zskaD1q9enUdEBBgjYcWQohyKS4uLklr7VvUbVYp6oCAAGJjY63x0EIIUS4ppU7d6jYZ+hBCCDsnRS2EEHZOiloIIeycVcaoi5KXl0diYiLZ2dm2ekpxh1xdXfHz88PZ2dnoKEIIbFjUiYmJeHp6EhAQgFLKVk8rSkhrTXJyMomJiTRo0MDoOEIIbDj0kZ2dTbVq1aSk7ZxSimrVqslvPkLYEZuOUUtJlw3ycxLCvsibiUIIURpOboBtn4LZXOoPLUVtY8uWLeOdd94plcf68MMPyczMvP55r169uHbtWqk8thCiBHLSYOnjsONzyC/9YUMp6ruQn59f4vv07duXF198sVSe/+aiXr58OVWrVi2VxxZClMCqVyHlDPT7FFwql/rDV6ii7tevH1FRUQQHBzN79mwAPDw8ePbZZ4mMjKRr165cvnwZgE6dOjF58mTatm1LSEgIO3bsAGDq1KmMHz+eHj16MGLECLKzsxk9ejTNmzcnIiKCtWvXAvD+++8zZswYAPbu3UtISAiZmZl89dVXPPHEEwCMGjWKiRMn0rlzZwIDA1m/fj1jxoyhadOmjBo16nruiRMnEh0dTXBwMK+99hoAH3/8MefOnaNz58507twZKDh1Pykp6frzh4SEEBISwocffghAQkICTZs2Zdy4cQQHB9OjRw+ysrKs+ZILUf4d+x3i5kLbJ8C/lVWewmbT82407af9HDiXWqqP2axOFV7rE3zbbebMmYOPjw9ZWVm0aNGCBx98kIyMDCIjI5kxYwavv/4606ZNY+bMmQBkZGSwZcsWNmzYwJgxY9i3bx8AcXFxbNq0CTc3N2bMmAEUlPGhQ4fo0aMHR44cYfLkyXTq1ImYmBjeeustPvvsMypX/uv/tFevXmXNmjUsW7aMPn36sHnzZr744gtatGhBfHw84eHhvPXWW/j4+GAymejatSt79uzh73//O++//z5r166levXq//OYcXFxzJ07l+3bt6O1plWrVnTs2BFvb2+OHj3KggUL+Pzzzxk0aBCLFy9m2LBhpfEjEKLiyU6BZU9C9SDS204hPzOXqpVdSv1pKtQR9ccff0xYWBitW7fmzJkzHD16FAcHBx5++GEAhg0bxqZNm65vP2TIEAA6dOhAamrq9fHfvn374ubmBsCmTZsYPnw4AE2aNKF+/focOXIEBwcHvvrqK4YPH07Hjh259957i8zUp08flFI0b96cmjVr0rx5cxwcHAgODiYhIQGA77//nsjISCIiIti/fz8HDhy47fe5adMm+vfvj7u7Ox4eHgwYMICNGzcC0KBBA8LDwwGIioq6/hxCiDuw8iVIOw/9PuWNlSfp9dFGMnNLPiRaHEOOqIs78rWGdevWsXr1arZu3UrlypXp1KlTkXOFb5yadvM0tT8/d3d3v/61210c+OjRo3h4eHDu3LlbblOpUiUAHBwcrn/85+f5+fmcPHmS6dOns3PnTry9vRk1alSxc5xvl+nG53B0dJShDyHu1OEVED8P2j3DqhQ/FsbGManTPVR2Kf1arTBH1CkpKXh7e1O5cmUOHTrEtm3bADCbzSxatAiAb7/9lnbt2l2/z8KFC4GCI1QvLy+8vLz+8rgdOnRg/vz5ABw5coTTp0/TuHFjUlJSeOqpp9iwYQPJycnXn6OkUlNTcXd3x8vLi4sXL7JixYrrt3l6epKWllZkpqVLl5KZmUlGRgYxMTG0b9/+jp5fCFGEjKSCIY+azUlq8Qz/WLKXZrWrMLlbkFWezpAjaiP07NmTWbNmERoaSuPGjWndujVQcHS8f/9+oqKi8PLyul7OAN7e3rRt25bU1FTmzJlT5ONOmjSJCRMm0Lx5c5ycnPjqq6+oVKkSEydOZNKkSQQFBfHll1/SuXNnOnToUOLcYWFhREREEBwcTGBg4P8MoYwfP5777ruP2rVrX38TEyAyMpJRo0bRsmVLAMaOHUtERIQMcwhRGrSGn5+GrGvo4TH848fDpOXks2BwOC5O1jn2Vbf7NflORUdH65svHHDw4EGaNm1a6s91tzw8PEhPT//L1zt16sT06dOJjo42IJXx7PXnJYTh9nwPS8ZBt6l8X+khXli8h1fub8rY9oF39bBKqTitdZGFU2z9K6UaK6Xib/iTqpSafFeJhBCiLLp2Bn55Duq14lTjR5n2037aBFZjzL3WXcCs2KEPrfVhIBxAKeUInAVirJrKhoo6moaCNx+FEOI6sxmWTgRtIv+BWUz+fi+ODooZg8JwcLDu+jglHaPuChzXWt/y2l5CCFEubZ0JCRvhgU+YGZ/PH6ev8e8hEdSp6mb1py7pyPdgYEFRNyilxiulYpVSsX+e3SeEEOXChX2w5g1o0ptdPr3495pjDIioS5+wOjZ5eouLWinlAvQFfijqdq31bK11tNY62te3yCueCyFE2ZOXVfDmoZs36T1m8NTCeGpVcWXqA7Y7H6QkQx/3Abu01hetFUYIIezOb6/BpQMwbDGv/naBs1ez+GFCG6q42u5SdSUZ+hjCLYY9Krp169axZcsWo2MIIUrbkVWw4zNoPYmYtCbE/HGWp7oGEVXfx6YxLCpqpVRloDuwxLpxyiYpaiHKofRL8OMkqBnCqYjneCVmHy0DfHiiS0ObR7GoqLXWmVrralrrFGsHsqZ58+bRsmVLwsPDeeyxxzh16hSNGjUiKSkJs9lM+/btWbVqFVD0kqgAK1euJDIykrCwMLp27UpCQgKzZs3igw8+IDw8/PriR0KIMsxshqWTICeNvH6z+fuiQzg4KD4YHI6jlafiFcWYU8hXvAgX9pbuY9ZqDvfd+sopBw8eZOHChWzevBlnZ2cmTZrE+vXrmTJlChMmTKBVq1Y0a9aMHj16AEUviWo2mxk3bhwbNmygQYMGXLlyBR8fHyZMmICHhwfPPfdc6X5PQghjbP8Ujv0GvaYz/Q8Hdp+5xn+GRlLXBlPxilJh1vr4/fffiYuLo0WLFgBkZWVRo0YNpk6dyg8//MCsWbOIj4+/vv3HH39MTEzBeT1/Lol6+fJlOnToQIMGBWch+fjYdpxKCGED5/4oeAOxSW/WVenLZ0tiGdrKn17NaxsWyZiivs2Rr7VorRk5ciRvv/32/3w9MzOTxMREoOAsRU9Pz1suiaq1lit0C1Ge5aTBojHgUYPLnafz7Ow9NKnlyau9mxkaq8Isc9q1a1cWLVrEpUuXALhy5QqnTp1iypQpDB06lNdff51x48YBt14StU2bNqxfv56TJ09efwy49XKjQogyZvnzcDUBU//ZPLXsFJm5JmY+EoGrs6OhsSpMUTdr1ow333yTHj16EBoaSvfu3UlISGDnzp3Xy9rFxYW5c+fSs2dP8vPzCQ0N5dVXX72+JKqvry+zZ89mwIABhIWFXb8yTJ8+fYiJiZE3E4Uoy+K/hd0LoMML/Pt4DbYcT2Za32Aa1vA0OpkscyqKJj8vUaFcPgyzO0HdKLa0/ZKhc2PpH1GXGQPDbDbceVfLnAohRLmWlwU/jAJnN5L+NpO/f7+Xe3w9eLNfiN28J1VhZn0IIUSRVkyBSwcwPbKIv/98gfScPOaPbWWVax/eKZseUVtjmEWUPvk5iQpj90LY9V+4dzIfJviz5XgybzwQQuNaxo9L38hmRe3q6kpycrKUgJ3TWpOcnIyrq6vRUYSwrkuH4OfJ4N+WdX6P8e81xxgU7cfA6HpGJ/sLmx3b+/n5kZiYiKxVbf9cXV3x8/MzOoYQ1pObAT+MBOfKnO/xHybP2UeTWp68/kCI0cmKZLOidnZ2vn5GnxBCGEZr+OVZuHyY3EcWM+HHc+SbNJ8OizJ8vvStyKwPIUTFEje3YL50pxeZtr8Gu89cY/rAUBpUdzc62S1JUQshKo6zcQWzPBp2Y5HHI8zffprHOgbSM8S4dTwsIUUthKgYMq/A9yPBoyaH2s7g5aX7aRNYjed7NDY6WbHsZ6KgEEJYi9kEi8dC+kXSHvmZ8YtO4l3ZhY+HRODkaP/Hq1LUQojyb93bcPx3zPd/wOPrFedTslj4WBt8PSsZncwi9v9fiRBC3I1Dv8CG9yBiGDOS27DhyGVefyCESH9vo5NZTIpaCFF+JR2FJY9BnQh+rf88n6w7wZCW/gxp6W90shKRohZClE85abBwGDi5cLzzLJ5ZcogI/6pM7WvsRQDuhIxRCyHKH7MZYiZA0hHSB/3A6JjzuFdy4tOhUVRyss+TWm5HjqiFEOXPxulw6GdM3d5gwmZPLqRkM2t4FLW8yuYaNhYVtVKqqlJqkVLqkFLqoFKqjbWDCSHEHTm0HNa+BaGD+deVzmw6lsSb/cvWm4c3s/SI+iNgpda6CRAGHLReJCGEuEOXD8OS8VA7nMV1n+PLzQmMahvAIDtcEa8kih2jVkpVAToAowC01rlArnVjCSFECWVegW8fBmc3drf7hBe/PUq7htV55f6yf0k5S46oA4HLwFyl1B9KqS+UUn9ZvUQpNV4pFauUipWlTIUQNmXKL7icVupZLt3/JWOWnMfPuzKfPBJZJs48LI4l34ETEAl8qrWOADKAF2/eSGs9W2sdrbWO9vX1LeWYQghxG6tehpPryek5nRGrINdk5vMR0XhVdjY6WamwpKgTgUSt9fbCzxdRUNxCCGG82DmwfRbmVpN4/EAzjlxMY+YjkTSs4WF0slJTbFFrrS8AZ5RSfy4x1RU4YNVUQghhiRPrYfnz0LA775iGsPrgJab2DaZjUPn6rd7SE16eBOYrpVyAE8Bo60USQggLJB+H70dAtYYsajCN2T8VzPAY0SbA6GSlzqKi1lrHA9FWziKEEJbJuloww8PBkZ1tPuXFRafoGORbLmZ4FEVOIRdClC35ubBwOFxN4Eyf7xiz9DL3+How85Gysbb0nZCiFkKUHVrDz09DwkZS7vuEwb864OoCc0a3wNO1fMzwKEr5/O9HCFE+bXof4ueR1+55RuwM4EpGLnNGtqBuVTejk1mVFLUQomzYtxh+fx1zyEAmJv6NPWdT+GhwOM39vIxOZnVS1EII+3dqC8RMQPu34Q2Hiaw+dImpfYLpEVzL6GQ2IUUthLBvSUdhwRCoWp+v/P/F3B0XGN8hkJFtA4xOZjNS1EII+5V+CeY9CA5OrIqYybTV5+kdWpsXezYxOplNyawPIYR9ykmH+QMh/RLxXefz+M9XaNXAh+kDw3BwUEansykpaiGE/THlFayGd2EPp7p/wdAVedzj68HsEdG4Ope9S2ndLSlqIYR90Rp+mgzHfiO583s8uMaLqpUd+Gp0S7zcyu9c6duRohZC2Jc1b0L8PDJaP8uAHY3IN+fx3ZhWZfZ6h6VB3kwUQtiP7Z/Bxunkho9g4KFOXErNYc6oFjSs4Wl0MkNJUQsh7MO+xbBiCqbGvRl58WGOXErn02GRZfqitKVFiloIYbzja2DJY2j/NjyR8zhbT6YwfWAYnRrXMDqZXZCiFkIYKzEWvhuG9g3iFbeXWXHoKtP6BtMvoq7RyeyGFLUQwjiXDsK8B9EeNZhR8x3m707hme5BFeqsQ0tIUQshjHE1Ab7pD85uzAn8gJk70ni0XQOe7NLQ6GR2R6bnCSFsL/U8fP0A5GXxQ+hs3tiQyaBoP17u1RSlKtZZh5aQI2ohhG1lJMM3/SAjieXhn/D8hnx6h9bm7QGhFe7UcEvJEbUQwnayU2DeALiawPoWnzJpraJrkxp88HA4jlLStyRFLYSwjT8XWbq4j22tZjJ6XSXubViNT4ZG4lxOr3VYWqSohRDWl5cFCwZD4k7+aPU+Q9d7ER3gzRcjWlTIRZZKyqKiVkolAGmACcjXWkdbM5QQohzJz4GFwyBhE/tavcegjTUJ8/NizqgWuLlISVuiJEfUnbXWSVZLIoQof/JzC5YrPbaaQy3+xYBNfjSt7clXY1riUUl+obeUDAwJIazDlAeLRsPh5RyJnkrfLYEE1fLgmzGtqOJaMZcrvVOWFrUGViml4pRS44vaQCk1XikVq5SKvXz5cuklFEKUPaZ8WPwoHPqZY1Gv0ntbExrW8GDeo63wqiwlXVKWFvW9WutI4D7gcaVUh5s30FrP1lpHa62jfX19SzWkEKIMMeXDknFw4EeORfyDXtuDucfXg/ljW1G1sovR6coki4paa32u8O9LQAzQ0pqhhBBllCkfloyF/Us4FvYCvXaE0aiGB9+ObYW3u5T0nSq2qJVS7kopzz8/BnoA+6wdTAhRxlwv6RiOhE6hV2wkjWt5Ml9K+q5Z8rZrTSCm8Px7J+BbrfVKq6YSQpQtpryCMekDP3Kw+Qv0iQ0nuI4nX4+RMenSUGxRa61PAGE2yCKEKIvyc+CH0XD4F/aGTKFfXDjh9aoyd3QLmd1RSmQioxDizuVlw/cj4Oiv/BHyMg/GBRMd4M2cUS1knnQpkldSCHFncjPhu0fgxFq2NXuVwbFNadewOrNHRFHZRaqlNMmrKYQouexU+PZhOLONtU2mMnpXEN2a1mDmI5GydocVSFELIUom8wrMfwh9fje/NHqDJ+Ib0Du0Nh88HC6r4FmJFLUQwnJpF+Gb/ujkoywIeJOX9vgzKNqPtweEynrSViRFLYSwzLXT8PUD6LQLzKrzL949UIdx7Rvwklw+y+qkqIUQxbt8BL7ph85N5+3q7zD7aHWe6xHE450bSknbgBS1EOL2zu6C+Q9hxoEp7m+zKKEqrz8QzIg2AUYnqzCkqIUQt3ZiPXz3CCZXb8bpV9h4oQofDw6nT1gdo5NVKFLUQoiiHVgGix8l16sBgzKf50imJ3NGRdG+kayOaWtS1EKIv9r5JSx/jnTfcHonPUma8uS78S0I9atqdLIKSYpaCPF/tIZ1b8P6d0mq05nuZ0bj6VmFxWNaElDd3eh0FZYUtRCigCkffnkGdv2X4379ue/EgwTVrsrcUS3x9axkdLoKTYpaCAG5GQUr4B39lW11RzH4WHc6BtXgP0MjcZfFlQwnPwEhKrr0S/DtIPT53Syq9QzPH4/m4eh6vNk/RE4JtxNS1EJUZJePFKzbkX6J97z/yX8SgnimexBPdpETWeyJFLUQFVXC5oI50sqJJ11eZ/XFenz4cCj9IuoanUzcRIpaiIpozw/w4ySyPOoxKO0ZTusafPNoFK0CqxmdTBRBilqIikRrWP8urHubpGot6HnhMTyqVmfJqBbc4+thdDpxC1LUQlQUedmw7AnY+wP7fHvT/8wgIhvUYNawKLlKuJ2TohaiIki/DAuHwpntxFQby9NnOjMwqh5v9W+Oi5PM7LB3UtRClHcX9sGCwZgzknjH40U+PxfKy72aMrZ9A5nZUUZYXNRKKUcgFjirte5tvUhCiFJz6BdYPI5cZ0/G6GnEpwUwZ2QEnZvUMDqZKIGSHFE/BRwEqlgpixCitGgNG6fDmrdI9gqmT9IkXKrWIWZENI1qehqdTpSQRYNTSik/4H7gC+vGEULctdwM+GEUrHmT3d7daXvxOe4JbMiPj7eTki6jLD2i/hB4AbjlT1kpNR4YD+Dv73/3yYQQJXf1FCwcir6wj3mej/Lq+S6Max/IlJ5NcJLTwcusYotaKdUbuKS1jlNKdbrVdlrr2cBsgOjoaF1qCYUQljmxDn4YjcmUz3NOL7PiWggfDQ7lgXA507Css+SI+l6gr1KqF+AKVFFKzdNaD7NuNCGERbSGrZ+gf3uVFPcGDEx9gqwqASweE0VwHS+j04lSUGxRa63/AfwDoPCI+jkpaSHsRE46LHsS9i9hn2cHBl8eSVSQPx89HC4nsZQjMo9aiLIq6RgsHIZOOsxXbqOYdrk7f+/SiKe6BeHoIPOjy5MSFbXWeh2wzipJhBCWO7AMfnycXBx5XL/E9qxQvhwZTtemNY1OJqxAjqiFKEtMebB6KmydyTn3ZjyUPAGfuoH8MjSKej6VjU4nrESKWoiyIvUcLBoDp7eywq0PTyU/xIMt7+G1Ps1wdXY0Op2wIilqIcqCY6thyXhMuVm8op5iWUZb3hvcXKbeVRBS1ELYM1M+rHsbvXEGl90CGZLxIs41m7BsaKSsH12BSFELYa9SzsLisXB6C7+79uCJq0N4sFUjXu0tQx0VjRS1EPbo8EpYOpH8vGxe4Ul+yWrPjEdCuT+0ttHJhAGkqIWwJ3nZsPo12D6LRNdGDM+YQNV6TVk+OEJmdVRgUtRC2IvLh2HRo3BxL4ude/NSykOM69SMp7o1wlkWVKrQpKiFMJrWEDcXvfIlslUlnsx7nn0ubZg7Noy291Q3Op2wA1LUQhgpI6lgrY7Dy9njEsnY1EdpGdqMlf1CqFpZ1uoQBaSohTDK0d/QPz6OOfMK0/VI5uXcx9SBzRkQWVeuZSj+hxS1ELaWmwGrXoXYLznvEsCYrDeoUj+c5YPC5A1DUSQpaiFsKTEWYh5DJx9nnurDuxkDebJnCGPbB8qKd+KWpKiFsIX8XFj/DnrTB1x18mVS7suk1mzNoofDaFJLrhctbk+KWghru7AXYibAxX385NCFVzOGMrJTc57o0ggXJ5l2J4onRS2EtZjyYOMM9Ib3SHOowtO5z5Lo24lvxoQS6lfV6HSiDJGiFsIaLuyFpRPhwl5Wqg68nDmcRzqG8Z+uDankJOt0iJKRohaiNOXnwIb30Js+IE158lzu05yu0YWvB4YRUlcuNCvujBS1EKUlMRb94+Ooy4f4mY5MyxvGiK4RzOx4j4xFi7siRS3E3cpJhzVvorfP4opjdZ7NfYFUv04seDCURjU9jU4nygEpaiHuxtHV6J8no1LO8K25Bx+Zh/B470iGt66Pg8yLFqVEilqIO5F2EVa+CPuXcMbBj6dzXsOnaQeW9g2mTlU3o9OJcqbYolZKuQIbgEqF2y/SWr9m7WBC2CWzGXb9F/3ba5hyM/k4/yGWuD3EK0Mj6BlSy+h0opyy5Ig6B+iitU5XSjkDm5RSK7TW26ycTQj7cmEv+uenUYk7iVPBTMkZTcc297KyRxAeleSXU2E9xe5dWmsNpBd+6lz4R1szlBB2JTsV1r+L3vYpacqDqbkTOF6nDx/1by5T7oRNWHQYoJRyBOKAhsAnWuvtRWwzHhgP4O/vX5oZhTCG1rBvMeZfX0KlX2KhqTP/cRrGxH4tmB5dT94sFDZjUVFrrU1AuFKqKhCjlArRWu+7aZvZwGyA6OhoOeIWZdulg+jlL6ASNnCYQF7KnUaT6M78+LcmeLvLgv7Ctko0sKa1vqaUWgf0BPYVs7kQZU92Cqx7F719Fum48W7eaA7UHsC0frI+hzCOJbM+fIG8wpJ2A7oB71o9mRC2ZDZD/DzMq6dBZjLf5Xfmy0rDeax/C16P9JNhDmEoS46oawP/LRyndgC+11r/bN1YQtjQ6e2Yl7+Aw4V4dusgpuU/Tat7u7K0S0M8XZ2NTieERbM+9gARNsgihG1dO43+7TXU/iUk48MbuY+T1bg/H9zfjAbV3Y1OJ8R1MvlTVDw5abDpA8xb/k2eSTErfwBrqg1hSp9I2jasbnQ6If5CilpUHKZ8+ONrTL+/hWNWEstMbfmi0kiG39+WJVH15JqFwm5JUYvyT2s48iv5q/6JU/Jhdpmb8P/007Tp0IOFHQJxl7MKhZ2TPVSUb2fjMP36TxxPbyJR1+bd/MlUjXyQT7oHUaOKq9HphLCIFLUon5KOYf79dRwO/kgKVXg/bzTJQYN59r5gGtaQNaJF2SJFLcqX1HOY1/8/2PU12dqZ2fkD2FVnGE/dH0lUfW+j0wlxR6SoRfmQeQW98X3M22djNpuYn9+VVdWGM65Xa54K8kUpeaNQlF1S1KJsy05Bb5mJacsnOORnstTUjsVVhjOkR3vmNa8tZxSKckGKWpRNOWno7bPJ3/QRzrkprDK1ZIH7UB7o3o2vw+vg5CgXkxXlhxS1KFty0tE7Pi8o6JyrbDBF8I3bMP7WrQdfRvrJ1b5FuSRFLcqGnHT0zi/I2/gRLjlX2GwK4xvXF+nU8z4+i/ajkpOj0QmFsBopamHfslMx7/ic/E0f45J7jW2m5sx3e54OXe/nP1FS0KJikKIW9inzCuZts8jf+ikuealsNoXxvfs/6NS1F/+OkCEOUbFIUQv7knaR/C0z0Tu/xDk/g3WmKJZWeYTu3Xry71B5k1BUTFLUwj5cOUnuxo9w2D0fZc5nhaklv1UbRq9u3ZnZrKZMsxMVmhS1MNb5PWSvm4HL4WUoFN/ndyDObwT9u3bgo4bV5EQVIZCiFkbQGk6sJWPtB7gnbiBfu/FfUy9ONRrJkK6teMTPy+iEQtgVKWphO/m56P1LyFj3ER5XD5Chq/KZHkx22CiGdQrDv1ploxMKYZekqIX1ZV0lb8dc8rZ8SuWcS5w312WB0ySqtx3G6DaN8HZ3MTqhEHZNilpYT9JRMjfOxHnvdzibs9lmCuHXKhMJ7fQgUyJkDrQQlpKiFqXLbEYf/53UdTPxOrsOJ+1EjPleDvgPpUeXrrwRKG8QClFSUtSidGSnkrfrW7I2z6JKxklydFU+YSA54aN4sEMEg6rJVb2FuFPFFrVSqh7wNVALMAOztdYfWTuYKCMuHSRt46e47P+BSuZM9pvvYYX7M/i3H8qoqAC5HqEQpcCSf0X5wLNa611KKU8gTin1m9b6gJWzCXuVn4vpwE+kbJyFz+UduGhnfja34VjAYNp3+hsvyvCGEKWq2KLWWp8Hzhd+nKaUOgjUBaSoK5qrCaRt+RKH+Hm4510h3ezLt87DcYwawYB2oTwoF4sVwipK9HupUioAiAC2F3HbeGA8gL+/fylEE3bBlEf+weVc3QBHmfwAAAqTSURBVPQ51S5sorKGNeZI9tR6huYdBzChaW1Zf0MIK7O4qJVSHsBiYLLWOvXm27XWs4HZANHR0brUEgpjJB/n2uY5OO39Fo+8K+RpH750fAgdMYL77o2mu4+cnCKErVhU1EopZwpKer7Weol1IwnD5GaSsyeGlC1zqHElFg/twDodzoHazxDcfgCj5ehZCENYMutDAV8CB7XW71s/krAprTGf2cnljXPwOv4jruZMMs01+dx1OC5Rw+jZJpxuMvYshKEsOaK+FxgO7FVKxRd+7SWt9XLrxRJWl3KWq1u/QcfPxyf7NJ66EitpzcV7HiKy/f2MDfCRmRtC2AlLZn1sAuRfbHmQk05GfAypO76hZvIOvNHsMDdhd7VnqNXmYf4W3gg3FzmtWwh7I2cjlHemfHKPrObylm+onrgad51NkrkG81wfxiFiMF3btqall5vRKYUQtyFFXR5pTf7p7VzYNA+vEz/habqGu3bnF8cOpAUNILrDfQyv4yVDG0KUEVLU5YXWmC7s48Kmebgd/RGf3PP4amfWqyjO+/cm6N4BPNCoDo5ySSshyhwp6jLOfPEQ57YsoNLhpfhmJ1BTO7CV5iTUGkGdNoPoGNJAlhMVooyToi6DTJcOc27LAlwO/UjN7BPU0YqdNGGV71P4thpEu7AmtHeRH60Q5YX8ay4LtCbv/D7ObV1IpaO/UCv7BHW1YhdBrKn+BN4tBtIuIoRWslKdEOWS/Mu2V2Yz2adjObf1ezxOrqRG7hnqaUUcTdhQ40m8ox+ibXgI0VLOQpR78q/cnuTnknpoLRd3xlA98Te8TUnU047EqmDW1xpIjZYP0ap5E1o4y5izEBWJFLXBdNZVLu5aTtruZdS5vJEqOgNn7cIOx3CS6k+kbsv+RDdpQBtZY0OICkuK2gC5l46TuH0J6sgK6qXFUwsTTroKmyu1ITOwJ43a9KGDf02Z5yyEAKSobcOUx9WD67kYtwyvxLXUzjtNIHBU+7HS6yGcmvSieetu9PDxMDqpEMIOSVFbSd7VRM7s/Im8Q79S78o2vMmisnYi3iGY2Dr98A7vS2R4BI1kGp0QohjSEqUlP4dL+9dzOX45Vc6up17uCQKBC9qHLZU7kRfYhcCWvWnpX1uGNIQQJSJFfae0Jv3sARJjf8HhxBrqpe6iBjlU1Y7sdWzK3lqTqBLyN0Kj2tLNzcXotEKIMkyKugRyrp3ndOwKco6soVbSVqqbk2gCnNK12FKlJ/qeLjRo0ZPIOvJGoBCi9EhR30Z+xlXO/LGatIO/431pK/XyEmgEXNPuHKgUzi6/cfiG9SQ4OJT6TjJ9TghhHVLUNzBlp3Emfg0pB9fgeX4r9XOO0EBpsrQLB5ybcajORDyadaNZRDvausvlqYQQtlGhizov8xqn49eSdmgdnhe3Uz/nMAGYydOOHHRszLqaI3EN6kxQVBeivKsYHVcIUUFVqKLOunaJhPjfyTqykaqXd1I/9yj3KE2uduSIYyM21hhKpYYduSeyC6HVqxkdVwghgPJc1FqTnHiYs3vWkn9yK75Xd1HPdIamQI525rBTYzbWHolrw44ERnYixMfH6MRCCFGkclPUptxsTu/fRvKhjTif24lf2h6qcZVqQIp255hrMCfq9MEjqB0NwzsS6ilnAQohyoayWdRak3zuJIn7NpB7chteyfEE5B6lgcqnAZBIDY56RHKgbit8mnakYXAUUc7ORqcWQog7UmxRK6XmAL2BS1rrEOtH+qus9BRO7dtM6rFtuFzYRd30/fhyhWoUDGMcd27EzpoDcarfkrrNO1O3XgB+Mo9ZCFFOWHJE/RUwE/jaulEK5OZkc/rADq4c2446+wfVU/fhbzpNE6UBSFS1OOkZwZHa0XgHtaFBcGuaubnZIpoQQhii2KLWWm9QSgVYO0huTjYJ77UnIO8EDVU+AFfw5LRrU7b79sAtoAV+Ie3wq1kXP2uHEUIIO1JqY9RKqfHAeAB/f/8S39+lkiuplesT594Sl3rR1GrShjoBQfg4yBl/QoiKTWmti9+o4Ij6Z0vHqKOjo3VsbOzdJRNCiApEKRWntY4u6jY5XBVCCDsnRS2EEHau2KJWSi0AtgKNlVKJSqlHrR9LCCHEnyyZ9THEFkGEEEIUTYY+hBDCzklRCyGEnZOiFkIIOydFLYQQds6iE15K/KBKXQZO3eHdqwNJpRintEiukpFcJSO5SqY85qqvtfYt6garFPXdUErF3ursHCNJrpKRXCUjuUqmouWSoQ8hhLBzUtRCCGHn7LGoZxsd4BYkV8lIrpKRXCVToXLZ3Ri1EEKI/2WPR9RCCCFuIEUthBB2zmZFrZTqqZQ6rJQ6ppR6sYjblVLq48Lb9yilIi29r5VzDS3Ms0cptUUpFXbDbQlKqb1KqXilVKleKcGCXJ2UUimFzx2vlPqnpfe1cq7nb8i0TyllUkr5FN5mzddrjlLqklJq3y1uN2r/Ki6XUftXcbmM2r+Ky2XU/lVPKbVWKXVQKbVfKfVUEdtYbx/TWlv9D+AIHAcCARdgN9Dspm16ASsABbQGtlt6Xyvnagt4F35835+5Cj9PAKob9Hp1ouCqOyW+rzVz3bR9H2CNtV+vwsfuAEQC+25xu833Lwtz2Xz/sjCXzfcvS3IZuH/VBiILP/YEjtiyw2x1RN0SOKa1PqG1zgW+Ax64aZsHgK91gW1AVaVUbQvva7VcWustWuurhZ9uA5tcW/duvmdDX6+bDAEWlNJz35bWegNw5TabGLF/FZvLoP3LktfrVgx9vW5iy/3rvNZ6V+HHacBBoO5Nm1ltH7NVUdcFztzweSJ//SZvtY0l97Vmrhs9SsH/mH/SwCqlVJwquLhvabE0Vxul1G6l1AqlVHAJ72vNXCilKgM9gcU3fNlar5cljNi/SspW+5elbL1/WczI/UsVXEM2Ath+001W28dK7SrkxVBFfO3meYG32saS+94pix9bKdWZgn9I7W748r1a63NKqRrAb0qpQ4VHBLbItYuCtQHSlVK9gKVAIwvva81cf+oDbNZa33h0ZK3XyxJG7F8Ws/H+ZQkj9q+SMGT/Ukp5UPCfw2StderNNxdxl1LZx2x1RJ0I1Lvhcz/gnIXbWHJfa+ZCKRUKfAE8oLVO/vPrWutzhX9fAmIo+BXHJrm01qla6/TCj5cDzkqp6pbc15q5bjCYm34tteLrZQkj9i+LGLB/Fcug/askbL5/KaWcKSjp+VrrJUVsYr19zBoD70UMxDsBJ4AG/N9gevBN29zP/w7E77D0vlbO5Q8cA9re9HV3wPOGj7cAPW2Yqxb/d8JSS+B04Wtn6OtVuJ0XBeOM7rZ4vW54jgBu/eaYzfcvC3PZfP+yMJfN9y9Lchm1fxV+718DH95mG6vtY6X24lrwjfai4J3S48DLhV+bAEy44YX4pPD2vUD07e5rw1xfAFeB+MI/sYVfDyx8wXcD+w3I9UTh8+6m4E2otre7r61yFX4+CvjupvtZ+/VaAJwH8ig4gnnUTvav4nIZtX8Vl8uo/eu2uQzcv9pRMFyx54afVS9b7WNyCrkQQtg5OTNRCCHsnBS1EELYOSlqIYSwc1LUQghh56SohRDCzklRCyGEnZOiFkIIO/f/AazctfHFIubeAAAAAElFTkSuQmCC\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"ez5m1YiPJJqP"},"source":["## 放射性物質の崩壊の時との違い\n","\n","`t` が大きくなると厳密解の方が少し大きくなります。\n","指数の肩の符号が違うだけなので放射性物質の崩壊で形式的に時間を負の方に伸ばせば同じようにズレが出てくるはずです。\n","\n","何はともあれ区間の分割数 `nt` を大きくしてみたのが次の結果です。\n","具体的には 101 から 1001 にしました。"]},{"cell_type":"code","metadata":{"scrolled":true,"id":"MKp0Vd2UJJqQ","outputId":"07e28c01-c276-4b74-beaf-5ac8677dff94"},"source":["import numpy as np\n","import matplotlib.pyplot as plt\n","\n","def malthus_euler(nt, init = 10):\n","    dt = 2 / (nt - 1)\n","    # 初期条件設定\n","    x = np.zeros(nt)\n","    x[0] = init\n","\n","    for i in range(1, nt):\n","        x[i] = x[i-1] + c * dt * x[i-1] # ここで方程式を解いている\n","\n","    return x\n","\n","# パラメータ設定\n","c = 1\n","init = 1\n","nt = 1001 # 変更\n","\n","# 近次解\n","x_approx = malthus_euler(nt, init)\n","plt.plot(np.linspace(0, 2, nt), x_approx)\n","\n","# 厳密解\n","t = np.linspace(0, 2, nt)\n","x_exact = init * np.exp(c * t)\n","plt.plot(t, x_exact)\n","\n","# 凡例\n","plt.legend(['approximation', 'exact'])\n","# 描画\n","plt.show()"],"execution_count":null,"outputs":[{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deVxVdeL/8dcHRFFAxX1BBHNDQBFx3zAnU0tLs23SMivTalpmmqn5NU3Ld/rWTMs0bWOOqWOWmpbtOpa7WS6YO24pCuKCiqAIyPL5/QH5dQzlovdyLvB+Ph4+gnvPuefN5fT2+LnnfI6x1iIiIt7Lx+kAIiJyaSpqEREvp6IWEfFyKmoRES+nohYR8XLVPPGiDRo0sGFhYZ54aRGRSikhIeGYtbZhSc95pKjDwsJYv369J15aRKRSMsbsv9hzGvoQEfFyKmoRES+nohYR8XIeGaMuSV5eHikpKeTk5JTXJuUy+fv7ExISgp+fn9NRRIRyLOqUlBSCgoIICwvDGFNem5UystZy/PhxUlJSCA8PdzqOiFCOQx85OTnUr19fJe3ljDHUr19f//IR8SLlOkatkq4Y9HsS8S76MFFExA2ydi7hzIo3obDQ7a+toi5nn3/+OS+99JJbXuv111/nzJkz574fOnQoJ0+edMtri0gZ5J4id95Eji19mzPZWW5/eRX1FcjPzy/zOsOHD+fJJ590y/YvLOqvv/6aunXruuW1RcR1Bz96nLpnj7Cyw3PUCghy++tXqaK+8cYb6dKlC5GRkUyePBmAwMBAfve73xEbG8vAgQNJS0sDID4+nkcffZRevXoRFRXF2rVrAXj22WcZP348gwYN4s477yQnJ4e7776b6OhoOnfuzNKlSwF47bXXGDduHABbtmwhKiqKM2fOMH36dB566CEAxo4dy8SJExkwYACtWrVi+fLljBs3joiICMaOHXsu98SJE4mLiyMyMpJnnnkGgDfeeIPU1FQGDBjAgAEDgKJL948dO3Zu+1FRUURFRfH6668DkJSUREREBPfddx+RkZEMGjSI7OxsT77lIpXeqa0LaP7TbD7xv5FRI0Z5ZBvldnre+Z77YhvbUzPd+podmtXmmWGRl1xm6tSp1KtXj+zsbLp27cpNN91EVlYWsbGxvPrqqzz//PM899xzvPXWWwBkZWWxevVqVqxYwbhx49i6dSsACQkJrFq1ipo1a/Lqq68CRWW8Y8cOBg0axK5du3j00UeJj49n/vz5vPDCC7z77rvUqlXrF5nS09NZsmQJn3/+OcOGDeO7775jypQpdO3alY0bNxITE8MLL7xAvXr1KCgoYODAgWzevJmHH36Y1157jaVLl9KgQYP/es2EhASmTZvGmjVrsNbSvXt3+vfvT3BwMLt372bWrFn861//4pZbbuHjjz9m9OjR7vgViFQ92ekUzH+IXTaEDnf8lRrVfD2ymSp1RP3GG2/QqVMnevToQXJyMrt378bHx4dbb70VgNGjR7Nq1apzy99+++0A9OvXj8zMzHPjv8OHD6dmzZoArFq1ijFjxgDQvn17WrZsya5du/Dx8WH69OmMGTOG/v3707t37xIzDRs2DGMM0dHRNG7cmOjoaHx8fIiMjCQpKQmAjz76iNjYWDp37sy2bdvYvn37JX/OVatWMWLECAICAggMDGTkyJGsXLkSgPDwcGJiYgDo0qXLuW2ISNmlfPAQAfnp/NjlJTqENvbYdhw5oi7tyNcTli1bxrfffsv3339PrVq1iI+PL/Fc4fNPTbvwNLWfvw8ICDj32KVuDrx7924CAwNJTU296DI1atQAwMfH59zXP3+fn5/Pvn37eOWVV1i3bh3BwcGMHTu21HOcL5Xp/G34+vpq6EPkMp1cP5eQlC+ZFTiam6+7zqPbqjJH1BkZGQQHB1OrVi127NjBDz/8AEBhYSHz5s0D4MMPP6RPnz7n1pkzZw5QdIRap04d6tSp84vX7devHx988AEAu3bt4sCBA7Rr146MjAweeeQRVqxYwfHjx89to6wyMzMJCAigTp06HDlyhAULFpx7LigoiFOnTpWY6dNPP+XMmTNkZWUxf/58+vbte1nbF5FfsqcO4/v1b9liW9F9zF+o5uvZKnXkiNoJgwcPZtKkSXTs2JF27drRo0cPoOjoeNu2bXTp0oU6deqcK2eA4OBgevXqRWZmJlOnTi3xdR944AEmTJhAdHQ01apVY/r06dSoUYOJEyfywAMP0LZtW9577z0GDBhAv379ypy7U6dOdO7cmcjISFq1avVfQyjjx49nyJAhNG3a9NyHmACxsbGMHTuWbt26AXDvvffSuXNnDXOIuIO1HHz/fhoWZLOn1yuMaBLs8U2aS/0z+XLFxcXZC28ckJiYSEREhNu3daUCAwM5ffr0Lx6Pj4/nlVdeIS4uzoFUzvPW35eI046vnEr9xY/xfp37ueORv+Lj454reY0xCdbaEgun1ON1Y0w7Y8zG8/5kGmMedUsyEZEKpODEfmoueYp1tgMDxz7jtpIuTalDH9banUAMgDHGFzgIzPdwrnJT0tE0FH34KCJyTmEhh2eMo05hIWm/+jtdgwNKX8dNyjoCPhD4yVp70Xt7iYhURocXvkzzk+v5pNGDDOnTvVy3Xdaivg2YVdITxpjxxpj1xpj1P1/dJyJSGeQc2ED9tX9jqenOsLueKPcZJl0uamNMdWA4MLek5621k621cdbauIYNS7zjuYhIxXM2i1Mf3MkxW5tao94hOLBG6eu4WVmOqIcAG6y1RzwVRkTE2yTPfoz6OSksifgfuke2diRDWYr6di4y7FHVLVu2jNWrVzsdQ0Tc7GTCJ7TYO4dPao7k5lG/diyHS0VtjKkFXAN84tk4FZOKWqTyKcxIxferh9lmw4m56xWqV3PuQm6XtmytPWOtrW+tzfB0IE+aOXMm3bp1IyYmhvvvv5/9+/fTpk0bjh07RmFhIX379mXRokVAyVOiAixcuJDY2Fg6derEwIEDSUpKYtKkSfz9738nJibm3ORHIlKBFRZyaPpd+BacZW+/f9C6aT1H4zhzCfmCJ+HwFve+ZpNoGHLxO6ckJiYyZ84cvvvuO/z8/HjggQdYvnw5TzzxBBMmTKB79+506NCBQYMGASVPiVpYWMh9993HihUrCA8P58SJE9SrV48JEyYQGBjI448/7t6fSUQccXjRqzRPX8uMhr9lzNVln/rB3arMXB+LFy8mISGBrl27ApCdnU2jRo149tlnmTt3LpMmTWLjxo3nln/jjTeYP7/oup6fp0RNS0ujX79+hIeHA1CvnrN/y4qI+2UfSKD+Dy+x1HTn+rFPesXNnp0p6ksc+XqKtZa77rqLF1988b8eP3PmDCkpKUDRVYpBQUEXnRLVWusVvzQR8ZCcTLJmjiHP1qbWzW9Tz4FT8UpSZaY5HThwIPPmzePo0aMAnDhxgv379/PEE09wxx138Pzzz3PfffcBF58StWfPnixfvpx9+/adew24+HSjIlKBWEvyjPHUzT3E8uiX6B7ZxulE51SZou7QoQN/+ctfGDRoEB07duSaa64hKSmJdevWnSvr6tWrM23aNAYPHkx+fj4dO3bk6aefPjclasOGDZk8eTIjR46kU6dO5+4MM2zYMObPn68PE0UqsKPLJ9MidQFza9/JqBE3Ox3nv1T5aU6lZPp9SVWSm7IZpgxkA+0Jf/Q/NKn7y/ubetqlpjmtMh8mioiUKPc0Ge+PBluLwhHvOlLSpakyQx8iIiVJnvkADXIOsDjiBXrHdHA6TonKtag9Mcwi7qffk1QVaSun0SL5M+YF3s6om527RLw05VbU/v7+HD9+XCXg5ay1HD9+HH9/f6ejiHhU7qHtBC1+gnVE0vuel/Hz8A1qr0S5jVGHhISQkpKC5qr2fv7+/oSEhDgdQ8Rzck9z8t+/ppqtTvbwSTSvF+h0oksqt6L28/M7d0WfiIhjrOXAjPsIyU5iTsQ/uL1LR6cTlcp7j/VFRDzg0DdvEnrwaz6qfSc33zzG6TguUVGLSJVxes9qGqx+jhUmjoH3/ZVqXjwufT6dRy0iVULhqaOcnTWGE7YetX89hYa1azodyWUV468TEZErUZBP6nu3Uys/gx97vkFM24r1eZmKWkQqvZRPniLk5HrmNvktw68d7HScMlNRi0ildmLDfEK2TeJLv2sZOe4PFXKqYhW1iFRauYd3UuOLB9lqWxFx9zsE1KiYH8upqEWkUrI5GZycOoqcQh+ODf0XVzVr4HSky6aiFpHKp7CQ5CmjqZd7kMXRLxPfvcTZQysMFbWIVDrJnzxF6LEVfNTgAUaNvM3pOFfMpaI2xtQ1xswzxuwwxiQaY3p6OpiIyOU4tmYOLba+wwK/axh+75/x8al4Hx5eyNWR9X8AC621o4wx1QHvm1lbRKq87AMbCVzwGzbSloh7JxNUs7rTkdyi1KI2xtQG+gFjAay1Z4Gzno0lIlI2NusYWTNuJd/WInvkdGIa13M6ktu4MvTRCkgDphljfjTGTDHGBFy4kDFmvDFmvTFmvaYyFZFyVZDHwcm3EpR3nO+7/oOenSKdTuRWrhR1NSAW+Ke1tjOQBTx54ULW2snW2jhrbVzDhg3dHFNE5OIOzHqMkIz1zGv2ODdeN9zpOG7nSlGnACnW2jXF38+jqLhFRByX+u1bhO55n0/9b+Smcb+vkFcelqbUorbWHgaSjTHtih8aCGz3aCoRERekb15Io1VP851PF3pNeAd/P1+nI3mEq2d9/Ab4oPiMj73A3Z6LJCJSupyD26g+/2722BDq3zmTRnV/8dFZpeFSUVtrNwIV+9IeEak0Ck+lcWraTVDoR9r1M+gb1szpSB6lKxNFpGLJyyH13ZEE5R1jdde36Nu1s9OJPE5FLSIVh7Xsn34PIac383HLpxl+3TCnE5ULFbWIVBjJnz1Py4Nf8lHtsdxy128q5RkeJVFRi0iFcGT1B7TY+BqLqsVz7YSX8asgN6Z1h6rzk4pIhXVy+xKCFz3MBiJoN34adWpVjjk8XKWiFhGvdiZlK9XmjuaAbYzf6Nm0bFR55vBwlYpaRLxWXnoK2dNuJKvQjyPDZhLdOszpSI5QUYuIV7LZJ0l7dzg18k+R0Odf9I6rujNXqKhFxPvknyV50igaZiexMPJlhl4zyOlEjlJRi4h3sZakqWMJzVjH3OZPcNPNY5xO5DgVtYh4lf0fPUFY6lfMrTOWm++pnLPhlZWKWkS8RvKC12iZ+C5fVx/MkImvVKlzpS9F74KIeIVDK6bTYs1zrPTtTtyD7xHo7+d0JK+hohYRx6UlfEbDJY+xzkQTdv9sGtUJdDqSV1FRi4ijTiYup/YX97KDcOrcPZcWVfCCltKoqEXEMaeSNuA35zYO2gYU3P4RbUObOh3JK6moRcQROYd3UzBjJBm2JkdvnE2ndq2djuS1VNQiUu7yTh7k1JTrKSjIZ+c1M+jRuZPTkbyailpEylXB6WOkvTOUmnknWdt7MgP69HE6ktdTUYtIuSnMSufwW4Opl3uQxTFvMGTQUKcjVQgqahEpFzYng4NvD6FB9j6+jnyVG0bc6nSkCkNFLSIeZ3NPk/zWMJpk7eKLdi8y4uY7nY5UoVRzZSFjTBJwCigA8q21cZ4MJSKVhz17hv1v30CLU5v55Kr/YdTt92n+jjJyqaiLDbDWHvNYEhGpfPJz2f/PmwjNSGBe6FPcPOYhlfRl0NCHiHhGQR77Jt1KWPpqPm72OKPuflwlfZlcLWoLLDLGJBhjxpe0gDFmvDFmvTFmfVpamvsSikjFU5DHvndvJ/zYUuY2epiR9z2Fj49K+nK5WtS9rbWxwBDgQWNMvwsXsNZOttbGWWvjGjZs6NaQIlKBFB9Jhx/9hnn1J3Dj/c/hq5K+Ii4VtbU2tfi/R4H5QDdPhhKRCir/LHv/eTPhaYuZ1+BBbnjgRc0p7QalvoPGmABjTNDPXwODgK2eDiYiFUz+Wfb+cxStji1lXqPfcOPEv6ik3cSVsz4aA/OLPwSoBnxorV3o0VQiUrHk57L3nZtodWIlcxs/wojxz1JNJe02pRa1tXYvoBlTRKRkeTlFJZ2+inlNHmPk+Gc0Ju1mZTmPWkTkv9i8bJLeHkmrk6uZ1/R3jLjvaZW0B6ioReSy2NzTJL0zgvCMtcxr9ntG3qtT8DxFRS0iZVZ4Jp3kt4cRenorc0Of4qa7f6+S9iAVtYiUSX7mUQ6/PZSmOXv5rM0LjLrjAV1x6GEqahFxWe6JZE78cygNzh7iP9F/Z8RNd6qky4GKWkRckn1kD6cmX0dgfgZLu05i2PWjnI5UZaioRaRUp5K3kjdtGH4FZ1nbbxpDBg5xOlKVoqIWkUtK37MWnw9uIr/Qh53XfMjAPv2djlTl6NIhEbmoIz8upMbMYZwurM6+4fPoq5J2hIpaREqUvOJ96n12Byk04vhtX9K9S1enI1VZKmoR+YV9X7xMiyUPsdW0pdo9C+gYEeF0pCpNY9Qi8n+sZfeHj9Nm9xRWVetJ64mzaVK/rtOpqjwVtYgUKchjz5SxtDn0JYtqXkf3B6dSJ9Df6VSCilpEKJq3Y+87o2id8T2fBY/l2omv4l9d9eAt9JsQqeJyT6ZyZNKNhGXvYn6L3zN83FOaAc/LqKhFqrDM/ZvInTGKBvkZfB35KjfePE6XhHshFbVIFXV4w9cEfn4PhbYGawfMZFj8IKcjyUWoqEWqoKRF/6T56qfYR3Oyb5lFfGSU05HkElTUIlVJYSG7Zj9B212TWevbmUb3zKZtsyZOp5JSqKhFqgibl82eyXfSNm0R39QcQtzE9wiuHeB0LHGBilqkCsg9mcqhd0fRJnsbnzWcwODxL1DDT//7VxT6TYlUcul71lHw4W00Lsjki/YvMfy2CTqzo4JxuaiNMb7AeuCgtfZ6z0USEXc5sGImjZY8xgkbROLA2QzrN9DpSHIZynJE/QiQCNT2UBYRcZfCQnbN+SNtd05ik2lP9Ts+oG+b1k6nksvk0ux5xpgQ4DpgimfjiMiVKszOZPebN9B25yQW+w+i2cPfEKGSrtBcPaJ+HfgDEHSxBYwx44HxAKGhoVeeTETK7PThPZx8bxStzibxebNHGDzuGar7+TodS65QqUfUxpjrgaPW2oRLLWetnWytjbPWxjVs2NBtAUXENYd+XEjBu/EEnT3K4rh/Mmz8cyrpSsKVI+rewHBjzFDAH6htjJlprR3t2Wgi4hJr2fXJ/3DV5tdIMs05ecN0BsXqbiyVSalH1NbaP1prQ6y1YcBtwBKVtIh3KMjOYNcbN9J2y6t8V6Mv/hOX0UUlXenoPGqRCupk0iayZ/6aVnmpfNnsN/zq7mc1h3QlVabfqrV2GbDMI0lExGVJy9+n0dLHybM1WNHzPa4fPNLpSOJB+utXpAKx+WfZOfO3tE96n82mPX63z+Dqdu2cjiUepqIWqSDOpO3nyNRf0z57K4sCb6Tr+LcJrh3odCwpBypqkQog+YdPqLPwYRraPL5q9wKDb3tQt8uqQlTUIl7M5ueyY+bjRCTNYCfhZN3wL67TWR1VjopaxEudPvwTadPvICInkW8Ch9P53rdoV7eO07HEASpqES+0f9Uc6n37GA1sIQsj/8qgUffjo6GOKktFLeJFbF42iTMeo0PyLLabq8i/aSqDo2OcjiUOU1GLeIkTe38ka9bddMjbxzdBI4m79w2C61x0HjSpQlTUIk4rLGTXF6/Q8se/UWhr8U3nN/nVDWN0FxY5R0Ut4qDsEwdJnjaWtqfWssavKw1+/S+uCQ93OpZ4GRW1iEP2fzeXOt/+lhaFOSwM+z0DRj+pG85KibRXiJSzgpzT7Pj3w0Qe+pidJpwzwycxuEsPp2OJF1NRi5Sjw1uWYj99gIj8Qyyqdxvdxr1K3SBdBi6XpqIWKQeFuWdInPUEEfveJ5UGrOw1lWsGjdAHhuISFbWIhx3Z/h35n9xPZH4yS4Kuo8Nd/6C/blcnZaCiFvEQm5fDttl/ImLPFNKox7Ju7zJg6K06ipYyU1GLeEDarjXkzB1PVF4SywKupc2dbxDfpInTsaSCUlGLuJHNy2HbnD/TbvcU0qnNki5vMWDYaB1FyxVRUYu4ycFNS+Dzh4kqSGZFrYG0Gv0mVzdv7nQsqQRU1CJXKC8rnZ0zf0fUoY85SEOWdZ1EvyG3abY7cRsVtcgV2LtyDrWXPElEYTqLg0fRcfTfiG9Q3+lYUsmoqEUuw5njKSS9/yAdTi5jt2nJT9dMZmCfa5yOJZVUqUVtjPEHVgA1ipefZ619xtPBRLxSYQE7vnqL5gkvcZXN45tm99Nj9LO0CajldDKpxFw5os4FrrbWnjbG+AGrjDELrLU/eDibiFc5vON7sj95hPZnd7LRN5pqN7zONR3jnI4lVUCpRW2ttcDp4m/9iv9YT4YS8Sa5p46z88MniEqdxwlqszjyBfqOmEh1P1+no0kV4dIYtTHGF0gAWgNvW2vXlLDMeGA8QGhoqDszijjDWnb+ZzKN1rxAZGEmy4NH0P72lxjYuLHTyaSKcamorbUFQIwxpi4w3xgTZa3desEyk4HJAHFxcTrilgot7acNZMx9mHY5W9jq046fhr7PgO79nY4lVVSZzvqw1p40xiwDBgNbS1lcpMLJPX2CHbP/RGTyLKoRwOK2f6LPLY9Sw8/P6WhShbly1kdDIK+4pGsCvwL+6vFkIuXIFuSz/as3ab7hNaLtKVbWHkrr2/7GwOYhTkcTcemIuinw7+Jxah/gI2vtl56NJVJ+9q9fgFn4RyLz97HJN4r8Qf9Lfw1ziBdx5ayPzUDncsgiUq7SU3aSOvdxIjNWFF363elV+gy7m2rVdDaHeBddmShVztmsDBI/+jMR+2cSZn1Z3Px+utz2J+Jr13Y6mkiJVNRSZdj8s2z78i2abfoHnexJVgZcQ8ioFxkY3sbpaCKXpKKWys9adi37kIBVLxBVcJDNvh3Ye/UU+vbW3BxSMaiopVI7sHExZxf8iba529lrWrAy7i16DbkDX18fp6OJuExFLZVS2t5NHJ3/RyJPfccRG8zS9k/TY8RvaOVfw+loImWmopZKJfPoAfbN+xNRRz6nJv4sCbmfzjf/PwbUret0NJHLpqKWSiHrRCq7P36eiIPziLCFrKo3ktY3PcvVIZp3Rio+FbVUaNkn09gx/wXa7/+QaHuW74OupdH1f6J/+2ino4m4jYpaKqScU+kkzn+RNntn0Mnm8EPAAIKHPk2fqFino4m4nYpaKpSzZzLZ/unLtNo1lc6cZo1/bwKufZpenXs6HU3EY1TUUiHkZqWT+NnfablrGjFksr56N/wGPkX37vFORxPxOBW1eLXsjOPs+OxvtNo7kxhOk+DXBeKfoEuvQRhjnI4nUi5U1OKVTqcfZtf8l2h7YDadyWZdjZ74DfgDsd0HqKClylFRi1fJPHKAPZ+9SETqPGJsHmsD+hMw8A907dLb6WgijlFRi1c4vj+R/V+/TOThz+lIAWuDBhJ87ZP0iNZdvkVU1OKo5C3LSf/mVaIyVhCIL+vqXEvjoX+kl86DFjlHRS3lzhYWsGvVx/Ddm7TL3UxtG8DKxqO56vrf0ie0ldPxRLyOilrKTcHZHLb9Zwp1N75Lu4IDHKI+y8IeI3r4b+hfr77T8US8lopaPC7rZBo7v3qT0D0z6GjT2eMTxqro/6XL0HuIr+nvdDwRr6eiFo85tGsDhxa9TsSxBcRylk1+Mezt/jJdBtxEa80HLeIyFbW4lS3IZ8eKubDmXSJyfiTY+rGh7iDqxP+GjjE9dA60yGVQUYtbZGeeIPHrt2m6830i7BEO0YAVoQ/SbuhD9GrSzOl4IhVaqUVtjGkBzACaAIXAZGvtPzwdTCqGQzsTOLT4Ldof/YpYctlarQP7Oj5B7LWj6VdDd1MRcQdXjqjzgd9ZazcYY4KABGPMN9ba7R7OJl4qL/cMid/OwH/Tv2l7djv1bDU21B5I7fiHiIztq+ENETcrtaittYeAQ8VfnzLGJALNARV1FXNk31aSF71Nm0Of05HTHDBNWRH2KO0Gj6dnk+ZOxxOptMo0Rm2MCQM6A2tKeG48MB4gNFS3P6osCvJySVw6G98N04jI+ZF61pcfA3rj2/UeOvUdRmg1X6cjilR6Lhe1MSYQ+Bh41FqbeeHz1trJwGSAuLg467aE4ojDezeTsmQK4SmfEcVJUmnIyhYTaTVoAt1ahDkdT6RKcamojTF+FJX0B9baTzwbSZySc/okid/OIGD7LNqe3U4D68Ommt3Z2+VuYuJH0szPz+mIIlWSK2d9GOA9INFa+5rnI0l5soWF7Fn/Dae+n0b79CV0JpckE8LKsEe46lfj6BIS5nREkSrPlSPq3sAYYIsxZmPxY//PWvu152KJpx0/uJe9i6fQbN/HtLGHOW1rsil4EAE97iKq60DCdOWgiNdw5ayPVYDOt6oEzmSeYOfSD6iROI/22Zvoaixb/DpyoMPDRP5qND2D6jgdUURKoCsTK7n83Gx2rPqE/I1ziMhcTWeTR7JpyuqQewiJv5voNlFORxSRUqioKyFbWMDu9d+QufZD2hxbTBSnOUFtEhoMp3b30XToEk8LDW2IVBgq6srCWg4kruHwdx8Smvo1bW0aZ2wNttbui2+nW4jqewO9amhKUZGKSEVdkVnLgcS1HF49i2ap/yG0MJVm1oet/l1IivgdHa6+jW61g51OKSJXSEVd0fxczt/PptnBhYQWptLcGrbViOHAVeNo3f92YpqEOJ1SRNxIRV0RFJfzoe/nFJfzweJy7nSunDuqnEUqLRW1lyrMz+enjUs5kfApzY8sOe/IuRPJV91N6/630bFJC6djikg5UFF7kbM5Z9j5/RfkbvmCVidW0IYMzlpfEv1jONBqnMpZpIpSUTssM/0oe1Z9jNn5Ne1OrSHa5HLa1iQxqAe7211Huz4j6BTcwOmYIuIgFXV5s5bUn7ZwcN1n1Er6lnY5m4k1haQRzKb6g/GPGk5Ez6F0rVnL6aQi4iVU1OUgNyeL3WsWcmbbAkLSVtLMHqYZkOTTgjXNxlAv9kbaxvanp6/mdhaRX1JRe8jR5F3s/+Ezqu/7lrZZG4gyZ8m21dlZszNJYeNo0e0Gwlq1J8zpoJpSt+QAAAc3SURBVCLi9VTUbpKbk8We9d9yetsiGh9ZQVjhARoBB01jNjYYhn+HIbTrMZiYgCCno4pIBaOivky2sICk7Ws5unEhASkraJ29hUiTx1nry84a0axuMYqmXW8grG1HmvtoXg0RuXwq6jJIO7iX/Wu/wuxbSljmesLJIJyiseYfG43Av/2vaN11ENG6bFtE3EhFfQmZJ4+xb/0icnYupunxHwgtTKEhcIy67A3qyp5WA2jZdShhIa001iwiHqOiPk9GehpJCd+Qs3s5DY6tIzx/L52MJdtWZ3fNjhxsfjMNYwbTqkM3GmiaUBEpJ1W6qDPS09i3/hty9iyn4XnFnGv92F2jAz80u4/aEfG07hxPx5oBTscVkSqqShX1iSMp7N+0lNyfVtHg2Dpa5e8l5hfFfDWtO/cjSsUsIl6i0ha1LSzgwK6NHNm6HJOyhqYZmwixh6gH5Fg/9tSIVDGLSIVQaYo658xp9m1aScauldQ8vJ6w7K20JIuWQDpBJNWKJrnJLdRt35fw6N5E6RJtEakgKmRR28JCDifv5tC27zi7fy3BxzcQnreHCFMAwH6fEHYGx0OL7jSOiie0dTSddS6ziFRQpRa1MWYqcD1w1FrryC2rM06kcWDzCrL2raXm0Y2EZCfSlAyaArnWj73V25LQ7Nf4X9WbsJgBtGzQhJZOBBUR8QBXjqinA28BMzwbpUhuzhmStv5A+u4fqHZoA41PbaOFTSUaKLSGZN8QfqrTiz3NYqnXrictI7oSoZu2ikglVmpRW2tXGGPCPB0kN+cM+1/pT1jeT7QrHsJII5iUWh1IaTySoKu6ExrVm5Z16+toWUSqFLeNURtjxgPjAUJDQ8u8fg3/WmTWaklCQA/8w+JoHtWPRs3DaeiugCIiFZSx1pa+UNER9ZeujlHHxcXZ9evXX1kyEZEqxBiTYK2NK+k5nQohIuLlVNQiIl6u1KI2xswCvgfaGWNSjDH3eD6WiIj8zJWzPm4vjyAiIlIyDX2IiHg5FbWIiJdTUYuIeDkVtYiIl3Ppgpcyv6gxacD+y1y9AXDMjXHcRbnKRrnKRrnKpjLmammtLfFibI8U9ZUwxqy/2NU5TlKuslGuslGusqlquTT0ISLi5VTUIiJezhuLerLTAS5CucpGucpGucqmSuXyujFqERH5b954RC0iIudRUYuIeLlyK2pjzGBjzE5jzB5jzJMlPG+MMW8UP7/ZGBPr6roeznVHcZ7NxpjVxphO5z2XZIzZYozZaIxx650SXMgVb4zJKN72RmPMn11d18O5fn9epq3GmAJjTL3i5zz5fk01xhw1xmy9yPNO7V+l5XJq/yotl1P7V2m5nNq/WhhjlhpjEo0x24wxj5SwjOf2MWutx/8AvsBPQCugOrAJ6HDBMkOBBYABegBrXF3Xw7l6AcHFXw/5OVfx90lAA4fer3iK7rpT5nU9meuC5YcBSzz9fhW/dj8gFth6kefLff9yMVe5718u5ir3/cuVXA7uX02B2OKvg4Bd5dlh5XVE3Q3YY63da609C8wGbrhgmRuAGbbID0BdY0xTF9f1WC5r7WprbXrxtz8AIW7a9hXl8tC67n7t24FZbtr2JVlrVwAnLrGIE/tXqbkc2r9ceb8uxtH36wLluX8dstZuKP76FJAINL9gMY/tY+VV1M2B5PO+T+GXP+TFlnFlXU/mOt89FP2N+TMLLDLGJJiim/u6i6u5ehpjNhljFhhjIsu4ridzYYypBQwGPj7vYU+9X65wYv8qq/Lav1xV3vuXy5zcv0zRPWQ7A2sueMpj+5jb7kJeClPCYxeeF3ixZVxZ93K5/NrGmAEU/Y/U57yHe1trU40xjYBvjDE7io8IyiPXBormBjhtjBkKfAq0cXFdT+b62TDgO2vt+UdHnnq/XOHE/uWyct6/XOHE/lUWjuxfxphAiv5yeNRam3nh0yWs4pZ9rLyOqFOAFud9HwKkuriMK+t6MhfGmI7AFOAGa+3xnx+31qYW//coMJ+if+KUSy5rbaa19nTx118DfsaYBq6s68lc57mNC/5Z6sH3yxVO7F8ucWD/KpVD+1dZlPv+ZYzxo6ikP7DWflLCIp7bxzwx8F7CQHw1YC8Qzv8NpkdesMx1/PdA/FpX1/VwrlBgD9DrgscDgKDzvl4NDC7HXE34vwuWugEHit87R9+v4uXqUDTOGFAe79d52wjj4h+Olfv+5WKuct+/XMxV7vuXK7mc2r+Kf/YZwOuXWMZj+5jb3lwXftChFH1S+hPwVPFjE4AJ570Rbxc/vwWIu9S65ZhrCpAObCz+s7748VbFb/gmYJsDuR4q3u4mij6E6nWpdcsrV/H3Y4HZF6zn6fdrFnAIyKPoCOYeL9m/Ssvl1P5VWi6n9q9L5nJw/+pD0XDF5vN+V0PLax/TJeQiIl5OVyaKiHg5FbWIiJdTUYuIeDkVtYiIl1NRi4h4ORW1iIiXU1GLiHi5/w+OQb9ScWGMMAAAAABJRU5ErkJggg==\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"stIRVhYaJJqT"},"source":["$nt = 101$ よりも近似の精度が良くなりました。\n","この辺は単純に振る舞ってくれるようです。"]},{"cell_type":"markdown","metadata":{"id":"PRtZT2Z5JJqT"},"source":["## 少し現実的な設定にしてみる\n","もう少し現実的な設定にしてみます。\n","個体数 $x (t)$ が大きくなると食糧が足りなかったり、\n","衛生状態の悪化で病気にかかりやすくなったりして増加率が減るはずです。\n","そこで増加率 $\\alpha$ を $\\alpha - \\beta x$ ($\\alpha, \\beta > 0$) に変えてみます。\n","\n","\\begin{align}\n"," \\frac{dx}{dt}\n"," =\n"," \\left(\\alpha - \\beta x \\right) x.\n","\\end{align}\n","\n","これには**ロジスティック方程式**という名前がついています。\n","悪化を $\\beta x$ と書くことに必然性はありません。\n","とりあえずやってみただけです。"]},{"cell_type":"markdown","metadata":{"id":"dCmkNNVRJJqU"},"source":["## ロジスティック方程式の厳密解\n","解き方はさておき解は次のように書けます。\n","あなたがもしこの関数を微分できるなら、\n","厳密解になることを計算して確認してみてください。\n","sympy を使ってみてもいいでしょう。\n","\n","\\begin{align}\n"," x (t)\n"," =\n"," \\frac{\\alpha}{\\beta} \\frac{1}{1 + e^{- \\alpha t}}.\n","\\end{align}\n","\n","定性的に言うと、増加率が減っていくことで人口が飽和していくことが分かります。\n","\n","さて、 もとの微分方程式を近似してみましょう。\n","次のように書けます。\n","\n","\\begin{align}\n"," \\frac{x_{n+1} - x_{n}}{h}\n"," =\n"," (\\alpha - \\beta x_{n}) x_{n}, \\quad\n"," x_{n+1}\n"," =\n"," x_n + h (\\alpha - \\beta x_{n}) x_{n}.\n","\\end{align}\n","\n","これも厳密解と比較しながら計算してみます。\n","数値的には $\\alpha = 2$, $\\beta = 1$ で計算しています。"]},{"cell_type":"code","metadata":{"scrolled":true,"id":"CsnaO6vfJJqU","outputId":"80c28c13-b3ad-47c5-f0a0-d8a7e7d15257"},"source":["import numpy as np\n","import matplotlib.pyplot as plt\n","\n","def logistics_euler(nt, init = 10):\n","    dt = T / (nt - 1)\n","    # 初期条件設定\n","    x = np.zeros(nt)\n","    x[0] = init\n","\n","    for i in range(1, nt):\n","        x[i] = x[i-1] + dt * (alpha - beta * x[i-1]) * x[i-1] # ここで方程式を解いている\n","\n","    return x\n","\n","# パラメータ設定\n","alpha = 2\n","beta = 1\n","init = 1\n","nt = 101\n","T = 5 # 時間変化を見る最大値\n","\n","# 近次解\n","x_approx = logistics_euler(nt, init)\n","plt.plot(np.linspace(0, T, nt), x_approx)\n","\n","# 厳密解\n","t = np.linspace(0, T, nt)\n","x_exact = (alpha / beta) / (1 + np.exp(- alpha * t))\n","plt.plot(t, x_exact)\n","\n","# 凡例\n","plt.legend(['approximation', 'exact'])\n","# 描画\n","plt.show()"],"execution_count":null,"outputs":[{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXyU9b328c93JhskYckCyo7KXsKqKAjiUgpWqFp7WqtW7ako2lb7HE/t8Xk8ak977Dm1tlVrKVbEFhFBxK3izr5Jwr6DrGExIYFANpKZ+T1/TKQogSQwyZ2ZXO/XK6/Mcs99XxPg4pff3Is55xARkejn8zqAiIhEhgpdRCRGqNBFRGKECl1EJEao0EVEYkScVxvOyMhwXbp08WrzIiJRKScn55BzLrO65zwr9C5dupCdne3V5kVEopKZ7T7dc5pyERGJESp0EZEYoUIXEYkRKnQRkRihQhcRiRE1FrqZdTSzuWa2ycw2mNn91SxjZva0mW03s7VmNrB+4oqIyOnUZrfFAPBvzrmVZpYK5JjZh865jSctMwboVvU1BPhz1XcREWkgNRa6c+4AcKDq9jEz2wS0B04u9G8Bf3Phc/EuM7NWZnZ+1WtFpB6EQo5AyBEMOQKhEIGgI+jC94MhRzAYIhCoJBSsJBQMhL8HArhQJaFgEBcM4EJBXChAKBSEULDqfhAXCoEL4oJBnAviQi58PxTCuSC48H1C7qT7IZwLYVXPOQe4EBA6cducw+EwFyJcF+HbOIcDDFf1mionTu8dXkf4+fDr/vn8qbeNk1cRCr+Ofy56yp0vnUa8ro/XzE5+T0By9+Fkjbi+TuuojTodWGRmXYABwPKvPNUe2HvS/dyqx75U6GY2HhgP0KlTp7olFWmEKoMhissDlFQEKKsIUlIRpPR4gLKKAMfLjhEoKyZwvJhgeTHueAmuogRXWQqVZVhlGRYox4Ll+IPl+ILH8Ycq8IfC3+NCFcS5CuJcJX5XSZyrJN5VEkeAOBcgjiAJBIi38O14AiQQJI4QcQTwm6510FgtC1WAl4VuZinALOAB59zRrz5dzUtO+dvknJsETAIYPHiw/rZJoxAIhigsreBwSSWFJRUUllRwpKyCY8eOETiWT7C4AFdaiK+8EH9FEQkVRSQGjpEYLCbZldCCUlKtjFRKOc/KSKacZMrx1bFQA/iptAQqLZGAJRDwxRP0JxC0eIK+BEK+JEK+FjhfHCFf+PlKXxzOnwC+OEK+eMwXh/PFgS8O88fjfHGYzw/+eLA4zO+v+h4HPj8+nx98fjA/5veHlzUf5ovDfL6q2/7wfTN8/qplzcAXXt4w8PvwmQ/44nHDsPDzZuHlzQdVt8182Bf3feHXGQY+C9/2hSsl/JjvxO3wegwj/Fqg6jFOuv3Fa//5EaFVrffEMl88bictYyfX2EnL+Kpfvvra48T2z/T4pdUvcc5qVehmFk+4zF92zr1ezSK5QMeT7ncA9p97PJGzFwo5DhUfZ9+RMg4WlXPwaDkHi8ooPpyPK9qHr3g/iWV5pFQeoi2HybAiMqyIvhSRZsdIsfLTrrvSEihPbEFlXCqBhFSCCW1xCamEElIoTkqlNCkVf9VXXFIK8c1SiE9KxZ+YDPHNIL551fdmEJcEcUnE+eOIA5o13I9IYkyNhW7h/7ZeADY55546zWJvAT82s+mEPwwt0vy5NITyyiB7CkvZkV/CroIS9hSWcuhQPhR8RrPivbRzn9PRPqeDHWKEHaKdFdDMKr68kjgoj29NRbNMQs3bYil9sRZtqGzRhvjUDGieDs3SoHkaNGsNSa2Ij08i3pu3LHJatRmhDwNuA9aZ2eqqxx4GOgE45yYC7wLXAtuBUuDOyEeVpqwiEOKz/GI2HTjKls+Psf3gMQrz9pJStI2LLJduto8Bvv3c5DtABkXhF1X97a5ITCPQoiP+1oNJSO8MLTtAi3aQ2g5anA/JbUiKSyDJu7cnEhG12ctlEaedLDqxjAPui1QoadoqgyE2HzjG6twjrM8tYn3uYSryt9PL7eBrvp2M9O1ign8vrdxRSAi/JpDYCjJ7EJc5GNIvgrQLIa0rtO5CQmLqF4uJxDTPTp8r8oXi4wFydh9m+Y4CsncfZnvuQXoHtzDYt4VxcZ/xn/YZyfHFAIT8iVib3th534K2X4O2vSGzF3HJGaf/MEqkiVChS4OrDIZYs/cIC7YdYuG2fLbm5jGQzQz3r+dXiVu5yP8ZPn8QZz5o0wvrcBO0HwTtBuDL7BneY0NETqFClwZRVFbJvC15fLQpj/lbPif9eC5X+1fxWLP19EncQJyrxPnisXaDofM46DwU63gJJKZ6HV0kaqjQpd4UlVby/saDvLvuAIu359Mr9Bk3JK3kkYRs2nxxHFqrnnDR3XDBlVjnyyAh2dvQIlFMhS4RVREIMW9LHq+v3Mcnm/PoEtrNbcnL+UPyElpVHMRZHNZhOPS8H7qNgtadvY4sEjNU6BIRewpKmfbpHmZm7yVQUsitzZfxWIuFnFe2DRf0Y12ugq89inUfHd6fW0QiToUuZ805x8Jth3hh0U7mb81jsH87z7VewGC3EH+oAlr1h5G/xfrcACnVXqRcRCJIhS51VhkMMXvVPv66cAc7Pj/CLckr+N+Mj2hbvAkqW8DgO2DgD+C8vl5HFWlSVOhSa8cDQV7LyeW5uZ9ReOQwP229lNtbvU3z8oOQ1AOu+B1kfQ8SU7yOKtIkqdClRsGQ4/WVufz+w60UFhXxi/SFfL/lbBLKDkOnoXD5M9Dt6zqwR8RjKnQ5o7mb8/jNnM3s+PwwP09fzO2tZpFQcgguvBqueAg66cJUIo2FCl2qtetQCb98ZyOfbP6cW1uu4/X0V0gu2Q1dhsOV/xc6X+Z1RBH5ChW6fEl5ZZBnP9nOpAU76BW3j6XnTeP8IzmQ2gOufw0uukZTKyKNlApdTsjeVcjPX1vLgUMF/Ln9h1x1eCZ2PBW++TsYeAf49ddFpDHTv1ChrCLI/7y3mZeW7mJcylbeTX+epIJ9MOBWuOaXkJzudUQRqQUVehO3cf9Rfjp9FQfy8pnR4W0uPvQGtLgI/mUOdB7qdTwRqQMVehPlnOPFxbv4zZzNXJ60g3fSnyPp0D647Mdw1f8LX+tSRKKKCr0JKj4e4N9nruG99ft58vy53HhkChbfHu6co71XRKKYCr2J2Z5XzD1TcziSv58F7abQsXAZ9L4exv4RmrXyOp6InAMVehMyd0seP5m2in7+Xfwj7Q8kFhXAdb+HQXdqV0SRGKBCbyL+vmw3j765nglpOTxY/ifMnw4/fA/aD/Q6mohEiAo9xoVCjifmbOL5hTt4tu0/uK5oGnQeBt95Sae0FYkxKvQYVhkM8e8z1/Du6t28ef7L9Dv8AQy4LTzNogsti8QcFXqMKq8M8uNpq1i+aSdzz/sz7Q/nhHdHHP6g5stFYpQKPQaVHA9w19+y2fTZTua1+SPpR7fDDZOg33e9jiYi9UiFHmPKKoL860sr2LnzM+ZnPkWL0n3wvWnQfZTX0USknqnQY0h5ZZDxf89mz85tfJz2v6QcPwy3vAZdh3sdTUQagAo9RlQEQkyYmsPW7dv4KO23pASOwm1vQMeLvY4mIg1EhR4DQiHHv81cw9ot2/kk/UlSKwtV5iJNkAo9Bvz63U0sWLOVT9J/R8vjB+HWWSpzkSZIhR7lnl+wg6mLtvBh+jOkle+B78+ALsO8jiUiHvDVtICZTTazPDNbf5rnW5rZ22a2xsw2mNmdkY8p1Xln7X6eeHcDr6Y/T8eS9diNk+DCK72OJSIeqbHQgSnA6DM8fx+w0TnXDxgJ/M7MEs49mpzJ+n1FPDhzNc+1foX+JYux0b+BPjd4HUtEPFRjoTvnFgCFZ1oESDUzA1Kqlg1EJp5UJ+9YOXf9LZv7EuYwuuwfMOx+uPQer2OJiMdqM0KvybNAL2A/sA643zkXqm5BMxtvZtlmlp2fnx+BTTc9xwNB7vl7Dv1Kl/Pj4N/D5zK/+jGvY4lIIxCJQv8GsBpoB/QHnjWzFtUt6Jyb5Jwb7JwbnJmpM/2djcff3sjRvet5NvFZ7Ly+cP2fwReJP0YRiXaRaII7gddd2HZgJ9AzAuuVr3hj1T7+sXwjM1s8TVxiMtz8CiQ09zqWiDQSkSj0PcDVAGbWFugB7IjAeuUk2/OO8X9nr+HFFs/TKpAH33sZWnbwOpaINCI17oduZq8Q3nslw8xygUeBeADn3ETgv4ApZrYOMOAh59yhekvcBJVWBJgwdSUT/G8xsGIFXPskdLzE61gi0sjUWOjOuZtreH4/oFP51aNfvr2RjEPLuS/xVfjat+HiH3kdSUQaIR0p2si9v+EgH69Yx9zUP2MtL4Kxf9QFKkSkWir0RizvWDkPz1rDpJTnSXZl8C9/g8RUr2OJSCOlQm+knHP8/LW13Fj5DoP8q8LXAW3Ty+tYItKIqdAbqWmf7uHA1hxeaDYdLhoDg3SKHBE5MxV6I7T/SBm/e3ctb6T8BV9CSxj3jObNRaRGKvRGxjnHw7PXca+bQafKnfCdmZCio2pFpGY6ZryReWP1Pgq3LuOHvndg4O26uLOI1JpG6I3IoeLj/Pdba5jV/K9Y8vkw6r+8jiQiUUSF3oj86p2N3B6YRSf/brhuJiS19DqSiEQRTbk0Est2FLBlzVImxL0BWd/TVIuI1JlG6I1AZTDEo2+s5fdJk/E1S4PRT3gdSUSikAq9EZiyeBcDC96md/w2GP08NE/zOpKIRCEVuscOFpXz0kfZvJ84AzoNh77f8TqSiEQpzaF77Ik5m/gZU2lOOXzzdzqASETOmgrdQ6v2HGbfmk/4tm8+NvQnkNnD60giEsVU6B5xzvHf76zn14l/J9SiA4x40OtIIhLlNIfukXfXHaTrvjfpEb8TRk2GhGSvI4lIlFOhe+B4IMgf56xkesJMXPshWJ8bvY4kIjFAUy4eeGnJLr51bDpp7gg25gl9ECoiEaERegMrKqvkjU+W8EbcHOh3M7Qf5HUkEYkRGqE3sL8u3MGE4FTi4uLg6v/0Oo6IxBAVegM6VHycpYs+Zqx/Gb6hP4YW7byOJCIxRIXegP40dzs/ZRrBpNYw9CdexxGRGKNCbyD7jpTx2fJ3GeFbh3/Egzo1rohEnAq9gTzz0VYe9L1CILU9XPwjr+OISAzSXi4NYG9hKcdWvU5W/Gdw1Z8gPsnrSCISgzRCbwAT523lgbjXqEzrFt5VUUSkHmiEXs/2HymjeOUsusXlwlUvgs/vdSQRiVEaodezv8zbyo99s6hM6w69r/c6jojEMI3Q69HnR8spynmNbv59VaNz/f8pIvWnxoYxs8lmlmdm68+wzEgzW21mG8xsfmQjRq9J87Zxr82iQqNzEWkAtRkyTgFGn+5JM2sFPAeMc871AXQNNeBwSQWHs2fS3bePhKv+Q6NzEal3NbaMc24BUHiGRb4PvO6c21O1fF6EskW1vy3ZyXhe53hrjc5FpGFEYtjYHWhtZvPMLMfMfnC6Bc1svJllm1l2fn5+BDbdOJVVBNm59HV6+vaSOPJBjc5FpEFEomnigEHAN4FvAI+YWffqFnTOTXLODXbODc7MzIzAphunmdl7uC0wi+MpHeBr3/Y6jog0EZEo9FzgPedciXPuELAA6BeB9UalQDDE8nnvMMi3jYTh94NfOxKJSMOIRKG/CQw3szgzaw4MATZFYL1R6d31B7mpbCbHE9Owgbd5HUdEmpAah49m9gowEsgws1zgUSAewDk30Tm3yczeA9YCIeCvzrnT7uIYy5xzfPjJhzzjX0No6CMQ38zrSCLShNRY6M65Gk8+4pz7LfDbiCSKYit2Hebrha9QmZBM/CU6o6KINCztfhFBb8xdyrX+5TD4TmjWyus4ItLEqNAjZG9hKV13vIwPI37oBK/jiEgTpEKPkOmLNvA9/ycc7zEWWnbwOo6INEHapy4Cio8HCK6cSqqVwfCfeh1HRJoojdAj4PXs3dwS+gfFbS+GDoO8jiMiTZQK/RyFQo7PFr5KR18+KSPv9zqOiDRhKvRztOSzAq4rnU1x847Q41qv44hIE6ZCP0fz5n/Ixb6tJA69W5eXExFPqdDPwYGiMrrtnk6FL4n4QTrMX0S8pUI/B7MXreNbvsVU9P6ODiQSEc9pt8WzVBkMUZnzd5KsEobrQCIR8Z5G6Gfpg/X7uSEwhyOZF0PbPl7HERFRoZ+tDfNfo5Mvn9Qr7vM6iogIoEI/KzvyixmSP4vihEz8va7zOo6ICKBCPyvvL1zGFf61MOgO8Md7HUdEBFCh11lFIETSuqmE8JFy6Z1exxEROUGFXkefbNjLdaGPKWw3Elq29zqOiMgJKvQ62rZwJpl2lNYj7vY6iojIl6jQ62BvYSkD8mZzNKEt/u5f9zqOiMiXqNDr4INFS7nctx438Ac6b4uINDoq9FoKhhzxa/5OEB8th/7Q6zgiIqdQodfSkq0HGBP4mEPnj4QW7byOIyJyChV6LW1eMItMO0raiLu8jiIiUi0Vei0UlVVyQe5sjsWlE999lNdxRESqpUKvhY9WrOUKW0VZ738Bv05QKSKNkwq9Fo4tn0qchcgcrg9DRaTxUqHXYPvnR7n82HscbNkfy+zudRwRkdNSoddgyfz3uci3n+ZDbvc6iojIGanQzyAYcrTYPJ1yS6LFoO94HUdE5IxU6GewbMterg4u5lCnMZCY6nUcEZEzUqGfwc5Fr5JqZfowVESiQo2FbmaTzSzPzNbXsNzFZhY0s5siF887pRUBOue+zeH480i84HKv44iI1Kg2I/QpwOgzLWBmfuB/gPcjkKlRWJizjqGspaTnTeDTLzIi0vjV2FTOuQVAYQ2L/QSYBeRFIlRjcHj5NPzmaDfiDq+jiIjUyjkPPc2sPXADMLEWy443s2wzy87Pzz/XTdeb/KPlDCicw76Ur+HL7OZ1HBGRWonEXMIfgIecc8GaFnTOTXLODXbODc7MzIzApuvH4sVz6eHbS/yA73sdRUSk1iJxYpLBwHQzA8gArjWzgHPujQis2xtrplNJHG0uu9nrJCIitXbOhe6c6/rFbTObArwTzWW+4/MjDCubS26bEXRtnuZ1HBGRWqux0M3sFWAkkGFmucCjQDyAc67GefNos2bBm9xgRcRfeqvXUURE6qTGQnfO1XrewTl3xzml8ZhzjuZbZlNiybTK+qbXcURE6kQ7WJ9k4548hlUu5UD7b0B8ktdxRETqRIV+kq0LZpBi5bQdqukWEYk+KvQqoZAjfedbHPGnk9pzpNdxRETqTIVeZfW2nVwazOFQl+vA5/c6johInanQq+xZNJ0EC9J++A+8jiIiclZU6EAgGKLD3rc5GN+RZp0HeR1HROSsqNCBles3MtBt4li36yF8xKuISNRRoQN5y17BZ45OI27zOoqIyFlr8oVeGQzR6cD77E3sRuJ5PbyOIyJy1pp8oa9cs5ostlHe/VteRxEROSdNvtAPLZsOQKcrdDCRiES3Jl3oFYEQF+S9z+5mvUnM6FrzC0REGrEmXegrV35KL3ZxvOf1XkcRETlnTbrQj6x4lRBGl+G3eB1FROScNdlCP14ZoFv+++xsnkVCWgev44iInLMmW+irc5ZyIfsI9LrB6ygiIhHRZAu9KOc1Qhhdh3/P6ygiIhHRJAu9IhDiwvwP2dG8Hwmtzvc6johIRDTJQl+zsmq6pcc4r6OIiERMkyz0ouyZhJzRZUStL5cqItLoNblCDwRDdM37kB3Ns0hq3c7rOCIiEdPkCn3tquVcSC4VPcZ6HUVEJKKaXKEfXjGDkDO6jvi+11FERCKqSRV6MOTo8vmH7GjWl2Zp7b2OIyISUU2q0NevWcGF7KWsu6ZbRCT2NKlCL1gxE4ALRuhgIhGJPU2m0J1ztD/wIZ8l9iI5o5PXcUREIq7JFPqWTevo4XZSfMG1XkcREakXTabQDyybAUCXy3UwkYjEpiZR6M452uR+wM74i2jZvpvXcURE6kWNhW5mk80sz8zWn+b5W8xsbdXXEjPrF/mY52bnjm30CW3hSOfRXkcREak3tRmhTwHO1IQ7gSucc1nAfwGTIpArovYuCU+3dBz2XY+TiIjUnxoL3Tm3ACg8w/NLnHOHq+4uAxrd5X9a736Pvf6OZHTN8jqKiEi9ifQc+r8Cc073pJmNN7NsM8vOz8+P8Karty93D30q15Pf8RsNsj0REa9ErNDN7ErChf7Q6ZZxzk1yzg12zg3OzMyM1KbPaMeiGfjNcd6l32mQ7YmIeCUuEisxsyzgr8AY51xBJNYZKck73+Ogrw3tegzxOoqISL065xG6mXUCXgduc85tPfdIkVNQcIg+5avYd941YOZ1HBGRelXjCN3MXgFGAhlmlgs8CsQDOOcmAv8JpAPPWbg0A865wfUVuC62LJrFUAuQNugGr6OIiNS7GgvdOXfGQyudcz8CfhSxRBEUt+UfFNKSLv2v8jqKiEi9i9kjRYtLiuldspydGSMxf0Q+KhARadRittA3LXqLFCsnud+3vI4iItIgYrbQAxvfpphmdBvyTa+jiIg0iJgs9IqKCnocWci2lsPwJyR5HUdEpEHEZKFvWv4BaXYMf29dak5Emo6YLPTitW9y3MXTfdj1XkcREWkwMbf7RygYomv+PLYkDyIrpZXXcUQalcrKSnJzcykvL/c6itQgKSmJDh06EB8fX+vXxFyhb1u3jB7kceCi+7yOItLo5ObmkpqaSpcuXTAdPd1oOecoKCggNzeXrl271vp1MTflcij7dULOuPDym7yOItLolJeXk56erjJv5MyM9PT0Ov8mFXOF3nb/R2xN7E2rNo3utOwijYLKPDqczZ9TTBX63h2buSi0k6OdR3kdRUSkwcVWoS+dCUDHyzTdIiKn99Zbb/Gb3/wmIuv6wx/+QGlp6Yn71157LUeOHInIuusqpgq9xe4P2OXrxPkXfM3rKCLSQAKBQJ1fM27cOH7xi19EZPtfLfR3332XVq282cMuZvZyKcg/QK/j68jueCddvA4jEgUef3sDG/cfjeg6e7drwaNj+5xxmeuvv569e/dSXl7O/fffz/jx40lJSeHuu+9m7ty5tG7dmunTp5OZmcnIkSPp378/n376KUePHmXy5MlccsklPPbYY+zfv59du3aRkZHB5MmTmTBhAtnZ2cTFxfHUU09x5ZVX8tRTT7F+/XomT57MunXruPnmm/n000+ZMWMG2dnZPPvss9xxxx00a9aMzZs3s3v3bl588UVeeuklli5dypAhQ5gyZQoAEyZMYMWKFZSVlXHTTTfx+OOP8/TTT7N//36uvPJKMjIymDt3Ll26dCE7O5uMjAyeeuopJk+eDMCPfvQjHnjgAXbt2sWYMWO4/PLLWbJkCe3bt+fNN9+kWbNm5/zzj5kR+vZFs/CbI32wzn0u0phNnjyZnJwcsrOzefrppykoKKCkpISBAweycuVKrrjiCh5//PETy5eUlLBkyRKee+45fvjDH554PCcnhzfffJNp06bxpz/9CYB169bxyiuvcPvtt1NeXs4DDzzA9u3bmT17NnfeeSd/+ctfaN68+SmZDh8+zCeffMLvf/97xo4dy89+9jM2bNjAunXrWL16NQC//vWvyc7OZu3atcyfP5+1a9fy05/+lHbt2jF37lzmzp37pXXm5OTw4osvsnz5cpYtW8bzzz/PqlWrANi2bRv33XcfGzZsoFWrVsyaNSsiP9uYGaHHb5tDHulcmHW511FEokJNI+n68vTTTzN79mwA9u7dy7Zt2/D5fHz3u98F4NZbb+XGG288sfzNN4cvyTBixAiOHj16Yn563LhxJ0a1ixYt4ic/+QkAPXv2pHPnzmzdupWsrCymTJlCVlYWd999N8OGDas209ixYzEz+vbtS9u2benbty8Affr0YdeuXfTv358ZM2YwadIkAoEABw4cYOPGjWRlZZ32fS5atIgbbriB5ORkAG688UYWLlzIuHHj6Nq1K/379wdg0KBB7Nq166x+ll8VE4VeVlJMz5IVrM/8Jm18MfNLh0jMmTdvHh999BFLly6lefPmjBw5stp9rU/eZe+ru+99cf+LooTwgTins23bNlJSUti/f/9pl0lMTATA5/OduP3F/UAgwM6dO3nyySdZsWIFrVu35o477qhxH/EzZTp5G36/n7KysjOuq7Ziov02L3mL5nac5lnjvI4iImdQVFRE69atad68OZs3b2bZsmUAhEIhXnvtNQCmTZvG5Zf/8zftV199FQiPeFu2bEnLli1PWe+IESN4+eWXAdi6dSt79uyhR48eFBUVcf/997NgwQIKCgpObKOujh49SnJyMi1btuTzzz9nzpw5J55LTU3l2LFj1WZ64403KC0tpaSkhNmzZzN8+PCz2n5txcQIvWLDOxyjGT0uHeN1FBE5g9GjRzNx4kSysrLo0aMHl156KRAebW/YsIFBgwbRsmXLEyUO0Lp1a4YOHXriQ9Hq3Hvvvdxzzz307duXuLg4pkyZQmJiIhMmTODee++le/fuvPDCC1x55ZWMGDGizrn79evHgAED6NOnDxdccMGXpm7Gjx/PmDFjOP/88780jz5w4EDuuOMOLrnkEiD8oeiAAQMiNr1SHTvTrwX1afDgwS47O/uc1xMMBDjyqwvYlTqIQf82OwLJRGLXpk2b6NWrl9cxTpGSkkJxcfEpj48cOZInn3ySwYMbxXXnG1x1f15mluOcq/YHEvVTLltzPiGdIuh5rddRREQ8FfVTLkdWvUGF89N92I01LywijVJ1o3MIf4gqtRfVI3TnHO0/n8uWpH6ktkr3Oo6IiKeiutD3bF1NJ7ef0gtGex1FRMRzUV3o+5eFd0HqMkwn4xIRiepCT9v7Edv8F9G2w4VeRxER8VzUFvqhg3voVrmFQx2u8TqKiDQC8+bNY8mSJV7H8FTUFvpni2fhM0fbi7V3i4io0CGKd1tM3P4e+60NXXtf7HUUkeg05xdwcF1k13leXxhz5gtHTJ06laeffpqKigqGDBnCww8/zDXXXMPSpUtJS0vjiiuu4JFHHmHUqFHVnmoX4L333uPhhx8mGAySkZHBCy+8wMSJE/H7/UydOpVnnnmm3g+zb4yistBLi4voWZrDmrbX05icc+wAAAYMSURBVE4n4xKJGps2beLVV19l8eLFxMfHc++99zJ//nweeugh7rnnHoYMGULv3r0ZNSp8GcnJkyeTlpZGWVkZF198Md/+9rcJhULcddddLFiwgK5du1JYWEhaWhr33HMPKSkpPPjggx6/S+9EZaFvXvwWA62SZJ2MS+Ts1TCSrg8ff/wxOTk5XHxx+DfrsrIy2rRpw2OPPcbMmTOZOHHiifOPQ/Wn2s3Pz2fEiBF07doVgLS0tAZ/H41VjYVuZpOB64A859wp13az8Lks/whcC5QCdzjnVkY66MmCG9/hKMn0uOQb9bkZEYkw5xy33347TzzxxJceLy0tJTc3FwgfNZqamnraU+065045pa6E1Wa+YgpwpiN3xgDdqr7GA38+91inF6is4KKixWxtcRnxCYk1v0BEGo2rr76a1157jby8PAAKCwvZvXs3Dz30ELfccgu//OUvueuuu4DTn2r3sssuY/78+ezcufPEOuD0p7FtSmosdOfcAqDwDIt8C/ibC1sGtDKz8yMV8Ku2Zn9Ma45hPa+rr02ISD3p3bs3v/rVrxg1ahRZWVl8/etfZ9euXaxYseJEqSckJPDiiy8yevRoAoEAWVlZPPLIIydOtZuZmcmkSZO48cYb6dev34krHY0dO5bZs2fTv39/Fi5c6OXb9EytTp9rZl2Ad04z5fIO8Bvn3KKq+x8DDznnTjk3rpmNJzyKp1OnToN2795d58Cbl39Axdz/5YJ7Z5LSonWdXy/SlDXW0+dK9ep6+txIfCha3WRWtf9LOOcmAZMgfD70s9lYzyGjYMios3mpiEhMi8Q+f7lAx5PudwBOf/E+ERGpF5Eo9LeAH1jYpUCRc+5ABNYrIvXAq6uUSd2czZ9TbXZbfAUYCWSYWS7wKBBftcGJwLuEd1ncTni3xTvrnEJEGkRSUhIFBQWkp6dr179GzDlHQUEBSUlJdXpdjYXunLu5hucdcF+dtioinujQoQO5ubnk5+d7HUVqkJSURIcOHer0mqg8UlREzk58fPyJIywl9uhEKCIiMUKFLiISI1ToIiIxolZHitbLhs3ygbofKhqWARyKYJxooPfcNOg9Nw3n8p47O+cyq3vCs0I/F2aWfbpDX2OV3nPToPfcNNTXe9aUi4hIjFChi4jEiGgt9EleB/CA3nPToPfcNNTLe47KOXQRETlVtI7QRUTkK1ToIiIxIuoK3cxGm9kWM9tuZr/wOk99M7PJZpZnZuu9ztJQzKyjmc01s01mtsHM7vc6U30zsyQz+9TM1lS958e9ztQQzMxvZquqrnwW88xsl5mtM7PVZnbKVd3Oef3RNIduZn5gK/B1whfWWAHc7Jzb6GmwemRmI4BiwtdtPeUSgLGo6pq05zvnVppZKpADXB/jf84GJDvnis0sHlgE3F91nd6YZWb/BxgMtHDOxfyFgs1sFzDYOVcvB1JF2wj9EmC7c26Hc64CmE74ItUxqxYX6Y45zrkDzrmVVbePAZuA9t6mql9VF1kvrrobX/UVPaOts2BmHYBvAn/1OkusiLZCbw/sPel+LjH+D72pq7pA+QBgubdJ6l/V9MNqIA/40DkX6+/5D8DPgZDXQRqQAz4wsxwzGx/plUdbodf6gtQS/cwsBZgFPOCcO+p1nvrmnAs65/oTvi7vJWYWs1NsZnYdkOecy/E6SwMb5pwbCIwB7quaUo2YaCt0XZC6iaiaR54FvOyce93rPA3JOXcEmAeM9jhKfRoGjKuaU54OXGVmU72NVP+cc/urvucBswlPI0dMtBX6CqCbmXU1swTge4QvUi0xpOoDwheATc65p7zO0xDMLNPMWlXdbgZcA2z2NlX9cc79h3Oug3OuC+F/x5845271OFa9MrPkqg/5MbNkYBQQ0b3XoqrQnXMB4MfA+4Q/KJvhnNvgbar6VXWR7qVADzPLNbN/9TpTAxgG3EZ41La66utar0PVs/OBuWa2lvDA5UPnXJPYla8JaQssMrM1wKfAP5xz70VyA1G126KIiJxeVI3QRUTk9FToIiIxQoUuIhIjVOgiIjFChS4iEiNU6CIiMUKFLiISI/4/btPQ7rBjcwEAAAAASUVORK5CYII=\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"S8HLShVYJJqW"},"source":["## 時間の刻みを変えてみる\n","近似解と厳密解にちょっとズレが見えます。\n","そこで `nt = 101` から `nt = 1001` にしてみましょう。"]},{"cell_type":"code","metadata":{"scrolled":true,"id":"rbme9MuAJJqX","outputId":"178c9818-9f9c-4edc-dfe0-bbf975c7393b"},"source":["import numpy as np\n","import matplotlib.pyplot as plt\n","\n","def logistics_euler(nt, init = 10):\n","    dt = T / (nt - 1)\n","    # 初期条件設定\n","    x = np.zeros(nt)\n","    x[0] = init\n","\n","    for i in range(1, nt):\n","        x[i] = x[i-1] + dt * (alpha - beta * x[i-1]) * x[i-1]\n","\n","    return x\n","\n","# パラメータ設定\n","alpha = 2\n","beta = 1\n","init = 1\n","nt = 1001\n","T = 5 # 時間変化を見る最大値\n","\n","# 近次解\n","x_approx = logistics_euler(nt, init)\n","plt.plot(np.linspace(0, T, nt), x_approx)\n","\n","# 厳密解\n","t = np.linspace(0, T, nt)\n","x_exact = (alpha / beta) / (1 + np.exp(- alpha * t))\n","plt.plot(t, x_exact)\n","\n","# 凡例\n","plt.legend(['approximation', 'exact'])\n","# 描画\n","plt.show()"],"execution_count":null,"outputs":[{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU9b3/8ddnshCyQAIJspOAQIASECO4IOCGoIKK9lpvW8VWqcttsff2d7U+6rV2+bW/LrZF28ulGtGr4o5Ui7YuKCprosgaFgmEEIQQloQsJJP5/v5IRJSEJDDJycy8n4/HONuZc96TyJvD92zmnENEREKfz+sAIiISHCp0EZEwoUIXEQkTKnQRkTChQhcRCRPRXi04NTXVpaene7V4EZGQlJeXt985l9bYe54Venp6Orm5uV4tXkQkJJnZzqbe05CLiEiYUKGLiIQJFbqISJhQoYuIhAkVuohImGi20M2sn5ktMbNNZrbBzGY3Mo2Z2Rwz22Zma81sTNvEFRGRprRkt0U/8B/OuY/MLAnIM7M3nXMbj5tmKjC44TYO+O+GexERaSfNFrpzbg+wp+FxuZltAvoAxxf61cCTrv5cvCvMLNnMejV8VkQ84Jyjrq4Ov7+Wutoa/LW11PlrCPhr8dfVUuevo67OT6DOT8A5XKCu4RbAuYb7QB3O1UEgQKCuDtznrwXg2LQNj93nj/045xpec/WfcQ4XcARwgKP+LoDh+OIM3o76J/X3DjDncDjs89cbvlf957547fjPcWx6js0DAl9ajh2/rIZ5HveD++Lhl/7TyKnG3QkPTnzuqF/ecRIGj2fkxBmN/dpOS6sOLDKzdOAsYOVX3uoD7DrueVHDa18qdDObBcwC6N+/f+uSioSQWn8dlVVVVFeUUV1RRk1lGUcry/FXl+OvOoL/aAV1tVW4mmqoO4rzV2P+o5i/Gqs7itUdxVdXQ1RdNVGBo0QHaogK1BDlavAF/ERRR5SrIwo/Uc5PFAGi8RPl6ojm85ufaHPeHT0oXxJwduzxykAteFnoZpYIvATc7Zwr++rbjXzkhL/OnHPzgHkA2dnZurKGdFjOOcrKj3BofzHlB/ZRXV5K7ZED1FUexFUdxKoP4Tt6mJiaMmJry4irK6dToJI4V01nV0081XS1Orq2crnVLoYai6WGGGot9tjN74vFb7HURHXB+aIJ+GJwFoXzReMsuv7eF9NwHw0Nz/FFY1FfPLeoGCwqGmt43eeLAvNhUVGY+eoff/7acff2pefHTe+Lwnf8+z4fPvPVz9/nw8zw+Qyz+tcxw4yGe8MwsPrb54+tfoL6ewO++jpgPt+xaepfB6z+NXxfvNbUtNYwrfmMz+vLfA0ZaIh07HHDo8+zH/PF+8c+dBLHb7A8rzX/U7RCiwrdzGKoL/OnnXMvNzJJEdDvuOd9geLTjycSXH5/Hfv2l3Bwz3bK9+2g5kARrnwvUVX7ia3eT3ztAbr4D5LiDtHVqpos5FoXxRFLoMKXRGVUEjWxyVTF9CEQnUAgJgEXmwCxCVhsIr5OiUTFJRIdl0R0fBKx8UnEdU4iNi6eTnGdiYnrTGyneCwqljgz4tr1JyLhpNlCt/q/fh4DNjnnHmpisr8B/2Zmz1K/MfSwxs/FC8459h04wN6CjZTvzsdf8im+siLiq/eQXLuXHoH99LYqen/lc2UkcjgqhYroFA7GZ1ISn4ZLSCMq6Qw6dUklLrE7cV26kZCcSkLXVGLikkgxI8WTbynSuJasoV8AfBtYZ2ZrGl67D+gP4JybCywGrgC2AZXALcGPKvIFf12AwqJd7N2aS3XRWqIPbKNL5U7OqN1NTzvAGcdNe4guHIzpwZHEdLYmno917Uds9/4knZFOt57pJHbrRZfoWLp49m1EgqMle7l8QONj5MdP44C7ghVK5HjVNX62bl5P6eYPcXvW0aVsM/1qtjPQDjGwYZpDdGF/p7581vVc9qQMJK7nEFL6DSdtwDCS4xJJ9vQbiLQPbQCXDqdoz14KPlnK0R0r6Vq6hkE1+Yy0cgCOEkNxTDp70i5gT48RdM0YTc/BZ5Oc3FOlLRFPhS6eK95/gG2571Cz7V16lq4kM7CNvhYg4IzimH4U9ZjEnv5j6Tn8ArqnjyIjSv/bijRGfzKk3fnrAqxfv4aSvEV02/0OX/NvorfVUksUhXGZbOj9XVKGTaLPiPH0jU+mr9eBRUKECl3aRXWNn7zlb1H5ySIyDixlNEUAFEUPYPOAb5A8/FL6jrqYQZ21aVLkVKnQpc3U1gX4KHc5h1ctIHP/P7nA9uInioKEUWwcdBP9z59B356DtQYuEiQqdAm67YW72Prmo6TvWsg4dlKH8WliNluz/p2M8TcwOEF7b4u0BRW6BEV1jZ8VS1/Hcuczruo9BlotOzsNZfOwn5Ax6VsMSe7ldUSRsKdCl9NysLySVX9/jAH5jzGJAirozKd9rqbXJXcwYFC21/FEIooKXU5J0WclrH31EUYVPc3lVkJxdD+2jvkFZ14ykxGdkryOJxKRVOjSKsX7D7Dmpd9ybvGTXGFH2J6QRdGE39B37DXg0xUNRbykQpcWKTlUzqqX/kB2YQ5X2EG2dhmHu+J+Bg670OtoItJAhS4ndbTWz5JF8xmx7jdcaXspSBjFviseY/DXLvE6moh8hQpdGuWcY9mKZcS++WOmBD5hd2w6e6Y8TcaYK5s9kb+IeEOFLifYU7KfT/73Hi45vJBqi2Pr2f/F4Ctmg86hItKh6U+oHBMION5a/AKZq3/CFNtLfp9rGfSN3zC4Sw+vo4lIC6jQBYBde/aS/+TdTK5azGcxvdk3/SUysy71OpaItIIKXXj3ndcZ9N5sLrZ9bB54M0O+8SssNsHrWCLSSir0CFZedZR3H7+fKXsf5VBUd0qve4WhIyZ5HUtETpEKPUJtK9jBwf+9iWmBT9iaegkZtzxKdGI3r2OJyGlQoUegD95/h4y3ZtHXDlNw3i8ZPPku7YooEgZU6BGkLuD4+zMPc9nWn1MR1YWKG18jY/A4r2OJSJCo0CNE1dFalvz395l+aAEFCVn0nvU8nXRKW5GwokKPAKWHy1n3l29zxdElbO5zHUNvmQvRsV7HEpEgU6GHucI9n1Hy168zKbCWLSN+yNDrH9B4uUiYUqGHsW0FBdQ9cQ1Z7GLH+N8y5NJZXkcSkTakQg9T+Vu3EPv0NfSlhJKr5pOePd3rSCLSxlToYWhDfj4Jz15DDw5y8NoF9B6lQ/hFIoEKPcxs3LSRpGevobuVUX79c/T62iSvI4lIO1Ghh5FtBQXEPzeDblZO1Tde5IzM8V5HEpF2pItAhonC4mLqnriWnpRy5LoFpKrMRSKOCj0MfLb/AAcfncFACim98jF6jpzkdSQR8UCzhW5mOWa2z8zWN/F+VzN71cw+MbMNZnZL8GNKU8orqyic+3VG1uVTfPEc+pyjvVlEIlVL1tDnA1NO8v5dwEbn3ChgEvB7M9NhiO3A769j5V9uY6w/l0/HPsiACd/yOpKIeKjZQnfOLQUOnGwSIMnMDEhsmNYfnHjSFOccb+T8lEuPvMqmjJkMvnK215FExGPBGEN/BBgGFAPrgNnOuUBjE5rZLDPLNbPckpKSICw6cv1j4eNcsfthtqRMZNi3/+B1HBHpAIJR6JcDa4DewGjgETPr0tiEzrl5zrls51x2WlpaEBYdmXJXfcCFn9zLrrjBnPm9Z8CnbdsiEpxCvwV42dXbBhQAmUGYrzSiaM8eeiz+Lkd98fSY9TK+uESvI4lIBxGMQi8ELgEwszOAocD2IMxXvqLqaC1FOTfRixJqZuTQuXs/ryOJSAfS7JGiZraA+r1XUs2sCHgAiAFwzs0Ffg7MN7N1gAH3OOf2t1niCOWc451H7+XK2lVsGfMThoy82OtIItLBNFvozrkbm3m/GJgctETSqHcXP8fUfY+R3+NyMqf9yOs4ItIBaWtaCCjYUcDXVv0nu2P6M+S7ObpAhYg0Sifn6uCqa/yUPHUrva0S378u0kZQEWmS1tA7uCVP/pyx/lx2nP1jug88y+s4ItKBqdA7sJUr3uPiXY+wucsFDL3q372OIyIdnAq9gzpUVkbqG3dR4UtiwHce17i5iDRLhd5BrX7iXgaxi/IpfyIu+Qyv44hICFChd0ArP3iLi/YvYH2P6QwYd7XXcUQkRKjQO5hDZeV0f+uHHPIlM+SmOV7HEZEQokLvYFY9cR9nUkjF5b8nNjHF6zgiEkJU6B3IR6uWcvH+p9iQdiUDzp3hdRwRCTEq9A6iuqaWTm/8iHJfEoO+raEWEWk9FXoH8f5zDzEisJl9595PXJdUr+OISAhSoXcAO3cVcs62P7G18yiGTr7V6zgiEqJU6B5zzrFzwX+QYNV0u+FhHUAkIqdMhe6x5UteY0LlP8nPuJnu6aO8jiMiIUyF7qHqozWkvn8/e31pDP/GL7yOIyIhToXuoWUvzWGIK+DQBT8hqlOC13FEJMSp0D1Ssn8/WZvnsC1uBEMvvtnrOCISBlToHln77AOk2mESpv9WG0JFJChU6B7I37Se8SXPsb77FHoNv8DrOCISJlTo7cw5x8FFP8aZMeCG33gdR0TCiAq9nX284h3Oq15KfsZMknoM8DqOiIQRFXo7CgQc9s7POEQSw66/z+s4IhJmVOjtaOU7Czmrdg2FI26nU4JOjSsiwaVCbye1/jqSl/1fSiyVEVf/h9dxRCQMqdDbyfLFTzIssJWSs+8mKraz13FEJAyp0NtB9dEa+n78e3ZH9WHY1Nu9jiMiYUqF3g5Wv/YoA90ujpx/DxYV43UcEQlTKvQ2drSmhr7r/0xhdDpDLvqW13FEJIyp0NvYyr/PJ8MVUXHuDzFflNdxRCSMNVvoZpZjZvvMbP1JpplkZmvMbIOZvRfciKGrptZPn7UPsyuqH5laOxeRNtaSNfT5wJSm3jSzZOAvwHTn3Ajg68GJFvpWvv4kg1wh5WNnY1HRXscRkTDXbKE755YCB04yyb8CLzvnChum3xekbCGt1l/HGR/PYbevN8Munel1HBGJAMEYQx8CpJjZu2aWZ2Y3NTWhmc0ys1wzyy0pKQnCojuuFf94pv7iFdk/0J4tItIuglHo0cDZwJXA5cD9ZjaksQmdc/Occ9nOuey0tLQgLLpjCtQFSMv7E3t8PRk++btexxGRCBGMQi8C3nDOVTjn9gNLgYi+2nHe0kVkBrayL+sOLDrW6zgiEiGCUeiLgAvNLNrM4oFxwKYgzDdkRS1/hAN0ZcTUWV5HEZEI0uyuF2a2AJgEpJpZEfAAEAPgnJvrnNtkZm8Aa4EA8KhzrsldHMPdho+XMaYml4/O/DfGdIr3Oo6IRJBmC905d2MLpvkt8NugJApxh99+iEo6kTntbq+jiEiE0ZGiQVSwfQvnlL9Dfq9riO8avht9RaRjUqEHUeHi32M4Mq76P15HEZEIpEIPkpL9+zi7ZBEbUy4ipc9gr+OISARSoQfJplcfJtGq6D5Za+ci4g0VehAcrTnK4J3PkB83ij7Dz/M6johEKBV6EHz85gJ6sZ/AWF2NSES8o0I/Tc454j9+lM+sB8Mm/ovXcUQkgqnQT9OmNcvJ8q9j95Bv6hS5IuIpFfppOvjun6kilsypd3kdRUQinAr9NOz9rJgxh/7JptSpJCTrQCIR8ZYK/TTkL/4zna2GXpfN9jqKiIgK/VQdrTnKkMJnyY8bRa+hZ3sdR0REhX6qPnn7eXqxn7pzdIpcEekYVOinKGbNE5RYN4ZNvMHrKCIigAr9lOzYtolR1bns6DcDX7SuFyoiHYMK/RTsensuABmX3+FxEhGRL6jQW6m6upqhexaxMWEsqX3O9DqOiMgxKvRWWvP2s/TgIL5zvuN1FBGRL1Ght1LcJ09SYt3JvPA6r6OIiHyJCr0Vtm/dQNbRj9g54DptDBWRDkeF3gq735qLAwZdfqfXUURETqBCb6Hq6mqG7V3ExsTzSOmV4XUcEZETqNBbaN2S50jlML7sW7yOIiLSKBV6C0WtXcB+Uhh24bVeRxERaZQKvQX2FheSVbmSgt5XaWOoiHRYKvQW2PbWY0RbgD4X3ep1FBGRJqnQm+ECAXoXvMzWmEx6Dx7tdRwRkSap0Jux8aOlZLhCyjJ1VkUR6dhU6M04vOxxql0Mwy6b6XUUEZGTUqGfRGXlEUaUvsmm5InEd+nmdRwRkZNSoZ/E2refpatVED/uJq+jiIg0q9lCN7McM9tnZuubme4cM6szs+uDF89bndYvYK+lMmTclV5HERFpVkvW0OcDU042gZlFAf8P+EcQMnUIxYWfklWdR2G/q7GoaK/jiIg0q9lCd84tBQ40M9n3gZeAfcEI1REUvJNDlDn6X3yb11FERFrktMfQzawPcC0wtwXTzjKzXDPLLSkpOd1FtxnnHL0KX2VzzHDOSB/mdRwRkRYJxkbRPwL3OOfqmpvQOTfPOZftnMtOS0sLwqLbxra1KxkY2En5EJ23RURCRzAGh7OBZ80MIBW4wsz8zrlXgjBvT5Qse5J0F8WQi7V3i4iEjtMudOfcsZODm9l84LVQLnO/38+gvW+wKeEcsrr39DqOiEiLNVvoZrYAmASkmlkR8AAQA+Cca3bcPNRsWP46oyhlz8j7vI4iItIqzRa6c+7Gls7MOTfztNJ0AFV5C6ggjsyJOneLiIQWHSl6nKrKCoYfXEJ+8iTi4pO8jiMi0ioq9OOsf/cFulgl8dkt/keJiEiHoUI/jm/9C5SSzNBzdai/iIQeFXqDA/v3MrJiBdt7Xq7LzIlISFKhN8hf8jSx5ift/G97HUVE5JSo0Bt02bKQIl8f0keO9zqKiMgpUaEDnxVtZ3jNOnb3uwrqj3gVEQk5KnSg4L2n8Zmj74Xf9DqKiMgpU6EDKQWvsT0qgz5njvI6iojIKYv4Qt9dsJlMfz4lA7SrooiEtogv9J3vPw1Afw23iEiIi/hCT925mG3RZ9IrY7jXUURETktEF/quTzcypG4rpelXeR1FROS0RXahNwy3ZEzUcIuIhL6ILvQehYvZEpNJj35DvI4iInLaIrbQd2z+hDMD2zk8UMMtIhIeIrbQd3/4DAAZEzTcIiLhISIL3TlHr6LXyY8dQWqfgV7HEREJiogs9O2b8hgY2En5oGleRxERCZqILPTPli0g4IxB2rtFRMJIxBW6c44+u98gP24k3Xr29zqOiEjQRFyh78z/iHRXRLn2bhGRMBNxhb5n+fMADJpwg8dJRESCK+IKvUfRP8iPGU5qr3Svo4iIBFVEFXrRpxsYFCjgcPpUr6OIiARdRBX6rg+fBaD/eA23iEj4iahC7174Blujz6TXgKFeRxERCbqIKfTPdm1jiH8Lpf0u9zqKiEibiJhCL3j/OQD6nP8Nj5OIiLSNiCn0rgWLKfANoN/gLK+jiIi0iWYL3cxyzGyfma1v4v1vmtnahtsyMxsV/Jinp+SzQjJrNrC3z2VeRxERaTMtWUOfD0w5yfsFwETnXBbwc2BeEHIF1adLn8dnjp7n/ovXUURE2kx0cxM455aaWfpJ3l923NMVQN/TjxVcCZ/+nSLrxYBh53gdRUSkzQR7DP27wOtNvWlms8ws18xyS0pKgrzoxh3av5fM6k8o6nUZ5ouYTQYiEoGC1nBmdhH1hX5PU9M45+Y557Kdc9lpaWnBWvRJbX7/BWKsjtRzvt4uyxMR8UqzQy4tYWZZwKPAVOdcaTDmGSyxW17lM0tj0KjxXkcREWlTp72Gbmb9gZeBbzvntpx+pOApO3yA4ZV57Ey7WMMtIhL2ml1DN7MFwCQg1cyKgAeAGADn3Fzgv4DuwF/MDMDvnMtuq8CtseX9l8i2WrpmX+d1FBGRNteSvVxubOb9W4Fbg5YoiHz5r1JKMkPGXOJ1FBGRNhe24xBHqysYWr6Sbd0m4osOyqYCEZEOLWwLPX/530mwauJGTvM6iohIuwjbQq9e9yoVLo7M8670OoqISLsIy0IP1NUx6MBSNieNo1NcvNdxRETaRVgW+paP3yWVQwSGau1cRCJHWBb6wbxXqHVRDB6v3RVFJHKE5e4fffa+TX5cFiNTUr2OItKh1NbWUlRURHV1tddRpBlxcXH07duXmJiYFn8m7Aq9cMsa+gd2s3LgTV5HEelwioqKSEpKIj09nYYDAaUDcs5RWlpKUVERGRkZLf5c2A25FK94CYAB51/vcRKRjqe6upru3burzDs4M6N79+6t/pdU2BV6cuGbbI06k579zvQ6ikiHpDIPDafyewqrQt+/ZydDavPZ3/dSr6OIiLS7sCr07R++WH+pubEzvI4iIh3Y3/72N379618HZV5//OMfqaysPPb8iiuu4NChQ0GZd2uFVaF3+vQNdtsZpOtScyIRw+/3t/oz06dP59577w3K8r9a6IsXLyY5OTko826tsNnL5UjZQYZVfkRez6/TR+c+F2nWg69uYGNxWVDnObx3Fx6YNuKk01xzzTXs2rWL6upqZs+ezaxZs0hMTOR73/seS5YsISUlhWeffZa0tDQmTZrE6NGjWbVqFWVlZeTk5DB27Fh++tOfUlxczI4dO0hNTSUnJ4c77riD3NxcoqOjeeihh7jooot46KGHWL9+PTk5Oaxbt44bb7yRVatW8fzzz5Obm8sjjzzCzJkz6dy5M/n5+ezcuZPHH3+cJ554guXLlzNu3Djmz58PwB133MHq1aupqqri+uuv58EHH2TOnDkUFxdz0UUXkZqaypIlS0hPTyc3N5fU1FQeeughcnJyALj11lu5++672bFjB1OnTmX8+PEsW7aMPn36sGjRIjp37nzaP/+wab4tH75CrPnpOvpqr6OIyEnk5OSQl5dHbm4uc+bMobS0lIqKCsaMGcNHH33ExIkTefDBB49NX1FRwbJly/jLX/7Cd77znWOv5+XlsWjRIp555hn+/Oc/A7Bu3ToWLFjAzTffTHV1NXfffTfbtm1j4cKF3HLLLfzP//wP8fEnng7k4MGDvPPOO/zhD39g2rRp/PCHP2TDhg2sW7eONWvWAPDLX/6S3Nxc1q5dy3vvvcfatWv5wQ9+QO/evVmyZAlLliz50jzz8vJ4/PHHWblyJStWrOCvf/0rH3/8MQBbt27lrrvuYsOGDSQnJ/PSSy8F5WcbNmvogU2vcZAuDMnWuc9FWqK5Nem2MmfOHBYuXAjArl272Lp1Kz6fjxtuuAGAb33rW8yY8cV2sBtvrL8kw4QJEygrKzs2Pj19+vRja7UffPAB3//+9wHIzMxkwIABbNmyhaysLObPn09WVhbf+973uOCCCxrNNG3aNMyMkSNHcsYZZzBy5EgARowYwY4dOxg9ejTPP/888+bNw+/3s2fPHjZu3EhWVlaT3/ODDz7g2muvJSEhAYAZM2bw/vvvM336dDIyMhg9ejQAZ599Njt27Diln+VXhUWh19YcZUjZcvKTJzI2JtbrOCLShHfffZe33nqL5cuXEx8fz6RJkxrd1/r4Xfa+uvve588/L0qoPxCnKVu3biUxMZHi4uImp+nUqRMAPp/v2OPPn/v9fgoKCvjd737H6tWrSUlJYebMmc3uI36yTMcvIyoqiqqqqpPOq6XCYshl88o36EIFMcN1Mi6Rjuzw4cOkpKQQHx9Pfn4+K1asACAQCPDiiy8C8MwzzzB+/BcXdX/uueeA+jXerl270rVr1xPmO2HCBJ5++mkAtmzZQmFhIUOHDuXw4cPMnj2bpUuXUlpaemwZrVVWVkZCQgJdu3Zl7969vP7668feS0pKory8vNFMr7zyCpWVlVRUVLBw4UIuvPDCU1p+S4XFGnrFJ4uocrFkXqDxc5GObMqUKcydO5esrCyGDh3KueeeC9SvbW/YsIGzzz6brl27HitxgJSUFM4///xjG0Ubc+edd3L77bczcuRIoqOjmT9/Pp06deKOO+7gzjvvZMiQITz22GNcdNFFTJgwodW5R40axVlnncWIESMYOHDgl4ZuZs2axdSpU+nVq9eXxtHHjBnDzJkzGTt2LFC/UfSss84K2vBKY+xk/yxoS9nZ2S43N/e05+MCAfb+bDB74ody1n8uDkIykfC1adMmhg0b5nWMEyQmJnLkyJETXp80aRK/+93vyM7uENedb3eN/b7MLM851+gPJOSHXD5dt4ye7Kd28FSvo4iIeCrkh1z2r36ZDGcMHq+TcYmEqsbWzqF+I6q0XMivofcofpvNnb5GSlovr6OIiHgqpAt99/ZNDAzsoGzAZK+jiIh4LqQLfdfyFwDod97XPU4iIuK9kC70pJ3/pMCXTp+BHW+rvYhIewvZQj9YsofMo+vZ21uH+otI/QbUZcuWeR3DUyFb6Ns+fJEoc3TP1rnPRUSFDiG822L0ltf5jFTOzDrf6ygioen1e+GzdcGdZ8+RMPXkF4546qmnmDNnDjU1NYwbN4777ruPSy+9lOXLl9OtWzcmTpzI/fffz+TJkxs91S7AG2+8wX333UddXR2pqak89thjzJ07l6ioKJ566ikefvjhNj/MviMKyUKvqigns2I1a9Om0VPnPhcJGZs2beK5557jww8/JCYmhjvvvJP33nuPe+65h9tvv51x48YxfPhwJk+u33MtJyeHbt26UVVVxTnnnMN1111HIBDgtttuY+nSpWRkZHDgwAG6devG7bffTmJiIj/60Y88/pbeCclC37zsVUZbDfFZOneLyClrZk26Lbz99tvk5eVxzjn1VxWrqqqiR48e/PSnP+WFF15g7ty5x84/Do2farekpIQJEyaQkZEBQLdu3dr9e3RUzRa6meUAVwH7nHNfa+R9A/4EXAFUAjOdcx8FO+jxaja8ShnxZJ47pS0XIyJB5pzj5ptv5le/+tWXXq+srKSoqAioP2o0KSmpyVPtOudOOKWu1GvJeMV84GTNORUY3HCbBfz36cdqWp3fz+BDH7Cly/nExHZq/gMi0mFccsklvPjii+zbtw+AAwcOsHPnTu655x6++c1v8rOf/YzbbrsNaPpUu+eddx7vvfceBQUFx+YBTZ/GNpI0W+jOuaXAgZNMcjXwpKu3Akg2szY7Dn9z7lukUIZlXtVWixCRNjJ8+HB+8YtfMHnyZLKysrjsssvYsWMHq1evPlbqsbGxPP74456f2EIAAARSSURBVEyZMgW/309WVhb333//sVPtpqWlMW/ePGbMmMGoUaOOXelo2rRpLFy4kNGjR/P+++97+TU906LT55pZOvBaE0MurwG/ds590PD8beAe59wJ58Y1s1nUr8XTv3//s3fu3NnqwPkr/0nNkt8w8M4XSOyS0urPi0Syjnr6XGlca0+fG4yNoo0NZjX6t4Rzbh4wD+rPh34qC8scNxnG6dwtIiJfFYx9/oqAfsc97ws0ffE+ERFpE8Eo9L8BN1m9c4HDzrk9QZiviLQBr65SJq1zKr+nluy2uACYBKSaWRHwABDTsMC5wGLqd1ncRv1ui7e0OoWItIu4uDhKS0vp3r27dv3rwJxzlJaWEhcX16rPNVvozrkbm3nfAXe1aqki4om+fftSVFRESUmJ11GkGXFxcfTt27dVnwnJI0VF5NTExMQcO8JSwo9OhCIiEiZU6CIiYUKFLiISJlp0pGibLNisBGj9oaL1UoH9QYwTCvSdI4O+c2Q4ne88wDmX1tgbnhX66TCz3KYOfQ1X+s6RQd85MrTVd9aQi4hImFChi4iEiVAt9HleB/CAvnNk0HeODG3ynUNyDF1ERE4UqmvoIiLyFSp0EZEwEXKFbmZTzGyzmW0zs3u9ztPWzCzHzPaZ2Xqvs7QXM+tnZkvMbJOZbTCz2V5namtmFmdmq8zsk4bv/KDXmdqDmUWZ2ccNVz4Le2a2w8zWmdkaMzvhqm6nPf9QGkM3syhgC3AZ9RfWWA3c6Jzb6GmwNmRmE4Aj1F+39YRLAIajhmvS9nLOfWRmSUAecE2Y/54NSHDOHTGzGOADYHbDdXrDlpn9O5ANdHHOhf2Fgs1sB5DtnGuTA6lCbQ19LLDNObfdOVcDPEv9RarDVgsu0h12nHN7nHMfNTwuBzYBfbxN1bYaLrJ+pOFpTMMtdNa2ToGZ9QWuBB71Oku4CLVC7wPsOu55EWH+Bz3SNVyg/CxgpbdJ2l7D8MMaYB/wpnMu3L/zH4H/BAJeB2lHDvinmeWZ2axgzzzUCr3FF6SW0GdmicBLwN3OuTKv87Q151ydc2409dflHWtmYTvEZmZXAfucc3leZ2lnFzjnxgBTgbsahlSDJtQKXRekjhAN48gvAU875172Ok97cs4dAt4FpngcpS1dAExvGFN+FrjYzJ7yNlLbc84VN9zvAxZSP4wcNKFW6KuBwWaWYWaxwDeov0i1hJGGDYSPAZuccw95nac9mFmamSU3PO4MXArke5uq7Tjnfuyc6+ucS6f+z/E7zrlveRyrTZlZQsNGfswsAZgMBHXvtZAqdOecH/g34B/Ubyh73jm3wdtUbavhIt3LgaFmVmRm3/U6Uzu4APg29WttaxpuV3gdqo31ApaY2VrqV1zedM5FxK58EeQM4AMz+wRYBfzdOfdGMBcQUrstiohI00JqDV1ERJqmQhcRCRMqdBGRMKFCFxEJEyp0EZEwoUIXEQkTKnQRkTDx/wGtyvQbLJkujgAAAABJRU5ErkJggg==\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"uVYJ89GEJJqZ"},"source":["こちらも時間の刻みを細かくすると近似がよくなりました。"]},{"cell_type":"markdown","metadata":{"id":"SZojZaIiJJqZ"},"source":["## 「近似の精度が気になります」\n","これもよく近似できていることは認めてもらえるのではないでしょうか?\n","しかし気になる人は「近似の精度はどう決めるのか」と考えているはずです。\n","\n","もちろんいたって真っ当な指摘です。\n","ただ、かなり難しい問題がたくさんあり、\n","いまの時点でなかなかスカっと綺麗に回答できません。\n","あとで少し考えることにしてここではこのまま進みます。"]},{"cell_type":"markdown","metadata":{"id":"RpIIKdRAJJqZ"},"source":["## 生物学として正しいの?\n","あなたは数値計算の精度とかそれ以前の問題として、\n","次のように思っているかもしれません。\n","\n","- **この微分方程式で実際の生物の増減にあてはめられるのか？**\n","- **もっと現実的なことを考える必要があるのではないか?**\n","\n","もちろんこれもいたって真っ当な指摘です。\n","\n","例えばマルサスモデルでもロジスティック方程式でも右辺が気になります。\n","実際には成長しないと子作りできません。\n","その成長が必要だという情報が方程式に盛り込まれていません。\n","\n","生物で考えるなら食物連鎖があるわけで、\n","食べられて個体数が減ることだってあります。\n","天敵が増えたらその個体の数はあおりを受けて激減するはずです。\n","\n","もちろんこんなことは折り込み済みで、\n","例えば**遅延型微分方程式**、\n","**ロトカ-ボルテラ方程式**と呼ばれる方程式が対応しています。\n","これを調べるのも意味があることで面白いことでもありますが、\n","今回はこのくらいにしておきましょう。"]},{"cell_type":"markdown","metadata":{"id":"jh6HQ6fiJJqa"},"source":["## 今回のポイント\n","簡単にまとめると次の通りです。\n","\n","- 実際にどのくらいよく本当の解を近似できているのかを調べた。\n","- 経済学や生物のように、 あまり数学との関係がなさそうな分野と数学の関係を紹介した。"]},{"cell_type":"markdown","metadata":{"id":"wotjsGh6JJqa"},"source":["## 数学的ポイント: 数列と漸化式\n","長くなってきましたが最後にちょっと大事な話。\n","\n","式 (1) で $\\alpha = 1$、 $h = 1$ とすると $x_{n+1} = 2 x_n$ になります。\n","この漸化式は高校でも出てくるやつで、 もちろん公比 $2$ の等比数列の漸化式です。\n","底が $2$ か $e$ かの違いはあっても指数関数で書けることは同じです。\n","これはたまたまではなくて数学的に意味があることです。\n","\n","$x_{n+1} - x_{n}$ のように数列の隣の項の差を差分と言うことがあります。\n","もちろん微分との関係を強く意識した言葉です。\n","数列の問題、 もっと言えば漸化式の問題は微分方程式とも深い関係があります。\n","\n","広く力学系と呼ばれる分野にまで広がります。\n","力学系まで行くととんでもなく難しくなるのでとてもここで紹介はできません。\n","でも幾何学とも関係してきたり整数論とも関係してきたり、\n","射程範囲は広いです。\n","\n","ここで必要以上にあまり深入りはしませんが、\n","数学や物理に興味があるなら覚えておくといろいろ楽しいことがあります。"]},{"cell_type":"markdown","metadata":{"id":"QlMaLdlDJJqb"},"source":["# アンケート\n","毎回アンケートを取っています。\n","質問や要望がある場合もこちらにどうぞ。\n","\n","- [アンケートへのリンク](https://goo.gl/forms/hn7bUP4sblqOkBcI3)\n","\n","アンケートは匿名なので気楽にコメントしてください。\n","直接返事してほしいことがあれば、\n","メールなど適当な手段で連絡してください。"]}]}