物体の放物運動【Pythonと物理#1】

今回のターゲット：物体の放物運動
Pythonで計算する：$x$と$y$を時間$t$で表す
実行結果：$x$と$y$を時間$t$で表す
Pythonで計算する：$y$を$x$で表す
実行結果：$y$を$x$で表す
最後に

今回のターゲット：物体の放物運動

　Pythonを使って物体の放物線運動を解いてみます。

下の図のように、原点Oにある質量mの物体を初速$v_{0}$で打ち出す場合を考えます。位置を $\boldsymbol{r} = (x, y)$ とすれば、運動方程式は、

$$m\frac{d^{2}}{dt^{2}}\boldsymbol{r} = m\boldsymbol{g}$$

となります。ただし、$\boldsymbol{g} = (0, -g)$とします。

成分に分けると、以下のようになります。

$$
\begin{eqnarray}
\left\{
\begin{array}{l}
m\cfrac{d^{2}x}{dt^{2}} = 0\\
m\cfrac{d^{2}y}{dt^{2}} = -mg
\end{array}
\right.
\end{eqnarray}
$$

1階微分方程式にしておきます。

$$
\begin{eqnarray}
\left\{
\begin{array}{l}
\cfrac{dx}{dt} = v_{0}\cos{\theta}\\
\cfrac{dy}{dt} = -gt+v_{0}\sin{\theta}
\end{array}
\right.
\end{eqnarray}
$$

Pythonで計算する：$x$と$y$を時間$t$で表す

それではPythonで計算していきます。コードは下のようにしました。

ただし、重力加速度$g$ = 9.8 m/s、初速$v_{0}$ = 30.0 m/s、$\theta$ = 45°、初期位置$\boldsymbol{r}_{0} = (0, 0)$としました。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint 
import matplotlib.pyplot as plt

#諸々の条件設定
g = 9.8 #重力加速度は9.8 m/s^2としておきます
x0 = 0.0 #初期位置 x = 0 mとしておきます
y0 = 0.0 #初期位置 y = 0 mとしておきます
v0 = 30.0 #初期速度 v0 = 30.0 m/sとしておきます

#θ=45°として、初速のx、y成分を求めておきます。
vx0 = v0/np.sqrt(2) #初期速度のx成分を計算しておきます
vy0 = v0/np.sqrt(2) #初期速度のy成分を計算しておきます
t = np.arange(0, 10, 0.1) #時間は100秒程度に設定しておきます

#計算に必要な関数(微分方程式)の定義
#x座標の計算（戻り値はx方向の速度）※わざわざ計算する意味はほとんどありません
def xcoordinate(x, t, g):
    return vx0
#y座標の計算（戻り値はy方向の速度）
def ycoordinate(y, t, g): #y座標の計算
    return -g * t + vy0

#x座標とy座標の計算　1階微分方程式を解く
#注意：odeintは微分方程式(おそらく1階微分方程式のみ?)を解くライブラリ
#引数は(解きたい微分方程式, 初期値, 変数, args = (パラメータ1, パラメータ2))
#なお、パラメータはタプルで渡して、1つしかない場合は、args = (パラメータ1, )とする
x = odeint(xcoordinate, x0, t, args = (g, ))
y = odeint(ycoordinate, y0, t, args = (g, ))

#以下は可視化のためにグラフを書くコード（詳細は割愛）
fig, ax = plt.subplots()

ax.plot(t, x, color='crimson', label='x')
ax.plot(t, y, color='cornflowerblue',  label='y')

ax.set_xlabel('Time [sec]')
ax.set_ylabel('x, y [m]')

ax.tick_params(direction='in', axis='x')
ax.tick_params(direction='in', axis='y')

ax.set_xlim(0, 5)
ax.set_ylim(0, 30)
ax.legend()

plt.show()

実行結果：$x$と$y$を時間$t$で表す

それでは実行結果を見てみましょう。下のようなグラフになると思います。

グラフを見ると、$x$方向は時間$t$に関する1次関数、$y$は時間$t$に関する2次関数になっています。これらの違いは力が働いているかどうか違いですね。
現実世界では$x$方向、$y$方向ともに空気抵抗が加わるので運動の様子がまた変わりますね。

Pythonで計算する：$y$を$x$で表す

続いて、$y$を$x$の式で表すと、以下の式になります。

$$
y = -\cfrac{g}{2v_{0}^{2}\cos{^{2}\theta}}x^{2} + \cfrac{\sin{\theta}}{\cos{\theta}}x
$$

この式から横軸$x$、縦軸$y$のグラフを描いてもいいのですが、先ほどのコードをほんの少しだけ変更するだけでできますね。以下のようなコードにしました。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint 
import matplotlib.pyplot as plt

#諸々の条件設定
g = 9.8 #重力加速度は9.8 m/s^2としておきます
x0 = 0.0 #初期位置 x = 0.0 mとしておきます
y0 = 0.0 #初期位置 y = 0.0 mとしておきます
v0 = 30.0 #初期速度 v0 = 30.0 m/sとしておきます

#θ=45°として、初速のx、y成分を求めておきます。
vx0 = v0/np.sqrt(2) #初期速度のx成分を計算しておきます
vy0 = v0/np.sqrt(2) #初期速度のy成分を計算しておきます
t = np.arange(0, 5, 0.1) #時間は5秒程度に設定しておきます

#計算に必要な関数(微分方程式)の定義
#x座標の計算（戻り値はx方向の速度）※わざわざ計算する意味はほとんどありません
def xcoordinate(x, t, g):
    return vx0
#y座標の計算（戻り値はy方向の速度）
def ycoordinate(y, t, g): #y座標の計算
    return -g * t + vy0

#x座標とy座標の計算　1階微分方程式を解く
#注意：odeintは微分方程式(おそらく1階微分方程式のみ?)を解くライブラリ
#引数は(解きたい微分方程式, 初期値, 変数, args = (パラメータ1, パラメータ2))
#なお、パラメータはタプルで渡して、1つしかない場合は、args = (パラメータ1, )とする
x = odeint(xcoordinate, x0, t, args = (g, ))
y = odeint(ycoordinate, y0, t, args = (g, ))

#以下は可視化のためにグラフを書くコード（詳細は割愛）
fig, ax = plt.subplots()

ax.plot(x, y, color='cornflowerblue') #←変更箇所

ax.set_xlabel('x [m]') #←変更箇所
ax.set_ylabel('y [m]') #←変更箇所

ax.tick_params(direction='in', axis='x')
ax.tick_params(direction='in', axis='y')

ax.set_xlim(0, 100) #←変更箇所
ax.set_ylim(0, 30)

plt.show()