Last active
August 29, 2015 14:21
-
-
Save kenjisato/08439241d90756a5248c to your computer and use it in GitHub Desktop.
2015 経済動学 講義資料 (2015/05/25)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "# ラムゼーモデルの局所動学\n", | |
| "\n", | |
| "### 経済動学 講義資料\n", | |
| "\n", | |
| "神戸大学 佐藤健治\n", | |
| "\n", | |
| "2015-5-25\n", | |
| "\n", | |
| "https://bitbucket.org/kenjisato/economic_dynamics\n", | |
| "\n", | |
| "See https://bitbucket.org/kenjisato/economic_dynamics/src/0230a3f1c5b7075b77e9a23611e56d743bbbc903/Codes/2015-05-25/?at=master" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": { | |
| "collapsed": true, | |
| "slideshow": { | |
| "slide_type": "skip" | |
| } | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "%matplotlib inline" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "## 目標\n", | |
| "\n", | |
| "1部門モデルを考える. \n", | |
| "\n", | |
| "生産関数と効用関数を特定化して、オイラー条件\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| " u_2(k_{t-1}, k_t) + \\rho u_1(k_t, k_{t+1}) = 0\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "によって定まる力学系を可視化したい\n", | |
| "\n", | |
| "3つのプログラミング・パラダイム\n", | |
| "\n", | |
| "1. 手続き型\n", | |
| "1. オブジェクト指向\n", | |
| "1. [Optional] 関数型" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "## 手続き型 (素朴な方法)\n", | |
| "\n", | |
| "### パラメータ" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "## Unicode identifiers are only valid in Python3\n", | |
| "## In Python2, use less readable notation such as\n", | |
| "# alpha = 0.4\n", | |
| "# rho = 0.9\n", | |
| "\n", | |
| "A = 1.1 \n", | |
| "α = 0.4\n", | |
| "ρ = 0.9 " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "collapsed": true, | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "### 生産関数 \n", | |
| "\\begin{align}\n", | |
| " f(k) = Ak^\\alpha\n", | |
| "\\end{align}" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "def f(k):\n", | |
| " return A * k ** α" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "### 消費の効用関数\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| " U(c) = \\log c\n", | |
| "\\end{align}" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "from math import log\n", | |
| "\n", | |
| "def U(c):\n", | |
| " return log(c)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "### 規約型効用関数\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| " u(x, y) = U(f(x) - y)\n", | |
| "\\end{align}" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "def u(x, y):\n", | |
| " return U(f(x) - y)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "### オイラー条件\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| " u_2(k_{t-1}, k_t) = \\frac{-1}{Ak_{t-1}^\\alpha - k_t}\\\\\n", | |
| " u_1(k_t, k_{t+1}) = \\frac{\\alpha A k_t^{\\alpha - 1}}{Ak_t^\\alpha - k_{t+1}}\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "より\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| " k_{t+1} = F(k_{t-1}, k_t) = (1+\\rho \\alpha) A k_t^\\alpha - \\rho \\alpha A^2 k_{t-1}^\\alpha k_t^{\\alpha-1}\n", | |
| "\\end{align}" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 6, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "def F(x, y):\n", | |
| " return ((1 + ρ * α) * A * y ** α \n", | |
| " - ρ * α * (A ** 2) * (x ** α) * (y ** (α - 1)))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "### 相空間上の力学系" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 7, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "def G(x):\n", | |
| " return [\n", | |
| " x[1],\n", | |
| " F(x[0], x[1])\n", | |
| " ]" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "### シミュレーション" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 8, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "[0.8, 0.43]\n", | |
| "[0.43, 0.4063168783383537]\n", | |
| "[0.4063168783383537, 0.509937247085211]\n", | |
| "[0.509937247085211, 0.6875981937381642]\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "duration = 4\n", | |
| "x0 = [0.8, 0.43]\n", | |
| "\n", | |
| "x = x0[:]\n", | |
| "for t in range(duration):\n", | |
| " print(x)\n", | |
| " x = G(x)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "以上をまとめたものが listing01.py" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 9, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "#!/usr/bin/env python3\r\n", | |
| "\r\n", | |
| "from math import log\r\n", | |
| "\r\n", | |
| "## Unicode identifiers are only valid in Python3\r\n", | |
| "## In Python2, use less readable notation such as\r\n", | |
| "# alpha = 0.4\r\n", | |
| "# rho = 0.9\r\n", | |
| "\r\n", | |
| "A = 1.1\r\n", | |
| "α = 0.4\r\n", | |
| "ρ = 0.9\r\n", | |
| "\r\n", | |
| "\r\n", | |
| "def f(k):\r\n", | |
| " \"\"\"Production function\"\"\"\r\n", | |
| " return A * k ** α\r\n", | |
| "\r\n", | |
| "\r\n", | |
| "def U(c):\r\n", | |
| " \"\"\"Utility function\"\"\"\r\n", | |
| " return log(c)\r\n", | |
| "\r\n", | |
| "def u(x, y):\r\n", | |
| " \"\"\"Reduced form utility function\"\"\"\r\n", | |
| " return U(f(x) - y)\r\n", | |
| "\r\n", | |
| "\r\n", | |
| "def F(x, y):\r\n", | |
| " \"\"\"Solution of Euler equation\"\"\"\r\n", | |
| " return ((1 + ρ * α) * A * y ** α\r\n", | |
| " - ρ * α * (A ** 2) * (x ** α) * (y ** (α - 1)))\r\n", | |
| "\r\n", | |
| "\r\n", | |
| "def G(x):\r\n", | |
| " \"\"\"Dynamical system\"\"\"\r\n", | |
| " return [\r\n", | |
| " x[1],\r\n", | |
| " F(x[0], x[1])\r\n", | |
| " ]\r\n", | |
| "\r\n", | |
| "if __name__ == \"__main__\":\r\n", | |
| "\r\n", | |
| " duration = 4\r\n", | |
| " x0 = [0.8, 0.43]\r\n", | |
| "\r\n", | |
| " x = x0[:]\r\n", | |
| " for t in range(duration):\r\n", | |
| " print(x)\r\n", | |
| " x = G(x)\r\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "%cat listing01.py" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "### 補足\n", | |
| "\n", | |
| " if __name__ == \"_main__\":\n", | |
| " [do something]\n", | |
| " \n", | |
| " \n", | |
| "これはPythonのイディオムで、 `import`文で呼び出された場合には実行されずに、`run` で呼び出された場合には実行されます. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "### 可視化の基本" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 10, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGgRJREFUeJzt3XuUXWV9xvHvkxBEsIgaiZgEg4JcxIRLSWhQGOQ2CeWq\nFQNIAREqjbW0BbysQlZpbbWLii6WGFNu1QqooHIzEQpTaZpQaEmAmmAiILlJBdFapTaRX/94d2QY\nZuZcZp+zL+f5rHWWc+Zs9vntlfDw+tvv+25FBGZmVh/jii7AzMzy5WA3M6sZB7uZWc042M3MasbB\nbmZWMw52M7OaaRjskq6R9LSkR0b4/DRJKyU9LGmppOn5l2lmZs1qZsR+LdA/yuePA4dGxHTgMuCL\neRRmZmbtaRjsEXEf8Nwony+LiJ9lb+8HpuRUm5mZtSHvHvsHgDtzPqeZmbVgm7xOJOlw4GzgkLzO\naWZmrcsl2LMbpouA/ogYtm0jyZvSmJm1ISLUyvFjDnZJuwK3AKdHxNrRjm21uCqRtCAiFhRdR6fU\n+frqfG3g66u6dgbFDYNd0g3AYcBESeuAS4EJABGxELgEeA1wlSSAzRExs9VCzMwMJMYD+wNHZq+W\nNQz2iJjX4PNzgHPa+XIzs14nIWB3Xgzyw4FNwN3A54AjWj1nbjdPjYGiC+iwgaIL6KCBogvosIGi\nC+iwgaILaJXEzqTA3hrm40lB/g1gfgSbBh3b+vm79aANSVHnHruZ2UgkXgm8Ezgqe00D/pkU5ncB\nj0UwbBi3k50OdjOznEmMA2YAR5OCfBawkhTidwP/FsHm5s7lYDczK4TEZFKIH01qrzwHfIcU5gMR\n/Hd753Wwm5l1hcT2wKGkID8a2AX4J1KQfyeCH+bzPQ52M7OOyGav7Asck70OBlYAS0hh/mAEv87/\nex3sZma5kXgdqb1yDGlU/itSkC8B7mm3vdJaDQ52M7O2ZYuDZpKCvB/YhzR7ZWuYrx1p9krnanKw\nm5m1RGIXUpDPId303AAszl5LI/hVgeU52M3MGpGYQOqPz8lebyJNQVwMLIlgQ4HlvYyD3cxsGBJv\nJLVW5pJWfD5BenbEYmB5BFsKLG9UDnYzM0BiG9KioLnZa+uo/E5gcQQ/KrC8ljjYzaxnSUwkjcqP\nJc1gWUcK8m8Dy8o8Kh+Ng93MekY2r3w/UpAfC+wN3EMW5mXrlbfLwW5mtSaxA2nmytYw/wVwR/a6\nr+gZLJ3gYDez2pF4E/C72esQ4AHgduCOCL5fZG3d4GA3s8rLFgnNIgX5ccAkUp/8NuCuCH5WYHld\n52A3s0qSeBXphufxpFksPyIF+e2kLW5z34OlKhzsZlYZElNII/LjSS2W5cCtwO0RPFlgaaXiYDez\n0ho0i+X47DWNNIPlVtKKz45vqFVFDnYzK5Vs+f5hwAmkMN8MfIsU5kurOre8m9rJTj/M2sxyJbEj\naaHQidn/fp8U5nOAVd3eHbEXecRuZmMm8QbSiPxE4B3AUrKReQQbi6yt6tyKMbOukXgLcFL22oe0\nodY3Sas+3S/PiYPdzDomu/k5gxfDfGfSqPwbwL11XPVZBg52M8uVxDjSYqF3Aydnv76FFObLe3l+\nebf45qmZjVm25e2hpDA/CXgOuDn7+WHf/Cw/B7uZIbEt8C7gPaSpiU+RwvzwCB4rsjZr3bhGB0i6\nRtLTkh4Z5ZjPSVojaaWk/fMt0cw6QWI7ieMkrict4b8EWAXMjODACD7pUK+mhsEOXEuaizosSXOB\n3SNiD+Bc4KqcajOznGVhfoLEl4FNwIXAg8D0CGZHcHkETxRbpY1Vw1ZMRNwnadoohxwPXJ8de7+k\nnSRNioin8ynRzMZCYjvgGOC9pA22VgJfBS6MYFORtVln5NFjn0x6BNVW64EpgIPdrCASr+DFMD8W\nWAF8DfjTKj3v09qT183ToVNxfNfcrMuyG6BHkcL8OOAR0sj8zxzmvSWPYN8ATB30fkr2u5eRtGDQ\n24GIGMjh+816VjY18XDgfaTl/KuBm4CPeSl/NUnqA/rGdI5mFihlPfbbIuLtw3w2F5gfEXMlHQxc\nEREHD3OcFyiZ5SBbNHQIKczfA/yQFOZfjXhJW9RqoCMLlCTdQNp2c6KkdcClwASAiFgYEXdKmitp\nLenBsme1XrqZjSZbzn8AMA84hbRo6CZgdgQ/KLI2Kx9vKWBWYhJvBU4lBfo2wA3AjRE8Wmhh1jXe\nUsCsBiTeSGqznEqadXYTcAbp2Z+emGANecRuVgISryZtsnUaqeXyTeArpF0TvdFWD/OI3axCsumJ\n/cDppDnn9wBfAO6I4Pkia7Nq84jdrIuym6AHk8L8vaS9Wf4R+FoEPymyNisnj9jNSkrizcD7SYH+\na+BLwEERPFlkXVZPDnazDpHYCfg90o3PPUk3QU8DHvBNUOskt2LMcpStBD0K+H1S//xu0uj82xH8\nX5G1WTW5FWNWEIl9SGH+ftKmeNcD57tvbkVwsJu1KWu1vI+02noqaWR+ZATfK7Qw63luxZi1INun\n5V3A2aS9zb8DXAd8J4ItBZZmNdVOdjrYzZogMQ04kzQ6f5b0ZLGvRPBsgWVZD3CP3SxH2cMqTgTO\nAfYn7dNyYgQPFVqYWQMOdrMhJPYlhflppMfIXQ18I4L/LbQwsyY52M0AiR1I2+F+ENiV1GqZFcHj\nhRZm1gb32K2nSewHnEua3fIvwCLSnHPfCLVScI/drAnZ6Px9wHnAG4C/B6ZHsL7Qwsxy4hG79QyJ\ntwN/QHpoxX3AQmCJt8W1MvOI3WwIie1IzwX9A2AaqdXi0bnVmoPdakliN1Kr5WxgBXA5cJt759YL\nHOxWG9mq0H7gfNKe59cDh0SwptDCzLrMwW6VJ/Fa0orQ84GfAlcC743gl4UWZlYQB7tVlsR04MOk\nHvodpAVF93uvc+t1DnarlGy/8+OBPwJ2Jz0jdK8Ini60MLMScbBbJWTtlnOAPwTWA58lLfPfXGhh\nZiXkYLdSk9gL+AhpQdFtwLsjeLDYqszKzcFupSMh4EjgAuAA0kKivSP4UaGFmVWEg91KI1tMdCop\n0AE+A5zsXRXNWuNgt8JJTAQ+RJquuIIU7P/k2S1m7RlXdAHWuyR2l/g8sAZ4E3BEBHMiuNuhbta+\nhsEuqV/SaklrJF08zOcTJS2WtELSo5LO7EilVhsSsyRuBpYDz5H65+f4IdBm+Rh1d0dJ44HHSDey\nNgAPAPMiYtWgYxYAr4iIj0mamB0/KSK2DDmXd3fsYdkN0TnAxaTR+WeAqyP4n0ILMyu5TuzuOBNY\nGxFPZl9wI3ACsGrQMZuA6dnPOwLPDg11613ZgqJTSIH+AvBp4KvejMuscxoF+2Rg3aD364FZQ45Z\nBNwjaSPwW8B78yvPqkrilaT9Wy4EfghcRNr73L1zsw5rFOzN/Ev4cWBFRPRJegtwl6QZEfHzoQdm\nbZutBiJioOlKrRIkdiTtfX4BqXV3agTLiq3KrDok9QF9YzlHo2DfAEwd9H4qvOwBBbOBvwKIiB9I\negLYE16+OjAiFrRdqZWaxOtI+7ecD9wFHB3BI8VWZVY92YB3YOt7SZe2eo5Gs2IeBPaQNE3StqRe\n6a1DjllNurmKpEmkUPeT3XuExM4SnyJNWZwMzI7gVIe6WXFGHbFHxBZJ84ElwHjg6ohYJem87POF\nwCeBayWtJP2H4qKI+EmH67aCSbyB1Dc/E7gB2C+CpwotyswAP8zaWiSxC2mGyxnAl4BPR7Ch2KrM\n6qud7PTKU2uKxC4SVwD/SZq2+LYIPuJQNysfB7uNSmKSxOWkQA9gnwj+JIJNBZdmZiNwsNuwJF6X\n3RRdBUwA9o3gAm+da1Z+DnZ7CYlXSywAvg/sBMyI4I8i2FhsZWbWLAe7AWmlqMSFpGmLuwEHRXBe\nxEtWHptZBXg/9h4nMQE4G7iEtNtin3dZNKs2B3uPkhgHvAf4S+Ap4MQIHii2KjPLg4O9B0kcTtpl\ncRzwhxHcVXBJZpYjB3sPkdgX+BSwN/AJ4KYIXii2KjPLm2+e9oBscdEi4B7SBl17R3CDQ92snhzs\nNSaxg8QlwCOkR9C9NYIrIvhVwaWZWQe5FVND2Y3R00nbKS8lTV18otiqzKxbHOw1IzEbuIK0/P+U\nCP614JLMrMsc7DUhMYV0Y/Qw4KPAV9xDN+tN7rFXnMR2Eh8HVgJPAHtF8GWHulnv8oi9wiSOBT4L\nPErqo/vJVWbmYK8iid1Igb4XaYHRkoJLMrMScSumQrK2y58DDwD/CrzdoW5mQ3nEXhESRwKfJz3w\n4sAIflhwSWZWUg72kpOYBPwdcAgwP4LbCy7JzErOrZiSkhgncQ5p1eh60jNGHepm1pBH7CUksSfw\nRWA74MgIHi64JDOrEI/YS0RiQjYnfSlwMzDboW5mrfKIvSQkDgCuAX6Eb46a2Rh4xF4wiVdIXAYs\nBj4DzHGom9lYeMReIIkDgeuAx4EZEWwqtiIzqwOP2AuQ9dIXAN8G/ob0vFGHupnlwiP2LpN4G/AP\nwNPAfhFsLLgkM6uZhiN2Sf2SVktaI+niEY7pk/SQpEclDeReZQ1k89L/GBgAvgAc61A3s05QRIz8\noTQeeAw4EthA2qNkXkSsGnTMTqTpecdExHpJEyPimWHOFRGhvC+gCiQmA9cDOwDvj2BtwSWZWUW0\nk52NRuwzgbUR8WREbAZuBE4YcsypwM0RsR5guFDvZRInAf8BfBd4p0PdzDqtUY99MrBu0Pv1wKwh\nx+wBTJB0L/BbwGcj4kv5lVhNEtuTpi8eRbo5uqzgksysRzQK9pH7NC+aABwAHAFsDyyTtDwi1gw9\nUNKCQW8HImKgyTorJbtBehPpqUb7RfDfBZdkZhUhqQ/oG8s5GgX7BmDqoPdTSaP2wdYBz0TE88Dz\nkr4LzABeFuwRsaD9UstPQsDZpCmMFwHXRTT1H0czMwCyAe/A1veSLm31HI167A8Ce0iaJmlb4BTg\n1iHHfAt4h6TxkrYntWq+12ohVSexA+kG6QXAoRFc61A3syKMOmKPiC2S5gNLgPHA1RGxStJ52ecL\nI2K1pMXAw8ALwKKI6Klgl9gH+DqwHJgZwS8LLsnMetio0x1z/aKaTneUOAW4Erg4gmuKrsfM6qWd\n7PTK0zZJTAA+DRwPHB3BQwWXZGYGONjbIvF64GvAL4HfjuC5gksyM/sNbwLWomzf9AeAfwGOc6ib\nWdl4xN6CQf30D0Xw9aLrMTMbjoO9CRLjgEuBM4GjIlhRbEVmZiNzsDcg8UrS/PQppKmMTxdckpnZ\nqNxjH4XEzsC9wGbgXQ51M6sCB/sIJPYiLThaApwewf8WXJKZWVPcihmGxCHALcBFEVxfdD1mZq1w\nsA8hcSLwRdIDMZYUXY+ZWavcihlE4lzg88Ach7qZVZVH7BmJjwLnknZm9FOOzKyyej7Ysz3U/wb4\nXdKj6zYUXJKZ2Zj0dLBnC4+uBH6bNFJ/tuCSzMzGrGeDXWI86SbpnsCRfnydmdVFTwZ7FurXkVaT\n9kfwP8VWZGaWn54L9izUrwV2AY71047MrG56Ktiznvrfkx7K7VA3s1rqmWDPZr9cBbyFNE/doW5m\ntdQTwZ6F+t8CM0jb7v6i4JLMzDqmJ4Id+DhwDHBYBD8vuhgzs06qfbBLnA+cRVp89JOi6zEz67Ra\nB7vEScAnSKG+qeh6zMy6obbBLjEbWEi6Ufp40fWYmXVLLXd3lNiTtJ/6GRH8e9H1mJl1U+2CXeK1\nwO3AJyJYXHQ9ZmbdpojozhdJERHq7HcwAVgMPBTBn3Xyu8zMuqGd7KxNsA9agDQFOCGCX3fqu8zM\nuqWd7GzYipHUL2m1pDWSLh7luIMkbZF0cisF5Og84B3AqQ51M+tlowa7pPGk/cr7gX2AeZL2HuG4\nT5HaIB1ttwxHYhbwF8BJ3n7XzHpdoxH7TGBtRDwZEZuBG4EThjnuw8DXgR/nXF9DEq8HvgZ8MII1\n3f5+M7OyaRTsk4F1g96vz373G5Imk8L+quxX3Wna85steG8AvhTBt7r1vWZmZdZogVIzIX0F8NGI\nCElilFaMpAWD3g5ExEAT5x/Nx0j/cbpkjOcxMysFSX1A35jOMdqsGEkHAwsioj97/zHghYj41KBj\nHufFMJ8I/BL4YETcOuRcuc6KyVaW3gIc6AdQm1ldtZOdjUbsDwJ7SJoGbAROAeYNPiAi3jyogGuB\n24aGet4kdgL+kdRXd6ibmQ0yarBHxBZJ84ElwHjg6ohYJem87POFXahxOFcBd0RwW0Hfb2ZWWpVb\noCTxe8BlwP4RPD/2yszMyqsTrZhSyaY2fo40X92hbmY2jEqN2CVuAp6K4MKcyjIzK7Vaj9gl3g3s\nB5xZcClmZqVWiWCX2JHUgjnFLRgzs9FVohUjcTnw2gjOyrksM7NSq2UrRmJf4AzgbUXXYmZWBaV+\nglK2x/qVwIII/qvoeszMqqDUwQ68D9gR+ELRhZiZVUVpe+wS2wKPAWdG8M+dq8zMrLw68gSlAp0L\nrHKom5m1ppQjdolXAWuAORGs6GxlZmblVacR+wXAvQ51M7PWlW7ELjGR1FufFcHazldmZlZedRmx\nfwS42aFuZtaeUo3YJXYAngRm+8HUZmb1GLGfBdznUDcza19pthSQ2Ab4E+C0omsxM6uyMo3YTwY2\nRrCs6ELMzKqsFMGe7QlzIfC3RddiZlZ1pQh2YCawE/jh1GZmY1WWYD8TuC6CF4ouxMys6gqf7iix\nHbAB2D+Cp7pSjJlZRVR1uuNxwEMOdTOzfJQh2H8fuL7oIszM6qLQVozEG4BVwJQIftGVQszMKqSK\nrZjTgG861M3M8lN0sJ+B2zBmZrkqLNgldgXeCHy3qBrMzOqoqWCX1C9ptaQ1ki4e5vPTJK2U9LCk\npZKmN3Hao4C7PXfdzCxfDYNd0njgSqAf2AeYJ2nvIYc9DhwaEdOBy4AvNvHdRwF3tVaumZk10syI\nfSawNiKejIjNwI3ACYMPiIhlEfGz7O39wJTRTigxDjgCB7uZWe6aCfbJwLpB79dnvxvJB4A7G5xz\nf+CZiJec18zMctDMfuxNT3SXdDhwNnBIg0PdhjEz65Bmgn0DMHXQ+6mkUftLZDdMFwH9EfHccCeS\ntCD9NP8MeNUi+OsWyzUzqzdJfUDfmM7RaOWppG2Ax0g98Y3AvwHzImLVoGN2Be4BTo+I5SOcJyJC\nEtsDTwNvjODnYynezKzu2ll52nDEHhFbJM0HlgDjgasjYpWk87LPFwKXAK8BrpIEsDkiZo5wykNJ\nm3451M3MOqDre8VIXA78NILLuvLFZmYVVpW9Yg4ClhbwvWZmPaGIYJ9E6tWbmVkHFBHsOwP/VcD3\nmpn1hK4Gu8S2wKuAn3bze83Mekm3R+yvB37sjb/MzDqn28HuNoyZWYd1O9gnkRYnmZlZh3jEbmZW\nMw52M7OacbCbmdWMg93MrGYc7GZmNeNgNzOrGU93NDOrmSJG7D/u8neamfWUbgf7ryJ4vsvfaWbW\nU7od7O6vm5l1mIPdzKxmHOxmZjXjYDczq5luB7unOpqZdZhH7GZmNeNgNzOrGQe7mVnNONjNzGrG\nwW5mVjPdDvafdPn7zMx6TleDPYJfd/P7zMx6UcNgl9QvabWkNZIuHuGYz2Wfr5S0f/5lmplZs0YN\ndknjgSuBfmAfYJ6kvYccMxfYPSL2AM4FrupQraUmqa/oGjqpztdX52sDX18vajRinwmsjYgnI2Iz\ncCNwwpBjjgeuB4iI+4GdJE3KvdLy6yu6gA7rK7qADuoruoAO6yu6gA7rK7qAsmkU7JOBdYPer89+\n1+iYKWMvzczM2tEo2KPJ86jNf87MzHK2TYPPNwBTB72fShqRj3bMlOx3LyOp1oEv6dKia+ikOl9f\nna8NfH29plGwPwjsIWkasBE4BZg35JhbgfnAjZIOBn4aES/bxTEiho7qzcysA0YN9ojYImk+sAQY\nD1wdEasknZd9vjAi7pQ0V9Ja4BfAWR2v2szMRqSIWndHzMx6Tu4rT+u+oKnR9Uk6LbuuhyUtlTS9\niDrb0cyfXXbcQZK2SDq5m/WNVZN/N/skPSTpUUkDXS5xTJr4uzlR0mJJK7LrO7OAMtsi6RpJT0t6\nZJRjqpwro15fy7kSEbm9SO2atcA0YAKwAth7yDFzgTuzn2cBy/OsoZOvJq/vd4BXZz/3V+X6mrm2\nQcfdA9wOvLvounP+s9sJ+E9gSvZ+YtF153x9C4C/3nptwLPANkXX3uT1vRPYH3hkhM8rmytNXl9L\nuZL3iL3uC5oaXl9ELIuIn2Vv76c6c/qb+bMD+DDwdeDH3SwuB81c36nAzRGxHiAinulyjWPRzPVt\nAnbMft4ReDYitnSxxrZFxH3Ac6McUuVcaXh9reZK3sFe9wVNzVzfYB8A7uxoRflpeG2SJpPCYuu2\nEVW6QdPMn90ewGsl3SvpQUnv71p1Y9fM9S0C3iZpI7AS+EiXauuGKudKqxrmSqPpjq2q+4KmpuuU\ndDhwNnBI58rJVTPXdgXw0YgISeLlf45l1sz1TQAOAI4AtgeWSVoeEWs6Wlk+mrm+jwMrIqJP0luA\nuyTNiIifd7i2bqlqrjSt2VzJO9hzXdBUQs1cH9mNjUVAf0SM9n8fy6SZazuQtF4BUo92jqTNEXFr\nd0ock2aubx3wTEQ8Dzwv6bvADKAKwd7M9c0G/gogIn4g6QlgT9J6laqrcq40pZVcybsV85sFTZK2\nJS1oGvov/a3AGVmhIy5oKqmG1ydpV+AW4PSIWFtAje1qeG0R8eaI2C0idiP12T9UkVCH5v5ufgt4\nh6TxkrYn3YT7XpfrbFcz17caOBIg6z/vCTze1So7p8q50lCruZLriD1qvqCpmesDLgFeA1yVjWw3\nR8TMompuVpPXVllN/t1cLWkx8DDwArAoIioR7E3++X0SuFbSStKg7qKIqMRTzSTdABwGTJS0DriU\n1DqrfK5A4+ujxVzxAiUzs5rp9jNPzcyswxzsZmY142A3M6sZB7uZWc042M3MasbBbmZWMw52M7Oa\ncbCbmdXM/wN2vHdz0f3yIAAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x108221a90>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "# listing02.py\n", | |
| "\n", | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as plt # プロッティング・ライブラリを読み込む\n", | |
| "from listing01 import f\n", | |
| "\n", | |
| "fig = plt.figure()\n", | |
| "ax = fig.add_subplot(111) # 作図領域の作成\n", | |
| "\n", | |
| "grids = np.linspace(0.0, 1.2, 120)\n", | |
| "ax.plot(grids, f(grids))\n", | |
| "\n", | |
| "plt.show() # 図の表示" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "### 軌道の可視化" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 11, | |
| "metadata": { | |
| "collapsed": false, | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHiVJREFUeJzt3X2YVXW99/H3R4QYJbKioQTMOWCmEKYmlHZsTC0Yj2jW\nyXzIo6Z5KqBjh7SH+1auOqcnqAS8jhmpWeekVvaARZDdNncdE7P7iIqhAWI8aE6a+dQmIb/3H7+F\nDsPM7D0ze++1Hz6v6+KSvfeatb/rAj/+/K7f77cUEZiZWePYI+8CzMysvBzsZmYNxsFuZtZgHOxm\nZg3GwW5m1mAc7GZmDaZosEu6WtIjku7p4/MzJN0l6W5Jt0qaWv4yzcysVKWM2K8BZvTz+QPA0REx\nFfg08NVyFGZmZoNTNNgj4pfA4/18fltEPJG9vB0YX6bazMxsEMrdY38fsLzM5zQzswHYs1wnknQM\ncC5wVLnOaWZmA1eWYM9umC4FZkREr20bSd6UxsxsECJCAzl+yMEuaT/ge8CZEbG+v2MHWlw9kTQ/\nIubnXUelNPL1NfK1ga+v3g1mUFw02CVdB7wFGCNpM3ApMBwgIq4ELgFeClwhCWB7REwbaCFmZgYS\nw4BDgeOyXwNWNNgj4rQin58HnDeYLzczawZqa+ugtXUuLS0jKRS20dW1ODZuXA4gIWASLwT5McDD\nwM+AxcCxA/2+st08NTrzLqDCOvMuoII68y6gwjrzLqDCOvMuoD9qa+tg8uRFzJs36fk3F3zpNXrZ\nZ7/D4x8fQwrzYaQg/z4wO4KHn//5QTSwVa0HbUiKRu6xm5lNkTomwNwXw8inYNtmWHzvtGlz+fzn\n377bwbO/0cW9V/8bcDNwfwS9hvFgstMjdjOzMpgidbwJFi1NbRUAzuKlb1i7fb8XPdfbD4x5YG0E\nSypRizcBMzMrg31puah7qAN8g8dfPmrbul5znUKhUKlaHOxmZoMgsZfEDIkvSazZzmG9Ls587ab1\nv2fhwl2ngi9YsIGuroqM1sGtGDOzPmmsOhjNXEYwkmfZxh9PWMYTP9oLeDvwRmA1sBJ430hu/RTw\ntp7nePkzz2zl3nuXcPHFc2hpaaFQKNDVtWTnrJiK1O2bp2Zmu9NYdTBuj8s56bm259/87qt3sP6i\nlWz74NeAWyJ4cudHvfXYz4MNq2DumohBh/hgstPBbmaWyRYHTQPezrij/oXzf/WS3Q76T1bEupjZ\n289PkTr2gzmjoOVpKGyCJUMJ9VSTZ8WYmQ2IxKtIrZWZpDnlW4EV7P27B4FDdvuBEbT0da4sxHPf\n4dbBbmZNRWI4qT8+M/v1atLioBXARyLYCqADHp1Kb8H+LBWbzVIuDnYza3gS+5KeBNdBWqK/kTSy\nngOsimDHbj/0JItZxkRmdZvCuIwNPFmZuefl5B67mTUciT2B6aQg7+CFUflyYEUEfyjpPGlWzBxG\n0MKzFHiSJfHI0HrmA+Wbp2bWtCTGkEblJ5CmHW4mBflPgNt6HZXXAd88NbOmke2K+Hq4+kK4dxbs\nPQqeeBQm3QBzpu7slTcjB7uZ1Q2JvUkzV05Iv34asGZv+NLOaYlj4fwOmLIS1jRtsHtLATOraRKv\nlviQxE9I+5TPBdYCx8CX18CX9tn1J5ZOgv3mVL/S2uERu5nVlGyR0HTgH4ATgbGkPvlVwHsieOKF\nY188svezjOpzrnkzcLCbWe4kRpFueM4izWL5A3AT8H7g1xH8rfeffGpb7+8/XfNzzSvJwW5muZAY\nTxqRzwKOAlYBy4D5ETxY2lk2L4bzJ6b2y07nbYBNNT/XvJI83dHMquKFWSzMyn7tT5qOuAxY2X1D\nrYGdd0pH6qmPakkj9U1LItbkvqy/XDyP3cxqSrZ8/y3ASaQw3w78kBTmt9br3PJq8jx2M8udxGjS\nQqGTs3/+jhTmM4G1fT3b08rHI3YzGzKJV5JG5CcDbwZuJRuZR/BQnrXVO7dizKxqJCYC78h+HUza\nHfEHwE8G2y+33TnYzaxispufh/BCmLeSRuXfB34ewV9zLK9hOdjNrKwk9iAtFnoncEr29vdIYb6q\n7/nlVi6+eWpmQ5ZteXs0KczfATwO3Jj9/m7f/Kx9DnYzQ2IE8FbgXaSpiZtIYX5MBPfnWZsNXNFN\nwCRdLekRSff0c8xiSesk3SXp0PKWaGaVIDFS4kSJa0lL+C8hba41LYLDI/hM91CX1CFphaTO7J8d\nedVu/StlxH4NsAT4Rm8fZn+4kyLiAEnTgStIzxM0sxojMZL04OZ/JG19ew/wHeCTEWzp++fUASyC\nbo+Jg4mSiKjuE4WsuKIj9oj4JanH1pdZwLXZsbcD+0gaW57yzGyospH5SRL/Rdr29kLgV8DBERwd\nwZL+Qj0zl11Dnex1U2+PW6vK0WMfR3oE1U5bgPHAI2U4t5kNgsSLSCPzd5NG5qtJI/N/LfV5nz30\nsT0uTb09bq0q183TnlNxfNfcrMqyG6DHk8L8RFKb5dvAvEGGeXd9bI9LU2+PW6vKEexbgQndXo/P\n3tuNpPndXnZGRGcZvt+saWVTE48B3kNazn8fcAPw8TIv5V8MTGTXdswG0v03KyNJ7UD7kM5RygIl\nSfsDN0XE63r5rAOYHREdkt4IXBYRu9089QIls/LIFg0dRQrzdwG/J4X5tyN2aYuW+XvVQeqpt5BG\n6kt847TyKrLyVNJ1pG03x5D65pcCwwEi4srsmMtJu7g9A5wTEf9TjuLMLMmW8x8GnAacSprQcANw\nfQQb8qzNKstbCpg1GInXAKeTAn1P4DpSmK/JtTCrGm8pYNYAJPYltVlOJ806uwE4i/TsT09MsKI8\nYjerARIvIW2ydQap5fID4FukXRO90VYT84jdrI5k0xNnAGeS5pzfAnwF+HGEpxHa4HnEblZF2U3Q\nN5LC/N2kvVn+C/hOBH/KszarTR6xm9Uoib8D3ksK9L8B3wSOiODBPOuyxuRgN6sQiX1Im22dBRxI\nugl6BnCHb4JaJbkVY1ZG2UrQ44F/IvXPf0Yanf8kgmfzrM3qk1sxZjmROJgU5u8lbYp3LfBB980t\nDw52s0HKWi3vAc4h7Zf0TeC4CH6ba2HW9NyKMRuAbJ+WtwLnAh3AT4GvAz+NYEeOpVmD8pYCZhUi\nsT9wNml0/hjpyWLfiuCxHMuyJuAeu1kZZQ+rOBk4DziUtE/LyRHcmWthZkU42M16kJhCCvMzgLuA\nq4DvR/T5sAmzmuJgNwMk9iZth3s+sB+p1TI9ggdyLcxsENxjt6Ym8Xrg/aTZLf8NLCXNOfeNUKsJ\n7rGblSAbnb8HuAB4JfA1YGoEW3ItzKxMPGK3piHxOuCfSQ+t+CVwJbDS2+JaLfOI3awHiZGk54L+\nM7A/qdVStdF5W5s6WluZ29LCyEKBbV1dLN640c8JtcpysFtDkmgjtVrOBVYDXwRuqmbvvK1NHZMn\ns2jePCbtfG/hQia2tQmHu1XSHnkXYFYuEntIdEj8CLiD9ND1oyJ4WwTfr/YN0dZW5nYPdYB585jU\n2sqcatZhzccjdqt7Ei8jrQj9IPBn4HLg3RH8Jc+6WloY2cf7LdWuxZqLg93qlsRUYA6ph/5j0oKi\n22tlr/NCofcFTYWCH3tnleVWjNUViT0lTpHoBJYDvwdeG8GZEayqlVAH6Opi8cKFrO/+3oIFbOjq\nYkleNVlz8HRHqwtZu+U84EPAFmARaZn/9lwLKyKbFTOnpYWWQoFCVxdLfOPUBsK7O1rDkXgt8GHS\ngqKbgMUR/Cbfqsyqx/PYrSFICDgOuBA4jLSQ6KAI/pBrYWZ1wsFuNSNbTHQ6KdABvgyc4l0VzQbG\nwW65kxgDfIA0XXE1Kdj/Ty3dCDWrJ54VY7mRmCTxH8A64NXAsRHMjOBnDnWzwSsa7JJmSLpP0jpJ\nF/fy+RhJKyStlrRG0tkVqdQahsR0iRuBVcDjpP75eX4ItFl59DsrRtIw4H7SjaytpGXap0XE2m7H\nzAdeFBEflzQmO35sROzocS7Pimli2Q3RmcDFpNH5l4GrIng618LMalwlZsVMA9ZHxIPZF1wPnASs\n7XbMw8DU7Pejgcd6hro1L4k9SU8muhh4DvgC8G0/yMKscooF+zhgc7fXW4DpPY5ZCtwi6SHgxcC7\ny1ee1SuJFtL+LR8lrQ69iLT3uXvnZhVWLNhL+ZfwE8DqiGiXNBG4WdIhEfFUzwOzts1OnRHRWXKl\nVhckRpP2Pr+Q1Lo7PYLb8q3KrH5Iagfah3KOYsG+FZjQ7fUE2O0BBUcC/w4QERskbQQOhN1XB0bE\n/EFXajVN4uXAXNKUxZuBt0VwT75VmdWfbMDbufO1pEsHeo5is2J+AxwgaX9JI0i90mU9jrmPdHMV\nSWNJoe4nuzcJiVaJz5OmLI4DjozgdIe6WX76HbFHxA5Js4GVwDDgqohYK+mC7PMrgc8A10i6i/Qf\niosi4k8VrttyJvFKUt/8bOA64PURbMq1KDMDvAmYDZDEq0gzXM4Cvgl8IYKt+VZl1rgGk51eeWol\nkXiVxGXAvaRpi5Mj+LBD3az2ONitXxJjJb5ICvQADo7gIxE8nHNpZtYHB7v1SuLl2U3RtaSHQk+J\n4EJvnWtW+xzstguJl0jMB34H7AMcEsHcCB7KtzIzK5WD3YC0UlTio6Rpi23AERFcELHLymMzqwPe\nj73JSQwHzgUuIe222O5dFs3qm4O9SUnsAbwL+DdgE3ByBHfkW5WZlYODvQlJHEPaZXEP4EMR3Jxz\nSWZWRg72JiIxBfg8cBDwSeCGCJ7LtyozKzffPG0C2eKipcAtpA26DorgOoe6WWNysDcwib0lLgHu\nIT2C7jURXBbBX3MuzcwqyK2YBpTdGD2TtJ3yraSpixvzrcrMqsXB3mAkjgQuIy3/PzWCX+VckplV\nmYO9QUiMJ90YfQvwMeBb9dxDb1NbRyutc1toGVmgsK2LrsUbY+PyvOsyqwcO9jonMRL4CPCvwBXA\nBRE8nW9VQ9Omto7JTF40j3mTdr63kIUT29SGw92sON88rWMSJwBrgGmkPvr/qvdQB2ildW73UAeY\nx7xJrbTOyasms3riEXsdkmgDFgGvJS0wWplzSWXVQsvIPt5vqXYtZvXII/Y6IjFS4n8DdwC/Al7X\naKEOUKCwrY/3C9WuxaweOdjrhMRxwN3AYcDhEXyuUeejd9G1eCEL13d/bwELNnTRtSSvmszqiZ95\nWuMkxgJfAo4CZkfwo5xLqopsVsycFlpaChQKXXQt8Y1Ta0aDyU4He43KFhmdC3wGuAb4VATP5FuV\nmVXbYLLTN09rkMSBwFeBkcBxEdydc0lmVkfcY68hEsMlPkHaBuBG4EiHupkNlEfsNULiMOBq4A+k\nm6O/z7kkM6tTHrHnTOJFEp8GVgBfBmY61M1sKDxiz5HE4cDXgQeAQyJ4ON+KzKwReMSeg6yXPh/4\nCfA50vNGHepmVhYesVeZxGTgG8AjwOsjeCjnksyswRQdsUuaIek+SeskXdzHMe2S7pS0RlJn2ats\nABJ7SPwL0Al8BTjBoW5mldDvAiVJw4D7geOAraQ9Sk6LiLXdjtmHND3v7RGxRdKYiHi0l3M17QIl\niXHAtcDewHsjWF/kR8zMgMFlZ7ER+zRgfUQ8GBHbgeuBk3occzpwY0RsAegt1JuZxDuA/wF+Afy9\nQ93MKq1Yj30csLnb6y3A9B7HHAAMl/Rz4MXAooj4ZvlKrE8Se5GmLx5Pujl6W84lmVmTKBbspWwk\nM5y04+CxwF7AbZJWRcS6ngdKmt/tZWdEdJZYZ13JbpDeANxFukH6ZM4lmVmdkNQOtA/lHMWCfSsw\nodvrCaRRe3ebgUcjogAUJP0COATYLdgjYv7gS619EiJt3PU54CLg6xEl/cfRzAyAbMDbufO1pEsH\neo5iPfbfAAdI2l/SCOBUYFmPY34IvFnSMEl7kVo1vx1oIfVOYm/SDdILgaMjuMahbmZ56HfEHhE7\nJM0GVgLDgKsiYq2kC7LPr4yI+yStID0E4jlgaUQ0VbBLHAx8F1gFTIvgLzmXZGZNzPuxD5HEqcDl\nwMURXJ13PWbWWLwfexVJDAe+AMwC3hbBnTmXZGYGONgHReIVwHeAvwBviODxnEsyM3ueNwEboGzf\n9DuA/wZOdKibWa3xiH0AuvXTPxDBd/Oux8ysNw72EmQPlr4UOBs4PoLV+VZkZtY3B3sREi2k+enj\nSVMZH8m5JDOzfrnH3g+JVuDnwHbgrQ51M6sHDvY+SLyWtOBoJXBmBNtyLsnMrCRuxfRC4ijge8BF\nEVybdz1mZgPhYO9B4mTgq6QHYqzMu55aoLa2Dlpb59LSMpJCYRtdXYtj48bleddlZr1zsHcj8X5g\nPjAzgv+Xczk1QW1tHUyevIh58yY9/+bChRPV1obD3aw2uceekfgY8DHSzowO9Z1aW+fuEuoA8+ZN\norV1Tk4VmVkRTT9iz/ZQ/xzwD6RH123NuaTa0tIyso/3W6pciZmVqKmDPVt4dDnwBtJI/bGcS6o9\nhULvs4EKhUKVKzGzEjVtK0ZiGLAUmAoc51DvQ1fXYhYu3PUB3AsWbKCra0lOFZlZEU25H3sW6l8n\nrSY9MYKn862otmWzYubQ0tJCoVCgq2uJb5yaVcdgsrPpgj0L9WuAfYFZftqRmdUyP2ijiKyn/jXS\nQ7lPcKibWSNqmmDPZr9cAUwkzVN3qJtZQ2qKYM9CfQFwCGnb3WdyLsnMrGKaItiBTwBvB94SwVN5\nF2NmVkkNH+wSHwTOIS0++lPe9ZiZVVpDB7vEO4BPkkL94bzrMTOrhoYNdokjgStJN0ofyLseM7Nq\nachglziQtJ/6Wd7Qy5qVxqqD0cxlBCN5lm08yeJ4JLywrAk0XLBLvAz4EfDJCFbkXY9ZHjRWHYxn\nEbN4YWfOZUzUWOFwb3wNtVeMxHDgO8API7gq73rMcjOaubuEOsAsJjEab7fcBBom2LO56kuAAnBx\nzuWY5WsEvW+3PAJvt9wEiga7pBmS7pO0TlKfgSnpCEk7JJ1S3hJLdgHwZuD0CP6WUw1NZ4rUMVNa\n8W6pc6a0YorUkXdNBjzbx8PXn8XbLTeBfnvskoaR9is/DtgK3CFpWUSs7eW4zwMrgKpv9CUxHfgU\ncFQET1b7+5vVFKnjTbBoKS/8L//5MHGKxJpwHzdXT7KYZUzs0WPfwJN4u+UmUOzm6TRgfUQ8CCDp\neuAkYG2P4+YA3wWOKHeBxUi8gtRXPz+CddX+/mY2AeZ2D3WApTCpI/19cLDnKB6J5Ror+E/mMIIW\nnqXAkyzxjdPmUCzYxwGbu73eAkzvfoCkcaSwfysp2KuzDzDPb8F7HfDNCH5Yre+15MX03scdhfu4\ntSALcQd5EyoW7KWE9GXAxyIiJIl+WjGS5nd72RkRnSWcvz8fJ90nuGSI57FBeIre+7hP4z6u2WBJ\nagfah3KOYsG+lbR3+U4TSKP27g4Hrk+ZzhhgpqTtEbGs58kiYv7gS91VtrJ0NnC4b5bmYzMsPh8m\ndm/HnAcbNuE+rtlgZQPezp2vJV060HP0+wQlSXsC9wPHAg8BvwZO63nztNvx1wA3RcT3evmsbE9Q\nktgHuBOYG8FN5TinDc4UqWM/mDMKWp6GwiZY4hunZuVT9icoRcQOSbOBlcAw4KqIWCvpguzzKwdd\n7dBcAfzYoZ6/LMQd5GY1pO6eeSrxj8CngUMj3Ms1s8bW8M88zaY2Lgbe4VA3M+tdXY3YJW4ANkXw\n0TKVZWZW0xp6xC7xTuD1wNk5l2JmVtPqItglRpNaMKe6BWNm1r+6aMVIfBF4WQTnlLksM7Oa1pCt\nGIkpwFnA5LxrMTOrBzW9H3u2x/rlwPwIuvKux8ysHtR0sAPvAUYDX8m7EDOzelGzPXaJEaTtDM6O\n4P9WrjIzs9o1mB57LY/Y3w+sdaibmQ1MTY7YJUYB64CZEayubGVmZrWrkUbsFwI/d6ibmQ1czY3Y\nJcaQeuvTI1hf+crMzGpXo4zYPwzc6FA3MxucmhqxS+wNPAgc6QdTm5k1xoj9HOCXDnUzs8GrmS0F\nJPYEPgKckXctZmb1rJZG7KcAD0VwW96FmJnVs5oI9mxPmI8CC/Kuxcys3tVEsAPTgH3AD6c2Mxuq\nWgn2s4GvR/Bc3oWYmdW73Kc7SowEtgKHRrCpKsWYmdWJep3ueCJwp0PdzKw8aiHY/wm4Nu8izMwa\nRa6tGIlXAmuB8RE8U5VCzMzqSD22Ys4AfuBQNzMrn7yD/SzchjEzK6vcgl1iP2Bf4Bd51WBm1ohK\nCnZJMyTdJ2mdpIt7+fwMSXdJulvSrZKmlnDa44Gfee66mVl5FQ12ScOAy4EZwMHAaZIO6nHYA8DR\nETEV+DTw1RK++3jg5oGVa2ZmxZQyYp8GrI+IByNiO3A9cFL3AyLitoh4Int5OzC+vxNK7AEci4Pd\nzKzsSgn2ccDmbq+3ZO/15X3A8iLnPBR4NGKX85qZWRmUsh97yRPdJR0DnAscVeRQt2HMzCqklGDf\nCkzo9noCadS+i+yG6VJgRkQ83tuJJM1Pv5t9FoxaCp8dYLlmZo1NUjvQPqRzFFt5KmlP4H5ST/wh\n4NfAaRGxttsx+wG3AGdGxKo+zhMRIYm9gEeAfSN4aijFm5k1usGsPC06Yo+IHZJmAyuBYcBVEbFW\n0gXZ51cClwAvBa6QBLA9Iqb1ccqjSZt+OdTNzCqg6nvFSHwR+HMEn67KF5uZ1bF62SvmCODWHL7X\nzKwp5BHsY0m9ejMzq4A8gr0V6Mrhe83MmkJVg11iBDAK+HM1v9fMrJlUe8T+CuCP3vjLzKxyqh3s\nbsOYmVVYtYN9LGlxkpmZVYhH7GZmDcbBbmbWYBzsZmYNxsFuZtZgHOxmZg3GwW5m1mA83dHMrMHk\nMWL/Y5W/08ysqVQ72P8aQaHK32lm1lSqHezur5uZVZiD3cyswTjYzcwajIPdzKzBVDvYPdXRzKzC\nPGI3M2swDnYzswbjYDczazAOdjOzBuNgNzNrMNUO9j9V+fvMzJpOVYM9gr9V8/vMzJpR0WCXNEPS\nfZLWSbq4j2MWZ5/fJenQ8pdpZmal6jfYJQ0DLgdmAAcDp0k6qMcxHcCkiDgAeD9wRYVqrWmS2vOu\noZIa+foa+drA19eMio3YpwHrI+LBiNgOXA+c1OOYWcC1ABFxO7CPpLFlr7T2teddQIW1511ABbXn\nXUCFteddQIW1511ArSkW7OOAzd1eb8neK3bM+KGXZmZmg1Es2KPE82iQP2dmZmW2Z5HPtwITur2e\nQBqR93fM+Oy93Uhq6MCXdGneNVRSI19fI18b+PqaTbFg/w1wgKT9gYeAU4HTehyzDJgNXC/pjcCf\nI2K3XRwjoueo3szMKqDfYI+IHZJmAyuBYcBVEbFW0gXZ51dGxHJJHZLWA88A51S8ajMz65MiGro7\nYmbWdMq+8rTRFzQVuz5JZ2TXdbekWyVNzaPOwSjlzy477ghJOySdUs36hqrEv5vtku6UtEZSZ5VL\nHJIS/m6OkbRC0urs+s7OocxBkXS1pEck3dPPMfWcK/1e34BzJSLK9ovUrlkP7A8MB1YDB/U4pgNY\nnv1+OrCqnDVU8leJ1/cm4CXZ72fUy/WVcm3djrsF+BHwzrzrLvOf3T7AvcD47PWYvOsu8/XNBz67\n89qAx4A98669xOv7e+BQ4J4+Pq/bXCnx+gaUK+UesTf6gqai1xcRt0XEE9nL26mfOf2l/NkBzAG+\nC/yxmsWVQSnXdzpwY0RsAYiIR6tc41CUcn0PA6Oz348GHouIHVWscdAi4pfA4/0cUs+5UvT6Bpor\n5Q72Rl/QVMr1dfc+YHlFKyqfotcmaRwpLHZuG1FPN2hK+bM7AHiZpJ9L+o2k91atuqEr5fqWApMl\nPQTcBXy4SrVVQz3nykAVzZVi0x0HqtEXNJVcp6RjgHOBoypXTlmVcm2XAR+LiJAkdv9zrGWlXN9w\n4DDgWGAv4DZJqyJiXUUrK49Sru8TwOqIaJc0EbhZ0iER8VSFa6uWes2VkpWaK+UO9rIuaKpBpVwf\n2Y2NpcCMiOjvfx9rSSnXdjhpvQKkHu1MSdsjYll1ShySUq5vM/BoRBSAgqRfAIcA9RDspVzfkcC/\nA0TEBkkbgQNJ61XqXT3nSkkGkivlbsU8v6BJ0gjSgqae/9IvA87KCu1zQVONKnp9kvYDvgecGRHr\nc6hxsIpeW0T8XUS0RUQbqc/+gToJdSjt7+YPgTdLGiZpL9JNuN9Wuc7BKuX67gOOA8j6zwcCD1S1\nysqp51wpaqC5UtYRezT4gqZSrg+4BHgpcEU2st0eEdPyqrlUJV5b3Srx7+Z9klYAdwPPAUsjoi6C\nvcQ/v88A10i6izSouygi6uKpZpKuA94CjJG0GbiU1Dqr+1yB4tfHAHPFC5TMzBpMtZ95amZmFeZg\nNzNrMA52M7MG42A3M2swDnYzswbjYDczazAOdjOzBuNgNzNrMP8f5DmIw0g6nr0AAAAASUVORK5C\nYII=\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x108250048>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "# listing03.py\n", | |
| "\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "from listing01 import f, G\n", | |
| "\n", | |
| "duration = 10\n", | |
| "x0 = [0.8, 0.43]\n", | |
| "\n", | |
| "fig = plt.figure()\n", | |
| "ax = fig.add_subplot(111)\n", | |
| "\n", | |
| "grids = np.linspace(0.0, 1.2, 120)\n", | |
| "ax.plot(grids, f(grids))\n", | |
| "\n", | |
| "x = x0[:]\n", | |
| "for t in range(duration):\n", | |
| " ax.plot(x[0], x[1], marker='o', linestyle='')\n", | |
| " x = G(x)\n", | |
| " \n", | |
| "plt.show()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "これでは動きが分からない .... " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 12, | |
| "metadata": { | |
| "collapsed": false, | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAAEACAYAAAC+rrMfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH55JREFUeJzt3XmUXFW5/vHv050ECBASCAQIIQkkMgcRCIMgYVDDIFzh\nJxBEmWRQUZxxJFFUFNHLYsFlBrnKoKJXAQOCYgSZFJFAZJAAQSAMEgaBkJB0P78/9mnodCpd1dVV\ndU51vZ+1enWfrl3nvOlUvbX3PnuQbUIIoae2vAMIIRRTJIcQQkmRHEIIJUVyCCGUFMkhhFBSJIcQ\nQkllk4OkSyQ9J+n+FTz+YUmzJd0n6TZJk2ofZgih0SqpOVwKTO3l8ceA99ieBJwKXFCLwEII+Sqb\nHGzfCrzUy+N32H4lO7wL2KBGsYUQclTrPodjgJk1PmcIIQeDanUiSbsDRwPvrtU5Qwj5qUlyyDoh\nLwSm2i7ZBJEUkzhCyIlt9fU5/U4OkjYEfgUcbntub2WrCbCeJM2wPSPvOLorYkxQzLgipspU+8Fc\nNjlIuhLYDRgp6UlgOjAYwPb5wCnACOBcSQBLbE+uJpgQQu1IbAEcW+3zyyYH29PKPP4x4GPVBhBC\nqB2JocDBpKQwjjQUoSqtPkJyVt4BlDAr7wBWYFbeAZQwK+8ASpjV6AtKGiJd/l7pX1eAnwQ+BHwf\nGGvz9arP26jFXiS5aH0OITQTSYLBp8BqW0HbeBi0AbxvHTgBGGt4/n9hm1Ns/tXjeVW992p2KzOE\nUD9S++dgrcPhzU1hj1VgF+Bw4C/A6cBt29n/vqeW14zkEEKBSWoD1oPVVoIDJsFx7TAWuATYDnju\nFVi0re1Ha33tSA4hFJCkY2H16TB+FHwSOKId5gh+BNwBrAUseAIW7Wj72XrEEMkhhIKRaIPfL4bX\n14d9Be3AVZ1w5Ouw2srw+svw/L9g8R7d5jXVXKvfrQihMCRGSHwWFj0GIy+D6wQfMJzZAR9bAIsn\nw8t/gEVbwuKd65kYIO5WhJCrdAfiZ++FSSfBuN3hxiXw/WFwO8BnYcgPYPAL8PqOtp9Q9kbq4zXi\nbkUIzUJa/Rj46DfhtvXTWKXfLoYPDoF5r8Kij9i+JpXT2vDmubafAuhrYuiPSA4hNJDEOsBxsOBE\nuGsdGCtYH3h9MDz5HCzazvb8rvK2v5ZXrNHnEEIDSGwrcRnwMLwxAXZ/Bt4j2MdwQSd87SV4fafu\niSFv0ecQQp1IDAIOBE4CNoA3zoeJO8DT+2dFfgQrfxLaXoeFO5Sb1Vx9HNHnEEIhSKxJmvj0SWAe\nvHEmjNgNFn8nK3Em+Iu2l0rt/4HOX9QrMfRH1BxCqBGJTUm1hEOBa2DRWTB0H/C3Uom2K6HzKNuL\nGxtX1BxCaDgJAXsCnwO2Bc6DhZvBqgcAd6dSbX+Czv3sjtdyC7QKkRxCqILESsBhwGeBNlh6Jmz1\ndXhoZ+CZVKrtYejc2e54Mb9IqxfNihD6QGIt4OPgT8KzT8Hp/4SfbAyvTYIRnfD8qtD+ACx5f9fY\nhLxV+96LW5khVEBigsQ5wFxgI9B7YfwcuO4QuGQH+NcqcNQQWPV+WLJdURJDf0RyCKEXEjtK/JI0\nFfJlYHObo23mwOLrYG47LAJmGs56EV7dy/Yb+UZdG9GsCKGHNCuSfYEvkXZw+2/gEpvX0uNaEzQP\nvHp6xvjF8NwSWDjZ9oP5RL1i0awIoZ8khkgcBcwBvgmcA0y0OcvmNUltUtsPgQVZYtgbGARPz4FF\n/1XExNAfUXMILU9ideA40p2HB0iLs95s47fLaDfeWjy27SLoPMF2R/bYKkVuSsQ4hxD6SGJt4NOk\nFVr/AOxvc8+yZTQC9BgwHDCwrt3xfPcyRU4M/RHNitByJDaUOAt4GFgb2Mnm0O6JIWtCnAG8CB4O\n7Gu7zfbzKzjtgBM1h9AyJDYBvgwcAFwEbGF3DVjqXk67ArekikLbpdB5bFcTopVEcggDnsQ7ga8C\nU4CzgQk2y41alDQ8a0KMyH61rt3xXMMCLZhoVoQBS2IHiWuBmcBdwEY23+qZGCRJajsdeAk8AtjP\ntmy3bGKAqDmEAUhiV+AbwCakOw8fsllUuqx2AW7NmhCXQecxrdiEKKVszUHSJZKek3R/L2XOkvSI\npNmStqltiCGUJyGJPSRmAT8GfkYao/A/pRKDpOFS2wvArdmv1rM7jozE8LZKmhWXAlNX9KCkfYAJ\ntieS7hWfW6PYQigrSwp7AbcA55G2gtrE5mKbN5cvL0lt3yM1IdYCPpA1IeqyMUwzK9ussH2rpHG9\nFNkfuCwre1fKyBrV6u21UF/d1lGYAYwETgWuslnhJ7+k/YHfZE2I/4XOo6OmsGK16HMYDTzZ7fgp\n0nj0SA6h5vqaFCRNgMGHwKpHAhOg/QXo2MruiJpCGbXqkOw5NLNha+uH1iGxG/AtYL3s+5WlkoKk\ndWDICbDKEbDGenBIG8wbDHfeDv/Z1XZno2NvRrVIDk8DY7odb5D9bjmSZnQ7nGV7Vg2uHwY4iZ1J\nNYRxpKRwuc3SXp6yEIYcA0dsmDaevQeY8iosPLgVEoOkKaQxHf1ju+wX6T/l/hU8tg8wM/t5R+DO\nFZRzJdeKr/jq+gJvC54JfgL8MfDgyp6HQDcChj8bNnwN2g/N+9+T398RV/O8sjUHSVcCuwEjJT0J\nTAcGZ1c83/ZMSftImgu8DhzV74wVWprE5qQawk7Ad4EP2lS0YrOkYcArqWWrW2HKu2HojfbSq+oX\n8cAUU7ZDYUiMI62jsDfwA+Acm4WVP1+7Azdnh1sB/wC+B3zP9ku1jbZ5VPvei+QQcicxCvg6aTXn\ns4Ef2vyn8udL0PZj6Pwo6AXwBm7w3hBFFitBhaYjsYbEqaQFVpYCm9pM72NiGA50Zonhy3bn2pEY\naiPmVoSGy/Z8+ARp+vT1wLtsnuj7ebQXcFN2uLndOaCWactbJIfQMNnCrYeRbkvOAfa0mdP380jQ\nfgVwKLQ9DZ0b2V5uqHTon0gOoSGy+Q8/ABYDR9jcUt15tCawADoAfd7u+FENwwzdRHIIdSWxFSkp\nbExacOVqu7oRtJL2Jq3NALCJ3fnP2kQZSokOyVAXEutJXAT8nvSG3sLmF9UkhjSTsv3qdJ62ecBK\ntiMx1Fkkh1BTEkMlTiH1KbxImj59lktMn67sfFqLdDfiINCn7Y7x0b/QGNGsCDXRrbPxNOB2YDub\nx/t3Tu0HXJsdTrA7H+1flKEvIjmEfssmRp1JGrN8iM3t/Tuf2qD9N8B+0PZP6NzS9pJaxBoqF8kh\nVE1iDGmNxl1JnY2X2/Rr1qOktYHns7sRH7c7zut/pKEa0ecQ+kxiFYlvAPeStqTf1OYnNUgMHwS6\nNo3ZyO6MxJCjqDmEimWrMH0Q+CHwN2Bbm3n9P6/aoP164H3Q9g/ofKft3tZrCA0QySFURGIz4CzS\nKkzH2G/NfuzneTUKeDZrRhxrd1xUi/OG/otmReiVxOoSPyCt7nwt8M4aJoYPAV1rOY61OyMxFEgk\nh1BStuT7ocCDpM1mt8zGK/S7up82qR10M/BzaP8bMNj2v/p73lBb0awIy5HYFDiHtLrzoTZ/rt25\ntR4wP2tGHGkvvaxW5w61FTWH8JbsLsS3gT+TmhDb1jgxTAPmZ4dj7M5IDAUWNYcAgMRUUm3hbmCS\n/dabuAbnVju03QFsD+13Qscujs1kCi+SQ4uTWA/4b2B74JM2N9T2/JoIg26BpesCH7aXXlHL84f6\niWZFi5JokzgeuA94DNiqlolB0hhpjZ9C+4PAWsAWtiMxNJGoObSgbMzChUA7sHs1qzGt+NwaBavN\ngKFHwAmDYH4n/PpC+9UHanWN0Bix+nQLkRhCWrfxU6S9Js/t75DnZc+voTBoDmw5Bq4flHLP2EXw\nxgTbJXdBC/VX7Xsvag4tQmIycDEwj7Sg65O9P6PvbC+U9H64959wqGGnDhj080gMzSlqDgOcxFDS\n7lGHA58l7Uhdl//0tPCrFoBHwEovQ+eqsGRT24/V43qhMlFzCMuR2JVUW7iH1OH47/pdS4K2x6Fz\nBLALLJ4HOiYSQ/OKmsMAlNUWvgscTLo9+X/1vZ4EbbOhcytgL9t/qOf1Qt/EjlcBeGtVpnuBdUi1\nhbomhqT9T1li2C8Sw8ARzYoBQmJlUt/CR4BPNCYpgDToOujYFTjY9m8bcc3QGGVrDpKmSnpI0iOS\nTi7x+EhJN0i6V9IcSUfWJdKwQhLbkIY9b0wa+tyoxHAldOwLOsL2LxpxzdA4vfY5pDHxPAzsBTwN\n/BWYZvvBbmVmkPYR+IqkkVn5UT1X8ok+h9qTaAdOBj4DfB74ab3uRCx/7UEXQMexoE/Ynec24pqh\nOvW6WzEZmGt7XnaRq4ADSHP8uzwDTMp+HgYsiCW+6k9iPPAT4E3S7Mmaj1tY8bXbz4DOY0FfjMQw\ncJVrVoyGZV50T2W/6+5CYAtJ84HZwEm1Cy/0lC3C8lHgL8CvgL0anBhmQOfnoe2bducZjbpuaLxy\nNYdKqqhfBe61PUXSxsBNkra2/WrPglkTpMss27MqjjQgMRw4F9iKtEP1fY29ftvnwdOh7Yd2x4xG\nXjtUTtIUYEp/z1MuOTwNjOl2PIZUe+huZ+A7ALYflfQ4sAmpg2wZtmdUHWmLk3g3cDlwHbC9zRuN\nvX7bx8FnQPuF9tIvNPLaoW+yD91ZXceSpldznnLNiruBiZLGSRoCHAJc06PMQ6QOy66VhDchTQEO\nNSDRLvF14JfAp2xOzCExfBT8P9B+lb30uEZeO+Sn15qD7aWSTgR+R5pid7HtByUdnz1+Pmkk3qWS\nZpOSzZdsv1jnuFuCxPrATwGROh0bPoFJ0v8DLoP2mfbSaY2+fshPDJ8uKIn3Az8m9TF8x6bhy6q9\nvZFt+2320l0aff1QG9W+9yI5FIzEIOCbwBHA4fbbbcfGxqE9gD9A2xzonORGvVBCzcWszAEgW8/x\nSmAJac2F58s8pU5xaGdSYngiEkPriolXBSGxG6kD+I/A1BwTw7bAbaCXoHN8JIbWFc2KnGWb036B\nNPz5ozY35heLtgTuBzpJu1DVbAm5kJ9oVjQhidWBS4CxwGSb3LaEk/QOUmIAGBKJIUSzIicS7wDu\nBF4Cds05MYwjTZiDVGOIDWdCJIc8SOxL2nLuTJvjbBbnF4tGA49nhyvFpLnQJZoVDZT1L3wF+ARw\ngM0d+cajdXh7OPwqtt/MM55QLJEcGiRb1/ESYDypf6Fme1FWF4/WBJ7LDle1vSjPeELxRLOiASRG\nA7eQxi/sVoDEMAxYkB0Os70wz3hCMUVyqDOJ7YC7gKtJtypz/YSWtCrwSnY4otTU+hAgmhV1JXEQ\ncB5wrM2v849HqwCvZYdr2345z3hCsUVyqIOs4/GLpD0p329zT84hIWkloKv5sJ7tF/KMJxRfJIca\nyyZOnQPsAOxkL7c4TsNJGgxvNWc2tP1snvGE5hDJoYYkVgN+TurL2dUm9/a8pEGkRWgBNrbdsPUm\nQ3OLDskakRhFWpprPvCBgiQGwaC52eHmsW9l6ItIDjUgMRG4HbiW1Pm4JOeQsqbE6lfA0rHAgd33\nGgmhEjErs58ktiUt+jrd5oK844GucQyrz4TNt4W/tcPSoTEsunXFRro5kNgduJ60N2VREsNoWO0e\nOHg7+NbKMOzhSAyhGpEcqiRxAPAz4JBG7U1ZTpp2vfLf4bBxcOFK8DfDoj/lHVdoTpEcqiDxYdLg\npn1s/ph3PN08Dp3fgwvaYazh+oWw8La8gwrNKZJDH0kcB3yftOPUchv35Mn2Enjz9nT0/B1w61DS\ntnkh9Fl0SPaBxEnAZ0mJ4dG84+kp3bqkE8C2JG1D2qow1oFsYbFMXJ1JfAk4jjSr8om84ylNn8+2\nNx0LYPvv+cYTmlnUHCog8RXgKGD3PHadqkQ2qWohtN1sd+yZdzyhOOJWZp1kieFIYEpRE0PS9vv0\nvXPffOMIA0U0K3qRNSWOJCWGZ3IOZ4UkbQTsDPq03RkrOoWaiGbFCkh8BjiR1MdQ4BpD+ttmP7ZF\n52PoqW7NCklTJT0k6RFJJ6+gzBRJf5c0R9KsvgZRNBLHA58B9miCxHBQ9uN2kRhCLfVac5DUTtrP\nYC/gaeCvwLTuk3gkDQduA95v+ylJI0stJNIsNQeJw4HTSE2Jwt2u7C77/1kKmm93js47nlBM9ao5\nTAbm2p6XBthwFXBAjzKHAb+0/RRAM68wJLE/cAZp9aZCJ4ZE56bv3iLfOMJAVC45jAa6Lw7yVPa7\n7iYCa0r6o6S7JX2klgE2SraR7UWktRgeyDuectLS8j4W2s6JtSBDPZS7W1FJG3Yw8C5gT2AocIek\nO20/0rOgpBndDmfZnlVhnHUl8U7gF8ChNn/NO57K6JH039N5Ut6RhGKRNAWY0t/zlEsOTwNjuh2P\ngeXWRHwSeMH2G8Abkm4BtgaWSw62Z1Qfan1IjAd+S5p2fXPe8VRC0g7AmsD+sa9l6Cn70J3VdSxp\nejXnKdesuBuYKGmcpCHAIcA1Pcr8BthFUrukoaSFVQtfLQeQWAu4ATjN5uq846lENn/iTgDb1+Yc\nThjAeq052F4q6UTgd0A7cLHtByUdnz1+vu2HJN0A3Eea9HOh7cInB4mVSYnu1zZn5x1P5ZadPxFC\nvbTkICiJNuCK7PAwO81kLLqYPxGqEbMy++ZUUv/Jns2SGJK2m1LlLOZPhPprueSQDXKaBuyQ976V\nfZHNn3g36KSYPxEaoaWaFRI7kvoZ9rCZk2csfRXzJ0K1Ysp2GRLrk3a6ProJE0PMnwgN1xI1B4mV\ngD8B19h8N48YqtVt/sQzduf6eccTmk/UHHp3Nmnw1ml5B9J3b82f2DzfOEKrGfAdkhJHA+8mdUA2\nVZU8zZ8gmz/REfMnQkMN6GaFxDbAjcB7bJpur0ipbQF4TWBQDJMO1YpxDj1IDCNNpvpUcyaGt+ZP\nHBCJIeRhQNYcJERae+IlmxMacc1a6rn/RM7hhCYXNYdlHQtsCuyYdyBVOjX7Pi7PIEJrG3A1B4nN\ngFuAXW0eqvf1aq3b/Ik/2h175B1PaH5xK5O3xjNcCXy1GRNDZuX0bdCN+YYRWt2AqjlIfJ+0bN1B\nzXbbskuaQ7HKQ8ASeOMdtgu9+nUovpavOUjsAnwEOL5ZE0NmPdhoIXxmMAw7L+9gQusaEMlBYjXg\nMuAEm3/nHU8/rQejBV8fDKvsIel9eQcUWtOASA6kYdG32sstYdeM1oUNh8AzwLShsNqlklbKO6jQ\nepr+VqbErsCBwJZ5x1IbbevDTwenVfIBhvwYWAuYn19MoRU1dYekxCrAbOCLNr+p5bnzImknYDHw\nN4hBUKH/WrVD8mvA7IGSGABs32H7HuD7eccSWlvT1hy6DXba2h54VW5JE0h7f6xr+7m84wnNq6Vq\nDtncifOA6QMxMWS69uo8LNcoQstqyuRAesOsCpyfdyD18vZycPpSvpGEVtV0yUFidVJ7/ESbgT6V\n+TrwunkHEVpT0yUH4BvA7+20JdwAdyaApJXzDiS0nqbqkJTYCPgrsIXNs7WJrLiy/UkXA3vbviHv\neEJzqva912zJ4WfA/TbfrlFYhZftV3GT7QE9jFpSG7A6MJw06GtD0noW47Ovrp9X6/a0bWzf29BA\nm9CAX+xFYidgZ+CovGNpLD0Bfm/eUVQia/6MyL7WY9k39XjQRuBR/bzKS+DHSTXIZp9HU2hlk4Ok\nqaS2bztwke2Sg3MkbQ/cARxs+1e1DDK7dXk68A2bhbU8d/H5dOAcZem/67eStgL2s12z5fazPTLW\nIL25R7LsJ3f3N3qVHyoGWAI83u1rXrfvLwAvA6/Eupn567VZkb1YHgb2Ap4mZetpth8sUe4mYCFw\nqe1fljhX1c0Kib2BM4BJLXCHYhmSRpI+ITcD5oE+BGt8ATrfAUvetBeu0a2sSLd4uz69x7D8p/d4\n8PB+RvUs+DGWf4M/A7wEvGw79vMsiHo1KyYDc23Pyy5yFXAALLea86dIW81t39cAyslqDd8GTmm1\nxABg+4X0nm+/EoZMgO2Az60GewLDVu62h2alZwR4lWXf1F1fTwILSJ/er8bWe62tXHIYTXrBdHkK\n2KF7AUmjSQljD1JyqPUL6sDse02bKs0iDaMe/ABoIrQNhg07UgXtrdbVFaRk3fUGf563q+ZLcgk6\nDAjlkkMlb/QzgS/bdlatXWH1RdKMboezbM/q7cRZreEU4GtNvrpT1WzPBbZIf9s3N4Er9oTrPwiv\n7wysAnzR9kAdQh6qIGkKMKXf5ynT57AjMMP21Oz4K0Bn905JSY/xdkIYSfpIO9b2NT3O1ed2j8QH\ngG8B72rV5LAikgYB7wIesv2fvOMJxVWXcQ7ZC/BhUgN3PvAXSnRIdit/KXBtqbsVfQ0wqzXcCfzA\n5upKnxdCWFZdOiRtL5V0IvA70q3Mi20/KOn47PF6TnzaExhGi/Y1hJC3wo6QlLgB+LnNJXUMK4QB\nb0CNkJTYAtiadBckhJCDos7K/Axwrs3ivAMJoVUVrlkhsTbwT2ATm+frH1kIA9tAWibuGOD/IjGE\nkK9C9TlItAHHAtPyjiWEVle0msOepHH/f807kBBaXdGSw3HA+TEaMoT8FaZDUmIUaTTmWJtXGhJU\nCC1gIHRIHgpcE4khhGIoUnI4DLg87yBCCEkhkoPERGAs8Ie8YwkhJIVIDqRaw89sluYdSAghyT05\nZFOzo0kRQsHknhyArYAhxNiGEAqlCMlhf+A3MbYhhGIpSnK4pmypEEJD5ToISmJ9YA4wyiZWSg6h\nDpp1ENR+wA2RGEIonryTQzQpQiio3JoVEiuR9kYcY/NyQ4IIoQU1Y7NiMvBwJIYQiinP5LAbMCvH\n64cQepFncpgC/CnH64cQepFLn4PEENJuztHfEEKdNVufw/bAI5EYQiiuvJJD9DeEUHB5Jofobwih\nwBqeHLIp2tsBdzX62iGEylWUHCRNlfSQpEcknVzi8Q9Lmi3pPkm3SZrUy+nWBzqA56qMOYTQAGWT\ng6R24GxgKrA5ME3SZj2KPQa8x/Yk4FTggl5OOQm4L6Zoh1BsldQcJgNzbc+zvQS4ih67X9u+w3bX\nqtF3ARv0cr5JwH3VBBtCaJxKksNo4Mlux09lv1uRY4CZvTweySGEJlBJcqi4+i9pd+BoYLl+iW62\nBmZXes4QQj4q2Uj3aWBMt+MxpNrDMrJOyAuBqbZfKnUiaaVT4eRN4PSDpMVr2J5VRcwhhF5ImkKa\nntC/85QbPi1pEGmbuj2B+cBfgGm2H+xWZkPgZuBw23eu4DwGbwNcbrNFfwMPIVSm2uHTZWsOtpdK\nOhH4HdAOXGz7QUnHZ4+fD5wCjADOlQSwxPbkEqeL/oYQmkRDJ16BfwgssDmtIRcNITTNxKsJpCZK\nCKHgGp0chkPMxAyhGTQ6OQwDXilbKoSQu0YnhzWI5BBCU4jkEEIoKY/k8J8GXzOEUIVGJ4cOm8UN\nvmYIoQqNTg7RpAihSURyCCGUFMkhhFBSJIcQQkmNTg5xpyKEJhE1hxBCSZEcQgglRXIIIZQUfQ4h\nhJKi5hBCKCmSQwihpEgOIYSSIjmEEEqKDskQQklRcwghlBTJIYRQUqOTw8IGXy+EUKWGJge78k15\nQwj5anTNIYTQJCI5hBBKiuQQQiipbHKQNFXSQ5IekXTyCsqclT0+W9I2tQ8zhNBovSYHSe3A2cBU\nYHNgmqTNepTZB5hgeyJwHHBunWKtOUlT8o6hpyLGBMWMK2Kqr3I1h8nAXNvzbC8BrgIO6FFmf+Ay\nANt3AcMljap5pPUxJe8ASpiSdwArMCXvAEqYkncAJUzJO4BaKZccRgNPdjt+KvtduTIb9D+0EEKe\nyiWHSsclqMrnhRAKalCZx58GxnQ7HkOqGfRWZoPsd8uRVLikIWl63jH0VMSYoJhxRUz1Uy453A1M\nlDQOmA8cAkzrUeYa4ETgKkk7Ai/bfq7niWz3rF2EEAqs1+Rge6mkE4HfAe3AxbYflHR89vj5tmdK\n2kfSXOB14Ki6Rx1CqDvZhavphxAKoOYjJIs4aKpcTJI+nMVyn6TbJE3KO6Zu5baXtFTSgUWISdIU\nSX+XNEfSrLxjkjRS0g2S7s1iOrIBMV0i6TlJ9/dSptGv8V5jquo1brtmX6Smx1xgHDAYuBfYrEeZ\nfYCZ2c87AHfWMoYqY9oJWCP7eWoRYupW7mbgOuCgvGMChgP/ADbIjkcWIKYZwGld8QALgEF1jmtX\nYBvg/hU83tDXeIUx9fk1XuuaQxEHTZWNyfYdtrsWormL+o/TqOTvBPAp4Grg33WOp9KYDgN+afsp\nANsvFCCmZ4Bh2c/DgAW2l9YzKNu3Ai/1UqThAwPLxVTNa7zWyaGIg6Yqiam7Y4CZdYwHKohJ0mjS\nG6FrOHq9O4cq+TtNBNaU9EdJd0v6SAFiuhDYQtJ8YDZwUp1jqkTRBwZW9Bovdyuzr4o4aKric0va\nHTgaeHf9wgEqi+lM4Mu2LUks/zfLI6bBwLuAPYGhwB2S7rT9SI4xfRW41/YUSRsDN0na2vardYqp\nUoUcGNiX13itk0NNB001MCayDpoLgam2e6syNiqmbUljRyC1pfeWtMT2NTnG9CTwgu03gDck3QJs\nDdQrOVQS087AdwBsPyrpcWAT0hidvDT6NV6RPr/Ga9wpMgh4lNSBNITyHZI7Uv/Ov0pi2pDU8bVj\nvTuOKo2pR/lLgQPzjgnYFPg9qaNwKHA/sHnOMf0ImJ79PIqUPNZswP/hOCrrkKz7a7zCmPr8Gq9p\nzcEFHDRVSUzAKcAI4Nzsk3qJ7ck5x9RQFf7fPSTpBuA+oBO40PYDecYEfBe4VNJsUh/al2y/WK+Y\nACRdCewGjJT0JDCd1OTK5TVeSUxU8RqPQVAhhJJimbgQQkmRHEIIJUVyCCGUFMkhhFBSJIcQQkmR\nHEIIJUVyCCGUFMkhhFDS/wfEYiTLwh5xBwAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x108416ba8>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "# listing04.py\n", | |
| "\n", | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "from listing01 import f, G\n", | |
| "\n", | |
| "duration = 10\n", | |
| "x0 = [0.8, 0.43]\n", | |
| "\n", | |
| "fig = plt.figure()\n", | |
| "ax = fig.add_subplot(111, aspect='equal')\n", | |
| "\n", | |
| "grids = np.linspace(0.0, 1.2, 120)\n", | |
| "ax.plot(grids, f(grids))\n", | |
| "\n", | |
| "\n", | |
| "x = x0[:]\n", | |
| "for t in range(duration):\n", | |
| " x1 = G(x)\n", | |
| " dx = [x1[0] - x[0], x1[1] - x[1]]\n", | |
| " ax.plot(x[0], x[1], marker='', linestyle='', color='black')\n", | |
| " ax.arrow(x[0], x[1], dx[0], dx[1], length_includes_head=True)\n", | |
| " x = x1\n", | |
| "plt.show()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "\n", | |
| "## OOP\n", | |
| "\n", | |
| "対象と動作をまとめてやる. たとえば次のような記法でシミュレーションができるようにしたい. \n", | |
| "\n", | |
| " >>> ramsey = Ramsey(A=1.1, α=0.4, ρ=0.9)\n", | |
| " >>> sim = Simulation(ramsey, x0=[0.8, 0.43], duration=5)\n", | |
| " >>> path = [x for x in sim]\n", | |
| " >>> path\n", | |
| " [[0.8, 0.43],\n", | |
| " [0.43, 0.4063168783383537],\n", | |
| " [0.4063168783383537, 0.5099372470852112],\n", | |
| " [0.5099372470852112, 0.6875981937381646],\n", | |
| " [0.6875981937381646, 0.8712750948974344]]\n", | |
| " " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 13, | |
| "metadata": { | |
| "collapsed": true, | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# listing05.py の一部\n", | |
| "\n", | |
| "class Ramsey:\n", | |
| " \"\"\"One-sector Ramsey model\"\"\"\n", | |
| "\n", | |
| " def __init__(self, A, α, ρ):\n", | |
| " self.A = A\n", | |
| " self.α = α\n", | |
| " self.ρ = ρ\n", | |
| "\n", | |
| " def f(self, x):\n", | |
| " return self.A * x ** self.α\n", | |
| "\n", | |
| " def U(self, x):\n", | |
| " return log(x)\n", | |
| "\n", | |
| " def u(self, x, y):\n", | |
| " return self.U(self.f(x) - y)\n", | |
| "\n", | |
| " def forward(self, x):\n", | |
| " \"\"\"1ステップの時間発展\"\"\"\n", | |
| " A, α, ρ = self.A, self.α, self.ρ\n", | |
| " return [\n", | |
| " x[1],\n", | |
| " (1 + ρ * α) * A * x[1] ** α - ρ * α * (A ** 2) * (x[1] ** (α - 1)) * (x[0] ** α)\n", | |
| " ]" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 14, | |
| "metadata": { | |
| "collapsed": true, | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# listing05.py の一部\n", | |
| "\n", | |
| "class Simulation:\n", | |
| " \"\"\"Simulation of a dynamical system\"\"\"\n", | |
| "\n", | |
| " def __init__(self, sys, x0, duration):\n", | |
| " self.sys = sys\n", | |
| " self.x0 = x0\n", | |
| " self.duration = duration\n", | |
| "\n", | |
| " def __iter__(self):\n", | |
| " x = self.x0[:]\n", | |
| " for _ in range(self.duration):\n", | |
| " yield x\n", | |
| " x = self.sys.forward(x)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 15, | |
| "metadata": { | |
| "collapsed": false, | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "[[0.8, 0.43],\n", | |
| " [0.43, 0.4063168783383537],\n", | |
| " [0.4063168783383537, 0.5099372470852112],\n", | |
| " [0.5099372470852112, 0.6875981937381646],\n", | |
| " [0.6875981937381646, 0.8712750948974344]]" | |
| ] | |
| }, | |
| "execution_count": 15, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "ramsey = Ramsey(A=1.1, α=0.4, ρ=0.9)\n", | |
| "sim = Simulation(ramsey, x0=[0.8, 0.43], duration=5)\n", | |
| "\n", | |
| "path = [x for x in sim]\n", | |
| "path" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 16, | |
| "metadata": { | |
| "collapsed": false, | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAAEACAYAAAC+rrMfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH41JREFUeJzt3XmYXFW57/HvrzsJECAkEAwQAgESmSSIQJilGTyGQbjC\nFQiiTDKoKA5HcSRRVFTQw8OFE2bkqoCKHgUMCIoRZFJAAhGCRAgCYZAwEwhJ9+/8sXZDp1Ppqq6u\nqr2r6/08Tz/du2vV3m86VW+ttfYaZJsQQuitLe8AQgjFFMkhhFBSJIcQQkmRHEIIJUVyCCGUFMkh\nhFBS2eQg6RJJz0i6fwWPf1jSbEn3SbpV0qTahxlCaLRKag6XAlP6ePwR4L22JwGnARfUIrAQQr7K\nJgfbtwAv9PH47bZfyg7vBNavUWwhhBzVus/hWGBmjc8ZQsjBkFqdSNIewDHALrU6ZwghPzVJDlkn\n5IXAFNslmyCSYhJHCDmxrf4+Z8DJQdIGwK+AI2zP66tsNQHWk6TptqfnHUdPRYwJihlXxFSZaj+Y\nyyYHSVcAuwOjJT0OTAOGAtg+HzgVGAXMkASwxPbkaoIJIdSOxJbAcdU+v2xysD21zOMfAz5WbQAh\nhNqRGA4cQkoK40lDEarS6iMkZ+UdQAmz8g5gBWblHUAJs/IOoIRZeVxUYiuJ/wc8DnwI+B6woc3X\nqj5noxZ7keSi9TmE0IyU2u9rweabwLRD4b0fgHVWAV0CXGTzr17lq3rvRXIIoQlI7Z+DUUfAkvVg\n3FrwccGh7XB3J7z4dTj0DJulpZ9b3XuvZuMcQgi1J6kNWBdWWwkOnAQntsMkYC6w+8vwQIftv9Xl\n2lFzCKF4JB0Hq0+D0WPgk8CR7TBHMAPoMly3EF7b1fZDFZwrag4hDAYSbfD7xfDaerCfoB24sguO\neg1WWxneeAZe28X2Y/WMo9XvVoRQGBKjJD4LbzwCoy+DawUfMJzVCR9bCIsnw4t/gNe2rXdigGhW\nhJCrdOfhZ++DSSfD+D3ghiXwvRFwG8BnYdgZMPQ5eG1H248peyP18xrRrAihWUirHwsf/Qbcul4a\nq/TbxfDBYTD/FXjjI7avTuW0Nrw5w/YTAP1NDAMRySGEBpJ4B3A8LDwJ7nwHbChYD3htKDz+DLyx\nne0F3eVtfzWvWKPPIYQGkNhW4jLgIXh9AuzxFLxXsK/hgi746gvw2k49E0Peos8hhDqRGAIcBJwM\nrA+vnw8Td4AnD8iK/BBW/iS0vQaLdig3q7n6OKLPIYRCkFiTNPHpk8B8eP0sGLU7LP52VuIs8Bds\nL5XaX4auX9QrMQxE1BxCqBGJzUi1hMOAq+GNs2H4vuBvphJtV0DX0bYXNzauqDmE0HASAvYCPgds\nC5wHizaHVQ8E7kql2v4EXfvbna/mFmgVIjmEUAWJlYDDgc8CbbD0LNjqazB3Z+CpVKrtIeja2e58\nPr9IqxfNihD6QWIt4OPgT8LTT8D3/wE/3gRenQSjuuDZVaH9AVjy/u6xCXmr9r0XtzJDqIDEBIlz\ngXnAxqD3wUZz4NpD4ZId4F+rwNHDYNX7Ycl2RUkMAxHJIYQ+SOwo8UvgduBFYAubY2zmwOJrYV47\nvAHMNJz9PLyyt+3X8426NqJZEUIvaVYk+wFfJO3g9l/AJTavpse1Jmg+ePX0jI0WwzNLYNFk2w/m\nE/WKRbMihAGSGCZxNDAH+AZwLjDR5mybVyW1SW0/ABZmiWEfYAg8OQfe+D9FTAwDETWH0PIkVgeO\nJ915eIC0OOtNNn67jHbnrcVj2y6CrhNtd2aPrVLkpkSMcwihnyTWBj4NnAj8ATjA5p5ly2gU6BFg\nJGBgHbvz2Z5lipwYBiKaFaHlSGwgcTbwELA2sJPNYT0TQ9aEOBN4HjwS2M92m+1nV3DaQSdqDqFl\nSGwKfAk4ELgI2NLuHrDUs5x2A25OFYW2S6HruO4mRCuJ5BAGPYl3A18BOoBzgAk2y41alDQya0KM\nyn61jt35TMMCLZhoVoRBS2IHiWuAmcCdwMY23+ydGCRJavs+8AJ4FLC/bdlu2cQAUXMIg5DEbsDX\ngU1Jdx4+ZPNG6bLaFbgla0JcBl3HtmITopSyNQdJl0h6RtL9fZQ5W9LDkmZL2qa2IYZQnoQk9pSY\nBfwI+BlpjMJ/l0oMkkZKbc8Bt2S/WtfuPCoSw9sqaVZcCkxZ0YOS9gUm2J5Iulc8o0axhVBWlhT2\nBm4GzgMuATa1udjmzeXLS1Lbd0lNiLWAD2RNiKcbG3nxlW1W2L5F0vg+ihwAXJaVvTNlZI1p9fZa\nqK8e6yhMB0YDpwFX2qzwk1/SAcBvsibE/4euY6KmsGK16HMYS9r2u9sTpPHokRxCzfU3KUiaAEMP\nhVWPAiZA+3PQuZXdGTWFMmrVIdl7aGbD1tYPrUNid+CbwLrZ9ytKJQVJ74BhJ8IqR8Ia68KhbTB/\nKNxxG7y8m+2uRsfejGqRHJ4ExvU4Xj/73XIkTe9xOMv2rBpcPwxyEjuTagjjSUnhpyvabj6zCIYd\nC0duAD8E7gE6XoFFh7RCYpDUQRrTMTC2y36R/lPuX8Fj+wIzs593BO5YQTlXcq34iq/uL/C24Jng\nx8AfAw+t7HkIdANg+LNhg1eh/bC8/z35/R1xNc8rW3OQdAWwOzBa0uPANGBodsXzbc+UtK+kecBr\nwNEDzlihpUlsQaoh7AR8B/igTUUrNksaAbyUWra6BTp2geE32EuvrF/Eg1NM2Q6FITGetI7CPsAZ\nwLk2iyp/vvYAbsoOtwL+DnwX+K7tF2obbfOo9r0XySHkTmIM8DXSas7nAD+webny50vQ9iPo+ijo\nOfD6bvDeEEUWK0GFpiOxhsRppAVWlgKb2UzrZ2IYCXRlieFLdtfakRhqI+ZWhIbL9nz4BGn69HXA\ne2we6/95tDdwY3a4hd01qJZpy1skh9Aw2cKth5NuS84B9rKZ0//zSNB+OXAYtD0JXRvbXm6odBiY\nSA6hIbL5D2cAi4EjbW6u7jxaE1gInYA+b3f+sIZhhh4iOYS6ktiKlBQ2IS24cpVd3QhaSfuQ1mYA\n2NTu+kdtogylRIdkqAuJdSUuAn5PekNvafOLahJDmknZflU6T9t8YCXbkRjqLJJDqCmJ4RKnkvoU\nnidNnz7bJaZPV3Y+rUW6G3Ew6NN250bRv9AY0awINdGjs/F04DZgO5tHB3ZO7Q9ckx1OsLv+ObAo\nQ39EcggDlk2MOos0ZvlQm9sGdj61QftvgP2h7R/Q9S7bS2oRa6hcJIdQNYlxpDUadyN1Nv7UZkCz\nHiWtDTyb3Y34uN153sAjDdWIPofQbxKrSHwduJe0Jf1mNj+uQWL4INC9aczGdlckhhxFzSFULFuF\n6YPAD4C7gW1t5g/8vGqD9uuA/4C2v0PXu233tV5DaIBIDqEiEpsDZ5NWYTrWfmv24wDPqzHA01kz\n4ji786JanDcMXDQrQp8kVpc4g7S68zXAu2uYGD4EdK/luKHdFYmhQCI5hJKyJd8PAx4kbTb7rmy8\nwoCr+2mT2iE3AT+H9ruBobb/NdDzhtqKZkVYjsRmwLmk1Z0Ps/lz7c6tdYEFWTPiKHvpZbU6d6it\nqDmEt2R3Ib4F/JnUhNi2xolhKrAgOxxnd0ViKLCoOQQAJKaQagt3AZPst97ENTi32qHtdmB7aL8D\nOnd1bCZTeJEcWpzEusB/AdsDn7S5vrbn10QYcjMsXQf4sL308lqeP9RPNCtalESbxAnAfcAjwFa1\nTAySxklr/ATaHwTWAra0HYmhiUTNoQVlYxYuBNqBPapZjWnF59YYWG06DD8SThwCC7rg1xfarzxQ\nq2uExojVp1uIxDDSuo2fIu01OWOgQ56XPb+Gw5A58K5xcN2QlHs2fANen2C75C5oof6qfe9FzaFF\nSEwGLgbmkxZ0fbzvZ/Sf7UWS3g/3/gMOM+zUCUN+HomhOUXNYZCTGE7aPeoI4LOkHanr8p+eFn7V\nQvAoWOlF6FoVlmxm+5F6XC9UJmoOYTkSu5FqC/eQOhz/Xb9rSdD2KHSNAnaFxfNBx0ZiaF5RcxiE\nstrCd4BDSLcn/6e+15OgbTZ0bQXsbfsP9bxe6J/Y8SoAb63KdC/wDlJtoa6JIWn/U5YY9o/EMHhE\ns2KQkFiZ1LfwEeATjUkKIA25Fjp3Aw6x/dtGXDM0Rtmag6QpkuZKeljSKSUeHy3pekn3Spoj6ai6\nRBpWSGIb0rDnTUhDnxuVGK6Azv1AR9r+RSOuGRqnzz6HNCaeh4C9gSeBvwJTbT/Yo8x00j4CX5Y0\nOis/pvdKPtHnUHsS7cApwGeAzwM/qdediOWvPeQC6DwO9Am7a0YjrhmqU6+7FZOBebbnZxe5EjiQ\nNMe/21PApOznEcDCWOKr/iQ2An4MvEmaPVnzcQsrvnb7mdB1HOgLkRgGr3LNirGwzIvuiex3PV0I\nbClpATAbOLl24YXeskVYPgr8BfgVsHeDE8N06Po8tH3D7jqzUdcNjVeu5lBJFfUrwL22OyRtAtwo\naWvbr/QumDVBus2yPaviSAMSI4EZwFakHarva+z12z4PngZtP7A7pzfy2qFykjqAjoGep1xyeBIY\n1+N4HKn20NPOwLcBbP9T0qPApqQOsmXYnl51pC1OYhfgp8C1wPY2rzf2+m0fB58J7RfaS/+zkdcO\n/ZN96M7qPpY0rZrzlGtW3AVMlDRe0jDgUODqXmXmkjosu1cS3pQ0BTjUgES7xNeAXwKfsjkph8Tw\nUfB/Q/uV9tLjG3ntkJ8+aw62l0o6CfgdaYrdxbYflHRC9vj5pJF4l0qaTUo2X7T9fJ3jbgkS6wE/\nAUTqdGz4BCZJ/xe4DNpn2kunNvr6IT8xfLqgJN4P/IjUx/Btm4Yvq/b2Rrbtt9pLd2309UNtVPve\ni+RQMBJDgG8ARwJH2G+3HRsbh/YE/gBtc6Brkhv1Qgk1F7MyB4FsPccrgCWkNReeLfOUOsWhnUmJ\n4bFIDK0rJl4VhMTupA7gPwJTckwM2wK3gl6Aro0iMbSuaFbkLNuc9j9Jw58/anNDfrHoXcD9QBdp\nF6qaLSEX8hPNiiYksTpwCbAhMNkmty3hJL2TlBgAhkViCNGsyInEO4E7gBeA3XJODONJE+Yg1Rhi\nw5kQySEPEvuRtpw7y+Z4m8X5xaKxwKPZ4UoxaS50i2ZFA2X9C18GPgEcaHN7vvHoHbw9HH4V22/m\nGU8olkgODZKt63gJsBGpf6Fme1FWF4/WBJ7JDle1/Uae8YTiiWZFA0iMBW4mjV/YvQCJYQSwMDsc\nYXtRnvGEYorkUGcS2wF3AleRblXm+gktaVXgpexwVKmp9SFANCvqSuJg4DzgOJtf5x+PVgFezQ7X\ntv1invGEYovkUAdZx+MXSHtSvt/mnpxDQtJKQHfzYV3bz+UZTyi+SA41lk2cOhfYAdjJXm5xnIaT\nNBTeas5sYPvpPOMJzSGSQw1JrAb8nNSXs5tN7u15SUNIi9ACbGK7YetNhuYWHZI1IjGGtDTXAuAD\nBUkMgiHzssMtYt/K0B+RHGpAYiJwG3ANqfNxSc4hZU2J1S+HpRsCB/XcaySESsSszAGS2Ja06Os0\nmwvyjge6xzGsPhO22Bbuboelw2NYdOuKjXRzILEHcB1pb8qiJIaxsNo9cMh28M2VYcRDkRhCNSI5\nVEniQOBnwKGN2puynDTteuW/weHj4cKV4G7DG3/KO67QnCI5VEHiw6TBTfva/DHveHp4FLq+Cxe0\nw4aG6xbBolvzDio0p0gO/SRxPPA90o5Ty23ckyfbS+DN29LRs7fDLcNJ2+aF0G/RIdkPEicDnyUl\nhn/mHU9v6dYlXQC2JWkb0laFsQ5kC4tl4upM4ovA8aRZlY/lHU9p+ny2vemGALb/lm88oZlFzaEC\nEl8Gjgb2yGPXqUpkk6oWQdtNdudeeccTiiNuZdZJlhiOAjqKmhiStt+n71375RtHGCyiWdGHrClx\nFCkxPJVzOCskaWNgZ9Cn7a5Y0SnURDQrVkDiM8BJpD6GAtcY0t82+7EtOh9Db3VrVkiaImmupIcl\nnbKCMh2S/iZpjqRZ/Q2iaCROAD4D7NkEieHg7MftIjGEWuqz5iCpnbSfwd7Ak8Bfgak9J/FIGgnc\nCrzf9hOSRpdaSKRZag4SRwCnk5oShbtd2VP2/7MUtMDuGpt3PKGY6lVzmAzMsz0/DbDhSuDAXmUO\nB35p+wmAZl5hSOIA4EzS6k2FTgyJZqTv3jLfOMJgVC45jAV6Lg7yRPa7niYCa0r6o6S7JH2klgE2\nSraR7UWktRgeyDuectLS8j4O2s6NtSBDPZS7W1FJG3Yo8B5gL2A4cLukO2w/3LugpOk9DmfZnlVh\nnHUl8W7gF8BhNn/NO57K6OH039N1ct6RhGKR1AF0DPQ85ZLDk8C4HsfjYLk1ER8HnrP9OvC6pJuB\nrYHlkoPt6dWHWh8SGwG/JU27vinveCohaQdgTeCA2Ncy9JZ96M7qPpY0rZrzlGtW3AVMlDRe0jDg\nUODqXmV+A+wqqV3ScNLCqoWvlgNIrAVcD5xuc1Xe8VQimz9xB4Dta3IOJwxifdYcbC+VdBLwO6Ad\nuNj2g5JOyB4/3/ZcSdcD95Em/Vxou/DJQWJlUqL7tc05ecdTuWXnT4RQLy05CEqiDbg8OzzcTjMZ\niy7mT4RqxKzM/jmN1H+yV7MkhqTtxlQ5i/kTof5aLjlkg5ymAjvkvW9lf2TzJ3YBnRzzJ0IjtFSz\nQmJHUj/DnjZz8oylv2L+RKhWTNkuQ2I90k7XxzRhYoj5E6HhWqLmILES8Cfgapvv5BFDtXrMn3jK\n7lov73hC84maQ9/OIQ3eOj3vQPrvrfkTW+QbR2g1g75DUuIYYBdSB2RTVcnT/Amy+ROdMX8iNNSg\nblZIbAPcALzXpun2ipTaFoLXBIbEMOlQrRjn0IvECNJkqk81Z2J4a/7EgZEYQh4GZc1BQqS1J16w\nObER16yl3vtP5BxOaHJRc1jWccBmwI55B1Kl07Lv4/MMIrS2QVdzkNgcuBnYzWZuva9Xaz3mT/zR\n7twz73hC84tbmbw1nuEK4CvNmBgyK6dvQ27IN4zQ6gZVzUHie6Rl6w5uttuW3dIcilXmAkvg9Xfa\nLvTq16H4Wr7mILEr8BHghGZNDJl1YeNF8JmhMOK8vIMJrWtQJAeJ1YDLgBNt/p13PAO0LowVfG0o\nrLKnpP/IO6DQmgZFciANi77FXm4Ju2a0DmwwDJ4Cpg6H1S6VtFLeQYXW0/S3MiV2Aw4C3pV3LLXR\nth78ZGhaJR9g2I+AtYAF+cUUWlFTd0hKrALMBr5g85tanjsvknYCFgN3QwyCCgPXqh2SXwVmD5bE\nAGD7dtv3AN/LO5bQ2pq25tBjsNPW9uCrckuaQNr7Yx3bz+QdT2heLVVzyOZOnAdMG4yJIdO9V+fh\nuUYRWlZTJgfSG2ZV4Py8A6mXt5eD0xfzjSS0qqZLDhKrk9rjJ9kM9qnM14LXyTuI0JqaLjkAXwd+\nb6ct4Qa5swAkrZx3IKH1NFWHpMTGwF+BLW2erk1kxZXtT7oY2Mf29XnHE5pTte+9ZksOPwPut/lW\njcIqvGy/ihttD+ph1JLagNWBkaRBXxuQ1rPYKPvq/nm1Hk/bxva9DQ20CQ36xV4kdgJ2Bo7OO5bG\n0mPg9+UdRSWy5s+o7Gtdln1TbwTaGDxmgFd5AfwoqQbZ7PNoCq1scpA0hdT2bQcusl1ycI6k7YHb\ngUNs/6qWQWa3Lr8PfN1mUS3PXXz+PnCusvTf/VtJWwH7267ZcvvZHhlrkN7co1n2k7vnG73KDxUD\nLAEe7fE1v8f354AXgZdi3cz89dmsyF4sDwF7A0+SsvVU2w+WKHcjsAi41PYvS5yr6maFxD7AmcCk\nFrhDsQxJo0mfkJsD80EfgjX+E7reCUvetBet0aOsSLd4uz+9x7H8p/dG4JEDjOpp8CMs/wZ/CngB\neNF27OdZEPVqVkwG5tmen13kSuBAWG4150+Rtprbvr8BlJPVGr4FnNpqiQHA9nPpPd9+BQybANsB\nn1sN9gJGrNxjD81KzwjwCsu+qbu/HgcWkj69X4mt91pbueQwlvSC6fYEsEPPApLGkhLGnqTkUOsX\n1EHZ95o2VZpFGkY99AHQRGgbCht0pgraW62ry0nJuvsN/ixvV82X5BJ0GBTKJYdK3uhnAV+y7axa\nu8Lqi6TpPQ5n2Z7V14mzWsOpwFebfHWnqtmeB2yZ/rZvbgqX7wXXfRBe2xlYBfiC7cE6hDxUQVIH\n0DHg85Tpc9gRmG57Snb8ZaCrZ6ekpEd4OyGMJn2kHWf76l7n6ne7R+IDwDeB97RqclgRSUOA9wBz\nbb+cdzyhuOoyziF7AT5EauAuAP5CiQ7JHuUvBa4pdbeivwFmtYY7gDNsrqr0eSGEZdWlQ9L2Ukkn\nAb8j3cq82PaDkk7IHq/nxKe9gBG0aF9DCHkr7AhJieuBn9tcUsewQhj0BtUISYktga1Jd0FCCDko\n6qzMzwAzbBbnHUgIrapwzQqJtYF/AJvaPFv/yEIY3AbTMnHHAv8TiSGEfBWqz0GiDTgOmJp3LCG0\nuqLVHPYijfv/a96BhNDqipYcjgfOj9GQIeSvMB2SEmNIozE3tHmpIUGF0AIGQ4fkYcDVkRhCKIYi\nJYfDgZ/mHUQIISlEcpCYCGwI/CHvWEIISSGSA6nW8DObpXkHEkJIck8O2dTsaFKEUDC5JwdgK2AY\nMbYhhEIpQnI4APhNjG0IoViKkhyuLlsqhNBQuQ6CklgPmAOMsYmVkkOog2YdBLU/cH0khhCKJ+/k\nEE2KEAoqt2aFxEqkvRHH2bzYkCBCaEHN2KyYDDwUiSGEYsozOewOzMrx+iGEPuSZHDqAP+V4/RBC\nH3Lpc5AYRtrNOfobQqizZutz2B54OBJDCMWVV3KI/oYQCi7P5BD9DSEUWMOTQzZFezvgzkZfO4RQ\nuYqSg6QpkuZKeljSKSUe/7Ck2ZLuk3SrpEl9nG49oBN4psqYQwgNUDY5SGoHzgGmAFsAUyVt3qvY\nI8B7bU8CTgMu6OOUk4D7Yop2CMVWSc1hMjDP9nzbS4Ar6bX7te3bbXevGn0nsH4f55sE3FdNsCGE\nxqkkOYwFHu9x/ET2uxU5FpjZx+ORHEJoApUkh4qr/5L2AI4BluuX6GFrYHal5wwh5KOSjXSfBMb1\nOB5Hqj0sI+uEvBCYYvuFUieSVjoNTtkUvn+wtHgN27OqiDmE0AdJHaTpCQM7T7nh05KGkLap2wtY\nAPwFmGr7wR5lNgBuAo6wfccKzmPwNsBPbbYcaOAhhMpUO3y6bM3B9lJJJwG/A9qBi20/KOmE7PHz\ngVOBUcAMSQBLbE8ucbrobwihSTR04hX4B8BCm9MbctEQQtNMvJpAaqKEEAqu0clhJMRMzBCaQaOT\nwwjgpbKlQgi5a3RyWINIDiE0hUgOIYSS8kgOLzf4miGEKjQ6OXTaLG7wNUMIVWh0cogmRQhNIpJD\nCKGkSA4hhJIiOYQQSmp0cog7FSE0iag5hBBKiuQQQigpkkMIoaTocwghlBQ1hxBCSZEcQgglRXII\nIZQUySGEUFJ0SIYQSoqaQwihpEgOIYSSGp0cFjX4eiGEKjU0OdiVb8obQshXo2sOIYQmEckhhFBS\nJIcQQkllk4OkKZLmSnpY0ikrKHN29vhsSdvUPswQQqP1mRwktQPnAFOALYCpkjbvVWZfYILticDx\nwIw6xVpzkjryjqG3IsYExYwrYqqvcjWHycA82/NtLwGuBA7sVeYA4DIA23cCIyWNqXmk9dGRdwAl\ndOQdwAp05B1ACR15B1BCR94B1Eq55DAWeLzH8RPZ78qVWX/goYUQ8lQuOVQ6LkFVPi+EUFBDyjz+\nJDCux/E4Us2grzLrZ79bjqTCJQ1J0/KOobcixgTFjCtiqp9yyeEuYKKk8cAC4FBgaq8yVwMnAVdK\n2hF40fYzvU9ku3ftIoRQYH0mB9tLJZ0E/A5oBy62/aCkE7LHz7c9U9K+kuYBrwFH1z3qEELdyS5c\nTT+EUAA1HyFZxEFT5WKS9OEslvsk3SppUt4x9Si3vaSlkg4qQkySOiT9TdIcSbPyjknSaEnXS7o3\ni+moBsR0iaRnJN3fR5lGv8b7jKmq17jtmn2Rmh7zgPHAUOBeYPNeZfYFZmY/7wDcUcsYqoxpJ2CN\n7OcpRYipR7mbgGuBg/OOCRgJ/B1YPzseXYCYpgOnd8cDLASG1Dmu3YBtgPtX8HhDX+MVxtTv13it\naw5FHDRVNibbt9vuXojmTuo/TqOSvxPAp4CrgH/XOZ5KYzoc+KXtJwBsP1eAmJ4CRmQ/jwAW2l5a\nz6Bs3wK80EeRhg8MLBdTNa/xWieHIg6aqiSmno4FZtYxHqggJkljSW+E7uHo9e4cquTvNBFYU9If\nJd0l6SMFiOlCYEtJC4DZwMl1jqkSRR8YWNFrvNytzP4q4qCpis8taQ/gGGCX+oUDVBbTWcCXbFuS\nWP5vlkdMQ4H3AHsBw4HbJd1h++EcY/oKcK/tDkmbADdK2tr2K3WKqVKFHBjYn9d4rZNDTQdNNTAm\nsg6aC4EptvuqMjYqpm1JY0cgtaX3kbTE9tU5xvQ48Jzt14HXJd0MbA3UKzlUEtPOwLcBbP9T0qPA\npqQxOnlp9Gu8Iv1+jde4U2QI8E9SB9IwyndI7kj9O/8qiWkDUsfXjvXuOKo0pl7lLwUOyjsmYDPg\n96SOwuHA/cAWOcf0Q2Ba9vMYUvJYswH/h+OprEOy7q/xCmPq92u8pjUHF3DQVCUxAacCo4AZ2Sf1\nEtuTc46poSr8v5sr6XrgPqALuND2A3nGBHwHuFTSbFIf2hdtP1+vmAAkXQHsDoyW9DgwjdTkyuU1\nXklMVPEaj0FQIYSSYpm4EEJJkRxCCCVFcgghlBTJIYRQUiSHEEJJkRxCCCVFcgghlBTJIYRQ0v8C\nWPY0HuAQTe0AAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x1085cbc88>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "# listing05.py の一部\n", | |
| "\n", | |
| "ramsey = Ramsey(A=1.1, α=0.4, ρ=0.9)\n", | |
| "sim = Simulation(ramsey, x0=[0.8, 0.43], duration=10)\n", | |
| "\n", | |
| "fig, ax = plt.subplots(subplot_kw={'aspect':'equal'})\n", | |
| "\n", | |
| "grids = np.linspace(0.0, 1.2, 120)\n", | |
| "ax.plot(grids, ramsey.f(grids))\n", | |
| "\n", | |
| "for i, x in enumerate(sim):\n", | |
| " if i == 0:\n", | |
| " ax.plot(x[0], x[1], marker='', linestyle='', color='black')\n", | |
| " else: \n", | |
| " dx = [x[0] - x0[0], x[1] - x0[1]]\n", | |
| " ax.plot(x[0], x[1], marker='', linestyle='', color='black')\n", | |
| " ax.arrow(x0[0], x0[1], dx[0], dx[1], length_includes_head=True) \n", | |
| " x0 = x[:]" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "初期値を変えてシミュレーションしたものを重ねることも簡単にできる. `Simulation` クラスを書き換える. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 17, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# listing06.py の一部\n", | |
| "\n", | |
| "class Simulation:\n", | |
| " \"\"\"__init__ と __iter__ は変更なし\"\"\"\n", | |
| "\n", | |
| " def reset(self, *, x0=None, duration=None):\n", | |
| " if x0 is not None:\n", | |
| " self.x0 = x0[:]\n", | |
| " if duration is not None:\n", | |
| " self.duration = duration\n", | |
| "\n", | |
| " def plot(self, ax):\n", | |
| " for i, x in enumerate(self):\n", | |
| " if i == 0:\n", | |
| " ax.plot(x[0], x[1], marker='', linestyle='', color='black')\n", | |
| " else:\n", | |
| " dx = [x[0] - x0[0], x[1] - x0[1]]\n", | |
| " ax.plot(x[0], x[1], marker='', linestyle='', color='black')\n", | |
| " ax.arrow(x0[0], x0[1], dx[0], dx[1], length_includes_head=True)\n", | |
| " x0 = x[:]" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 18, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAAEACAYAAAC+rrMfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYXFXRh99f98wkmawkgUBCIEAiECBsIWwGgqCGsCoi\nm+wIiEFB+ED2VUVc4BMwYoSAyqIsYkAEQQwCQpRPCCBhCZtJQJAEyEK26a7vj3Mn00x6unt6+t7u\nman3efqZvn3PPae653Z1nTp1qmRmOI7jtCZVbQEcx6lNXDk4jpMXVw6O4+TFlYPjOHlx5eA4Tl5c\nOTiOk5eiykHSjZLelfR8G+ePkDRL0nOSnpA0pvJiOo6TNKVYDtOAiQXOvw7sZmZjgMuAn1dCMMdx\nqktR5WBmjwEfFDj/pJl9FB3OBNavkGyO41SRSvscjgfur3CfjuNUgbpKdSRpD+A4YNdK9ek4TvWo\niHKInJBTgYlmlncKIsk3cThOlTAztfeaDisHSRsAdwNfMbM5hdqWI2CcSLrYzC6uthy51KJMUJty\nuUylUe4Pc1HlIOk2YHdgsKS5wEVAPYCZXQ9cCKwFTJEEsMrMxpUjjOM4lUNiC+Cr5V5fVDmY2WFF\nzp8AnFCuAI7jVA6JRuDLBKUwghCKUBbdPUJyRrUFyMOMagvQBjOqLUAeZlRbgDzMqMagEltJXAPM\nBQ4Gvg9saMb5ZfeZVLIXSVZrPgfH6UxIagDWBYYBQ2HwhnDGnnDgONhA0PhT4Bdm/LvVdWV991w5\nOE4nQFJPSL8Njb1grMHhPeDgFDwBPHkHHHi42fZNbVxb1nevYnEOjuNUFkkjIX0m9NkShmwAn+0H\nx6Vhc+BXwB7L4cVvmy3/31jGd8vBcWoPSX0g9TaM7AvXAeOBRcBk4GPg0eWw/ASzlbeU0JdbDo7T\nFZBIge0Cjy+B0X3h18AuwAHACwtgWX/IHmSWiXWrglsOjlNFJKWATYBRsMkY+NaBsM+2sKgBrgFu\nBVYAnwLemgFLJwNrmdnj7RjDHZKO05mQNALYBsZNg9P7w34KK5EPAmcBKy+APkNh1YlQ90dYeoCZ\nZcsYx5WD43QWJOrhimmw2+EwUvATQjqUh4AvrISPP2NmT0gaBmwK/MXK/LK6z8FxOgES6wAnAifD\nMR/CZME9QAZYCVy3WjEAmNl8YH5VZHXLwXHiR2J74BvA/vDyk3Dw3tCceXHoKliehUVpaNqjPf6E\n0sYu77vX3cOnHSc2JOokvizxBHA3PLEcBg6AzZoVww+BnrBgASy7EprGVVoxdAS3HBynwkgMJGx8\n+jrYm/DzmfD1M8PUAYDrgVOanYuSRhZLd9Axedwh6ThVRZq7FdgFMGQfeHIeXJGCv2wEq9JgENYl\njzazvGHO8cnlysFxEkdCwJ7At2DZLnBlf3hrGYzuBSOBU4H3noeVY81sZXVkdJ+D4ySGRA+JY6Hp\nBVg4Db7fCBu/CxcDqV5BKVyZgQ9vgZVbV0sxdAS3HBynHUgMAr5G8CfMglvvh2O/C0N6w7yclvs1\nwYyHYfG+ZpZpo7tE8GmF48SIxEjgdOBwWPgnOGgTeHI09DFY0BhapYGxy2D2B2Bvw+LxZra8imID\nPq1wnFiQ2EniLuBJ4ENgNAw6FJ42aOjVohh6AhtnYPY/YdF4WLxXLSiGjuDKwXFaIZGS2E/iMeA2\n4FFgIzPOAy2A1GOwZCwsJoQ8D1gGK06Ht16CRZPM7PWcKnCdFp9WOE6ERANwBPA/wHLgSuBOM5ok\nCXQe2GUtVzTcD6v2hoYzzJZfJanRzD6uivAF8L0VjlMmEn0J+x1OB14kLDU8YhaCEyTtCTwcYhUE\n2G3AEbCyD/ANs+VXAdSiYugIbjk43RaJtQn7HU4G/gxcacY/W85rY+C1nCs+BNvQzBYlLGqHcIek\n45SIxAYSPwFeBtYGdjbj0GbFIKmvlH6RTygGdjPLrtXZFENHcOXgdBskNpWYBjxL8ClsYcbJZswJ\n51UnpaYAiyC7eXTZNWYmM3usSmJXDfc5OF0eiW2Ac4EJwLXASDMWfrKNjgVubPErMA9sVGdfjuwI\nbjk4XRaJHSXuBe4HZgIbm3Fps2KQNETSV6SGd4EbI6UA2HZm2eHdWTGAWw5OF0RiPHABIb3a94GD\nzVguqafEntBrH+hxIPQaBj0aYFl0pV1mZhdWTfAao5Qq2zcC+wDvmdlWbbT5CbA3IaH+MWb2TEWl\ndJwiRLsj9yBUfR8OfBf4lRm5G56+AkyF07OwTyrc/hOAlW8Am5rZqoTFrmlKmVZMAya2dVLSJGCk\nmY0irBVPqZBsjlMUCUnsBfwV+BlwI7CpGTe0UgwA08Ofu1OwGbA/sOJ0s+zGrhjWpKhyiLy0HxRo\nsj9wc9R2JjBA0pDKiOc4+clRCo8RnIw/AzY345dmNLW00wgpfbrU7wXg3fDqS8DnMrD0p2aZq5OX\nvnNQCZ/DMEKy/WbmAeuz+h/hOJUjJ7nKxcBgwuaG281YY1t0CHmuuwl6jIfF0Q/h0AwsnAmzG+Dj\n0xITvBNSKYdk6+irZMIunW6FxO7ApcB60d/b8imF0Fa9gaOgaXdoItyiVxmcPxeW7wX08qlEYSqh\nHOYTHEDNrE8befYlXZxzOMPMZlRgfKeLI7ELwUIYQVAKt+ROHVraqR+wD/Q/BtJ7Qibdcrb/Cjgv\nC0v3N7NltCxRdDkkTSB4WjvWTyl7K6KyXffmW62IHJKTzWySpJ2Aq81spzztfG+F0y6iWg+XAVtE\nf282o81fe6nuD9B/L1jYEF7ZAnh9CSwbCL2nQNPzcZWrr2Vi25Up6TZgd2CwpLnARUA9gJldb2b3\nS5okaQ6wFDi2vUI4Ti4SowkWws6EJckvmLGi8DUS6PUWxTAN+PoSWLahma2S9HVYY/XCKYDvynRq\nBokRwCWEmJkfANeZUXQbtKT1We0UF2ERbsAy+OAAs8xDccnbWfB8Dk6nRWIIcD5wOC17H4rufoys\nhcuA83Je/i30eANWDHDF0DFcOThVQ6I/cCZwCvBLYDMz/lvatRpG2BxFlIAF4NNm2Sck1eH3dofx\njVdO4kQ1H04HXiGsdG1nxumlKAZJklKX8Ik88LYM6J1Tmbqpu2+aqgSuXZ3EkEgRpg6XAS8Ae5rx\nQunXaygwv5W1cIGZXR6DuN0eVw5OIkShzj8AVgBHm/HX0q+VIHURYaUsl43N7I0Kiunk4MrBiRWJ\nrQhKYRNCwpU7mxO3lna91gPehiwt1kLqWciOrXYlqa6O+xycWJBYT+IXwMOEZCtbmHFH+xRD6nzg\n7ZZXDOAQs8y2rhjixy0Hp6JINBJWIL4J3EDYPv1h+/rQusA70RE5W3XWMrN29eWUj1sOTkWIqkR9\nhZDReQtgrBlntV8xpM5ltWKAyPl4K5ByxZAsbjk4HSbaGHU14Zt8iBl/a38fGgL8Jxyl3oHsetGp\nHc2yf6+QqE47cMvBKRuJ4RK3Ar8BriHUfyhDMaTOZrVigKAY9CHQw8xcMVQJVw5Ou5HoJXEBof7D\nHEJk46/MyLavH60jycCugFSuErggKiDjG6WqiE8rnJKJsjB9AfgR8H/A9ma8WV5fOouQGRrQS5Ad\nFzkfR5jZWxUR2OkQvivTKQmJzYGfELIwfcOMR8rrRyOgbiY0rQOp+yE7KTozK9SLsHZZH05xvFam\nEwsSfSV+QMjufC+wTXsVg6SUpM9La/055FxQI3BNi2LgCLPsNq4YagufVjh5iaYQhwA/JAQybWnW\nvqTBkgZB/fHQ53RYrw+c3gcuWw7vrgBODb9N2QFm9lHl34HTUdxycNZAYjOCQjgHONSMY9qrGCI2\nhex34Efrwst9YHkWFvWE7CDQ78wycsVQu7jPwVmNRC9C4pSTgcuBa/Mlcm1fn6nHwXaFpwnZBpcC\n7G5mJW+8cjqGZ4JyOoTEROA6wrd4jFnunoZy+pMgPRNsh2Cg7gjUGbCemXlNk06AK4dujsR6wFXA\nDsDXzXig432qJ7CM1SUlssCwDCy4yhVD58F9Dt2UaC/EScBzwOvAVhVSDANZoybEo8DCZbDck7J0\nItxy6IZEMQtTgTSwR3uyMRXuV5sBs3NeIWy3uGw52JXufOxcuOXQjZBokLiQELNwG7BrBRXD3nxC\nMUBQDI3fg78uh+VXVWIcJznccugmSIwj5Fd4k5DQdW7hK0rtVwJ+ARwHPbKQETQ1e8br4OMscIuZ\nLanEeE5yuOXQxZFolPghMJ1QPWr/CiqGPaDhI+A4+Cxwewp6iuBzSJtZxgL/qsR4TrK45dCFkRhP\nsBb+SXA4llQTorS+lQbdBSv7hkDK24BNgaVLgP4eCt35ceXQBYlStX0X+DJhefJ3le1ffYDFLenb\nfgOsDbyTAetnSUXWObHi04ouRpSV6VlgHYK1UGnFMBxYvOaZ6w2WftkVQ9fBlUMXQaKnxJXAXcDZ\nZhxuxoLKjqFdgH+veSb1AqwaDXZPJcdzqktR5SBpoqSXJL0q6ew85wdLekDSs5JekHRMLJI6bSKx\nLSHseRNC6HNFrYUwRupY4AlQK6sh/Q/IjjGzl9zP0LUoqByC04lrgYnAaOAwSZu3ajYZeMbMtgEm\nAD+KCpk6MSORljgXeJCQVelLlXQ6toyTvhbsRkjfD9a75Uz6r2ZN43wq0TUp9iUeB8wxszcBJN0O\nHMAng13eAcZEz/sBC8ysQzv5nOJIbAT8ClhJSNdWkeXJT44hQfoZyG4NOh8yF7H6ByX9oFnTxEqP\n6dQOxaYVw+ATN9286LVcpgJbSHobmEUoZuLEhIQkjgL+DtwN7BWTYugBZCGzNXAg2PlAfTib/r0r\nhq5PMcuhFHPxXOBZM5sgaRPgIUlbm9kaHm1JF+cczjCzGSVL6iAxAJgCbEWoUP1cPONobeC96HAM\n8BTQMxymf2PWdGgc4zqVQdIEwhS/QxRTDvOB4TnHwwnWQy67AN8BMLPXJL1BiIZ5unVnZnZx2ZJ2\ncyR2BW4B7gN2MGu987FS42grWK101gG9BtYYDtPTzJqOi2Ncp3JEP7ozmo8lta5OXhLFphVPA6Mk\njZDUQAiFm96qzUvAXpEQQwiK4fVyhHHWJHI6nk9YojzVjMkxKoYDaFEMvUBzwPqGw9QUVwzdi4KW\ng5k1SZpM8IangRvMbLakk6Lz1xMi8aZJmkVQNmeZ2cKY5e4WSAwFfk3Y+7y9GfPjGyt9AXBpcEBm\ntgctAOsXzqavMmv6VlxjO7WJ55CsUSQ+D9xE8DF8x4xYSs5HKxJ/gMzekLoWst8AvQc2OLRIXWGW\nOSeOsZ1k8BySXQSJOuAS4GjgMLOWuWPlx1IdaCFk+oKOhezNoP/kKIaLzTKXxDW+U9u4cqghonyO\ntwGrCDkX3itySQfGUj/go2hBahewp0DzwdaJWnzbLPP9uMZ3ah/fW1EjSOxOcAD/BZgYs2LYCGhO\n2bYB8BSk5oJFZe91ulnWFUM3xy2HKhNVljoTOAM4yow/xTuedqdlmasvsBRSb0E2Cm7TKWbZKXHK\n4HQOXDlUEYm+wI3AhsA4s3w7His5XupkYErkcBwKZCH1OmSjWBadYJa9IU4ZnM6DTyuqhMSnCJGH\nHwDj41cMdb8AmwLpu8DWJSiGOZAdEbU40hWDk4srhyogsQ/wOHC1GSeasSK+sZSS0i9D5njQ/5g1\nfSmcSb0M2Y2jZoeYZX8dlwxO58SnFQkS+RfOAU4BDjDjyXjHUy/g41Bxiklm2T+GuIbUbMiOipp9\nwcyTtDhr4sohIaK8jjcCGxH8Cx2qRVl8PK1L2E4PsLmZvRQphhcgu2n0+r5m9oc45XA6Lz6tSACJ\nYYRCMquA3RNQDNvSohgG5SiGWZAdHb3+OVcMTiFcOcSMxFhgJnAnYalyebzj6cuEVPQAPcxsYaQY\nnoHsVtHre5jZQ3HK4XR+XDnEiMRBwB+ByWZcYVZSfowOjJe6CvgNpJ8EUma2MlIMT4dsTgB82vNo\nOKXgPocYiByP/wOcCnzebPUveUzjaRj0+ylof9ANZk0nRK8rJIDNbBc13cnMZsYpi9N18F2ZFSba\nOHUdsCOwr9kayXEqOJYaoMcZkDofxveAv/0XlgxrzgIt1c2EzLio+fZmFquScmqTcr97Pq2oIBJ9\nCMlwNiQENsWpGPaCPnNg1/PguUYYuBJWXJ2jGP6RoxiOccXgtBdXDhVCYghhz8LbwH5m+apCVWos\nCep3BNaFM3vDWsA9KVh1Qzhf9zxkxsI+GahfSahX5zjtwpVDBZAYBfwNuBf4qhmr4hwv1IlY9T1Y\nUg+TgJMNGv4ALJTq34DMlnA+cHQa+v7dzGJdIXG6Jq4cOojE9oQYhu+bcUncKxItpF6OnhwN9yyC\nRT8MzsemEbCLwXnAvcvggzuSkcfparhy6AASexCWKk8x4+fJjVv3C8iOBCaa2S+haUtIXRqtSpwD\nfxP0AqYL7IGk5HK6Fq4cykTiAMJc/pA4alO2PW7qxGgT1Zlm9mB4NX03ZPcCjgC+F16rnwzL3gRe\nTUo2p2vhcQ5lIHEE8ENgktma9TniG1fjgeshfbdZ049CLdOeP4DMDqCjgY2itG9bm618TtKNXsfS\nKRePc2gnEicCFwKfM+PF5MbVhsCboPejPI/rQL/fwYqdoOeD8NHxwHxI3WaWOTwpuZzax+McEkDi\nm4Tyf7snrBj6AG+GI1sP2B0aZ8PXx8JuS+Gj26C5pkX2K0nJ5XRtfFpRIhJnAScSFMNbyY2rNKyO\nmegPvc6HhrPgt71COcQBWaA5PHqz5iAox+kobjmUgMQ5wAkkrxgEat7evRGwGJregY/TcK6FFdQe\nc4FvQurnZvZy2705Tvtwn0MRIsVwDLBH3HkY1hy77reQORjY3cz+Gl5THTQHWfXJBEWRhbAL052P\nzhq4zyEGoqnEMcCE5BVD6vSgGPS1ZsUQSDfnYVgXlv49UgybuGJwKo37HNpA4jRafAzvFGtf2bH1\nWeDHUcn7n+W8vgcwAXRCqH5tO0Pqx2YZr2ruVJyi0wpJE4GrCVW2f2Fma1RCkjQBuAqoB943swl5\n2nSaaYXESYREsLvFnTJ+zbE1CngFUm+YZTbOeT1KFqsPwQbB6sK6Pp1wChJLId3IU34tsBdhqewf\nkqab2eycNgMI+Qs+b2bzJA1urxC1hMRXCLuWJlRBMfQHXglH2ZGtzs4JAU62PqRvjHTDcFcMTlwU\n8zmMA+aY2Ztmtgq4HTigVZvDgbvMbB6Amb1feTGTQWJ/QuTj5814LdmxVQ98GB32zl2SlHREVKFq\nb2AEZI4GXdr8mTtOHBRTDsOAuTnH86LXchkFDJT0F0lPSzqykgImRVTI9heEXAyJBTiFsSVaCtuu\nb2Yf55xbC/h12HHJQ8AL4YxdnKSMTvejmEOyFJO1nhCEsyfQCDwp6SkzW2PDj6SLcw5n1EqiU4lt\ngDuAQ834R/ISpO+HTC9Cjsf5rU4uDH8yu0Lqjmh1Yl2fTjhtEfkAJ3S0n2LKYT4wPOd4OKyR+mwu\nwQm5DFgm6a/A1uTZDWhWe792EhsBfyBsu34k+fHTF0B2Iuhos+zMT55LnRM93R7YCrJfCCXtsu8m\nLafTeYh+dGc0H0u6qNyO2nwQlMdrwAigAXiWUD0pt81mwMOE1YxG4HlgdJ6+rNBY1XgEr7+9DDa5\nOuOzH2CQ/kmec8PDudSvo/+D1eJn6I/af5R73xS0HMysSdJk4MHoy3+Dmc2WdFJ0/noL1ZQeAJ4j\n2LxTzSzROXs5SPQkJIO9x4xrkx9fo8P46VlmTd9odU7QvFKSPQrSD0arE516JcjpXHTL8GmJFHBr\ndHi4GYluVpI0CGhe1Ulbq81SUvoXkD2esJ9iKPAE6BSz7JQk5XS6BuV+97qrcvgOwWGzp8Vcnm7N\nsdUArIgOe1mr5K+SxgCzQJeCXQ6sBKiVz87pfLhyKFkOvgJcCuxoxn+THVuC1VbKEDN7r9X51Zuq\nzExS3VOQ2REYYGYf4ThlEEuEZFdDYifgx8BnklYMgfTjke9g29aKITr/52bfgqTPEKpmHeOKwakG\n3UY5SAwlVLo+zqw5kCjJ8dNXQHYX4Mtm9uya57UnsBvoeLAlwJ9Bi82yNyctq+NAN5lWSPQAHgWm\nm/Hd5MfXIcDtkPqeWebcPOebN1UtNMsOktKzIbsZ0NfMliQtr9O1cJ9DwbGZSqgZd7BZUkVnmsfW\ntsA/If24WdP4/G1S74CtC/QG9gDuI1gYXpDG6TDuc2gDieOAXQkOyKQVwxAgKmCb2a2NNkcC6wIT\nCYFO94Hmm2VdMThVpUtbDhLbAn8i5GWYXax9ZcdWT2BZdNjDzFbmabMWsBDST5k17Syl5oENAxot\nhKM7Todxy6EVEv0Im6lOrYJiSNGiGAblUwwRzZuqdpP0JcKO131dMTi1QJfMISkhYCrwsBm3Jy9B\nalb0ZLSZLczXQkqdFz3dDugB3AGpl83sD0lI6DjF6KqWw1cJG8J2SnpgKX0dZLckWAB5LRZJGwCX\nQ+qXZplnpNSS4G7Ibp2osI5TgC6nHCQ2B74DjDcjUfNcSh0DdgqkzjfL5LUAoijJqPZF9hhJxxBW\nKfY0sxX5rnGcatClHJJRPMNM4DozpsY51ppjayfgSUjfb9a0T9vt6m6EzLGEbfDzoc9zYGuZLVkv\nKVmd7oXXrQhcCrxOSPeWGJLWB54ElkBm3wLtxgTFkLrYzN6CnmfAuhuD6pOT1nFKo8tYDhKfBn4L\nbJ3kvglJvYHmKMZ6M2tqo12rTVX6FDQ+C4/1gp2Xm63olZDITjejW1sOEn2Am4GTE1YMaVoUw4C2\nFEMgPSN6Migsdfa7DS7vAdsATT2ivhynZugSygH4HvCYGdOTGjB8wdNPRYcjC+2clLRXSBDLMdHS\n5tdg6bbwYips+ei5EhiQgNiOUzKdflohMZ5QT2NLMz6odP/5x9Qw6HcXLNoROMnMfl6gbfOmqv+a\nZdeJXlsH+AxwW07TkWaWaK0Mp3vQLacVEr2AGwiZoxNSDOkDofFFOGgs9FwK3FjkiihTtI1ofiXK\n5fC36HAEsAlUI7+E47RNp1YOwHnALDN+H/dAknpLfW+GdW6Bh/vB1mnocW8hP4OU+hKoL3Cw5RSq\niTgdwMzeMrPXzWxRnPI7TnvptEFQUbDTSYQaGUkwEZYdAWekYWfgG4vgo1vbaiyl94Yet0DTSsje\nlafJaZBsYlun+xDt7xlGuFnLolNaDtHeiZ8BF5nxdhJjmtldkDkTLgHGr4TnGwj1OvLIl94TGu+E\nqxugz7vWyrETRUlCqOTtOGUjaZCkvST9WEq9KskkGSHf4L+B35Tbd2e1HA4nhBxfn9SAkvoCV4E+\ngFn3Qv0gs+VrhGdLGg+9p8MfGuEDIP1Wnu6arZ0b4pTZ6RpIagQ2Bz4P7Eub1oABegvsXkIVt38S\nfFllWaidTjlI9AW+D3zJLGRjTYjIJ2DrmS1aIWmNpUdJa0Pd72FsD+gFPAOsfClPX+cCmNmCGOV1\nOhFRkNxGwO4EBbAvoZBUPpYQsoXdR3Bs/9vM2vwutBiq7aPTKQfgAsJW7KeKtqwQUury6OnY5s1R\nZvZh63Zm9l9JW8Gj82AcMLwJluRTDgeDWtccdbo40XRyPcLNsS9oP7B1ClzyCEEBPAy82rrGSdx0\nKuUgsTFwPLBFcmNqNHAepH5qlvm/Ei5ZFf09B/7TD+z5Vv1Fn7mtkWjW6RpI6g+MAfYG7Qu2Vf6W\nBug5sPuAB4Dn8/3oVItOFQQl8RvgeTMuL9q4AkiqJ6o4BaRaOxbzX5N+PsrnkLd9iJbkIaCPmS2t\nqMBOYkRpAD8F7EWYAuxRoPW7kR/gPuDvwH9KuZfKlEvABsA+oCPBdoLyKqZ1GstBYmdgF+DY5EZN\nPRf5cgaWphjUH9gSdLZZto32ugQMVwy1T7TfZQNCguJmP0DvNppnCF/+e4G/Am8U3mtTEfkagR0I\nDvojCY6uylGsDDchK/JLwKvA2QXa7QA0AV9s47wVG6vtvk1gj4EdU24f7R+Towh234GlX5P6S3RN\nqkC/Bnoqqffhj6L/ZwHrRPf5NcHbjxV4PAF8m5DerzFBGTcETiGkBigkX+vHY+V+94oJlQbmEEJ8\n64Fngc3baNfsPDmojb7KEjBca3uD/QssndA/Y0j4YNOPt+OaxuifcUUJbT6XxPvwxyc++z6E8oIX\nEkz7Al8ovUoom7gnIUFwknI2ElYsfgYsbaci+BXwOUIxpNw+rRxZik0rxgFzzOxNAEm3AwfAGtmc\nTyWUmtuhSH/tJgp4uhy40BJYuozmbP8JR5ndS78ydVM0BbmgQKMDo7+PlCObU5jIRzSSMP/fF9i7\nQOsPc/wATwLzzSyxiNUc38AkwpQgT+yCIG+pFc0Hm0bIrv6ixTR9KaYchgFzc47nEbTvasIORQ4g\n7DLcgfzvpiN8Mfp7d4X7bYP0rVEx202swNpxLpJ6EJYnbzDLrirQ8ruRvyHWuWhXJgoLHkpIHty8\nHDiwwCUPEhTAI4QfurbKBMRG5BsYS4tvoLG0Kw1CHMNNwB8JCiyxwkzFlEMpglwNfNsspDciqLu8\nSLo453CGmc0o1HFkNVwInGcWf7UqSbsBh4L+xyz7ejuu/HFk2U0u3M42JDEl17mRNJAQSTopUgCb\n5m9pgP4vZznwX2a2ODFBcyjNGijIrYSpwRMdeQ+SJgATyr1+NUXmPzsBD+Qcn0MrpyQhZ+Mb0WMx\n8C6wf56+2j3vAdsP7BkwJTDX60240xa387q6cF3qd0XaDYz63y7JOWwtPwje9W2Bswge/kJ+gLnA\nT4F9CD6h2O+JIrI3ArsBUwgRi/lkzrTx+tuEDOlbA3UJyGplXVek0zrgNYJDsoE2HJI57adRodWK\naIViJtiXEvpnr4z+eb3aed1l0XW9i7Q7M2rXn+AdH06YH48GhiXxHqvxiO6hkYQl6LtyPud8j48J\neUCPAjYGEnFAF5G/2Ro4mbBS0R4HoYGeJOweHl4thVaucig4rTCzJkmTCfO2NHCDmc2WdFJ0Ps6N\nT3sC/UjADJfS5xFWY3a1dpSii+a/50PqcbNMm3ELUeanUyMjYwE0roQ6g4YsLOgFTSfRiTdhReb0\nEILPaZ/B7f8eAAARB0lEQVRoGjC0wCWPEvwADwGvtOczj5PIN7A9Lb6BfDENKwk/lPm4HfglYVrQ\n+fNz1Kr2AnsA7LgE5BpFWLacVsa1p4drGVBCW0H9GdD4MUzNQtZgcfMvzD5U2Uwu8f32IwQEXUaw\nIgv9Yv4LuIKwLLdWtWVf839RijWgZW28/l703rYhZByv+nsq8n6tnOtqMnxaYgvCZpMRZsRWBSo3\nXTwlhkfnXCsgG+pbZjZrx3VbQt97YLehsG8v+Fru2UVgJwJ3m1mBVY/4iFZeRhIst32BzxZo/T4t\n24OfAt5uz2eYBFEOz7HAYQRroE+eZstoM7pQfwe7ifAe59ba+yuFsrcu1KL2ApsKdmH8MqWbg2HW\nLuP9HB1dO6SMa3tCn+tABvXnAWsDP2LNX6mLKcEqyem3jjC3H1mkXYrgRzoMuAX4KM/YzY8sMB04\nkbCXoCZ/KWmxBk6isG9gQYFzvyFYcf2r/X4q/NlYWdfVmoBga4N9ALZOzPJ8KbohDi33/YDe76AM\nnwU+1eq1nsBxhFD03Bv3TmBUG/3UgY6CPu9AfVP03gQMjsa4GjSnwJfCCL/85xP8Bn2Sui/K/Nx6\nAZ8GriOskOV7P4tDoFPeacH7wA8I/oWGar+fBD4vK+u6WhMQ7NtgN8Ysy6Bwk6SeKfP6A6IbbcOY\n5RSwG2hWq5v7lehLXwc6HPrOg7GL4S8GO7WxfIaBXguKgs9FiqOm/RytrIHHCyi2eYTItXzv+Wng\n64S9CTX9fmP8HK2c62rK5yCRImzwOsyMv8ckR+QrAAqUryvSR/MnHmth4DzjjiDUAz3yk2f2JuSr\n3Yrw3b8S4AiCeT3XEgwLLpfIN7A9Lb6BvnmaLSLk3tuwjW7uIlQ+e8xqKC9CtSnX51BryuGzBHNv\nWzNiEUxKT4XsCcBmZvZy+6/X3sD9wBZm9mLFBSxdjqHAIcAS6DMWeuwASz8FmZ7Q43mzxdtWS7ZC\nRMp5fYJGO5IwPcjHHFBfsCF5evmAsLfgdmCWVSEkujPRVZTDHcAjZkyJSYYdgadClevMJWX2URWr\noRSiL966wMZm9kQNyNOLsLX5MEJgUz5r4ENCluTR5A3n1z+j1YJ7gbcsqRu2C9HplYPEEOBlYEMz\n2qw72YHxexKWrMr+YkvaDNIvQuZYM7u5ogJ2YtphDTxP8CNs2cb5ewjTgkfNLJEKZt2BcpVDLWWC\nOhSYHodiCOjt4KPK++tV/GpJ0O8GWA5kYpKx9mllDRxJCIxqzULgX6ANwTbIeT3KpaiPomnBbYRp\nQWyxLE751JJyOJywA7PiSKnTwNYC9jCzJeX1kjoIBm8N41fAzYVCg7sEOdbARIISGN9G06cJjsJP\n0xJWPDC0NwgJVG8kxEq86dOCzkNNKAeJUQQP9J8r37dGAFdB6i6zzIwy++gDjT+Dab3D5sH64RUU\nserkWAOHEhRB/zzNFhCSogwk5PJsZmzO8+mEacEMM1sYj7ROUtSEciBYDb8xo6JJUKIEoW+Eo+zB\n5ffUeCmM7xN+POcAfTapgHiJk2MNfJ6gBHZro+nfgHmg7cBGRq8NIoRTQ7AUbiJMC57xaUHXpOrK\nIUro0rwLrsKkH4qyOg0t15yV1AC9hsNjyyHVAzZbDqkNil8ZL5LGAeubWd5dq5E1sC3BGjiKtq2B\nhwj3wSRaMhRFloERNlDZNIKz8HWfFnQjqh2lBTYG7A0qnNCF8CtnoGMr1N93Qn8cApyW1OeWR46N\nof/0ECY98E8E7//6wAmErdBtRRE+ClxDUAZttbmXEHqdaFJVf8R+z1g511V9KVPifGCwGadVcKwB\nwAeQmmOWGVWhPpcBPfO9hySQNAj6XAbZY+FL9TAoDVfla7oA+D2wBLQX2Og8bZYQpgW3EqYFiZZZ\nc5KlMy9l7k+oA1ARwjQgPTtMJ7KVLJvXE/h1BfsrSLSdfCzwDUgdCHVpqK+HkQp+v9UV1n4MDAAO\nJizTDiJs3CKaFszOmRa8Zp0glNqpDaqqHCSGEnIHPFaZ/jQS+k6HxesCe1uFwmolNWcEuq4S/bUx\nxkDgC6AzwVrlh8gugZ5vwqpPwcHpYPlv3HzyWzkN7yesFjxiZu/HJavTPai25bAv8IAZHU5sItUd\nCY1T4OBGuONdWPxgBeRrpjnhydOlyaJUoV/oKL3cloRMRF/75FkD9ATYFGBT0NmQaoD6zeA7Cjuq\nV8cdnQXcYVFdEcepJNVWDvvTQVNdUt8Qubj+PvD7Rrg7A5k7rELOFEmDQd+Lois/I2k+MB/4qPUY\nktaF3tcQsgbd9EkZmQQ6Ayxf4Z9rgOuB/wIXEFLc7xpOWQYyV8LyBrhoS8iOgqVDCdOcV10xOLFR\nLY8pWA+wxWAlZzpqo9+jIJWBHxiYwacWAbtVVvbU9yHdBANXwvqLoOcKqF8JvecTVi/qoOF06LUE\nxmWg8T5CjsHWCVsM9DxwDFFeRWAk6J5WbeYRIhPz1tykJZHL4CS93v7onI/W372Sr6uWgGDjwZ6u\nUN8XhS/VZ1ZAzyXEUAsAmACNS+B3kRKaZUFBMAn6zoGdl8CLBre3Xh6cRpgLpKN+BOwYdht+QiE8\nRchn0C0Tkvgjvke5yqFqS5nREuYAM87sYL/Ny5avQe+ZkE2bLTm0o/K2Mdb20PhnuLYf/C8wK3o/\nw7JwfCpUD1wbODgLmUEWJRyJIjUngW4IafBWcyehSFA7qms5TvvodFu2JR4G/teMezvQZ25Wpx5m\ntlJSfzOLbddk2Lbd+zFYNQBWXkLIRdAfUkdD77nQsA4sHwRLTyGsxLSu7XEV8F3z1QQnITpVnINE\nA6EgbweXMENhWmAbi5Yt41QMUf8vBQtCJ5vZ5ZK2gj7TIbU+LL4DrJGwuzR3E9kZwPVm1mbhG8ep\nNapiOUjsClxjxnYd6G9z4EVIXWuWObVScrZjfEHDZKi/Aq7tBW8pZJIHQtKH44A7rUr1JxynmU41\nrZA4lxAy/a0il7XVV9nFaMocb6iZvZ1zvB30+TH03hl+3RAynC0CNgcYaJ7FyKkhOtW0glAi7afl\nX55+PNptOTgBxfBpSN8vaRvCSsMrQApWfhQyGn0hDdk6sBTUC1btSXA0Ok6nJnHlEG3RHgvMLO96\nHUDwVxxhZgsqKVuesXaBxgdgfC948LWcU/uZrbgvzrEdp9qkSmkkaaKklyS9KunsPOePkDRL0nOS\nnpA0pkB3Qwk/+++2V1hJ/YF7IPWKmd3a3uvbOdZO0PgnuKs3/Kz5czrTzGRmrhicLk9R5RCt0V9L\niNgbDRwWOQNzeZ0QlTiGUIH55wW6HAM8Z1ZWXYqoUEl2q8LN2o+k3lGClEgx9PwzXNY7VFy/Exi5\nCnqMLdyL43QdSplWjAPmWBTDL+l2Qjm42c0NzOzJnPYzCclH2mIM8Fx7BZVSl0dPt7VYipjocGg8\nTdKDoNNgBXDJu1A/B5bPgqUvAC9VflzHqU1KUQ7DgLk5x/MIc/62OJ6wdbgtxhCyEZWMpE2B8yD1\nM7PMs+25tnT6T4bPbAZ3jwb+A7ah2UdeScnptpTicyjZ/Je0B2F9fw2/RA5bA7Pa0Wcdq3+xs6eU\nel17CHkgsqPg9hRMBeoHg46NYyzH6SyUYjnMB3JTsQ8nWA+fIHJCTgUmtrXOL/W4DM7eFK48SFrR\n38xmFB8+/Wi0bLlOfMuWPY+DI9PBdfLzpdDjDVhZ9XJyjlMOkiYAEzrcUQk7uuqA14ARhKIlzwKb\nt2qzASFn+06FdoaBbQP2r3bsJtsvXMeRMe5YS0Hv/8LBGei1FBpOI9pB6Q9/dIUHZe7KLGo5mFmT\npMnAg0AauMHMZks6KTp/PWEvwVrAlLAXilVmNi5PdyU7IyX1A6ZD6nWzzK9KuaY8dBosHQwPPgbL\njjSzt+Iby3E6D4mGT4P9CFhgxvdKaw9Euy3jlSv9EmRGW1IfhuMkSLnh0yUFQVWQkYRK2gWRUpdE\nT7ePWTFElasyh7licJxPkrRyGMDqQKb8SBoFdiGkpprZP+MVJ3U7gJnFtDzqOJ2XpPdW9APazLcQ\nRWO+Eo6yJ8UpSJQBeiwd2gDmOF2XpC2H/hRQDpB+JHoyJAEz/8vR33NiHsdxOiVJWw5tKoewn4Hd\nQm3L7HsJyHIbgJktSmAsx+l0VMNyaOvL+BpwnFn2priFiJLSAhwU91iO01lJeilzhRk9Exkwvww7\nAgtA3wY7ngSySDlOteksmaBiTf5anJ6HwopvgglYCjRGfx3HaUXS04oqK4eGAfAjwYnA6DTUvZ4z\nxXAcJ4duphzq+oUo74uA/xhkv2pR4RnHcT5JN5tWqG/YHjJxKSy70iwzvbryOE7tkrRyqPKyofrA\nhRlY8GdYdll1ZXGc2qabTSsyjbBgDiz2vRSOU4Ruphw+mA+LP2tmH1dXDsepfbqZcuAUM5tbvJnj\nOEkrh6r6HDyRi+OUTnezHBzHKRFXDo7j5MWVg+M4eXHl4DhOXrqVQ9JxnNJxy8FxnLy4cnAcJy9J\nKwePTHScTkKiysGs9KK8juNUl6QtB8dxOgmuHBzHyYsrB8dx8lJUOUiaKOklSa9KOruNNj+Jzs+S\ntG3lxXQcJ2kKKoeoPN21wERgNHCYpM1btZkEjDSzUYTMrVNikrXiSJpQbRlaU4syQW3K5TLFSzHL\nYRwwx8zeNLNVwO3AAa3a7A/cDGBmM4EBkoZUXNJ4mFBtAfIwodoCtMGEaguQhwnVFiAPE6otQKUo\nphyGAbnJUeZFrxVrs37HRXMcp5oUUw6lxiW0rqbj8QyO08kpln16PjA853g4wTIo1Gb96LU1CCXx\nagtJF1VbhtbUokxQm3K5TPFRTDk8DYySNAJ4GzgEOKxVm+nAZOD2qFL2h2b2buuOyqnV5zhO9Sio\nHMysSdJk4EFCNZgbzGy2pJOi89eb2f2SJkmaQ6g7eWzsUjuOEzuJVdl2HKdzUfEIyVoMmiomk6Qj\nIlmek/SEpDHVlimn3Q6SmiR9sRZkkjRB0jOSXpA0o9oySRos6QFJz0YyHZOATDdKelfS8wXaJH2P\nF5SprHvczCr2IEw95gAjgHrgWWDzVm0mAfdHz3cEnqqkDGXKtDPQP3o+sRZkymn3CHAfcFC1ZQIG\nAP8C1o+OB9eATBcD32uWB1gA1MUs13hgW+D5Ns4neo+XKFO77/FKWw61GDRVVCYze9LMmhPRzCT+\nOI1SPieAU4E7gf/GLE+pMh0O3GVm8wDM7P0akOkdoF/0vB+wwMya4hTKzB4DPijQJPHAwGIylXOP\nV1o51GLQVCky5XI8cH+M8kAJMkkaRvgiNIejx+0cKuVzGgUMlPQXSU9LOrIGZJoKbCHpbWAW8M2Y\nZSqFWg8MLOker3SV7VoMmiq5b0l7AMcBu8YnDlCaTFcD3zYzkyTW/MyqIVM9sB2wJ9AIPCnpKTN7\ntYoynQs8a2YTJG0CPCRpazNbHJNMpVKTgYHtuccrrRwqGjSVoExEDpqpwEQzK2QyJiXT9oTYEQhz\n6b0lrTKz6VWUaS7wvpktA5ZJ+iuwNRCXcihFpl2A7wCY2WuS3gA2JcToVIuk7/GSaPc9XmGnSB3w\nGsGB1EBxh+ROxO/8K0WmDQiOr53idhyVKlOr9tOAL1ZbJmAz4GGCo7AReB4YXWWZfgxcFD0fQlAe\nAxP4H46gNIdk7Pd4iTK1+x6vqOVgNRg0VYpMwIXAWsCU6Jd6lZmNq7JMiVLi/+4lSQ8AzwFZYKqZ\nvVhNmYDvAtMkzSL40M4ys4VxyQQg6TZgd2CwpLnARYQpV1Xu8VJkoox73IOgHMfJi6eJcxwnL64c\nHMfJiysHx3Hy4srBcZy8uHJwHCcvrhwcx8mLKwfHcfLiysFxnLz8PxbIkZvZJzhDAAAAAElFTkSu\nQmCC\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x10850a278>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%run listing06" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "初期値を決めたときに、軌道を計算・プロットする方法の説明は終わり.\n", | |
| "\n", | |
| "残る問題は安定多様体上の軌道をどのように計算するか？\n", | |
| "\n", | |
| "長期均衡に収束するような初期値をなんとか探せばよい \n", | |
| "\n", | |
| "あるいは、このモデルなら**長期均衡の近傍から出発する逆向きの軌道を計算する**方がスマートにできる. 生産関数と効用関数の特定化次第では逆向きの時間発展が計算できる. ヤコビアンの安定固有ベクトル方向に摂動を加えて、逆向きに時間発展させれば安定多様体上の軌道が計算できるはず " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "celltoolbar": "Slideshow", | |
| "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.4.3" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 0 | |
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/usr/bin/env python3 | |
| from math import log | |
| ## Unicode identifiers are only valid in Python3 | |
| ## In Python2, use less readable notation such as | |
| # alpha = 0.4 | |
| # rho = 0.9 | |
| A = 1.1 | |
| α = 0.4 | |
| ρ = 0.9 | |
| def f(k): | |
| """Production function""" | |
| return A * k ** α | |
| def U(c): | |
| """Utility function""" | |
| return log(c) | |
| def u(x, y): | |
| """Reduced form utility function""" | |
| return U(f(x) - y) | |
| def F(x, y): | |
| """Solution of Euler equation""" | |
| return ((1 + ρ * α) * A * y ** α | |
| - ρ * α * (A ** 2) * (x ** α) * (y ** (α - 1))) | |
| def G(x): | |
| """Dynamical system""" | |
| return [ | |
| x[1], | |
| F(x[0], x[1]) | |
| ] | |
| if __name__ == "__main__": | |
| duration = 4 | |
| x0 = [0.8, 0.43] | |
| x = x0[:] | |
| for t in range(duration): | |
| print(x) | |
| x = G(x) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/usr/bin/env python3 | |
| # -*- coding: utf-8 -* | |
| import numpy as np | |
| import matplotlib.pyplot as plt # プロッティング・ライブラリを読み込む | |
| from listing01 import f | |
| fig = plt.figure() | |
| ax = fig.add_subplot(111) # 作図領域の作成 | |
| grids = np.linspace(0.0, 1.2, 120) | |
| ax.plot(grids, f(grids)) | |
| plt.show() # 図の表示 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/usr/bin/env python3 | |
| import matplotlib.pyplot as plt | |
| from listing01 import f, G | |
| duration = 10 | |
| x0 = [0.8, 0.43] | |
| fig = plt.figure() | |
| ax = fig.add_subplot(111) | |
| grids = np.linspace(0.0, 1.2, 120) | |
| ax.plot(grids, f(grids)) | |
| x = x0[:] | |
| for t in range(duration): | |
| ax.plot(x[0], x[1], marker='o', linestyle=' ') | |
| x = G(x) | |
| plt.show() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/usr/bin/env python3 | |
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| from listing01 import f, G | |
| duration = 10 | |
| x0 = [0.8, 0.43] | |
| fig = plt.figure() | |
| ax = fig.add_subplot(111, aspect='equal') | |
| grids = np.linspace(0.0, 1.2, 120) | |
| ax.plot(grids, f(grids)) | |
| x = x0[:] | |
| for t in range(duration): | |
| x1 = G(x) | |
| dx = [x1[0] - x[0], x1[1] - x[1]] | |
| ax.plot(x[0], x[1], marker='', linestyle=' ', color='black') | |
| ax.arrow(x[0], x[1], dx[0], dx[1], length_includes_head=True) | |
| x = x1 | |
| plt.show() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/usr/bin/env python3 | |
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| class Ramsey: | |
| """One-sector Ramsey model""" | |
| def __init__(self, A, α, ρ): | |
| self.A = A | |
| self.α = α | |
| self.ρ = ρ | |
| def f(self, x): | |
| return self.A * x ** self.α | |
| def U(self, x): | |
| return log(x) | |
| def u(self, x, y): | |
| return self.U(self.f(x) - y) | |
| def forward(self, x): | |
| """1ステップの時間発展""" | |
| A, α, ρ = self.A, self.α, self.ρ | |
| return [ | |
| x[1], | |
| (1 + ρ * α) * A * x[1] ** α - ρ * α * (A ** 2) * (x[1] ** (α - 1)) * (x[0] ** α) | |
| ] | |
| class Simulation: | |
| """Simulation of a dynamical system""" | |
| def __init__(self, sys, x0, duration): | |
| self.sys = sys | |
| self.x0 = x0 | |
| self.duration = duration | |
| def __iter__(self): | |
| x = self.x0[:] | |
| for _ in range(self.duration): | |
| yield x | |
| x = self.sys.forward(x) | |
| if __name__ == "__main__": | |
| ramsey = Ramsey(A=1.1, α=0.4, ρ=0.9) | |
| sim = Simulation(ramsey, x0=[0.8, 0.43], duration=10) | |
| fig, ax = plt.subplots(subplot_kw={'aspect':'equal'}) | |
| grids = np.linspace(0.0, 1.2, 120) | |
| ax.plot(grids, ramsey.f(grids)) | |
| for i, x in enumerate(sim): | |
| if i == 0: | |
| ax.plot(x[0], x[1], marker='', linestyle='', color='black') | |
| else: | |
| dx = [x[0] - x0[0], x[1] - x0[1]] | |
| ax.plot(x[0], x[1], marker='', linestyle='', color='black') | |
| ax.arrow(x0[0], x0[1], dx[0], dx[1], length_includes_head=True) | |
| x0 = x[:] | |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| from collections import namedtuple | |
| from listing05 import Ramsey | |
| class Simulation: | |
| """Simulation of a dynamical system""" | |
| def __init__(self, sys, x0=None, duration=None): | |
| self.sys = sys | |
| self.x0 = x0 | |
| self.duration = duration | |
| def __iter__(self): | |
| x = self.x0[:] | |
| for _ in range(self.duration): | |
| yield x | |
| x = self.sys.forward(x) | |
| def reset(self, *, x0=None, duration=None): | |
| if x0 is not None: | |
| self.x0 = x0[:] | |
| if duration is not None: | |
| self.duration = duration | |
| def plot(self, ax): | |
| for i, x in enumerate(self): | |
| if i == 0: | |
| ax.plot(x[0], x[1], marker='', linestyle='', color='black') | |
| else: | |
| dx = [x[0] - x0[0], x[1] - x0[1]] | |
| ax.plot(x[0], x[1], marker='', linestyle='', color='black') | |
| ax.arrow(x0[0], x0[1], dx[0], dx[1], length_includes_head=True) | |
| x0 = x[:] | |
| if __name__ == "__main__": | |
| ramsey = Ramsey(A=1.1, α=0.4, ρ=0.9) | |
| sim = Simulation(ramsey) | |
| SimParam = namedtuple('SimParam', ('x0', 'duration')) | |
| params = [ | |
| SimParam((1.2, 0.460), 10), | |
| SimParam((1.2, 0.430), 10), | |
| SimParam((1.2, 0.425), 6), | |
| SimParam((0.05, 0.15), 10) | |
| ] | |
| fig, ax = plt.subplots(subplot_kw={'aspect':'equal'}) | |
| grids = np.linspace(0.0, 1.2, 120) | |
| ax.plot(grids, ramsey.f(grids)) | |
| for param in params: | |
| sim.reset(x0=param.x0, duration=param.duration) | |
| sim.plot(ax) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment