{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pythonの応用2 〜運動方程式の数値解〜\n",
    "\n",
    "前回の演習では一階の連立微分方程式の数値シミュレーションを取り扱った.\n",
    "今回はこれを応用して, 二階の微分方程式である運動方程式の解法を学び,\n",
    "スポーツの一場面を力学的に解析してみよう."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 物体の自由落下のシミュレーション\n",
    "\n",
    "ここでは, 以下のような質点の自由落下をシミュレートする.\n",
    "\n",
    "<img src=\"figs/freefall.png\" width=\"80%\">\n",
    "\n",
    "自由落下運動の模式図. 時刻 $t = 0$ に $x(0)=h$ で鉛直方向に速度 $u_0$ で\n",
    "運動している物体（図の赤丸）の軌道をシミュレートする.\n",
    "位置 $x$ および速度 $u$ の模式図は中・右図のようになる."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "質点の質量を $m$ 、重力加速度を $g$ 、初期位置を $h$ 、初速を $u_0$ とする。\n",
    "鉛直座標として $x$ をとると, 運動方程式は\n",
    "\n",
    "$$\n",
    "m\\frac{d^2x}{dt^2}=-mg, \\ \\ \\ x(0)=h, \\ \\ \\ \\frac{dx}{dt}(0)=u_0\n",
    "$$\n",
    "\n",
    "で与えられる.\n",
    "その解は\n",
    "\n",
    "$$\n",
    "x(t)=h+u_0t-\\frac{1}{2}gt^2.\n",
    "$$\n",
    "\n",
    "であるが, これをコンピュータを用いて近似的に求めてみよう.\n",
    "\n",
    "上式を以下のように変形することで一階の連立微分方程式が得られる.\n",
    "\n",
    "$$\n",
    "\\frac{dx}{dt}=u, \\ \\ \\ x(0)=h, \\\\[3pt]\n",
    "\\frac{du}{dt}=-g, \\ \\ \\ u(0)=u_0.\n",
    "$$\n",
    "\n",
    "これを前回学習したEuler法を用いて数値的に解けばよい. つまり，以下のように近似を行う．\n",
    "\n",
    "$$\n",
    "x_{n+1} = x_{n} + u_n \\times \\Delta t \\\\[3pt]\n",
    "u_{n+1} = u_n - g \\times \\Delta t\n",
    "$$\n",
    "\n",
    "プログラムの例を以下に示す.\n",
    "計算した $x, u$ の時間変化をグラフにし,　厳密解と比較してみよ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def next_x(x, u, dt):\n",
    "    \"\"\"\n",
    "    現在の x, u の値から、dt 後のxを求める関数。\n",
    "    args:\n",
    "    + x: 時刻 t での空間座標\n",
    "    + u: 時刻 t での速度\n",
    "    + dt: 線形近似する時間幅\n",
    "    return:\n",
    "     時刻t+dtでのxの値\n",
    "    \"\"\"\n",
    "    return x + u * dt\n",
    "\n",
    "def next_u_vacuum(u, dt, g):\n",
    "    \"\"\"\n",
    "    現在の u の値から、dt 後のxを求める関数。空気抵抗は考慮しない。\n",
    "    args:\n",
    "    + x: 時刻 t での空間座標\n",
    "    + u: 時刻 t での速度\n",
    "    + dt: 線形近似する時間幅\n",
    "    + g: 重力加速度の値\n",
    "    return:\n",
    "     時刻t+dt での u の値\n",
    "    \"\"\"\n",
    "    return u - g * dt\n",
    "\n",
    "def solve_freefall_vacuum(x0, u0, dt, tmax):\n",
    "    \"\"\"\n",
    "    t=0 に位置h、速度u0の物体の自由落下をシミュレートする関数。\n",
    "    args:\n",
    "    + x0: t=0での座標 [m]\n",
    "    + u0: t=0での速度 [m/s]\n",
    "    + dt: オイラー法を用いるための時間軸の刻み幅 [s]\n",
    "    + tmax: 追跡を行う最大時間 [s]\n",
    "    \"\"\"\n",
    "    # 計算結果を代入するリスト\n",
    "    x_list = []\n",
    "    u_list = []\n",
    "    t_list = []\n",
    "    \n",
    "    g = 9.81 # 重力加速度 [m/s^2]\n",
    "    \n",
    "    # 初期値の代入\n",
    "    x = x0\n",
    "    u = u0\n",
    "    t = 0.0\n",
    "    \n",
    "    while(t < tmax):\n",
    "        x_list.append(x)\n",
    "        u_list.append(u)\n",
    "        t_list.append(t)\n",
    "        \n",
    "        # x, u, tの値の更新\n",
    "        x = next_x(x,u,dt)\n",
    "        u = next_u_vacuum(u,dt,g)\n",
    "        t = t + dt\n",
    "    return x_list, u_list, t_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAADQCAYAAAAK/RswAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VGX6xvHvk0ZICITeQgfpECD0YpcmIogCKoKoWFHU1VW36K7rrmWVLogVUVFUFBUBBUWaCKH3jlTpvYbw/v7I+NssGyBAZs7M5P5c17mSOWXOzXjy+sw5532POecQERERCScRXgcQERERyWkqcERERCTsqMARERGRsKMCR0RERMKOChwREREJOypwREREJOyowBEREZGwowJHREREwo4KHBEREQk7UV4HuBBFihRx5cuX9zqGiPjMmzdvt3OuqNc5/EntjkhwyW67E1IFTvny5UlNTfU6hoj4mNmvXmfwN7U7IsElu+2OLlGJiIhI2FGBIyIiImFHBY6IhAUza2Nmq8xsrZk9lcVyM7NBvuWLzay+FzlFJDBU4IhIyDOzSGAo0BaoAXQ3sxpnrNYWqOKb+gDDAhpSRAIqpG4yzo75m/Yxdv4WoiIiiI40oiIjiI7I+JknKoKE2GjyxUaREBtFQp4oEmKjKRgfTeH4PERGmNfxReTiNALWOufWA5jZx0BHYHmmdToC7zvnHDDbzBLNrKRzbvul7vzHlTuZv2kfD19dhehIfW8UCQZhV+Bs2XeMb5f8Rlr6aU6lO06dPk1aujvvdhEGhfPloVhCxlQ0IQ9JBeMoWyiOsoUzfhaOj8FMRZBIECoNbM70egvQOBvrlAb+p8Axsz5knOWhbNmy5935z+v3MGLaeqau2kX/rslULpbvwtKLSI4LuwLnhrqluKFuqf+a55wj/bTj+KnTHD5+ikPH0zh04hSHfL/vO3KSnYdOsOvQif//uWzbQXYeOvFf7xMfE0nZwvFcVjwflxVPoFqJBC4rnkBSwbwqfETCiHNuBDACICUl5bzfkJ5pV516ZRJ55osltB80nWfaVeeOpuXULoh4KOwKnKyYGVGRRr7ICPLliaJEgdhsbXc8LZ0t+47y656jbNqb8XPjniOkbtzHuIXb/n+9fHmiqFoigbpJidQtU4DkMomULRSnxk0kcLYCZTK9TvLNu9B1Llrb2iVpUK4gT36+mGe/WsaUlTt5pUsdiufPXnsjIjkrVxQ4Fys2OpLKxRKoXCzhf5YdPJ7Gmh2HWPXbYVb9dpBl2w7y0ZxfeWfmaQAS46Kpm5RI/bIFaVqpMHXLFCBPVGSg/wkiucVcoIqZVSCjaOkG3HrGOl8BD/nuz2kMHMiJ+28yK5Y/lnd7NeSDXzbxwvjltB4wjX92qk272iVzcjcikg0qcC5S/thoGpQrRINyhf5/3qn006zacYhFmw+waPN+Fm7ez4Apq+k/GfJERdCgXEGaVCycUfAkJRITpZsRRXKCc+6UmT0ETAIigXecc8vM7D7f8uHAt0A7YC1wFLjTH1nMjB5NytGsUmEe+2QhD3w4n871SvNcx5rkj432xy5FJAuW0aEgNKSkpLhQGzJ9/9GTzNmwl9nr9zJ7/R5W/HYQ5zIua7WoXISrqhXjiqpFKabT2BKCzGyecy7F6xz+dCntTlr6aQb/sJahP66lRP5YXr2lLk0qFs7hhCK5S3bbHZ3B8bPEuBiuq1mC62qWADIKntnr9/LT6l1MXbWTict+A6BmqfxcWbUY19UsTu3SBXT/jkgYiI6M4LFrL+OKqkV57JOFdH9zNn1aVuSx6y7TJWsRP9MZHA8551i14xA/rNzJ1JW7mLdpH+mnHUkF89Kudkna1ipBcplEFTsStHQGJ/uOnDjFC9+u4KNfNlGtRAIDuiVTrUT+HEgokrtkt91RgRNE9h89yffLd/Dtku3MWLubtHRHqQKxtK1dko7JpXRmR4KOCpwL98PKHTz52WIOHjvFE62rcleLCkRokFGRbFOBE+IOHEtjyoodfLvkN6at3sXJ9NNUKZaPzvWT6FSvdLa7uov4kwqci7Pn8AmeGruE75fvoEnFQrx6SzKlE/Pm6D5EwpUKnDBy4Fga3y7ZzufztpD66z7MoEXlInSuX5q2tUoSG61r+eINFTgXzznHp6lb+NvXy4iIMJ7vWIuOyaV0llbkPFTghKmNu48wdsFWxs7fwpZ9xyiQN5ouDZK4rXFZKhbV8PASWCpwLt2mPUd5bMxCUn/dR/s6JXnhxlokxsX4bX8ioU4FTpg7fdoxe8MePvxlE5OW/sap045mlQpze5NyXFujuB74JwGhAidnpJ92DP9pHf2/X03hfDH8++a6tKxS1K/7FAlVKnBykZ2HjvNp6hY++mUTW/cfo1hCHno0KcftTcpRMF7fBMV/VODkrKVbD9Dvk4Ws3XmYXs3K81TbaroELXIGFTi5UPppx9RVOxn5869MW72L2OgIujRI4q4WFalQJN7reBKGVODkvONp6bw4YSXvzdpIpaLxDOxWj1qlCwRs/yLBLrvtjqfXMcysjZmtMrO1ZvaUl1nCQWSEcXX14rzfuxHfPdqKG+qWYszcLVz16lTueT+VORv2EkoFrUhuFBsdyXM31GTUXY04fOIUNw6dydAf15J+Wn+7IhfCszM4ZhYJrAauBbaQ8bC87s655WfbRmdwLtzOQ8cZ9fOvfDD7V/YdTaNh+YL0vaoKLasUUW8NuWQ6g+Nf+4+e5E9fLmX84u00KFeQ/rckU7ZwnCdZRIJFKJzBaQSsdc6td86dBD4GOnqYJywVS4jl8euqMuupq/nbDTXZsu8Yd7wzhxtfn8WUFTt0RkckiCXGxTCkez0Gdktm9Y5DtB04jU/mbtLfrUg2eFnglAY2Z3q9xTfvv5hZHzNLNbPUXbt2BSxcuMkbE0nPZuWZ+sQVvNCpFnsOn+CukalcP3gGE5du57ROf4sEJTOjY3JpJvZrRZ2kRP74+RL6jJrH7sMnvI4mEtSCvi+xc26Ecy7FOZdStKi6TV6qPFGR3Na4HD/+4Qpe7lKHIydOcd8H87lh6Ax+Wr1L3wxFglTpxLx8eHdj/ty+Oj+t2kWbAdOYsmKH17FEgpaXBc5WoEym10m+eRIA0ZER3JJShsmPXc6rN9dl/9E0er4zh+5vzmb+pn1exxORLEREGHe3rMhXfZtTJF8e7hqZytNjl3DkxCmvo4kEHS8LnLlAFTOrYGYxQDfgKw/z5EpRkRHc1CCJKY9fznMdarB252E6vz6Le95PZfWOQ17HE5EsVCuRn3EPNefeyyvy8dxNtB80XV9MRM7gWYHjnDsFPARMAlYAY5xzy7zKk9vliYqkV/MK/PTElTx+7WXMXreH1gOm8dTni9l1SNf6RYJNnqhInm5bndH3NCEt3dFl2Cxe+24VaemnvY4mEhQ00J9kad+Rkwz5cS0jZ20kNjqSB6+szJ3Ny2tUVfkv6iYeHA4eT+O5r5Yxdv5W6iQVoH/XZCrp2XQSpkKhm7gEsYLxMfzl+hp892grmlQszEsTV3Jt/5+YsGS7bkQWCTL5Y6N57ZZkXr+tPpv2HqX9oOmM+nmj/lYlV1OBI+dUsWg+3uqZwgd3NSYuOor7P5xP1xGzWb7toNfRRAAws1fMbKWZLTazL8ws8SzrbTSzJWa20MyC+5TMRWpXuyST+rWiUYXC/GXcMnq9O5edB497HUvEEypwJFtaVCnC+Idb8EKnWqzdeZgOQ2bw/DfLOazeG+K974Fazrk6ZIyO/vQ51r3SOZcczpfViuePZeSdDfl7x5r8siHjXroJS7Z7HUsk4FTgSLZFRUZwW+Ny/PD45XRtWIZ3Zm7g6lenMn6xLluJd5xz3/k6LQDMJmPIiVzNzLijaXm+6duSMoXiuP/D+Tw2ZiEHj6d5HU0kYFTgyAVLjIvhn51qM/b+ZhTJl4cHP5pPz3fnsnH3Ea+jifQGJpxlmQMmm9k8M+sTwEyeqVwsH5/f34yHr6rMlwu20nbAdH5Zv8frWCIBoQJHLlq9sgUZ92Bznu1Qg/m/7uO6AdMYNGUNJ0+pm6rkLDObbGZLs5g6ZlrnT8Ap4MOzvE0L51wy0BZ40MxanWN/YfOImOjICB67riqf3teMqEij25uz+deEFZw4le51NBG/UjdxyRE7Dx7n798s55vF26lWIoGXu9ShTlKW93pKGAmWbuJm1gu4F7jaOXc0G+s/Bxx2zv37fOuGU7tz5MQp/jF+BaPnbKJ6yfwM6JpM1RIJXscSuSDqJi4BVSx/LENurc9bd6Sw7+hJbhw6kxcnrOR4mr4lin+ZWRvgSeCGsxU3ZhZvZgm//w5cBywNXMrgEJ8nin91rs1bd6Sw69BxOgyewVvT1+thuxKWVOBIjrqmRnG+e/Rybm5QhuE/raPdwOnM3bjX61gS3oYACcD3vi7gwwHMrJSZfetbpzgww8wWAXOA8c65id7E9d41NYozsV8rWl1WlH+MX8Ftb/3Ctv3HvI4lkqN0iUr8Zsaa3Tw1djFb9x+jZ9Py/LFNNfLGaCTkcBIsl6j8KZzbHeccn8zdzN+/WU5khPGPG2vRMbm017FEzkmXqMRzLaoUYVK/VvRsWp73Zm2k/aDpLNq83+tYIuJjZnRrVJYJj7SkSrF8PPLxQvqOXsCBo+pOLqFPBY74VXyeKJ67oSYf3d2YY2npdB42i4GT13BKDwQUCRrlCscz5t6m/OG6y5iwZDutB0xjxprdXscSuSQqcCQgmlUuwsR+rehQpyT9J6/mpuE/s37XYa9jiYhPVGQED11VhS8eaE58nkhuf/sX/vb1MnUUkJClAkcCpkDeaAZ0q8fg7vXYuPsI7QfNYNTsXzUKskgQqZ1UgG/6tqRXs/K8O3MjHQbPYOnWA17HErlgKnAk4DrULcWkfq1IKV+Qv3y5lD6j5rH/6EmvY4mIT96YSJ67oSYjezfiwLE0Or0+k6E/riVd3cklhKjAEU+UKBDLyDsb8ef21Zm6aiftBk4nVd3JRYLK5ZcVZVK/VlxXowSvTFpF1zd+ZvPe846jKBIUVOCIZyIijLtbVuSz+5oRFRlB1xGzGfrjWg06JhJECsbHMOTWevTvWpdVvx2izYBpjEndrEvLEvRU4Ijn6pZJ5JuHW9C2Vsa3xJ7vzmHXoRNexxIRHzOjU70kJj7aitpJBXjys8XcO2oeew7r71SClwocCQr5Y6MZ3L0e/+pcmzkb9tJ24HRmrlU3VZFgUjoxLx/d3YRn2lVj6qpdtB4wnR9W7vA6lkiWVOBI0DAzujcqy1cPtaBgXDQ93v5Fl6xEgkxEhNGnVSXGPdScIvli6P1eKn/6YglHT57yOprIf1GBI0GnaokEvnywOe3rlOKVSau494N5HDyukVVFgkn1kvkZ91Bz+rSqyEdzNtF+0AwWbNrndSyR/6cCR4JSfJ4oBnVL5tkONfhx5U5uGDyDlb8d9DqWiGSSJyqSZ9pV56O7m3Dy1Gm6DP+Z/t+vJk0jlUsQUIEjQcvMuLN5BUb3acKRk+l0GjqLcQu3eh1LRM7QtFJhJvRrSce6pRg4ZQ1dhs3SSOXiORU4EvQali/E+L4tqFU6P498vJDnvlqmb4giQSZ/bDSvdU1m6K312bjnKO0GTddI5eIpFTgSEorlj+Wje5rQu3kF3pu1kTvensO+Ixr9WCTYtK9Tkkn9WtGwfCH+8uVS7nxvLjsPHfc6luRCKnAkZERHRvDXDjV49ea6zPt1Hx2HzmT1jkNexxKRM5QoEMv7vRvxtxtq8vO6PbTuP42JS7d7HUtyGU8KHDN7xcxWmtliM/vCzBK9yCGh6aYGSXx8bxOOnkyn8+uzmLJC43CIBBszo2ez8ox/uCVJBeO474P5/OHTRRxSj0gJEK/O4HwP1HLO1QFWA097lENCVP2yBfm6b3PKF4nj7vdTGTZ1na71iwShysXy8fn9zXjoysqMnb+FtgOnM2eDnjsn/udJgeOc+8459/uoULOBJC9ySGgrWSAvn97bjPa1S/LSxJU8NmYRx9PSvY4lImeIiYrgD62r8ul9TYkwo+uIn3lp4kpOnlJnAfGfYLgHpzcw4WwLzayPmaWaWequXbsCGEtCQd6YSAZ3r8fj117GFwu20v3N2Xo+jkiQalCuEN8+0pKuKWUYNnUdN+o+OvEjvxU4ZjbZzJZmMXXMtM6fgFPAh2d7H+fcCOdcinMupWjRov6KKyHMzOh7dRWG316f5dsO0ul1jcGRm5jZc2a21cwW+qZ2Z1mvjZmtMrO1ZvZUoHNKhnx5onjxpjq8eUcKOw4e5/rBM3h7xgY9kkVynN8KHOfcNc65WllM4wDMrBdwPXCb080TkgPa1CqZMSjgiVN0HjaLuRt1nT8X6e+cS/ZN35650MwigaFAW6AG0N3MagQ6pPzHtTWKM7FfK1pWLsLz3yynxzu/sP3AMa9jSRjxqhdVG+BJ4Abn3FEvMkh4ql+2IGMfaEahuBhue/MXvl60zetIEhwaAWudc+udcyeBj4GO59lG/KxoQh7e6pnCvzrXZsGm/bTuP02jlUuO8eoenCFAAvC975TycI9ySBgqVziesQ80I7lMIn1HL1APqwAzs1gz62JmA83sUzN738yeNLOaftxtX9+wE++YWcEslpcGNmd6vcU3L0u69y9wzIzujcry7cMtqVQsH498vJCHRy/gwFF1J5dL41UvqsrOuTKZTinf50UOCV+JcTGMursRN9QtxUsTV/LMF0s5pcc7+J2Z/Q2YCTQFfgHeAMaQca/di2b2vZnVuYj3Pdc9fcOAikAysB149VL/Hbr3L/DKF4nn03ub8vi1l/Htku20GTiNmWt3ex1LQliU1wFE/CVPVCQDuiZTplBehv64jl2HTjDk1nrERkd6HS2czXHOPXuWZa+ZWTGg7IW+qXPumuysZ2ZvAt9ksWgrUCbT6yTfPAkiUZER9L26CpdXLUq/TxZy21u/0Lt5BZ5sU1V/t3LBgqGbuIjfREQYT7SuxvMdazJl5Q56vP0LB47p1Le/OOfGnznPzCLMLL9v+U7nXGpO7tPMSmZ62QlYmsVqc4EqZlbBzGKAbsBXOZlDck6dpETG923JHU3L8c7MDXQYPINl2w54HUtCjAocyRV6NC3P4O71WLh5P13f+JmdB/XwP38ys4/MLL+ZxZNRcCw3syf8tLuXzWyJmS0GrgQe9WUoZWbfAvgGFn0ImASsAMY455b5KY/kgLwxkfy9Yy3eu7MhB46lcePQmQybuo50dSeXbLLs3HxpZilAS6AUcIyMBut759w+/8b7bykpKS41NUe//EkuM2PNbu4dlUrB+BhG3dWYCkXivY4U0sxsnnMuJYv5C51zyWZ2G1AfeAqY53s8S0hRu+O9fUdO8swXS5iw9DcalS/Eq7fUpUyhOK9jiUfO1u6c6ZxncMzsTjObT8azovICq4CdQAtgspmNNLMLvp4u4pUWVYowuk/Ggzq7DJvF0q067e0n0WYWDdwIfOWcSwP01VsuSsH4GF6/rT6v3VKXFdsP0nbgdD5N3azekXJO57tEFQc0d87d5Jz7p3PuLefcEOfcw865BkB/oIr/Y4rknDpJiXx2X1NioyPpNmI2s9RTwx/eADYC8cA0MysHHPQ0kYQ0M6Nz/SQm9GtJjVL5eeKzxdz3wTz2HjnpdTQJUucscJxzQ51zZx1a0jm30Dk3JedjifhXxaL5GPtAM0on5qXXu3P5fvkOryOFBTNrambmnBvknCvtnGvnG6l8Exn3x4hckqSCcYy+pwlPt63GDyt3cl3/afy4cqfXsSQIZesmY1/Pg9fMbKyZffX75O9wIv5UPH8sY+5tSvVS+bn/g3ka9Thn3AHMM7OPzayXmZUAcBlOeZxNwkRkhHHv5ZUY92ALCsfHcOd7c/nzl0s4elKHmPxHdsfB+RJ4G/ga0GhpEjYKxEXzwV2NuOu9VB75eAEnTp2mS4Mkr2OFLOfc/QBmVo2M5z69Z2YFgB+BicBM51y6hxEljNQolZ9xDzXn1e9W8daMDcxau4fXuiaTXCbR62gSBLLbTfy475Tzj865n36f/JpMJEASYqMZ2bsRzSsX4Q+fLmLUzxu9jhTynHMrnXP9nXNtgKuAGcDNZIxuLJJjYqMj+VP7Gnx4d2OOp6Vz07BZDJi8WiOXS7YLnIFm9qzv+nr93ye/JhMJoLwxkbx5RwrXVC/GX8YtY8S0dV5HCnlmVtD3WIbqwG/Au9np2ilyMZpVKsKEfq3oUKckAyav4abhP7Nh9xGvY4mHsnuJqjbQg4xvYr+Xxc73WiQsxEZHMuz2BvT7ZCH//HYlx06e5uGrK2NmXkcLOWb2PNALWI/aDAmQAnmjGdCtHldXL86fv1xKu4HT+fP11bm1UVn9HedC2S1wbgYqOufUH0/CWnRkBIO61SNvdCT9J6/mWFo6f2xTVY3jhbsFqKQ2Q7zQoW4pGpYvxBOfLeJPXyxlyoqdvHhTbYolxHodTQIou5eolgK6a0tyhcgI4+Wb6nBb47IM/2kdL01cpQHFLpzaDPFUiQKxjLyzEc91qMHMtbtpM2A6k5b95nUsCaDsnsFJBFaa2VzgxO8znXM3+CWViMciIoznO9YCYPhP6zCDJ1vrTM4F+BewwMyWojZDPBIRYfRqXoHmlYvQ75OF3DtqHjc3SOLZG2qSL092//cnoSq7/4Wf9WsKkSD0e5HjgGFT12HAEypysmsk8BKwBA0tIR6rUjyBLx5ozsApqxk2dR2zN+zhtVuSaVi+kNfRxI/OWeD4RiR15+oS/vs6OR9NxHsREcY/OtbCOXh9asaZnD9cpyInG4465wZ5HULkdzFRETzRuhpXVi3Go2MW0vWNn7nv8kr0u+YyYqKye7eGhJLz/Vf90cz6nvlATTOLMbOrzGwk0NN/8US8FxFhvHBjLbo3KsvQH9fx7+90T042TDezf2loCQk2KeULMeGRVtzcoAyvT11Hp9dnsmbHIa9jiR+c7xJVG6A3MNrMKgD7yXiqeATwHTDAObfAvxFFvPd7kQOOoT+uwzAev+4ynck5u3q+n00yzVM3cQkK+fJE8VKXOlxVvRhPj11C+8EzeKpNNXo1K09EhP6mw8U5Cxzn3HHgdeB1M4sGigDHnHP7AxFOJJhkFDm1cQ6G/LiWiAjjsWsv8zpWUHLO6cGaEvRa1yxB/bIF+ePni/n7N8v5YeVOXrm5DiUL5PU6muSAbF94dM6lOee2q7iR3Cwiwvhnp9rckpLEoClrGP6TRjzOzMxuN7OztitmVsnMWgQyk8i5FE3Iw9s9U/hnp9rM+3UfrftP04N3w4T6yYlcoIgI41+d63As7TQvTlhJfEwkPZqW9zpWsChMRvfwecA8YBcQC1QGLgd2A095F0/kf5kZtzYuS9NKhXn0k4X0Hb2AySt28PeOtSiQN9rreHKRVOCIXITICOO1W+py7GQ6fxm3jLwxUXoKOeCcG2hmQ8i416Y5UAc4BqwAejjnNnmZT+RcKhSJ57P7mvL61HUMnLKGORv28urNdWlWuYjX0eQiZKvAMbMazrnlZ8y7wjk31S+pREJAdGQEQ26tx90jU3nys0XExUTSrnZJr2N5zjmXDnzvm/zOzD4BqvpeJgL7nXPJWay3ETgEpAOn9OBPyUpUZAQPX12Fyy8ryqOfLOTWt37hrhYVeKJ1VWKjI72OJxcgu/fgjDGzP1qGvGY2mIyRSkVytdjoSEbc0YD6ZQvy8OgF/LByh9eRch3nXFfnXLKvqPkcGHuO1a/0raviRs6pbplExj/ckh5NyvH2jA3cMGQGy7cd9DqWXIDsFjiNgTLALGAusI2M088iuV5cTBTv3NmQ6iXzc98H85m1drfXkXIly+izfwsw2ussEh7yxkTy/I21ePfOhuw7mkbHoTMY/tM60k9rHKxQkN0CJ42M6+h5ybhhcINz7pKHXzezx83MmZkucEpIyx8bzfu9G1GhcDx3v5/K/E37vI6UG7UEdjjn1pxluQMmm9k8M+sTwFwS4q6sWoxJ/VpxdbXivDhhJd3fnM3mvUe9jiXnkd0CZy4ZBU5DMhqR7mb26aXs2MzKANcBuulQwkLB+BhG3d2IYgl5uPPdubl6dFQz+2tW0yW832QzW5rF1DHTat0599mbFr7LWG2BB82s1Tn218fMUs0sddeuXRcbW8JIofgYht1en3/fXJfl2w7SduB0Pp+3RaOaB7HsFjh3Oef+mmksnI7AV5e47/7Ak2R8qxIJC8USYhl1V2NioiK44505bN1/zOtIXjmSaUono6gof7Fv5py7xjlXK4tpHICZRQGdgU/O8R5bfT93Al8Ajc6x7gjnXIpzLqVo0aIXG1vCjJnRpUESEx5pSY2S+Xn800U88OF89h456XU0yUK2ChznXGoW80Zd7E5937q2OucWZWNdfZOSkFKmUBzv927E4ROnuOPtX3Jl4+ecezXT9AJwBVDRj7u8BljpnNuS1UIzizezhN9/J+Ps8VI/5pEwVqZQHKP7NOGpttWYvGIHrQdMY+qqnV7HkjP47RGq5zml/AyQrdPV+iYloah6yfy83bMhW/Ydo/d7czly4pTXkbwWB/hzoKBunHF5ysxKmdm3vpfFgRlmtgiYA4x3zk30Yx4Jc5ERxn2XV+LLB5tTMC6aXu/O5S9fLuXYyXSvo4mPBfr6oZnVBqYAv9+hlURGr6xGzrnfzrVtSkqKS039n5NJIkHr++U7uHdUKi2qFOWtO1KIifLbdwpPmNm8rLpcm9kS/nP5ORIoCvzdOTckkPlygtodOZ/jaen8e9Iq3pqxgYpF4unfNZm6ZRK9jhW2ztbunCngra1zbolzrphzrrxzrjywBah/vuJGJBRdW6M4L3auw7TVu3jis0Wczj3dS68HOvim64BSoVjciGRHbHQkf76+Bh/d3Zhjael0HjaLgZPXcCr9kjsbyyUIr6+TIkHoloZleLJNVcYt3Mbfv1meK3pdOOd+zTRtdc7l+mt0Ev6aVS7CxEdacX2dkvSfvJouw39mw+4jXsfKtTwvcHxncjQymoS1+y+vRO/mFXhv1kaG6QnkImGrQFw0A7vVY1D3eqzfdZh2A6fz0S+bcsUXm2DjeYEjkhuYGX9uX52OyaV4eeIqxi0atHwsAAATJElEQVTc6nUkEfGjG+qWYtKjrahfLpFnvljC3SNT2XXohNexchUVOCIBEhFhvNylDo0rFOKJTxcze/0eryOJiB+VLJCXUb0b89frazB97W7aDJjGd8t0u2mgqMARCaA8UZGM6JFC2cJx9Hk/NVePdiySG0REGL1bVGB83xYUzx9Ln1Hz+ONnizmsoSP8TgWOSIAViIvm3V4NiYmKpNe7c9l58LjXkUTEz6oUT+DLB5vzwBWVGDNvM+0GTmfer3u9jhXWVOCIeKBMoTje7dWQvUdO0nukBgIUyQ1ioiJ4sk01xtzblNPOcfPwn3ll0kpOnlJ3cn9QgSPikdpJBRh6Wz2WbztI39ELNGaGSC7RsHwhJjzSki4Nkhj64zo6D5vJ2p26XJ3TVOCIeOiqasV5/sZa/LByJ89+tUxdSUVyiYTYaF7uUpfhtzdg2/7jtB80g3dnbshNg4H6nQocEY/d1rgc911eiQ9/2cQb09Z7HUdEAqhNrRJM7NeSZpUK87evl9Pz3Tn8dkD35eUEFTgiQeDJ1lW5vk5JXpq4kknqRiqSqxRLiOWdXg35x421SN24j9YDpvHN4m1exwp5KnBEgkBEhPHvm+tSJymRfh8vZOnWA15HEpEAMjNub1KO8Q+3oHyReB76aAGPfrKQA8fSvI4WslTgiASJ2OhI3uzRgMS4aO55P1Xdx0VyoYpF8/H5fU3pd00Vvlq0jbYDpvHzOg0KejFU4IgEkWL5Y3mrZwoHjqVxz/upHE9L9zqSiARYVGQE/a65jM/vb0ae6EhufWs2L4xfrvbgAqnAEQkyNUsVYEDXZBZvPcDjny5SrwqRXCq5TCLjH27BbY3L8ub0Ddw4dCYrth/0OlbIUIEjEoSuq1mCP7apxvjF2xkwZY3XcUTEI3ExUfzjxtq826shuw+fpOOQmbzx0zrS9cXnvFTgiASpe1tV5OYGSQyaskZPHxfJ5a6sVoxJ/VpyZbWi/GvCSm59czZb9h31OlZQU4EjEqTMjBc61aZR+UI88dliFmza53UkEfFQ4Xx5GH57A17pUodl2w7SdsB0xs7fogFCz0IFjkgQi4mKYHiPBpTIH8s978/TAGAiuZyZcXNKGSY80pJqJRN4bMwiHvxoPvuOnPQ6WtBRgSMS5ArFx/B2zxSOnTzFvR/My7U9KczsZjNbZmanzSzljGVPm9laM1tlZq3Psn0hM/vezNb4fhYMTHKRnFemUBwf92nKk22q8v3yHbQeMI2fVu/yOlZQUYEjEgKqFE/gta7JLNq8nz9/uTS3npJeCnQGpmWeaWY1gG5ATaAN8LqZRWax/VPAFOdcFWCK77VIyIqMMB64ojJfPNCcAnmj6fnOHJ4dt5RjJ3Pnl6AzqcARCRGta5bgkaur8Nm8LYyctdHrOAHnnFvhnFuVxaKOwMfOuRPOuQ3AWqDRWdYb6ft9JHCjf5KKBFat0gX4um8LejevwMiff+X6wdNZvGW/17E8pwJHJIQ8cnUVrq1RnOfHr9Dopv9RGtic6fUW37wzFXfObff9/htQ/GxvaGZ9zCzVzFJ37dJpfwl+sdGR/LVDDT68uzFHTqTT+fVZDJ6yhlPpp72O5hkVOCIhJCLCeO2WulQoEs+DH80Pu26iZjbZzJZmMXXMyf24jGt8Z73O55wb4ZxLcc6lFC1aNCd3LeJXzSsXYVK/VrStXZJXv1/NLW/8zK97jngdyxMqcERCTEJsNCN6NCAt/TR93p8XVtfbnXPXOOdqZTGNO8dmW4EymV4n+eadaYeZlQTw/dyZc8lFgkeBuGgGd6/HwG7JrNl5mLYDpzN6zqZcd++eChyREFSxaD4GdavHit8O8uTni3Ndw3WGr4BuZpbHzCoAVYA5Z1mvp+/3nsC5iiaRkNcxuTST+rUiuUwiT49dwj3vz2P34RNexwoYFTgiIerKasV4onVVvl60jRHT1nsdx+/MrJOZbQGaAuPNbBKAc24ZMAZYDkwEHnTOpfu2eStTl/IXgWvNbA1wje+1SFgrlZiXD+5qzF+ur8G0Nbto3X8a3y/f4XWsgLBQ+uaXkpLiUlNTvY4hEjScczw0egETlmzn/d6NaVGlSED3b2bznHMp518zdKndkXCx6rdD9PtkISu2H6RbwzL85foaxOeJ8jrWBctuu+PZGRwz62tmK30Dd73sVQ6RUGZmvNKlDpWL5ePhjxewbf8xryOJSJCqWiKBLx9sxn2XV+KT1M20HTideb+G7yNgPClwzOxKMsakqOucqwn824scIuEgLiaKYbc34OSp0zzw4XxOnAqfm45FJGfliYrkqbbV+KRPU9JPO24ePotXv1tFWhh2J/fqDM79wIvOuRMAzjn1ZhC5BJWK5uOVLnVYuHk///hmhddxRCTINapQiIn9WtK5fhKDf1hL59dnsXbnYa9j5SivCpzLgJZm9ouZ/WRmDc+2ogbcEsmetrVL0qdVRUbN/pUvFmzxOo6IBLmE2Gj+fXNdht9eny37jtJ+0HRGztoYNr0y/VbgnGfAriigENAEeAIYY2aW1ftowC2R7HuydVUaVSjE02OXsPK3g17HEZEQ0KZWSSb1a0XTSoV59qtl3PHOHHYcPO51rEvmtwLnPAN2bQHGugxzgNNAYLt/iIShqMgIhtxaj/yx0dw3ah4Hj6d5HUlEQkCx/LG826shz99Yi7kb99J6wDS+XbL9/BsGMa8uUX0JXAlgZpcBMcBuj7KIhJViCbEMva0+W/Yd4w9jFoXN6WYR8S8zo0eTcox/uCXlCsXxwIfzeeyThSH7RcmrAucdoKKZLQU+Bno6tcIiOaZh+UI83a463y3fwRu5YBBAEck5lYrm47P7m/HI1VUYt2gbbQdMZ/b60Hu4rycFjnPupHPudt8lq/rOuR+8yCESzno3L0/7OiV5eeJKZq3TCVIRyb7oyAgevfYyPruvKdGRRvc3Z/PPb1eE1DAUelSDSJgyM166qQ7li8TzyMcL2XUo9zyDRkRyRr2yBfn2kZZ0b1SWEdPW03HIzJDpwKACRySM5csTxeu31efgsTQeG7OQ06d1JVhELkxcTBT/7FSbd3qlsPvwCW4YPJM3p60P+vZEBY5ImKtWIj9/u6Em09fs5vWpa72OIyIh6qpqxZnUrxVXVC3KC9+u4Na3ZrM1iB8PowJHJBfo2rAMHZNL8dr3q0PyZkERCQ6F8+XhjR4NeLlLHZZsOUCb/tP4YsGWoOytqQJHJBcwM17oVJvyheN55OMF7Dms+3FE5OKYGbeklGHCI62oWiKBRz9ZxEMfLWD/0ZNeR/svKnBEcol8eaIYcmt99h1N49Exi4L++rmIBLeyheP45N6mPNG6KpOW/UbrAdOYviZ4HqmkAkckF6lRKj/PdqjBtNW7GPbTOq/jiEiIi4wwHryyMl8+2JyE2Gh6vD2H575axvE077uTq8ARyWVubVSW6+uU5LXvVzN3416v44hIGKhVugDf9G3Bnc3L896sjbQfNJ0lWw54mkkFjkguY2b8q3NtyhTMS9+PFrD3SHBdNxeR0BQbHcmzHWoy6q5GHD5xik6vz2TID2s4lX7akzwqcERyoYTYaIbcWp+9R05qfBwRyVEtqxRlUr9WtKlVgn9/t5quI2bz654jAc+hAkckl6pVugB/vr46U1ft4t1ZG72Oc15mdrOZLTOz02aWkmn+tWY2z8yW+H5edZbtnzOzrWa20De1C1x6kdwlMS6Gwd3rMbBbMqt3HKLdwOl8MndTQLuTq8ARycV6NCnHNdWL8+KEFSzd6u318mxYCnQGpp0xfzfQwTlXG+gJjDrHe/R3ziX7pm/9lFNEyLgc3jG5NBP7taJOUiJ//HwJfUbNY3eAhqlQgSOSi5kZr3SpQ+H4PDw8egFHTpzyOtJZOedWOOdWZTF/gXNum+/lMiCvmeUJbDoROZvSiXn58O7G/Ll9dX5avYs2A6YxefkOv+9XBY5ILlcwPobXutZlw54j/O3rZV7HuVQ3AfOdc2f7itjXzBab2TtmVjCQwURys4gI4+6WFfn6oRYUyZeHu99P5emxi/36pUoFjojQrFIRHryiMmNSt/D1om3n38BPzGyymS3NYuqYjW1rAi8B955llWFARSAZ2A68eo736mNmqWaWumtX8AxcJhLqqpZIYNxDzbn38op8PHcz7QdNZ/6mfX7ZlwocEQHgkWuqUL9sIs+MXcLmvUc9yeCcu8Y5VyuLady5tjOzJOAL4A7nXJYjGDrndjjn0p1zp4E3gUbnyDHCOZfinEspWrTopfyTROQMeaIiebptdUbf04S0dEeXYbN47btVpOVwd3IVOCICQHRkBAO71QMIqUtVZpYIjAeecs7NPMd6JTO97ETGTcsi4pEmFQszoV9LOtVLYtAPaxk+NWdHV4/K0XcTkZBWplAcb/RoQOVi+byO8j/MrBMwGCgKjDezhc651sBDQGXgr2b2V9/q1znndprZW8Bw51wq8LKZJQMO2MjZL2WJSIDkj43m1Vvq0qZWCZpWKpyj723B+Ijzs0lJSXGpqalexxARHzOb55xLOf+aoUvtjkhwyW67o0tUIiIiEnZU4IiIiEjYUYEjIiIiYUcFjoiIiIQdFTgiIiISdkKqF5WZ7QJ+zcaqRch4AF+wCvZ8EPwZle/S5FS+cs65sB4JT+1OwCjfpQv2jAFtd0KqwMkuM0sN5q6rwZ4Pgj+j8l2aYM8XioL9M1W+SxPs+SD4MwY6ny5RiYiISNhRgSMiIiJhJ1wLnBFeBziPYM8HwZ9R+S5NsOcLRcH+mSrfpQn2fBD8GQOaLyzvwREREZHcLVzP4IiIiEgupgJHREREwk7IFThm1sbMVpnZWjN7KovlZmaDfMsXm1n97G4boHy3+XItMbNZZlY307KNvvkLzcwvjy/ORr4rzOyAL8NCM/trdrcNUL4nMmVbambpZlbItywQn987ZrbTzJaeZbnXx9/58nl6/IUqtTt+z6d25+zZgrrNyWZGb44/51zITEAksA6oCMQAi4AaZ6zTDpgAGNAE+CW72wYoXzOgoO/3tr/n873eCBTx+PO7AvjmYrYNRL4z1u8A/BCoz8+3j1ZAfWDpWZZ7dvxlM59nx1+oTmp3ApJP7c7Z9xfUbU42M3py/IXaGZxGwFrn3Hrn3EngY6DjGet0BN53GWYDiWZWMpvb+j2fc26Wc26f7+VsICmHM1xSPj9t66983YHROZzhnJxz04C951jFy+PvvPk8Pv5CldodP+fz07b+yhfQdifY25zsZPTq+Au1Aqc0sDnT6y2+edlZJzvbBiJfZneRUXn/zgGTzWyemfXJ4WwXkq+Z73TiBDOreYHbBiIfZhYHtAE+zzTb359fdnh5/F2oQB9/oUrtTmDyqd25OKHU5kAAj7+onHwzyT4zu5KM/9AtMs1u4ZzbambFgO/NbKWvMg6k+UBZ59xhM2sHfAlUCXCG7OgAzHTOZf7WEAyfX0gI4uNP/CiI/7ur3ckFAn38hdoZnK1AmUyvk3zzsrNOdrYNRD7MrA7wFtDRObfn9/nOua2+nzuBL8g4xRjQfM65g865w77fvwWizaxIdrYNRL5MunHGaeIAfH7Z4eXxly0eHn+hSu2On/Op3bkkQd/mgEfHnz9u7PHXRMYZp/VABf5z01TNM9Zpz3/fcDUnu9sGKF9ZYC3Q7Iz58UBCpt9nAW08yFeC/wwA2QjY5Pssg+Lz861XgIzrvfGB/Pwy7as8Z7+ZzrPjL5v5PDv+QnVSuxOQfGp3zp0xqNucbGT05Pjzyz/UnxMZd4yvJuPu8D/55t0H3Of73YChvuVLgJRzbetBvreAfcBC35Tqm1/RdwAuApZ5mO8h3/4XkXEzWLNzbRvofL7XvYCPz9guUJ/faGA7kEbGNe27guz4O18+T4+/UJ3U7vg9n9qds2cL6jYnmxk9Of70qAYREREJO6F2D46IiIjIeanAERERkbCjAkdERETCjgocERERCTsqcERERCTsqMCRHGVmiWb2wDmW5zWzn8ws8hzrTDazgv5JKCLhRG2OnI0KHMlpicBZGxugNzDWOZd+jnVGnec9RER+pzZHsqQCR3Lai0AlM1toZq9ksfw2YByAmZU0s2m+dZeaWUvfOl+R8cReEZHzUZsjWdJAf5KjzKw88I1zrlYWy2KATc65Er7XjwOxzrkXfKeP45xzh3zL1gBNXKZnloiInEltjpyNniYugVQE2J/p9VzgHTOLBr50zi3MtGwnUApQYyMiF0ttTi6mS1QSSMeA2N9fOOemAa3IeMLte2Z2R6Z1Y33ri4hcLLU5uZgKHMlph4CErBY45/YBkWYWC2Bm5YAdzrk3yXgYW33ffCPj6cIbAxFYREKa2hzJkgocyVG+69czfTfwZXXD33dAC9/vVwCLzGwB0BUY6JvfAJjtnDvl77wiEtrU5sjZ6CZjCSgzqw886pzrcY51BgJfOeemBC6ZiIQjtTm5l87gSEA55+YDP55r0C1gqRoaEckJanNyL53BERERkbCjMzgiIiISdlTgiIiISNhRgSMiIiJhRwWOiIiIhB0VOCIiIhJ2/g8QAGd4lgxsYAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3ef404bac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 数値解\n",
    "x_list, u_list, t_list = solve_freefall_vacuum(1.0, 0.0, 0.001, 1.3)\n",
    "\n",
    "plt.figure(figsize=(8,3))\n",
    "# t vs x のグラフ\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(t_list, x_list)\n",
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('x (m)')\n",
    "# t vs u のグラフ\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(t_list, u_list)\n",
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('u (m/s)')\n",
    "# 2つのグラフを綺麗に配置する。\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 空気抵抗下での自由落下のシミュレーション\n",
    "\n",
    "上の例では空気抵抗の影響を無視したが、実際の系では多少なりとも空気抵抗が影響を及ぼす。\n",
    "空気抵抗の大きさは物体の速度の大きさの二乗に比例し、その向きは物体の運動方向と逆向きとなる。\n",
    "具体的には、空気抵抗力の大きさ $D$ は物体が球形の場合、\n",
    "\n",
    "$$\n",
    "D=\\frac{1}{2}\\rho |\\mathbf{u}|^2 \\pi a^2 C_D,\n",
    "$$\n",
    "\n",
    "と表せる。\n",
    "なお、 $|\\mathbf{u}|$ は速度ベクトル $\\mathbf{u}$ の絶対値を表す。\n",
    "ここで、 $\\rho$ は空気の密度, $a$ は球の半径である。\n",
    "また、抗力係数 $C_D$ はある条件のもと（注）では約0.44で一定であることが知られている。\n",
    "\n",
    "解くべき微分方程式は、\n",
    "\n",
    "$$\n",
    "\\frac{dx}{dt}=u, \\ \\ \\ x(0)=h, \\\\[3pt]\n",
    "\\frac{du}{dt}=-g-\\frac{D}{m}\\frac{u}{|u|}, \\ \\ \\ u(0)=u_0.\n",
    "$$\n",
    "\n",
    "となる.\n",
    "\n",
    "空気抵抗の影響を含めた落下運動のシミュレーションを行うソースコードは以下のように書ける.\n",
    "計算結果を, 空気抵抗のない場合と比較してみよ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def next_u_air(u, dt, g, rho, a, m):\n",
    "    \"\"\"\n",
    "    現在の x, u の値から、dt 後のxを求める関数。\n",
    "    空気抵抗を考慮する。\n",
    "    args:\n",
    "    + u: 時刻 t での速度 [m/s]\n",
    "    + dt: 線形近似する時間幅 [s]\n",
    "    + g: 重力加速度の値 [m/s^2]\n",
    "    + rho: 空気の密度 [kg/m^3]\n",
    "    + a: 物体の半径 [m]\n",
    "    + m: 物体の質量 [kg]\n",
    "    return:\n",
    "     時刻t+dt での u の値\n",
    "    \"\"\"\n",
    "    # 効力係数\n",
    "    CD = 0.44 \n",
    "    # 空気抵抗の大きさ\n",
    "    D = 0.5 * rho * u**2.0 * np.pi * a**2.0 * CD\n",
    "    return -g*dt -D/m * np.sign(u) * dt + u \n",
    "    # np.sign(x) 関数はxの正の場合は1を負の場合は-1を返す関数である。\n",
    "    # https://docs.scipy.org/doc/numpy/reference/generated/numpy.sign.html\n",
    "    # を参照のこと。\n",
    "\n",
    "\n",
    "def solve_freefall_air(x0, u0, dt, tmax, a, m):\n",
    "    \"\"\"\n",
    "    t=0 に位置h、速度u0の物体の自由落下をシミュレートする関数。\n",
    "    args:\n",
    "    + x0: t=0での座標 [m]\n",
    "    + u0: t=0での速度 [m/s]\n",
    "    + dt: オイラー法を用いるための時間軸の刻み幅 [s]\n",
    "    + tmax: 追跡を行う最大時間 [s]\n",
    "    + rho: 空気の密度 [kg/m^3]\n",
    "    + a: 物体の半径 [m]\n",
    "    + m: 物体の質量 [kg]\n",
    "    \"\"\"\n",
    "    # 計算結果を代入するリスト\n",
    "    x_list = []\n",
    "    u_list = []\n",
    "    t_list = []\n",
    "    \n",
    "    g = 9.81 # 重力加速度 [m/s^2]\n",
    "    rho = 1.261 # kg/m^3\n",
    "    \n",
    "    # 初期値の代入\n",
    "    x = x0\n",
    "    u = u0\n",
    "    t = 0.0\n",
    "    \n",
    "    while(t < tmax):\n",
    "        x_list.append(x)\n",
    "        u_list.append(u)\n",
    "        t_list.append(t)\n",
    "        \n",
    "        # x, u, tの値の更新\n",
    "        x = next_x(x,u,dt)\n",
    "        u = next_u_air(u,dt,g,rho,a,m)\n",
    "        t = t + dt\n",
    "    # 以下のように、Pythonでは関数の戻り値を複数にすることができる。\n",
    "    return x_list, u_list, t_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAADQCAYAAAAK/RswAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VFX+x/H3d1IJpAAJvTchIDWEgFRFf+DSFMVGExCx67q6uq7uurvuWtYuCoiAINJsuAoWBEFASkINiPTekV5DOL8/EndZlhIgyZ2ZfF7PMw+ZmXvvfBgmh++ce+455pxDREREJJj4vA4gIiIikttU4IiIiEjQUYEjIiIiQUcFjoiIiAQdFTgiIiISdFTgiIiISNBRgSMiIiJBRwWOiIiIBB0VOCIiIhJ0Qr0OcDHi4+NdpUqVvI4hIueQlpa22zmX4HWOy6W2RsR/5bSdCagCp1KlSqSmpnodQ0TOwcw2eJ0hN6itEfFfOW1ndIpKREREgo4KHBEREQk6KnBEJOiZWTsz+9nMVpvZE17nEZG8pwJHRIKamYUAA4H2QCJwm5kleptKRPJaQA0yzomNe44w5Ic1hPp8hIUYoSE+wkJ8hPmMqIhQYiJDiSkURkxkGDGFQikaFU5CdARhIar1RIJUMrDaObcWwMzGAp2B5Zdz0CnLd7B0y34eblsdM8uFmCKSm4KuwNl9+DiTlm4nI/MUJzNd1p+n3Hn3MYPihSMoFRtBqZhISsVGUrFYYSrHF6ZSfGEqFIsiPFQFkEiAKgtsOu3+ZqDJmRuZWX+gP0CFChUueNDpK3cxas4GMjJP8dj/XaEiR8TPBF2B07BCURY8fe1/PeacIyPTcfREJgeOZbD/aAYHjmVw4OhJfjl8gh0HjrHjwDG2HzjG5r1Hmb9+L/uPZvx7f59BuaJRVC9RhNplYkgsE0Ni6VjKFyukRk0kSDjnhgBDAJKSks7/rQh4tlNtMp3j7e/XkHnK8UT7mmoPRPxI0BU4Z2NmhIca4aE+YqPCKJ+DffYdOcG63YdZt/sw63cfZs3uw6zcfpBpP+/k1w6h6IhQEsvE0KhiUZIqFaVhhaLERYXn6d9FRC7aFvivX/ty2Y9dFp/PeK5LHUJ9xuAZa8nIdDzdoZaKHBE/USAKnEsRFxVOgwrhNKhQ9L8eP5aRyc/bD7Js6wGWb9vP0s37GTJjLW9/n1X1VC9RhEYVi9K0anGaV4uneJEIL+KLyH/MB6qbWWWyCptbgdtz48BmxrOdahPiM4bNWkfmqVP8uVNtFTkifkAFzkWKDAuhXvk46pWP+/djR09ksnjzPtI27CV1/S9MWrqNsfOzTvnXLhNDi+oJtKgeT6OKRYkMC/EqukiB5Jw7aWb3A18DIcAw59yy3Dq+mfFMh0RCzBg6cx2ZzvGXTnXw+VTkiHhJBU4uKBQeQkqV4qRUKQ5A5ilH+pb9/LBqFz+s2s17M9cyaPoaosJDaFUjgetql+TqK0oSGxXmcXKRgsE5NwmYlFfHNzOe+k0tQkKMwdPXknnK8VyXK1XkiHhIBU4eCPHZv3t57r+6OoePn2Tuuj1MXbGTb5btYHL6dkJ9RkqV4lxXuyTtapeiREyk17FF5DKYGU+0q0mozxg4bQ0nMx3Pd61LiIocEU+owMkHhSNCubpmSa6uWZK/dKrD4s37+Gb5Dr5Ztp1nJi7jz58vo1nVeDrXL0O7OqWIjlTPjkggMjN+d90VhPh8vPHdKjKd46Wb6qnIEfGACpx85vMZDSoUpUGFovy+XU1W7zzI54u28tmirTz20RL++Fk6bWuVpHP9MrSpWUITEIoEGDPjt9fWINRnvPLtSjJPOV6+uR6h+l0WyVcqcDxWrUQ0v73uCh65tgYLN+1j4sItfLFkG18u3UZCdAQ3NyrHrY0rUKF4lNdRReQiPHhNdUJ8xktf/0zmKcdrt9RXkSOSj1Tg+Akzo2GFrLl0/tghkRkrdzFm3iYGTV/D29+voXm1eG5LrsC1iSU1q7JIgLivTTVCfcY/Jq8g85TjjdsaqFdWJJ+owPFDYSE+rqlVkmtqlWTb/qNMSN3MuPmbuO/DBcQXiaB7SgW6p1QkXnPsiPi9u1tVJcRn/O3Ln7hv9ALeur2hvqSI5AP9lvm50rGFePCa6sx4vA3D72zMlWVjeG3KKpr9YyqPTVjM8q0HvI4oIhfQr0UV/twxkW+W7+De0WkcP5npdSSRoKcenAAR4jPaXFGCNleUYM2uQ4yYtZ6P0jYzIW0zTasU566WlWlzRQnNoCrip3pfVZkQn/H0xGUMGJXGO90baeJPkTykHpwAVDWhCH/tUoc5T17DE+1rsn7PYfqMSOX6N2byxZKtZF5g9XQR8UaPppX4+w1XMu3nXfQflcaxDPXkiOQVFTgBLDYqjAGtqjLj8Tb88+Z6HD+Zyf0fLuTaV6YzPnUTGZmnvI4oIme4vUkFXuxalx9W7aLf+6kcPaEiRyQvqMAJAmEhPm5qVI5vH2nFwNsbEhkWwuMfLaH1S98zas4GTpxUoSPiT7o1Ls9LN9Vj1prd9BkxnyMnTnodSSToqMAJIiE+4zd1S/Plg80Z3rsxJWMiePqzdK555Xs+TtusU1cifuSmRuV4pVs95q7bQ+/h8zl0XEWOSG7ytMAxs3Zm9rOZrTazJ7zMEkzMjDY1S/DxPc0Y3rsxMZFhPDphMf/32gwmL92Gcyp0RPzBDQ3K8dqtDUjbsJfew+Zx8FiG15FEgoZnBY6ZhQADgfZAInCbmSV6lScY/Vro/Ov+5rx9R0Occ9wzegEd35rJjJW7vI4nIkCnemV449YGLNq0j57D5nFARY5IrvCyBycZWO2cW+ucOwGMBTp7mCdo+XzG9VeW5ptHWvHPm+ux70gGPYfNo/fweazacdDreCIF3m/qluat2xuydPN+egydy/6jKnJELpeXBU5ZYNNp9zdnP/ZfzKy/maWaWequXep1uBwhPuOmRuX47tFWPHV9LdI27KXd6z/w1KdL2X3ouNfxRAq0dnVK8U73RizfdoDuQ+ey78gJryOJBDS/H2TsnBvinEtyziUlJCR4HScoRISGcFfLKkx/rA3dm1Rg7PxNtH7pe975fo3m5RDx0LWJJRncoxE/bz/I7e/OZe9hFTkil8rLAmcLUP60++WyH5N8UqxwOM92rsPXD7ckpUoxXvhqBde+Op0py3d4HU2kwLq6ZkmG9GzE6l2HuO3dOexR76rIJfGywJkPVDezymYWDtwKfO5hngKrWokiDO3VmNH9mhAZGkK/kan0HTGfjXuOeB1NpEBqfUUJ3uuVxLrdh7llyBx2HjjmdSSRgONZgeOcOwncD3wN/ASMd84t8yqPwFXV4pn0UAv+cH1N5qzdw7WvTuf1Kat02krEAy2qJzDizmS27jtKt8E/smXfUa8jiQQUT8fgOOcmOedqOOeqOuee8zKLZAkL8dG/ZVW+e7Q11yaW5NUpK7nu1RlMW7HT62giBU7TqsUZ1bcJew6foNugH1m/+7DXkUQCht8PMhZvlIqN5K3bGzK6XxPCQow7R8znvg8XsOugxgOI5KdGFYsy5q4Ujpw4SbfBP7J6p6Z2EMkJFThyXldVi2fyQy159NoafLtsB21fmc6E1E2aDVkkH9UpG8vY/k055eCWwXNYvvWA15FE/J4KHLmg8FAfD1xTnUkPtaBGySI89tESerw3T4OQRfLRFaWiGX93CuGhPm4d8iOLNu3zOpKIX1OBIzlWrUQRxvVvyl+71GHRpn1c99p03p2xlpOZWq1c/JOZvWRmK8xsiZl9amZxXme6HFUSijD+7qbERoXRfehc5q37xetIIn5LBY5cFJ/P6JFSkW8eaclVVeN5btJPdB30I2t2HfI6msjZfAvUcc7VBVYCT3qc57KVLxbFhLubUSImgp7D5jJz1W6vI4n4JRU4cknKxBViaK8kXr+1Phv2HOb6139g2Mx1nDqlsTniP5xz32RPSQEwh6wJRQNeqdhIxvVvSqXihenz/ny++0mTc4qcSQWOXDIzo3P9snzzcEuaVS3OX75Yzh1D57J5r8bmiF/qA0w+15OBtu5dQnQEY/unULNUNHePSuPLJdu8jiTiV1TgyGUrERPJsN6NeaHrlSzZvI92r/3A+Pm60kryh5lNMbP0s9w6n7bNU8BJYPS5jhOI697FRYXzQb8m1C8fxwNjFvDJgs1eRxLxG6FeB5DgYGbc0rgCzarG89hHi3n84yV8vWw7z3etS0J0hNfxJIg559qe73kz6w10AK5xQVh1x0SGMbJvMv3eT+XRCYs5lnGK25tU8DqWiOfUgyO5qnyxKD7sl8IzHRKZuXo37V+fwfSV/t/dL8HJzNoBjwOdnHNBe+40KjyUYb0b07pGAn/4dClDf1jrdSQRz6nAkVzn8xl9mlfm8/ubU7xwBL2GzeNvXyzn+EmtaSX57i0gGvjWzBaZ2SCvA+WVyLAQBvdI4vorS/G3L3/ilW9+1mliKdB0ikryzBWlopl4/1X8Y9JPDJ25jh/X7uGN2xpQNaGI19GkgHDOVfM6Q34KD/Xxxq0NKBKxlDemrmb/0Qz+1LE2Pp95HU0k36kHR/JUZFgIz3auw7s9k9i67ygd3pjJuPkb9c1SJI+Ehvh4oWtd7mpRmfd/3MCjExaTock4pQBSgSP54trEkkx+qCUNKsTx+4+Xcv+HCzlwLMPrWCJBycz4w/W1+N11Nfh04Rbu+SCNYxk6RSwFiwocyTelYiP5oG8Tft+uJl8v207HN2eybOt+r2OJBCUz4/6rq/PXzrWZ8tNOeg+fx0F9qZACRAWO5Cufz7indVXG9k/heMYpbnh7NmPm6ZSVSF7p0bQSr91Sn/nr93LH0Ln8cviE15FE8oUKHPFEUqVifPlgc5pULsaTnyzl0fGLOXLi5IV3FJGL1qVBWYb0aMTP2w/SbfCPbNt/1OtIInlOBY54pniRCEbcmczDbavz6aItdBk4i9U7D3odSyQoXVOrJO/3SWb7/mPc9M6PrN992OtIInlKBY54KsRnPNy2BiP7JLP70Ak6vTWLiYu2eB1LJCilVCnOmLtSOJqRyU2DfmT51gNeRxLJMypwxC+0qJ7ApAdbkFg6hofGLuLZfy3Tpa0ieeDKcrGMv7spYSHGrUN+JG3DL15HEskTKnDEb5SKjWRM/xR6N6vE8Fnr6fHeXPYcOu51LJGgU61EESYMaErxIhHcMXQu3/20w+tIIrlOBY74lbAQH3/uVJuXb67Hgo376PTWLNK36FJykdxWrmgUEwY0pXqJaPqPSmN86iavI4nkKhU44pe6NirHxwOa4Zyj6zuz+XThZq8jiQSd+CIRjOmfQrOqxXn8oyUMnLZaUzZI0FCBI37rynKxfP5Ac+qXj+ORcYv5y7+Wc1LjckRyVZGIUN7r1ZjO9cvw0tc/8+y/lnPqlIocCXwqcMSvxReJ4IN+TejdrBLDZq2ju8bliOS68FAfr3arT9/mlRkxez0PjF3I8ZNa2kECmwoc8Xtnjsvp8vYsVu7QfDkiucnnM57ukMgfrq/Jl0u2cefw+VraQQKaChwJGF0blWNc/xSOZZzixrdnM+3nnV5HEgk6/VtW5ZVu9Zi37hduHTKHnQePeR1J5JKowJGA0qBCUSbedxUVikXRd8R8hs1cp0GRIrnsxoblGNoribW7DtP1ndma9VgCkgocCThl4grx0T1NaVurJH/5Yjl/+DRdkwKK5LLWV5RgTP8UDh/PpOs7s1myeZ/XkUQuigocCUhR4aEM6t6Ie1tXZcy8jfR8bx77jmiVZJHcVL98HB8NaEpkWAi3DpnDtBU6LSyBQwWOBCyfz3i8XU1evrkeaRv2csPbs1mz65DXseQSmFmkmd1kZq+b2QQzG2lmj5tZba+zFXRVEorw6b3NqJJQmL7vz+eDORu8jiSSIypwJOB1bVSOD+9qwoGjGdwwcBZz1+7xOpJcBDN7FpgFNAXmAoOB8cBJ4Hkz+9bM6noYscArERPJuP5NaVUjgT9+ls4/Jv+kuXLE76nAkaCQVKkYn913FQnREfR4bx6fL97qdSTJuXnOuUbOuUedcx8656Y4575wzr3inOsI3AGEex2yoCscEcq7PZO4o0kFBk9fy4NjF3IsQ3PliP9SgSNBo3yxKD6+pxn1y8fx4JiFDJq+RldYBQDn3JdnPmZmPjOLyX5+p3MuNf+TyZlCQ3z8rUsdnmhfky+WbKPHe3PZe1hj38Q/eVLgmNlLZrbCzJaY2admFudFDgk+cVHhjOybTMd6ZXh+8gqenpiu5R0ChJl9aGYxZlYYSAeWm9ljXueS/2ZmDGhVlbdub8Dizfvp+s5sNuzRZeTif3JU4JhZkpk9kl2Y/MXMuplZ0ct43W+BOs65usBK4MnLOJbIf4kMC+H1W+ozoFVVPpizkbtHpXHkxEmvY8mFJTrnDgBdgMlAZaCHt5HkXDrULcPofk345cgJbnx7Ngs37vU6ksh/OW+BY2Z3mtkCsgqQQsDPwE6gOTDFzN43swoX+6LOuW+cc7/+jzMHKHexxxA5H5/PeKJ9Tf7auTbTft7JrUPmsOug1rDyc2FmFkZWgfO5cy4DyJVzjGb2qJk5M4vPjeNJlsaVivHJPc0oHBHKbe/O4av07V5HEvm3C/XgRAFXOee6Ouf+7pwb6px7yzn3oHOuEfAqUP0yM/Qh69vaWZlZfzNLNbPUXbt2XeZLSUHTo2klhvRIYtWOQ9zw9ixW79Rl5H5sMLAeKAzMMLOKwIHLPaiZlQeuAzZe7rHkf/16GXmt0jHcMzqNd2es1dg38QvnLXCccwOdc0fP8/wi59x3Z3vOzKaYWfpZbp1P2+Ypsi4FHX2e1xjinEtyziUlJCRc+G8kcoa2iSUZ2z+FYxmZ3DRIXen+xsyampk5595wzpV1zl3vsv6H3Ai0yYWXeBV4nFzqDZL/VbxIBGPuSqF9nVI8N+knnvh4KSdOauybeCunY3Aqm9krZvaJmX3+6+18+zjn2jrn6pzlNjH7mL2BDsAdTuW+5LF65eP45J6riC0Uxu3vzuV7LdTpT3oCaWY21sx6m1kpAJflsgZPZX+h2uKcW5yDbdVbfBkiw0J467aGPHh1NcalbtIVVuI5y0ltYWaLgfeApcC/y3Ln3PRLelGzdsArQCvnXI5bkqSkJJeaqqtF5dLtOnicXsPmsXLHQV7uVo/O9ct6HSmomFmacy7pEvetCbQH/g+IBaYBXwGznHPnnHDFzKYApc7y1FPAH4DrnHP7zWw9kOSc232hLGprLs/ERVt47KMllI6N5L1eSVQrEe11JAkiOW1nclrgzHXONcmVZFnHWw1EAL9OOTvHOTfgQvup0ZHccOBYBv1HpjJn7S880yGRPs0rex0paFxOgXPGcQqRdXqqPdD0Uo5pZlcC3wFHsh8qB2wFkp1z5x0Nq7bm8i3YuJf+I9M4fjKTgbc3pGUNDTGQ3JHbBc7tZA0m/gb496UozrkFlxPyYqnRkdxyLCOTh8cu4qtl27mvTVV+d90VmJnXsQLe5RY42dNPlAdCf30st9oZ9eDkvy37jtJ3xHxW7TzEMx0S6dWskteRJAjktJ0JvdAG2a4kaz6Kq/nPKSqXfV8k4ESGhTDwjob88bN0Bk5bw55DJ/hblzqEhmhyb6+Y2V+B3sBa1M4EhbJxhfj4nmY8NHYRf/p8Gat3HuKZjomE6fdM8kFOC5ybgSrOOY0Yk6AR4jP+fkMd4ouE8+bU1fxy+ARv3NaAyLAQr6MVVN2AqnnVzjjnKuXFceX8CkeEMrhHI178agWDZ6xl3e7DDLy9IbFRYV5HkyCX0zI6HdByChJ0zIxHr7uCP3dM5JvlO+g1bB4Hj2V4HaugUjsTpEJ8xpPX1+LFm+oyd90eOg2cycodB72OJUEupwVOHLDCzL7O6WXiIoGk91WVef3W+qRt2Ev3oXPZd0SdlR74B7BQ7Uzw6pZUnrH9UzhyIpMuA2fxVfo2ryNJEMvpKao/5WkKET/QuX5ZCoeHcu+HC7hl8BxG9UumRHSk17EKkveBFzhjOgoJLo0qFuOLB5pz96g0BnywgAeursYjbWvg82mQv+SuC61FZZA1383ZbqdvIxIM2iaWZHjvxmz85Qi3DJ7Dln3nnMhbct+R7NmMp53ZzkhwKRkTybi7U+iWVI43p67mrpGpHNCpYcllFzpFNc3MHjhzQU0zCzezq83sfaBX3sUTyX9XVYvng37J7D50nG6DfmTd7sNeRyoofjCzf2Qv3dDw15vXoSRvRISG8ELXuvylc22mr9xFl4FaK05y14UKnHZAJjDGzLaa2XIzWwesAm4DXnPOjcjjjCL5rlHFYoy5K4WjGZncPOhHVmy/7DUf5cIaACnA34GXs2//9DSR5Ckzo2fTSozu14T9RzLoMnAWU5bv8DqWBIkcTfQHYGZhQDxw1Dm3L09TnYMm35L8tmrHQbq/N5fjJ0/x/p3J1Cuvi3zOJ7dmMvaa2pr8t3XfUe4elcbSLft5pG0NHri6msblyFnltJ3J8WxLzrkM59w2r4obES9ULxnNhLubUSQilDuGzmXu2j0X3kkuipl1N7NztkVmVtXMmudnJsl/ZeIKMWFAU25sUJZXp6ykz/vztVinXBZNJylyARWKR/HRgGaUjImg57B5zFiplaZzWXGyLg8fZmb3mVk3M+tpZn8xs+nAi4DOWxQAkWEhvNytHn/tUofZq/fQ4c2ZLNms79RyaVTgiORAqdhIxt3dlCoJReg3MpVpP+/0OlLQcM69DjQExgAJwDXZ97cAPZxzXZ1zqzyMKPnIzOiRUpEJA5oCcNM7P/LBnA3kdDiFyK9yVOCYWeJZHmud62lE/Fh8kQg+7NeE6iWKcPfINA2GzEXOuUzn3LfOuT875+52zj3snBvsnNvodTbxRr3ycXzxQHOaVi3OHz9L57fjF3PkxEmvY0kAyWkPzngz+71lKWRmb5I166hIgVK0cDgf9kuhZulo7hmdxtfLtnsdSSRoFS0czvDejfnttTX4bNEWugycxZpdupRccianBU4ToDwwG5gPbAWuyqtQIv4sNiqMUX2bULtMLPeNXsDkpZpuXiSv+HzGg9dUZ2SfZHYdPE7nt2YxSb9zkgM5LXAygKNAISASWOec01TqUmDFFgpjVN+sy8bvH7OQfy3e6nUkkaDWonoCXz7Yguoli3Dv6AU8MzGdYxmZXscSP5bTtajmAxOBxmTNhTPIzLo6527Os2Qifi46Moz3+yRz5/B5PDR2Iaeco3P9sl7HClhm9szZHnfO/SW/s4h/KhNXiHH9m/LiVysYOnMd89fv5a3bG1A1oYjX0cQP5bQHp69z7pnT5sLpDGiVXynwikSEMuLOZJIrF+ORcYv4OG2z15EC2eHTbplAe6CSl4HE/4SH+vhjh0SG9U5i+/6jdHxzJh/p907OIkcFjnPuf6b0dM6Nyv04IoGncEQow3sn07RqcX730WLGp27yOlJAcs69fNrtOaA1UMXjWOKnrq5ZkskPteTKsrH8bsJifjtuEYeO6yor+Q/NgyOSCwqFh/Ber8Y0rxbP4x8tUZGTO6KAcl6HEP9VKjaSD+9K4ZG2WVdZdXxzJulb9nsdS/yEChyRXBIZFsK7PZNoUT2e33+8hE8WqNv8YpjZUjNbkn1bBvwMvOZ1LvFvIT7jobbV+fCuFI6eyOTGt2czfNY6TQwoKnBEclNkWAhDeiTRtEpxfjdhMRMXbfE6UiDpAHTMvl0HlHHOveVtJAkUKVWKM+mhFrSoHs+z/1pOnxHz2XnwmNexxEMqcERyWaHwEIb2SqJxpayBx18u0ZwdOeGc23DabYtzTgMq5KIUKxzO0F5J/LljIrPX7KHdaz/wjSbjLLBU4IjkgajwUIb1bkyjikV5cOxCvkpXkSOSH8yM3ldV5osHmlM6NpL+o9L4/UdLOKwByAWOChyRPFI4IpThdyZTr1ws93+4kG+1dpVIvqleMppP772Ke1tXZXzaJq5/4wfSNuz1OpbkIxU4InmoSEQoI/okU7tsLPeOTmPqChU5+c3MHjCzFWa2zMxe9DqP5J/wUB+Pt6vJuP5NOZnpuHnQbF75diUZmZqIvyBQgSOSx2IiwxjZJ5mapWIYMGoB01fu8jpSgWFmbYDOQD3nXG3gnx5HEg8kVy7G5Idb0KVBWd74bhU3DfqR1Tu1aGewU4Ejkg9+XbuqWoki9B+ZysxVu72OVFDcAzzvnDsO4Jzb6XEe8UhMZBivdKvPwNsbsmHPYa5/4weGzFhD5ildTh6sVOCI5JO4qHA+6NeEyvGF6TdyPnPW7vE6UkFQA2hhZnPNbLqZNfY6kHjrN3VL880jLWlVI4G/T1rBzYNms2aXenOCkQockXxUrHA4o/s1oXzRKPqOmM/CjRr0eLnMbIqZpZ/l1pmsBYWLASnAY8B4M7NzHKe/maWaWequXTqNGMxKREcypEcjXr+1Pmt2Heb613/g3Rlr1ZsTZFTgiOSz4kUi+KBfE+KjI+g1bB7Ltx7wOlJAc861dc7VOcttIrAZ+MRlmQecAuLPcZwhzrkk51xSQkJCfv4VxANmRuf6Zfn2ty1pWSOB5yb9pN6cIKMCR8QDJWMi+aBvEwpHhNLjvbka8Jh3PgPaAJhZDSAc0AAo+bdfe3Neu+W/e3NO6kqrgKcCR8Qj5YtFMbpfE8yg+9C5bPrliNeRgtEwoIqZpQNjgV5OixTJGcyMLg3K8u0jLWlRPas354a3Z2vhzgDnaYFjZo+amTOzs3YZiwS7KglFGNW3CUczMrlj6Fy279faObnJOXfCOdc9+5RVQ+fcVK8zif8qERPJuz0b8dbtDdi2/xidB87iH5N+4uiJTK+jySXwrMAxs/JkLai30asMIv6gVukY3u+TzJ5Dx7lj6Bz2HDrudSSRAsvM6FC3DN/9thXdksoxeMZarnttuuavCkBe9uC8CjwOqLtYCrz65eMY1rsxW/Ydpcd789h/NMPrSCIFWmxUGP+4sS7j+qcQFuKj17B5PDx2Ibv1BSRgeFLgZF++ucU5t9iL1xfxR02qFGdwjyRW7TxI7+HztDigiB9oUqU4kx9qwUPXVOfLpdto+8p0xqdu4pQuKfd7eVbgXGBuij8Az+TwOJpgSxwgAAAQ00lEQVSbQgqMVjUSePO2hizZvJ9+76dyLEPn/kW8FhEawiPX1mDSgy2oXqIIj3+0hJsH/8iyrRqE7M/yrMA519wUwFqgMrDYzNYD5YAFZlbqHMfR3BRSoLSrU4p/3lyXOev2cO/oBVoYUMRPVC8Zzbj+TXnpprqs332Yjm/O5JmJ6ew/olPK/ijfT1E555Y650o45yo55yqRNRFXQ+fc9vzOIuKvbmhQjr91qcPUFTv53YTF6g4X8RM+n3FzUnmmPtqaHikV+WDOBq5++XudtvJDmgdHxE/d0aQij7e7gomLtvLsv5ah6VtE/EdsVBjPdq7Dvx5oTqX4wjz+0RJuGqS5c/yJ5wVOdk+OZhYVOYt7WlXlrhaVef/HDbzx3Wqv44jIGWqXiWXC3VmnrTbsOUKnt2byx8+WaroHPxDqdQAROTcz4w/X12LvkQxenbKSuKgwejWr5HUsETnNr6etrqtdile/XcmoORuYuGgrD1xdjV7NKhERGuJ1xALJ8x4cETk/M+P5G6/k2sSS/OnzZUxctMXrSCJyFrGFwvhzp9p8/XALkioW5e+TVnDdqzP4Kn27TjF7QAWOSAAIDfHx5m0NaFK5GI+OX8y0FTu9jiQi51CtRDTD70xmZJ9kIkJ9DPggjVuHzNH4nHymAkckQESGhTC0VxI1S0dzz+g0Utf/4nUkETmPljUSmPRgC/7WpQ6rdh6i41szeWzCYnYc0Jpz+UEFjkgAiY4MY8SdyZSJLUSfEfP5adsBryOJyHmEhvjonlKR7x9rTf8WVZi4aCutXprG85NXaP6cPKYCRyTAxBeJYGTfZKLCQ+k5bB4b9hz2OpKIXEBMZBhPXl+L7x5tRfs6pRk8Yw0tXpzKO9+v0WrleUQFjkgAKlc0ilF9k8nIPEWP9+axU13eIgGhfLEoXr2lPpMebEGjikV54asVtP7nNMbM28hJzVqeq1TgiASo6iWjGXFnMrsPHafnsHnq7hYJILVKxzD8zmTG9U+hbFwhnvxkKde9OoNJS7fpiqtcogJHJIDVLx/HkB5JrNl1iLtGanFOkUDTpEpxPr6nGe/2TCI0xLh39AJ+88ZMvl6mS8svlwockQDXvHo8L3erz7z1v/Dw2EVkaj0ckYBiZlybWJLJD7XklW71OHLiJHePSqPDmzP5RoXOJVOBIxIEOtUrw9MdEvlq2XatWyUSoEJ8xo0NyzHlt614+eZ6HD5+kv7Zhc63y3fo9/oiaakGkSDRt3lldh44xuAZaykZE8l9bap5HUlELkFoiI+ujcrRuX4ZPlu0lTenruKukanUKRvDw9fU4JpaJTAzr2P6PRU4IkHk9+1qsuPAMV76+mcSoiPollTe60gicolCQ3zc1KgcXeqX4dOFW3hz6mr6jUylZqlo7mldld9cWZrQEJ2IORe9MyJBxOczXrypHi2qx/PkJ0uZumKH15FE5DKFhvi4Oak83z2adeoq85TjobGLaPPy94yas0EXF5yDChyRIBMe6uOd7o2oVTqa+0YvZOHGvV5HEpFcEJZ96urrh1sypEcjiheO4OnP0mn+wjTe/n41B45pqojTqcARCUJFIkIZ3juZhOgI+oyYz9pdh7yOJCK5xOczrqtdik/vbcaYu1JILBPDi1/9zFX/mMrzk1ewfb8m/gQVOCJBKyE6gpF9kvGZ0XNYwZzt2Mzqm9kcM1tkZqlmlux1JpHcYmY0rVqckX2S+eKB5rS8IoHBM9bQ/IWpPDR2IUs27/M6oqdU4IgEsUrxhRl+Z2N+OXyCXsPnc7DgdWG/CDzrnKsPPJN9XyTo1Ckby8DbGzL9d23o2bQS3/20k05vzeLmQbP5Kn1bgZwfSwWOSJCrWy6Ot+9oyKodBxnwQRrHTxaoAYkOiMn+ORbY6mEWkTxXoXgUz3RM5Mcnr+bpDolsP3CMAR8soNVL0xj6w9oC9SXHAmnioKSkJJeamup1DJGA9MmCzfx2/GI61ivD67fUx+fL/Xk0zCzNOZeU6we+RGZWC/gaMLK+0DVzzm240H5qayRYZJ5yfLt8O8Nmrmfe+l8oEhHKDQ3K0j2lIleUivY63iXJaTujeXBECogbG5Zjx4HjvPDVCsrERfJk+1peR8oVZjYFKHWWp54CrgEecc59bGbdgPeAtuc4Tn+gP0CFChXyKK1I/grxGe3qlKZdndIs3byf4bPWMS51E6PmbKBxpaJ0T6lIuzqliAgN8TpqrlMPjkgB4pzj6YnpfDBnI3/tUoceKRVz9fh+2IOzH4hzzjnLmvp1v3Mu5kL7qa2RYLb38AkmpG1i9NyNbNhzhOKFw+nWuDy3J1egfLEor+NdkHpwROR/mBl/7libbfuO8aeJ6ZSJjeSaWiW9jpWXtgKtgO+Bq4FVnqYR8QNFC4fTv2VV+jWvwszVu/lgzgYGT1/DoOlraF0jgVuTK3B1zRKEBfgsyerBESmAjpw4yS2D57B65yHG3Z1C3XJxuXJcP+zBaQ68TtaXuWPAvc65tAvtp7ZGCpqt+44ydt5Gxs7fxM6Dx4kvEk6X+mXp1rg8NUr611idnLYzKnBECqidB49x49uzOZaRyaf3XpUrXdP+VuBcKrU1UlCdzDzFjFW7GD9/M9+t2EFGpqNe+Ti6JZWjY70yxESGeR0xx+1MYPc/icglKxEdyYg7G5OR6eg9fB77jpzwOpKIeCw0xMfVNUsyqEcj5jx5DU93SOR4RiZPfZpO479N4eGxC/lh1a6AmFdHPTgiBdzctXvo8d486leIY1Tf5Mu6mkI9OCLBxznH0i37mZC6mc8WbeHgsZMkREfQoW5putQvS91ysWSN4c8fOkUlIjn2+eKtPDhm4WXPkaMCRyS4HcvIZNqKnUxctJWpK3ZyIvMUleML06leGbo0KEvl+MJ5nkFXUYlIjnWqV4Yte4/ywlcrKBtXiCfa1/Q6koj4ociwENpfWZr2V5Zm/9EMvkrfxsRFW3lj6ipe/24V9crF0ql+WdrXKUWZuEKeZlWBIyIADGhVhc17jzBo+hrKFS1E91yeI0dEgktsoTBuaVyBWxpXYPv+Y/xr8VYmLt7CX79Yzl+/WE6DCnFcX6c07eqU8mR+HRU4IgJkzZHzbKfabNt/jGcmplM6+OfIEZFcUio2krtaVuGullVYu+sQk9O381X6dp6b9BPPTfqJK8vG0v7KUrSvUzpfTmOBxuCIyBkOHz/JrUMubY4cjcERkdNt+uUIk9O3MWnpdhZt2gdAzVLRtKtTimsTS5JYOuaiByhrkLGIXLKdB49xw8DZHD95ik/vbZbj7mUVOCJyLlv2HeWr9O1MXrqNtI17cQ7KxhXi2sSStK1VkiZViuVo9mTNgyMil6xEdCTv92nMiZOZPP7REq/jiEgQKBtXiL7NK/PRPc2Y94e2vND1SmqVjmbMvI10f28uU5bvyNXX82wMjpk9ANwHZAJfOuce9yqLiPyvaiWiGda7MSVjIr2OIiJBJiE64t8DlI+eyGTm6t00q1o8V1/DkwLHzNoAnYF6zrnjZlbCixwicn5JlYp5HUFEglyh8BCuTcz9Cxq8OkV1D/C8c+44gHNup0c5REREJAh5VeDUAFqY2Vwzm25mjc+1oZn1N7NUM0vdtWtXPkYUERGRQJVnp6jMbApQ6ixPPZX9usWAFKAxMN7MqrizXNLlnBsCDIGsKxvyKq+IiIgEjzwrcJxzbc/1nJndA3ySXdDMM7NTQDygLhoRERG5bF6dovoMaANgZjWAcGC3R1lEREQkyHgy0Z+ZhQPDgPrACeB3zrmpOdhvF7Ahhy8Tj38XTf6eD/w/o7/ng4KXsaJzLiGXjuWZi2hrCtq/b17w93ygjLkh39uZgJrJ+GKYWao/z6jq7/nA/zP6ez5QxmAXCO+dv2f093ygjLnBi3yayVhERESCjgocERERCTrBXOAM8TrABfh7PvD/jP6eD5Qx2AXCe+fvGf09Hyhjbsj3fEE7BkdEREQKrmDuwREREZECSgWOiIiIBJ2AK3DMrJ2Z/Wxmq83sibM8b2b2RvbzS8ysYU73zceMd2RnW2pms82s3mnPrc9+fJGZpXqUr7WZ7c/OsMjMnsnpvvmY8bHT8qWbWaaZFct+Lj/ew2FmttPM0s/xvD98Di+U0dPPob/z97bG39uZHGZUW3PhfH7d1vh1O+OcC5gbEAKsAaqQNfvxYiDxjG2uByYDRtZaV3Nzum8+ZmwGFM3+uf2vGbPvrwfiPX4PWwNfXMq++ZXxjO07AlPz6z3Mfo2WQEMg/RzPe/o5zGFGzz6H/n7z97bG39uZi8iotubCGf26rfHndibQenCSgdXOubXOuRPAWKDzGdt0Bka6LHOAODMrncN98yWjc262c25v9t05QLk8yHHJ+fJo37zMeBswJg9ynJNzbgbwy3k28fpzeMGMHn8O/Z2/tzX+3s7kKGMe7ZuXGdXWXGQ+Lz+HgVbglAU2nXZ/c/ZjOdkmJ/vmV8bT9SWr+v6VA6aYWZqZ9fcwX7PsbsXJZlb7IvfNr4yYWRTQDvj4tIfz+j3MCa8/hxcrvz+H/s7f2xp/b2dAbU1+CaS2Jl8/h3m2mrhcmJm1IesfvPlpDzd3zm0xsxLAt2a2IrtCzk8LgArOuUNmdj1Zi6NWz+cMOdURmOWcO/0bhD+8hwHDjz+Hkgv8/N9XbU0B4cXnMNB6cLYA5U+7Xy77sZxsk5N98ysjZlYXGAp0ds7t+fVx59yW7D93Ap+S1c2Yr/mccwecc4eyf54EhJlZfE72za+Mp7mVM7qM8+E9zAmvP4c54uHn0N/5e1vj7+1MjjKqrckVft/WePY5zKvBPXlxI6vHaS1Qmf8Mmqp9xja/4b8HXM3L6b75mLECsBpodsbjhYHo036eDbTzIF8p/jMJZDKwMfv99Jv3MHu7WLLO/RbOz/fwtNeqxLkH1nn6OcxhRs8+h/5+8/e2xt/bmYvIqLYmZzn9uq3x13Ym1/+ieX0ja8T4SrJGhz+V/dgAYED2zwYMzH5+KZB0vn09yjgU2Assyr6lZj9eJftDuBhYllcZc5Dv/uzXX0zWoLBm59vXi4zZ93sDY8/YL7/ewzHANiCDrHPbff3wc3ihjJ5+Dv395u9tjb+3MznMqLbmwvn8uq3x53ZGSzWIiIhI0Am0MTgiIiIiF6QCR0RERIKOChwREREJOipwREREJOiowBEREZGgowJHcp2ZxZnZved5vpCZTTezkPNsM8XMiuZNQhEJBmpr5HxU4EheiAPO2egAfYBPnHOZ59lm1AWOISKitkbOSQWO5IXngapmtsjMXjrL83cAEwHMrLSZzcjeNt3MWmRv8zlZK/eKiJyL2ho5J030J7nOzCoBXzjn6pzluXBgo3OuVPb9R4FI59xz2d3IUc65g9nPrQJS3Glrl4iI/EptjZyPVhOX/BYP7Dvt/nxgmJmFAZ855xad9txOoAygRkdELpbamgJOp6gkvx0FIn+945ybAbQka5XbEWbW87RtI7O3FxG5WGprCjgVOJIXDgLRZ3vCObcXCDGzSAAzqwjscM69S9aibA2zHzeyVhpenx+BRSQgqa2Rc1KBI7ku+zz2rOyBfGcb+PcN0Dz759bAYjNbCNwCvJ79eCNgjnPuZF7nFZHApLZGzkeDjCXfmVlD4BHnXI/zbPM68Llz7rv8SyYiwURtTcGmHhzJd865BcC0802+BaSrwRGRy6G2pmBTD46IiIgEHfXgiIiISNBRgSMiIiJBRwWOiIiIBB0VOCIiIhJ0VOCIiIhI0Pl/mh4RoxC0Y20AAAAASUVORK5CYII=",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3ef404b978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 数値解。　複数の戻り値を受け取るには以下のようにする。\n",
    "x_list, u_list, t_list = solve_freefall_air(1.0, 0.0, 0.001, 1.3, 0.1, 0.1)\n",
    "\n",
    "plt.figure(figsize=(8,3))\n",
    "# t vs x のグラフ\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(t_list, x_list)\n",
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('x (m)')\n",
    "# t vs u のグラフ\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(t_list, u_list)\n",
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('u (m/s)')\n",
    "# 2つのグラフを綺麗に配置する。\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 課題\n",
    "上記の関数を用いて、以下を数値的に調べよ。結果をグラフに描画せよ。\n",
    "\n",
    "+ 大きいリンゴと小さいリンゴはどちらが早く落ちるか? リンゴの密度は等しいものとする.\n",
    "+ 同じ大きさのリンゴと鉄球はどちらが早く落ちるか?"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 空気抵抗下での野球ボール軌道のシミュレーション\n",
    "\n",
    "ここでは、以下の図に示すような野球ボールの軌道をシミュレートする。\n",
    "ここで野球ボールは、以下のように時刻 $t=0$ で 原点($x=0, z=0$)から,\n",
    "斜め方向に $x, z$ 方向の速度成分 $u_x, u_z$ で打ち出される.\n",
    "\n",
    "<img src=\"figs/shahou.png\" width=\"80%\">\n",
    "\n",
    "野球ボール軌道の模式図. 時刻$t = 0$ に $x(0)=0, z(0)=0$ で\n",
    "斜め方向に速度$u_x(0), u_z(0)$で打ち出された球状の物体（図の赤丸）が,\n",
    "重力および空気抵抗の影響を受けて運動する様子をシミュレートする.\n",
    "位置$z$および$z$方向速度$u_z$の模式図は中・右図のようになる."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "この場合, 空気抵抗の$x, z$方向の成分はそれぞれ\n",
    "\n",
    "$$ \n",
    "\\begin{aligned}\n",
    "D_x &= -D u_x / |\\mathbf{u}| \\\\[3pt]\n",
    "D_z &= -D u_z / |\\mathbf{u}| \\\\[3pt]\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "と表される.\n",
    "\n",
    "運動方程式は, 速度に関する微分方程式\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{\\mathrm{d}u_x}{\\mathrm{d}t} &= \\frac{D_x}{m},\\\\\n",
    "\\frac{\\mathrm{d}u_z}{\\mathrm{d}t} &= \\frac{D_z}{m} - g,\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "と, 位置に関する微分方程式\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{\\mathrm{d}x}{\\mathrm{d}t} &= u_x,\\\\\n",
    "\\frac{\\mathrm{d}z}{\\mathrm{d}t} &= u_z,\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "と表される.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 課題\n",
    "\n",
    "大谷翔平は野球ボールを最大球速 165km/h で投げることができる。\n",
    "彼は地面に対しどの方向(角度 $\\theta$)に対してもこの初速で野球ボールを投げ出すことができ、\n",
    "ボールの回転は考慮しないものとする。\n",
    "また、彼の身長は193cmで腕の長さは85cmとする。\n",
    "\n",
    "これら仮定のもとで、以下の状況をシミュレーションせよ。\n",
    "\n",
    "+ 空気抵抗がない場合($C_D=0$), 適当な角度$\\theta$に対するシミュレーションをおこない, 解析解と比較せよ.\n",
    "ボールのリリース点の高さは適当に設定すること.\n",
    "+ 空気抵抗がある場合のシミュレーションをおこない, 空気抵抗がない場合と比較せよ.\n",
    "\n",
    "\n",
    "各物理量の具体的な値については, 以下の表を参考にするとよい.\n",
    "\n",
    "|   パラメータ |      値      |    備考   |\n",
    "|:----------:|:------------:|:---------:|\n",
    "|   $m$ [kg] | 0.145 | 野球ボール |\n",
    "|  半径$a$ [m]|0.036 | 野球ボール |\n",
    "|   重力加速度$g$ [m/s $^2$ ] | 9.80665 | - |\n",
    "|空気の密度 $\\rho$[kg/m $^3$ ] | 1.261 | 280Kの場合 |\n",
    "|                     　      | 1.176 | 300Kの場合 |\n",
    "|空気の動粘性係数 $\\nu \\mathrm([m^2/s])$| 1.395 $\\times 10^{-5}$ | 280Kの場合 |\n",
    "|                                     | 1.579 $\\times 10^{-5}$ | 300Kの場合 |\n",
    "|抗力係数$C_D$[-] | 0.44 | $5 \\times 10^2 < Re < 1 \\times 10^5$のとき |\n",
    "|Reynolds数$Re$[-] | - | $2aU/\\nu$ |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ヒント\n",
    "\n",
    "以下のような入力・出力を持つ関数を作成するとシンプルなコードで問題を解くことができる。\n",
    "\n",
    "```python\n",
    "def next_uxuz_air(ux, uz, dt, g, rho, a, m):\n",
    "    \"\"\"\n",
    "    現在の (ux, uz) の値から、dt 後の (ux, uz) を求める関数。\n",
    "    空気抵抗を考慮する\n",
    "    args:\n",
    "    + ux, uz: 時刻 t での(x-z)方向速度 [m/s]\n",
    "    + dt: 線形近似する時間幅 [s]\n",
    "    + g: 重力加速度の値 [m/s^2]\n",
    "    + rho: 空気の密度 [kg/m^3]\n",
    "    + a: 物体の半径 [m]\n",
    "    + m: 物体の質量 [kg]\n",
    "    return:\n",
    "     時刻t+dt での u の値\n",
    "    \"\"\"\n",
    "    \n",
    "    #  -----------------------------\n",
    "    #  なにか ux, uz, dt, g, rho, a, m から\n",
    "    #  ux_next, uz_next を求める計算\n",
    "    #  -----------------------------\n",
    "       \n",
    "    return ux_next, uz_next\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 課題\n",
    "\n",
    "空気抵抗のもとでの放物運動が関係するスポーツシーンを何でもよいので一つ取り上げ, 適当な問題設定に対するシミュレーションを実施せよ.\n",
    "\n",
    "［例］\n",
    "\n",
    "+ 野球ボールを最も遠くまで飛ばすためには,\n",
    "大谷翔平は地面に対しどの角度で野球ボールを投げる必要があるか.\n",
    "+ 大谷翔平が投げたボールがストライクになるためには, どのような角度で投げればよいか.\n",
    "初速, ピッチャーとキャッチャーの距離やストライクゾーンの広さは適当に設定せよ.\n",
    "+ サッカーのフリーキックで壁を越えてゴールに入れるにはどのような初速、角度で蹴ればよいか。\n",
    "ゴールまでの距離や壁の高さは自由に設定せよ。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注：\n",
    "<font color=\"Gray\">\n",
    "Reynolds数$Re=2aU/\\nu$ ($\\nu$は空気の動粘性係数, $U$は系の代表速度)と呼ばれる無次元数が\n",
    "$5 \\times 10^2 < Re < 1 \\times 10^5$のときにこの近似が有効である.\n",
    "</font>\n"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}