{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Python の応用1 〜数値シミュレーション〜\n", "\n", "この章では, これまで学習した知識を用いて初歩的な数値シミュレーションを行う. \n", "本章と次章を学習することによって, 最終的には, 野球ボールの軌道を計算できるようになることを目標とする. \n", "\n", "多くの数値シミュレーションは, \n", "解析的に解くことのできない微分方程式を数値的に解く問題に帰着する. \n", "まずは微分方程式を数値的に解く最も簡単なオイラー法について少し学習し, \n", "その後, 実際の問題に適用する. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 常微分方程式とオイラー近似\n", "\n", "一階の微分方程式\n", "\n", "$$\n", "\\frac{dx}{dt}=f(x,t)\n", "$$\n", "\n", "を数値的に解くことを考える。\n", "計算機は離散的な値しか扱うことができないので、微分方程式を差分方程式に近似する。\n", "その最も単純な近似が次のオイラー法である。\n", "\n", "$$\n", "\\frac{x_{n+1}-x_{n}}{\\Delta t}=f(x_n,t_n).\n", "$$\n", "\n", "ここで、$\\Delta t$を十分に小さい時間刻み幅として、\n", "第 $n$ ステップにおける$t, x$をそれぞれ $t_n(=n \\Delta t), x_n$ としている。\n", "初期条件 $x(0)=x_0$ を与えた上で、 $n=0, 1, 2, \\cdots$に対して\n", "\n", "$$\n", "x_{n+1}=x_{n}+f(x_n,t_n)\\Delta t\n", "$$ \n", "\n", "として逐次 $x_n$ を求めていく。\n", "\n", "以下の例では, 常微分方程式\n", "$$\n", "\\frac{dx}{dt}=ax, \\ \\ \\ x(0)=1\n", "$$\n", "\n", "を数値的に解いている. \n", "なお, この方程式の解析解は\n", "$$\n", "x(t)=\\exp(at)\n", "$$\n", "\n", "で与えられる." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $\\frac{dx}{dt} = ax$ のシミュレーション" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "n = 1000 # ステップ数\n", "dt = 0.01 # tの刻み幅\n", "a = 1.0 # a の値\n", "\n", "# 初期条件\n", "x0 = 1.0\n", "# t, x の値を格納するリスト\n", "t_list = []\n", "x_list = []\n", "x_exact = [] # 解析解を格納するリスト\n", "\n", "# t=0 では x=x0\n", "t = 0.0\n", "x = x0\n", "\n", "# 値を格納する。\n", "t_list.append(t)\n", "x_list.append(x)\n", "x_exact.append(x)\n", "\n", "for i in range(n):\n", " f = a * x # 右辺を予め計算する。\n", " x = x + f * dt # Euler 法による時間発展\n", " t = t + dt # 時刻をdtだけ進める。\n", "\n", " # 計算した x, t の格納\n", " t_list.append(t)\n", " x_list.append(x)\n", " x_exact.append(np.exp(a*t))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAERCAYAAACO6FuTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FGXXx/HvHSBIFQgC0gKiKEpTHxELGsRXsSDPIypFioBiQQi9RIEASpHeBIXQRKQpvYpItYCAUkOoCSU0gRBqQnK/fyRAKIGQNrvZ3+e6crEzuzN7Mgl7crczxlqLiIjIzXg5HYCIiLguJQkREUmUkoSIiCRKSUJERBKlJCEiIolSkhARkUQpSYiISKKUJEREJFEulySMMc8bY1YaY0YaY55zOh4REU/mckkCsEAkkBU44HAsIiIeLc2ThDEmyBhzxBiz6br91Y0xwcaYEGNMx8v7rbUrrbWvAZ2AHmkdn4iIJC49WhLjgJcT7jDGeAHD4/c/AtQ1xjx03XGnAO90iE9ERBKROa3fwFq72hjje93uSsBOa20ogDFmClATCDbG/I+45HE3cYlEREQckuZJIhFFgP0Jtg8Qlziw1s4EZt7qYGOMSteKiCSDtdbcyetdceA6Say1+rKWbt26OR6Dq3zpWuha6Frc+is5nEoSB4HiCbaLxu9LssDAQJYvX56aMYmIZEjLly8nMDAwWcemV5Iw8V+XrQPuN8b4GmO8gTrAnDs5YWBgIH5+fqkXoYhIBuXn5+e6ScIYMxn4DShtjAkzxjS21sYALYAlwFZgirV2+52cVy2JOEqUV+laXKVrcZWuRcpaEia5/VROMsZYd4xbRMRJxhjsHQ5cOzW7KU2UKFGC0NBQp8OQJPD19WXfvn1OhyEit+G2SeLymETCpmRoaGiyR/AlfRlzR3/MiEgKLF++PNnd8xmquym+KeVARHKn9LMSSX/J6W5y23USIiKS9tw2SWh2k4hI0mh209X9btmF0bt3b/bu3cu3336b6ueuWrUqDRo0oEmTJnd87P79+3nkkUeIiIhI9TEEd/1Zibgzj5/d5K46d+7sdAgAlCxZkqCgIF544QUAihUrxunTpx2OSkScpO4mEZEMzh3KcqQ6dy3L0bdvX4oWLUru3LkpU6YMv/76K927d6dBgwZA3DReLy8vxo8fT/HixfHx8eGbb77hr7/+okKFCuTLl48WLVpcOV/CYxMeHxsbe8N779mzh2rVqpE/f34KFChA/fr1r7QUGjZsSFhYGDVq1CB37tz079//hnOFh4dTs2ZNfHx8KF26NGPGjLkmjtq1a9OoUSNy585NuXLl2LBhQ5pcQxG5My5dlkOuCgkJYcSIEaxfv57Tp0+zePFiSpQoAdy4bmDt2rXs2rWLqVOn0qpVK3r16sWyZcvYsmUL06ZNY9WqVVdee/2xiY0fWGsJCAjg8OHDbN++nQMHDlz5xZk4cSLFixdn3rx5nD59mnbt2t1wrtq1a1O8eHEOHz7M9OnTCQgIuKY1N3fuXOrVq0dERAQ1atSgefPmyb1UIuIiPC5JGJM6X8mRKVMmoqKi2LJlC5cuXaJ48eKULFnyJjEaunbtire3Ny+++CI5cuSgbt26+Pj4ULhwYapUqcLGjRvv+P1LlSpFtWrVyJw5Mz4+PrRu3ZoVK1Zc85rEBpP379/P77//Tt++fcmSJQsVKlTg/fffZ+LEiVde8+yzz/Lyyy9jjKFBgwZs2rTppucSEffhtkkiuWMS1qbOV3KUKlWKwYMHExgYSIECBahXrx7h4eE3fW2BAgWuPM6WLRsFCxa8ZvvMmTN3/P5Hjx6lbt26FC1alDx58lC/fn2OHz+epGPDw8PJly8f2bNnv7LP19eXgwevVngvVKjQlcfZs2fnwoULN+32EpH0pTEJN1KnTh1WrVpFWFgYAB07dkzR+XLkyMG5c+eubCeWdAACAgLw8vJi69atnDp1ikmTJl3TcrjVNNfChQtz4sQJzp49e2VfWFgYRYoUSVH8IpL2NCbhJkJCQvj111+JiorC29ubbNmykSlTphtedyfrBypWrMjKlSvZv38/ERER9OnTJ9HXRkZGkjNnTnLlysXBgwfp16/fNc8XKlSIPXv23DSWokWL8vTTT9O5c2cuXrzIpk2bCAoKumbQPCXfh4i4JiWJdHTx4kU6derEPffcQ+HChTl27Bi9e/e+4XW3G4hOuP3iiy9Su3ZtypcvzxNPPEGNGjUSfW23bt1Yv349efLkoUaNGtSqVeua13bq1ImePXuSL18+Bg4ceMPxP/zwA3v37qVw4cLUqlWLnj17UrVq1US/XxXxE3F/WnEtjtDPSiT9qcCfiIikKrdNElpxLSKSNCrwd3W/ujDchH5WIulP3U0iIpKqlCRERCRRShIiIpIoJQkREQ+w/tD6ZB2nmw6JiGRgR84c4bNlnzFz6/xkHa+WRAawYsUKihUrluzjc+XKxb59+1IUQ+PGjenatWuKziEiqScqJor+v/Xn4eFlWb00L3eN3pGs87htktA6iWsltQRG1apVGTt27DX7IiMjr9zXQkTcm7WW+SHzeWREWb5dsoKY0Wt4MuI1GtUZmKzzuW13U3IXhoiIZFTbj22nzeI2bNq/j6g5Q6nmW52+y8DXtzTgR+/e3e/4nG7bknBXffv25f777yd37tyULVuWWbNmATBhwgSqVKlC+/btyZcvH6VKlWLRokVXjhs/fjwPP/wwuXPn5v777+fbb7+96fn79+/PW2+9dc0+f39/Wrduzeeff86qVav49NNPyZ07Ny1btgTAy8vrSvXXCxcu0LZtW0qUKEHevHl57rnnuHjxIgDvvPMO9957L3nz5sXPz49t27al+vURkTt36sIpWi9qzdNjniN4/ssU/HETP31VnSlTwNc3hSe31rrdV1zYN0psvyuZMWOGPXz4sLXW2mnTptmcOXPaw4cP2/Hjx9ssWbLYoKAgGxsba0eOHGkLFy585bgFCxbYvXv3WmutXblypc2ePbvduHGjtdba5cuX22LFillrrQ0PD7c5c+a0ERER1lprL126ZAsUKHDltX5+fjYoKOiamLy8vOzu3buttdZ+8skntmrVqjY8PNzGxsba33//3UZFRVlrrR03bpw9e/asjYqKsq1bt7YVK1a8co733nvPdunSJcnXwR1+ViKu7lLMJTtq3Sh7T9+CtnTbZrZAiaM2KMjaS5du/vr4/3d39Hnrtt1NyWW6p075atsteSUlEpbnfvvtt+nVqxdr164FoESJEjRp0gSARo0a0bx5c44ePUqBAgV45ZVXrhxXpUoVXnrpJVatWkXFihWvOX+hQoV47rnnmD59Ok2bNmXhwoXcc889N7zumu8lvjyGtZZx48axdu3aK3eZq1y58pXXvffee1ced+3alcGDBxMZGUmuXLmSdS1EJPlW7FtBy4X+RB67m6jJi/hvzYp89g/kzp267+NxSSK5H+6pZeLEiQwaNOjKbKKzZ89y/PhxvLy8rrn9Z7Zs2bDWcubMGQoUKMDChQvp0aMHISEhxMbGcv78ecqXL3/T92jYsCGjRo2iadOmfP/997e8MVBCx48f5+LFi9x33303PBcbG0tAQAAzZszg+PHjGGMwxnD8+HElCZF0FHoqlPY/t+fXnX/itbQ/lXO/xYB5hvvvT5v305hEOgoLC6NZs2Z8/fXXnDx5kpMnT/LII4/cttBdVFQUb731Fh06dODYsWOcPHmSV155JdHj/vvf/7Jp0ya2bt3KvHnzePfdd688d6tZUPnz5+euu+5i9+7dNzw3efJk5s6dy7Jlyzh16hT79u1L2P0nImnsbNRZuv3ajYojH2fdgrLcM3U73we8zexZaZcgQEkiXZ09exYvLy/y589PbGws48aNY8uWLbc9LioqiqioKPLnz4+XlxcLFy5kyZIlib4+a9as1KpVi3r16vHkk09StGjRK88VLFjwhluUXmaMoXHjxrRp04bw8HBiY2P5448/iIqKIjIykqxZs5I3b17Onj1L586ddec5kXRgreWHzT/w4LAyTF68E6/RG2n3n65sWp+dF19M+/d3ySRhjMlujFlnjHnV6VhSU5kyZWjbti2VK1emUKFCbN26lWeffTbR11/+EM6ZMydDhw7l7bffJl++fEyZMoWaNWve8r0aNWrE5s2badiw4TX7/f39mT59Oj4+PrRq1eqa94G42VHlypXjiSeewMfHh06dOmGtpWHDhhQvXpwiRYpQtmxZnn766eReBhFJovWH1vPs2Cp0mNWfMxMm88q5yez8qxjNm0PmdBoscMn7SRhjugORwDZr7YKbPG9vFrfuUXDV/v37KVOmDIcPHyZnzpxOh3MD/axEEne5lMZPW+aTdc0XlLv0HoMHZeLhh1N2Xpe8n4QxJsgYc8QYs+m6/dWNMcHGmBBjTMcE+18EtgHHAPVnJENsbCwDBgygTp06LpkgROTmomKiGPDbAMoML8uyBXnJ+30woz9uyuJFKU8QyZUeDZZxwDBg4uUdxhgvYDhQDTgErDPGzLbWBgN+QHbgEeAckLyqVB7q3LlzFCxYkJIlS7Jw4UKnwxGRJLDWsmDnAvwXtib2WGnslDU0/7g0LYaCt7ezsaV5krDWrjbGXL/mrxKw01obCmCMmQLUBIKttZ/H72sIHE/r+DKa7NmzExkZ6XQYIpJEwceDabWoNX/v28vFOUN4u+IrfPEbFCjgdGRxnFonUQTYn2D7AHGJ4wpr7URuIWHtJj8/P/z8/FIvOhGRNHbqwim6L+/OuA2TyL4+gAf+/ZSh32Th0UdT7z2WL1+e4kKo6TJwHd+SmGutLR+/XQt42VrbLH67PlDJWtsyiefTwLWb089KPFVMbAxBG4P4bGlXch6syaUlPRnYswBvvQVpPas8OQPXTrUkDgLFE2wXjd+XZIGBgWpBiIhbWRm6khbz/Tl5JBdRUxfStP6jtN0I2bKl7fumpEWRXi2JEsS1JMrFb2cCdhA3cB0OrAXqWmu3J/F8N21JlChRgtDQ0FSKWtKSr69vim90JOIu4kppdGDZjj+IXdyPV0u8Td8+hiJF0jcOl2xJGGMmEzdjyccYEwZ0s9aOM8a0AJYQNw03KKkJ4lb0oSMiruRc9Dn6ru7LkN9HkGNLS0rsHcfwQdlJUDfT5aXH7KZ6iexfCCR7jqa6m0TEVVlrmbp1Km0XdSDL4We5a85G+nQpxrvvgpcDdS5cvrsptSXW3SQi4rQN4Rv4dH5L9h08x5npQ2n532fp1AlcYV2rS3Y3pRW1JETElRw9e5SAXz7jx03z8Fr+BVXzvke/uZkoWdLpyNSSEBFxTFRMFMP+HMYXK/qQY2cj8m7qwvABd/P8805HdiOPakmIiDhtfsh8Wi5oTfTh0mSatZqubR6k6WjIlMnpyFKP2yYJdTeJiFOCjwfTamEbNuzdw4XZQ/jA7xW6/AF58jgd2c2pu0lEJB3EldLowdj13+H9RwBP0JxB/b158EGnI0sadTeJiKSBmNgYxm4cS+efu3BX6BsU+H0rw/oUoHp1pyNLe0oSIiK3sDJ0JZ/O8+ff8FzE/LiQDh8/ysdDIUsWpyNLH26bJDQmISJpKSwijHaL27N0xx9cWvgV9R99hx6rDPnzOx3ZndOYhIhIKjkXfY6v1nzFoDXDyfp3Cx4+1Z5hA7NTrpzTkaWcxiRERJLpcimNNgs7YA48Q55fNzCkR3Fq1kz7Et6uTElCRDzehvANfDrPn937z3Jh1vcEvFuFVusga1anI3Oe2yYJjUmISEodPXuUgKWfMWPTPOyynvyvZGP6LMlEoUJOR5a6NCYhInIHLpfS6Lm8D1m3N6JEWBdGDLib//zH6cjSlsYkRERuY8HOBbSY15rzB+8n2+LVDPrsQWrX9uxxh1tRkhARjxB8PBj/BW34a89uouYMol3NV2n/J2TP7nRkrk1JQkQytFMXTtFjRU/GrJ2I12+deTnvLPrP9qZYMacjcw9KEiKSIV0updFpSRcy736Dktu28nW/AjzzjNORuRe3TRKa3SQiiVkVuopP5vpz9GAOvOYvoE/rx2j0rTO3DnUFmt0kIkJcKY22izrwc/BvxCzqxyfPv8NnAYbcuZ2OzDVodpOIeKS4Uhr9GLh6GJnWt6CKGcvgSdkpVcrpyNyfkoSIuC1rLdO2TqP1gg5cCn2Ke9dt4OvexalWzenIMg4lCRFxSxvDN/LJXH92hZ3h0rxJ9PqwCh8Mgcz6VEtVupwi4laOnj1K56WfM/2fudhfetL40cZ0/zUTefM6HVnGpCQhIm4hKiaK4WuH031ZbzJtbsgTkdsZPioPZco4HVnG5rZJQlNgRTzHgp0LaD63NWf3lyLfqtUM7/4gr76qUhpJpSmwIpIh7Ti+gxbz27Bu9y5iFwyiW71X+fRT8PZ2OjL3pCmwIpIhRFyIIHB5D8asnQirO1O75Ex6LfCmQAGnI/M8ShIi4jJiYmMY9/c4Oi7qgg15nXIHt/B1v4JUrOh0ZJ5LSUJEXMKq0FV8PMefw2E5yLpsPsMCHuPNNzXu4DQlCRFxVFhEGG0WdowvpfEVnWu8Q9s1hrvucjoyASUJEXHIuehzfLW6HwNWD4N1n/JGviD6/ZSdwoWdjkwSUpIQkXRlrWX6tum0nNeeqN1PcX/IBkb1Lc6TTzodmdyMkoSIpJuN4Rv5aLY/O0MjyfTzJIa2qkK9UZ5bwtsduFySMMY8BPgDPsAya+0oh0MSkRQ6evYonZZ8zvRNc4j9pSetnm9C5+WZyJnT6cjkdlx2MZ0xxgATrLUNb/KcFtOJuIGomCiG/RlXSoN/GvBCpq4M7pOHEiWcjswzueRiOmNMEPA6cMRaWz7B/urAYMALCLLW9k3wXA3gI+C7tI5PRNLGwp0L+Xh2a06HlqTI36v45suHeO45p6OSO5XmLQljzLPAGWDi5SRhjPECQoBqwCFgHVDHWht83bHzrLWv3+ScakmIuKiQf0P4ZE5r1u3ehdeSQXz1was0aQKZMjkdmbhkS8Jau9oY43vd7krATmttKIAxZgpQEwg2xjwPvAlkBeandXwikjoiLkQQ+GtPRq+dgF3diQ8rzqTbz97cfbfTkUlKODVwXQTYn2D7AHGJA2vtCmDF7U4QGBh45bGqwYo4JyY2hrEb40ppXNr+Os+c28LwrwvywANORyYpqf56WboMXMe3JOYm6G6qBbxsrW0Wv10fqGStbZnE86m7ScQFrA5bTbOZ/hwKzUa+tUMY1fVxXnrJ6agkMS7Z3ZSIg0DxBNtF4/clme4nIeKc/RH78Z/fgZ+D12CWfsUXdWrzcV9DlixORyY34/L3kzDGlCCuJVEufjsTsIO4getwYC1Q11q7PYnnU0tCxAHnos/Rd1V/Bqweiv3zU+qX7ECv7tnx8XE6MkkKl2xJGGMmA36AjzEmDOhmrR1njGkBLOHqFNgkJYjL1JIQST+XS2m0mNOe87sqU/HIekb19aVsWacjk6Rw+ZZEalNLQiT9bAzfSLOZ/oSERpJz1RBGdnqOGjVUwtsduWRLQkTc07Gzx2i/6DOmb5qDWd6DLq81pdWvmcia1enIJD25bZJQd5NI2ogrpTGCwF96EbuxAW/6BNNvah4KFXI6MkkudTeJSKpYtGsRH85sxcm9JXlg9yBG936Ixx5zOipJLepuEpFkCfk3hI9ntWHt7hCyrRjEty1epfYIo3EHcd8koe4mkZSLuBBB1196MnrdeFjdmXZVfqLTz95kz+50ZJKa1N0kInckrpTGeNov/Jzora/xcpYvGfJlQYoVczoySUvqbhKR21odtpr3f/Tn4L5sFNs6jzE9Hufpp52OSlyV2yYJdTeJ3Jn9EftpObcjS4JX472yL8Oa1KHhIKNbh3oAdTeJSKLOR5+n14p+DFwzFLu2OZ9U6EC3zjnIlcvpyCS9qbtJRK6IK6Uxg09nt+dsyJM8d2E9I4b5ct99Tkcm7kRJQiQD+vvw3zSd4c+OfREU2DCBqV2ep2pVp6MSd6TeSJEM5NjZYzSc/iFPf12dnT++y1el1hOyRAlCks9tWxIauBa5KjommiG/jyBw2ZfEbKzPeyWD6fVTHvLmdToycQUauBbxYAt3LqLZT605sceXx48N4tteZXjoIaejElekgWsRDxLybwjNfmzDuj07yLt2ENPavsZrr6mOhqSu245JGGMevsk+vzSJRkRuK+JCBC3mtqfC0KdZN8OP7gW3snfx60oQkiaS0pKYZoz5DvgKuCv+3/8AT6VlYCJyrVgby5j142i/8HOitr7KOz5b6f99Qe65x+nIJCNLSpJ4EugL/AbkAr4HnknLoJJCA9fiSdaEraHJdH8O7MtKmdC5jP3iP5Qv73RU4i7SdODaGOMNfAn8H5AT+NxaOyVZ75ZKNHAtnuLA6QN8MrMDS0NWkfOPrxj5SR3efFMlvCV5kjNwnZR1EuuA88ATQBWgrjFmejLiE5EkOh99ns9/7knpgRVZOu0BOuYIJmx+XWrVUoKQ9JWU7qam1tq/4h+HAzWNMQ3SMCYRj2WtZdrWGTSf1Z6zIZV4zfsvho4uQeHCTkcmnkrrJERcxN+H/6bxVH92hEZQcscQxnV7nkqVnI5KMpK06m4SkTR07Owx6k/5iKdGVGffvHp8+/h6tsxTghDXoMV0Ig6Jjolm0JoRdF/2JTF/16dlhe10m5WXHDmcjkzkKiUJEQcs2rmYpjNacWJPcapeXMnIIWXw9XU6KpEbuW2S0DoJcUc7/91Jk2lt+GtfMIU3DWJx59d47jlNV5K0pQJ/Ii7u9MXTdF74BUEbxuK9tiNfvdWSDxpnJVMmpyMTT6ICfyIuJtbGMnrd+PhSGq/QpOQW+kwrRO7cTkcmkjRKEiJpZE3YbzSa0pIDoVl58sQcxnz5Hx54wOmoRO6MkoRIKjtw+gAfzujIL7tWkn9jX+a0rstLL2ncQdyT1kmIpJLz0efpvLAnDwyoyPKfStGrcDD75tZTghC3ppaESApZa5myOa6Uxrmdlajj8xcDJ5YgXz6nIxNJOSUJkRT45/A/NJjsT0jYScoeGM+E7n488ojTUYmkHpdMEsaYmsBrxN2/Yqy19meHQxK5xrGzx2gxqwszt8/k7g09mPLR+9SskUkVWiXDcckkYa2dDcw2xuQB+gFKEuISomOiGbDya3os/wK76V0+ezqYjjPzkjWr05GJpI10SRLGmCDgdeCItbZ8gv3VgcHEDaAHWWv7Xnfo58CI9IhR5HYWhCym6bRWnNhbnBreKxgx6mEKFnQ6KpG0lS4rro0xzwJngImXk4QxxgsIAaoBh4i7uVEda21w/PN9gCXW2mU3OZ9WXEu62fnvTt77oS1/hW7n/t0D+e7z13nsMfUriftx2VLh1trVwMnrdlcCdlprQ6210cAUoCaAMaYFccnjLWNMs/SIUeR6py+e5sMZHSg7+Cm2LarC2Me3sOXHGkoQ4lGcHJMoAuxPsH2AuMSBtXYYMMyJoERibSzf/DmBDos+I2pbdVqV20Lg9EJky+Z0ZCLpzyUHrpMiMDDwymNVg5XUsjr0NxpMbsmh/d68EDWb0QOfoGhRp6MSSZ6UVH+9LN2qwBpjfIG5CcYkKgOB1trq8dudAHuTweubnUtjEpKqDpw+QNMpHVm+ZyVFd/RlUoe6PPWUupUkY3HZMYl4Jv7rsnXA/cYYX2OMN1AHmJPUkwUGBqY4Q4qcjz5P+3k9KdW/Ar/Nv4/hDwWz88d6ShCSoSxfvvya3pc7kV6zmyYDfoAPcAToZq0dZ4x5hWunwPZJ4vnUkpAUsdbywz8/0nx2O87ufIIPfPvRt3MJcuZ0OjKRtOOy95Ow1tZLZP9CYGFyzqk700ly/R3+D/UntSLkwAmeOjWe8T39KFnS6ahE0o7uTCeSBMfPHefj6V2YHfITBbZ2Z4L/+1Sr6rZzN0TumKuPSaQqjUlIUkXHRPPF0iEU61OGBXOy0rdYMKE/fqQEIR7D5cckUptaEpJU84OX0GRaK06GFqVu3sEM+fxh8uRxOioRZ7jsmIRIegs5vpMGk9qy8cA2yocPYkWX13noIc1YErlTbpskNHAtN3P64mlaz/yS77YEkWdrB2Z8MJ03XlWJVvFsGrgWjxdrYxmxZgKdf/6MS8HV+axyLzp9WogsWZyOTMR1qLtJPNLKfb/R4Ht/Dh3ITM0ssxk14gny53c6KpGMQUlC3NaB0wdo9F1HVoatoHRYX9Z3rkf58hp3EElNmgIrbud89Hlaz/yCUv0qsmHZfUysFMyWye8qQYgkQlNgxSNYa5m04SdazGnH+T2P4/9wP3q0LslddzkdmYh70JiEZFh/h2+izgR/dh/6lxejxzL2q6rce6/TUYlkfEoS4tKOnzvO+5O7MH/PTxTf3Z1Vbd+nciX92oqkF41JiEuKjomm+6KhFO39ML8s8ebrh7eza8pHShAiyaAxCclQ5mxdQtPprYjYX5QPig2mX4eHyZ7d6ahE3J/GJMSt7fx3F3XGtWVT+FaeOj2Q77rWwNdXM5ZEnKQkIY6LvBhJ82lf8ENwEAV2dWBx82m88JxKaYi4AiUJcUysjWXoiokE/BIAu1+mb9XNtOp5L15uO1ImkvEoSYgjlu/+nfrft+TIoczUyzubYaOeIHdup6MSkeu5bZJQFVj3tD/iAPXHd2LNweVUONaXZQF1Kf2Amg4iaUlVYMXlnY8+T4dZA/nmn0Hk2fURY5t04vWXcjodlohH0ewmcTnWWsav/Qn/ee2ICn2cz59YR0BgSTLrN0/ELei/qqSZ9Qc3UWdcK/YeOc4b3kGMGfwC+fI5HZWI3AklCUl1x88d570JXVkU9iOlDwWyof0HlC+rXzURd6T/uZJqomOiCVwwkv5/fkG23XWYVGc7td/Ih9F6OBG3pSQhqWLmPz/z/o+tiDxYmFYP/coXEx7B29vpqEQkpZQkJEV2HNvFO0Ft2Xp0K9UuDWRirxoULKimg0hG4bZJQusknBV5MZJm33/J9F1jKH6gPb/5T6PS4yqlIeKKtE5C0k2sjWXA0ol0XRFA5n0vM/C1Xrxf516NO4i4Aa2TkDS1dMfvNJjckuNHM/F+kVkMHF2JbNmcjkpE0pKShNzW/lMHqT2mE396kVPcAAAL0UlEQVQe+ZWnzvZh7Wf1KFZUpTREPIGShCTqwqUL+E8dwNhtgygQ9hFLPwqm6jMqpSHiSZQk5AbWWsasmUmbRe24dOBR+jy3ltY971MJbxEPpCQh11gbuok641sRdvwYdfKMYdTwF8ipxoOIx1KSEACOn/2Xd4O6svTQDCqc7MaSjs24/z79eoh4OpfrQDDGlDTGjDHGTHM6Fk8QHRNNhx+HUbhXGTb8lYlZL25nwzefKEGICODC6ySMMdOste8k8pzWSaSC6euX0mymP+eOFCbg0cF8/uEjZMrkdFQiklZccp2EMSYIeB04Yq0tn2B/dWAwca2ZIGtt37SOReJsO7yLd8a0Y/uJzdTwHsj4fm+QJ49Ww4nIjdKju2kc8HLCHcYYL2B4/P5HgLrGmIeuO06fWqns9IVI3hrViXJDKhMTWpl/PtjGrD41lSBEJFFpniSstauBk9ftrgTstNaGWmujgSlATQBjTD5jzEigojGmY1rH5wlibSxfzh/PPT0fZOkfh/nu6U1sH92JsmVUa0lEbs2p0ckiwP4E2weISxxYa08AH9/uBIGBgVceq9Bf4hZv+4OGk1ty4l8v/EvPondgJbJkcToqEUkPKSnsd1m6DFwbY3yBuZfHJIwxtYCXrbXN4rfrA5WstS2TeD4NXN9G6ImDvP1NJ9af+JUXYnszufO73JPf5SaziUg6csmB60QcBIon2C4avy/JVCr85i5cusAn3w1k4s6BFDvyIb+3CKZSRa2GE/FkLl8q3BhTgriWRLn47UzADqAaEA6sBepaa7cn8XxqSVzHWsuIX2fScWk7zOFHGVC9H83evk8lvEXkCpdsSRhjJgN+gI8xJgzoZq0dZ4xpASzh6hTYJCWIy9SSuOr3PZupM74VB08epUnh0QzrVo2sGpMWkXgu35JIbWpJxDka+S+1v+nKymMzePJCN6Z1bEbRwlopLSI355ItCUl9l2Iv0W7qKEZs7sE9R2rz8/vbeeGpfE6HJSIZkNsmCU/tbvr+96V8MqcVF08Uovczy2j7ZVmNO4jILam7yQNsPribt0a3Zdfpzbx99wCC2tUkRw5lBxFJOnU3ZUCnL0Ty7jdfsuDwGMpFtiO47RQeKHmX02GJiIdw2ySR0bubYm0s3Wd+R+91AeQ+9iI/1dtEzRcKOx2WiLghdTdlMHM3/kGTaf5EREDHCkPp/sGTunWoiKSYupvc3J5jh6g1shP/RP7Ca959mNTnXe7OrewgIs5RknAB56Mv0HTMQKaGDeT+iA/5x38H5R5UKQ0RcZ7bJomMMCZhrWXwolkErGiL94mKTKi5lvqv3ed0WCKSwWhMwg2tCN5MvYmtOHLmKM1LDWZA82pkdtuULSLuQGMSbiD81L+8/XU3fouYhh/d+Lvbh9zjox+DiLgmfTqlk0uxl2gxYRSjQ3pQ5FRt1ny0nacq+DgdlojILbltknCnMYmxy5fiv7AVMacLMrTaMj55q6zTIYmIB9GYhIvasHc3b49py77zm2lUaAAj/WuSNatKaYiIMzQm4SJOno2k7te9WPLvaJ6Ibsvy9lModq9KaYiI+1GSSEWxNpaAKZMY8E9nfCJeZHGjTfxfZZXSEBH3pSSRSn7680/e/7El585Bjyd/pFP9yirhLSJuT0kihULCD/HmyE5sO/8Lte7uzYTu9cmeTaU0RCRjcNsk4fTspnNRF2g4chAzDw/g4fPNCG4VTOkSuRyJRUTkVjS7KR1Za+kzexaBv7Ulx5kKBL0zgP/5qZSGiLg+zW5KY0s3baH+pFb8e/Ew7cp+S6/3X9S4g4hkaEoSSXDg3xPUGt6Vdeem8XK2bkzp8iF359KlE5GMT590txAdc4kPR3/DhH3due/8O2z4ZDsVH1QpDRHxHEoSiRi5+BfaLfXH63xBgl77hfdeKed0SCIi6U5J4jp/7NhN7bHtOHjpHz4sOYAhH/2XzJk18CAinsltk0RqT4E9fjqSd4b1ZnnktzybqS1rO/xAQR+V0hAR96cpsCkQExtL2wmTGB7cmXvPV2Ny095UqVAkVc4tIuJKNAX2Dk1e+Scfz/YnKtoy0O9HWr5Z2emQRERcikcmiS2hh3hrVGd2XlpKvcK9CWpZH+8sKqUhInI9j0oSkecvUG/YIOafGMCj9gN2twmmxL0qpSEikhiPSBLWWgKnzKb3hrbkuVie+e/+yStPlnI6LBERl5fhk8S8tVtoPLUVp2MP8/nj39DlXZXSEBFJqgybJPYdOcH/hnbjn0tTecOnK9+3/ogc2TLstysikiYy3KfmxehLNB7xDVMPd+ehmHfY+ul2yviqlIaISHK4XJIwxmQHvgYuAiustZOTeuygWb8QsMqfrJcKMPl/v1DbT6U0RERSwhXnfb4JTLfWfgi8kZQDVm7eQ9G2b9Jh9Qc0L9OTE4M8J0EkdxVlRqRrcZWuxVW6FimT5knCGBNkjDlijNl03f7qxphgY0yIMaZjgqeKAvvjH8fc6tyHT5zhmW4B+H3/BA/d/R+OdttG//f/h5eX54xM6z/AVboWV+laXKVrkTLp0ZIYB7yccIcxxgsYHr//EaCuMeah+Kf3E5coABL9tP9o5HcU6f0Q4WcO8GfjTSztGkDeXKq1JCKSmtJ8TMJau9oY43vd7krATmttKIAxZgpQEwgGZgLDjTGvAXMTO++kkGGMfGk6zV55Ko0iFxGRdCnwF58k5lpry8dv1wJettY2i9+uD1Sy1rZM4vncryqhiIgL8IgCf3f6TYqISPI4NbvpIFA8wXbR+H0iIuJC0itJGK4dhF4H3G+M8TXGeAN1gDnpFIuIiCRRekyBnQz8BpQ2xoQZYxpba2OAFsASYCswxVq7Pa1jERGRO+N2d6YzxlQHBhOX4IKstX0dDskRxpiiwESgIBALjLbWDnU2KmfFT63+CzhgrU3SQsyMyBhzNzAGKEvc70YTa+2fzkblDGNMa6ApcddhM9DYWhvlbFTpwxgTBLwOHEkwaSgvMBXwBfYB71hrI251HldccZ2o26yv8DSXgDbW2keAp4DmHnwtLvMHtjkdhAsYAiyw1pYBKgAe2Uo3xhQmrsfisfgPyczEdW17ihvWqAGdgKXW2geBZUDn253ErZIECdZXWGujgcvrKzyOtfawtfbv+MdniPsg8Nibc8e3rF4l7i9oj2WMyQ1UsdaOA7DWXrLWnnY4LCdlAnIYYzID2YFDDseTbqy1q4GT1+2uCUyIfzwB+O/tzuNuSaIIV0t2ABzAgz8YLzPGlAAqAh7ZpRBvENAecK/+09RXEjhujBlnjNlgjPnWGJPN6aCcYK09BAwAwoibPXnKWrvU2agcV8BaewTi/tAECtzuAHdLEnIdY0xOYAbgH9+i8Djxq/OPxLesrp9J52kyA48BI6y1jwHniOti8DjGmDzE/eXsCxQGchpj6jkblcu57R9V7pYktL4igfgm9AzgO2vtbKfjcdAzwBvGmD3AD0BVY8xEh2NyygFgv7X2r/jtGcQlDU/0IrDHWnsifkblT8DTDsfktCPGmIIAxphCwNHbHeBuSULrK641FthmrR3idCBOstYGWGuLW2vvI+53Ypm1tqHTcTkhvithvzGmdPyuanjuYH4YUNkYc5cxxhB3LTxtEP/6lvUc4L34x42A2/5x6VZlOay1McaYT4lbX3F5Cqyn/dABMMY8A7wLbDbGbCSu2RhgrV3kbGTiAloC3xtjsgB7gMYOx+MIa+1aY8wMYCMQHf/vt85GlX7i16j5AT7GmDCgG9AHmG6MaQKEAu/c9jzutk5CRETSj7t1N4mISDpSkhARkUQpSYiISKKUJEREJFFKEiIikiglCRERSZSShEgqMcbcbYz52Ok4RFKTkoRI6skLfOJ0ECKpSUlCJPX0Bu6Lr77qkTfDkoxHK65FUokxxheYe/kuYCIZgVoSIiKSKCUJERFJlJKESOqJBHI5HYRIalKSEEkl1toTwBpjzCYNXEtGoYFrERFJlFoSIiKSKCUJERFJlJKEiIgkSklCREQSpSQhIiKJUpIQEZFEKUmIiEii/h9Bb34badyvMgAAAABJRU5ErkJggg==", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# グラフに結果を描画する。\n", "plt.plot(t_list, x_list, label='simulation')\n", "plt.plot(t_list, x_exact, label='analytical')\n", "plt.xlabel('t')\n", "plt.ylabel('x')\n", "plt.yscale('log')\n", "# 凡例を表示する。\n", "plt.legend(loc='best')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 課題 1\n", "\n", "上記の例について, $\\Delta t$ = 0.01, 0.1 および $\\Delta t$ = 1.0 についてそれぞれ計算し, その結果を描画することで\n", "$\\Delta t$ の大きさが微分方程式の数値解にどのような影響を及ぼすか考察せよ." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 課題 2\n", "\n", "あらかじめ解析解の分かっている常微分方程式を何でもよいので一つ取り上げ, これを数値的に解いた場合の誤差を評価せよ.\n", "\n", "例えば、減衰振動の微分方程式は以下のように表される。 \n", "$$\n", "\\frac{d^2x}{dt^2} + 2\\zeta\\omega_0\\frac{dx}{dt} + \\omega_0^2 x = 0\n", "$$\n", "\n", "この常微分方程式の解析解は\n", "0 < ζ < 1のとき以下のように表される.\n", "$$\n", "x(t)=Ce^{{-\\zeta \\omega _{0}t}}\\cos \\left(\\omega _{0}{\\sqrt {1-\\zeta ^{2}}}t-\\alpha \\right)\n", "$$\n", "\n", "ただし、\n", "$$\n", "C=x_{0}{\\sqrt {1+\\left({\\frac {\\sigma +\\zeta }{{\\sqrt {1-\\zeta ^{2}}}}}\\right)^{2}}}\n", "$$\n", "および\n", "$$\n", "\\alpha =\\tan ^{{-1}}\\left(-{\\frac {\\sigma +\\zeta }{{\\sqrt {1-\\zeta ^{2}}}}}\\right)\n", "$$\n", "である。\n", "\n", "[https://ja.wikipedia.org/wiki/%E6%B8%9B%E8%A1%B0%E6%8C%AF%E5%8B%95]" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Lorenz 方程式\n", "\n", "気象学者のLorenzは1963年に熱対流の近似モデルとして\n", "以下の方程式を提案した(Lorenz方程式)\n", "> Lorenz, E. N., 'Deterministic nonperiodic flow', *J. Atms. Sci.* **20**, pp.130-141, (1963).\n", "\n", "$$\n", "\\begin{array}{ccc}\n", "\\dfrac{dx}{dt}&=&-ax+ay \\\\[3pt]\n", "\\dfrac{dy}{dt}&=&\\mu x-y-xz \\\\[3pt]\n", "\\dfrac{dz}{dt}&=&-bz+xy \\\\[3pt]\n", "\\end{array}\n", "$$\n", "\n", "Lorenzは数値計算により, この方程式の解が\n", "不規則で周期性をもたない振動をすることを発見した.\n", "決定論的な微分方程式の解がこのような予測不可能は振舞いを示すことは驚きをもって受け止められ,\n", "現在ではこのような非周期運動はカオスと呼ばれている.\n", "\n", "なお, Lorenz方程式は自明解 $x=y=z=0$ の他に,\n", "$\\mu >1$ のとき, 定常解 $x=y=\\pm \\sqrt{b(\\mu-1)}, z=\\mu-1$ をもつ.\n", "\n", "\n", "\n", "### 課題 3\n", "+ $a=10, b=8/3$ と固定した上で, Lorenz方程式のシミュレーションを実施せよ. \n", "$\\mu$ を $0, 1, 3, 10, 30,...$ の条件でシミュレーションを行い, どのような解が現れるか調べよ. \n", "$x, y, z$ の時間変化を示すグラフや上のような $xyz$ 空間における軌跡などをプロットしてみるとよい.\n", "+ 初期値がわずかに(例えば1％程度)異なる2ケースのシミュレーション結果を比較せよ. \n", "\n", "なお、3次元グラフは\n", "```python\n", "from mpl_toolkits.mplot3d import Axes3D\n", "\n", "ax = plt.subplot(111, projection='3d')\n", "ax.plot(x, y, z)\n", "```\n", "のようにして表示することができる。\n", "\n", "より大きく描画するためには\n", "\n", "```python\n", "from mpl_toolkits.mplot3d import Axes3D\n", "\n", "plt.figure(figsize=(10, 10))\n", "ax = plt.subplot(111, projection='3d')\n", "ax.plot(x, y, z)\n", "```\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] } ], "metadata": { "anaconda-cloud": {}, "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.6.3" } }, "nbformat": 4, "nbformat_minor": 0 }