葛西 真寿 弘前大学大学院理工学研究科
鉛直下向きの一様重力場中で，ある初速度をもって物体を空中に投げ出す。これを業界用語で「斜方投射」（Wikipedia 英語版タイトルでは Projectile motion）という。
空気抵抗が無視できる場合には，斜方投射された物体は放物線を描いて運動し，初速度一定の条件で水平な地上面から投射角 $\theta$ で投射された物体の水平到達距離が最大となるのは，$\theta = 45$° であることが知られている。
では，空気抵抗がある場合はどうなるか？
空気抵抗がある場合，初速度一定の条件で投射された物体の水平到達距離が最大となる投射角 $\theta$ を，ここでは Python（主に SymPy）を使って求めてみよう。
まずは予想をたててみよう。
空気抵抗がある場合の斜方投射の問題。初速度一定の条件で投射された物体の水平到達距離が最大となる投射角 $\theta$ は...
（以下では3次元ベクトルを \boldsymbol{} を使って $\boldsymbol{r}$ と書きたいところだが，なんとなく Jupyter Notebook では普通の文字 $r$ と区別がつきにくい懸念があるので，3次元ベクトルは \vec{} を使って $\vec{r}$ などと書くようにする。）
鉛直上向きを $+y$ 方向として，質量 $m$ の質点（物体）の運動方程式は，空気抵抗がない場合，以下のように書ける。
$$ m \frac{d^2\vec{r}}{dt^2} = - m\vec{g} $$ここで，$\vec{g} = (0, g, 0)$ は重力加速度ベクトルであり，$g$ は重力加速度の大きさ。
この式は直ちに解くことができて（両辺を時間で2回積分するだけ），初期条件 $t = 0$ で $\vec{r}(t=0) = \vec{r}_0, \vec{v}(t=0) = \vec{v}_0$ を使って以下のように書ける。
$$ \vec{r}(t) = \vec{r}_0 + \vec{v}_0\, t - \frac{1}{2} \vec{g}\, t^2 $$空気抵抗がある場合は，質点の速さがそれほど大きくない場合には，速度に比例するとしてよいだろうから，比例定数を $- m \beta$ として，運動方程式は以下のようになる。（マイナスがつくのは，抵抗というのは運動を妨げる向きに働くから。）
$$ m \frac{d^2\vec{r}}{dt^2} = - m\vec{g} - m \beta \frac{d\vec{r}}{dt} $$この運動方程式の解は
$$ \vec{r}(t) = \vec{r}_0 + \frac{1}{\beta} \vec{v}_0\,(1 - e^{-\beta t}) - \frac{1}{\beta^2} \vec{g}\,\left\{\beta t - \left(1 - e^{-\beta t}\right)\right\} $$となる。以下で詳細を説明する。
上の運動方程式の両辺を $m$ でわり，$\vec{v} = \frac{d\vec{r}}{dt}$ を使って書くと...
\begin{eqnarray} \frac{d\vec{v}}{dt} + \beta \vec{v} &=& - \vec{g} \\ e^{-\beta t} \frac{d}{dt} \left(\vec{v}\, e^{\beta t}\right) &=& - \vec{g} \\ \frac{d}{dt} \left(\vec{v}\, e^{\beta t}\right) &=& - \vec{g}\,e^{\beta t} \end{eqnarray}となる。両辺を $t$ で積分すると，
\begin{eqnarray} e^{\beta t} \,\vec{v} &=& -\frac{1}{\beta} \vec{g}\, e^{\beta t} + \vec{C}\\ \therefore\ \vec{v} &=& -\frac{1}{\beta} \vec{g} + \vec{C}\, e^{-\beta t} \end{eqnarray}積分定ベクトル $\vec{C}$ は初期条件を使って書ける。$t = 0$ で
\begin{eqnarray} \vec{v}_0 &=& -\frac{1}{\beta} \vec{g} + \vec{C} \\ \therefore\ \vec{C} &=& \vec{v}_0 + \frac{1}{\beta} \vec{g} \end{eqnarray}この $\vec{C}$ を代入すると，
\begin{eqnarray} \vec{v} = \frac{d\vec{r}}{dt} &=& -\frac{1}{\beta} \vec{g} + \left( \vec{v}_0 + \frac{1}{\beta} \vec{g} \right) e^{-\beta t} \\ &=& \vec{v}_0\, e^{-\beta t} - \frac{1}{\beta} \vec{g}\, (1 - e^{-\beta t}) \end{eqnarray}両辺をもう一度，$t$ で積分して...
\begin{eqnarray} \vec{r}(t) &=& -\frac{1}{\beta} \vec{v}_0\, e^{-\beta t} -\frac{1}{\beta}\vec{g}\,\left(t + \frac{1}{\beta} e^{-\beta t}\right) + \vec{C}' \\ &=& -\frac{1}{\beta} \vec{v}_0\, e^{-\beta t} -\frac{1}{\beta^2}\vec{g}\,\left(\beta t + \ e^{-\beta t}\right) + \vec{C}' \end{eqnarray}積分定ベクトル $\vec{C}'$ は初期条件を使って書ける。$t = 0$ で
\begin{eqnarray} \vec{r}_0 &=& -\frac{1}{\beta} \vec{v}_0 -\frac{1}{\beta^2}\vec{g} + \vec{C}' \\ \therefore\ \vec{C}' &=& \vec{r}_0 + \frac{1}{\beta} \vec{v}_0 + \frac{1}{\beta^2}\vec{g} \end{eqnarray}この $\vec{C}'$ を代入すると最終的な解が得られる。
空気抵抗がある場合の解
$$ \vec{r}(t) = \vec{r}_0 + \frac{1}{\beta} \vec{v}_0\,(1 - e^{-\beta t}) - \frac{1}{\beta^2} \vec{g}\,\left\{\beta t - \left(1 - e^{-\beta t}\right)\right\} $$が，$\beta \rightarrow 0$ で空気抵抗がない場合の解になることを示せ。単純に $\beta=0$ とすると分母がゼロになって発散してしまいそうだけど...
# SymPy をつかうときのおまじない
from sympy import *
from sympy.abc import *
from sympy import pi
# グラフを SVG でインライン表示するときのおまじない
from IPython.display import set_matplotlib_formats
%matplotlib inline
set_matplotlib_formats('svg')
ヒント：
SymPy を import すると，極限 $\displaystyle \lim_{\beta\rightarrow 0} $ が使えます。以下では $y$ 成分について $\beta\rightarrow 0$ の極限をとっています。（Python では文末が ; で終わる場合は結果が表示されません。）
# sympy.abc を import すると，文字変数は全て Symbol となるようですが，
# y0 などのような変数は以下のように宣言する必要あり。
var('y0 v0')
ybeta = y0 + 1/beta * v0 * (1-exp(-beta*t)) \
- 1/beta**2 * g * (beta*t - (1 - exp(-beta*t)))
limit(ybeta, beta, 0);
斜方投射の問題は，しばしば運動が $xy$ 平面内に限られることを自明であるかのごとく扱われるが，運動が（変直方法に平行な）平面内に限られることは以下のように説明することができる。
上記の解から直ちにわかるが，空気抵抗があろうがなかろうが，$\vec{r}_0 = \vec{0}$ とすると，位置ベクトル $\vec{r}(t)$ は2つの定ベクトル $\vec{v}_0$ と $\vec{g}$ の線形結合で書かれている。
つまり，位置ベクトル $\vec{r}(t)$ は2つのベクトルで張られる平面内に限られることになり，斜方投射された物体の運動は，この平面内に限られることになる。
これが斜方投射された物体の運動が平面内に限られる理由である。
以下では簡単のために，その平面を $xy$ 平面とし，
$$\vec{r} = (x, y, 0), \quad \vec{r}_0 = (0, 0, 0), \quad \vec{v}_0 = (v_0 \cos\theta, v_0 \sin\theta, 0)$$とする。
空気抵抗がない場合の解が放物線になること，つまり，$y$ が $x$ の2次関数として書けることを示せ。
ヒント：
$x, y$ それぞれが $t$ の関数として表されているので，$t$ を消去して，$y$ を $x$ の関数として表す，ということ。$y$ が $x$ の2次関数として書けるということは，$y$ の極値（斜方投射の場合は最大値）と，その値をとるときの $x$ や $t$ の値もすぐわかる。
すでに述べたように，斜方投射の運動は鉛直上向きを $+y$ 方向とした $xy$ 平面内に限られるとして話を続ける。また初期条件は，
$$\vec{r} = (x, y, 0), \quad \vec{r}_0 = (0, 0, 0), \quad \vec{v}_0 = (v_0 \cos\theta, v_0 \sin\theta, 0)$$解は
\begin{eqnarray} x &=& v_0\,t\, \cos\theta\\ y &=& v_0\,t\, \sin\theta - \frac{1}{2} g\,t^2 \end{eqnarray}$t = 0$ で投射された物体がふたたび地上に戻る時間 $t_1$ は $ y = 0$ を解いて
$$ t_1 = \frac{2\,v_0\,\sin\theta}{g} $$この時間を $x$ の式に代入すると，水平到達距離 $x_1$ は
$$ x_1 = v_0\, t_1\, \cos\theta = \frac{v_0^2 \sin 2\theta}{g} $$$x_1$ が最大となるのは $\sin 2\theta = 1$ すなわち $\displaystyle \theta = \frac{\pi}{4}$（45°）のときで，そのときの最大水平到達距離 $x_m$ は $$ x_m = \frac{v_0^2}{g} $$
このままでもいいですが，初速度の大きさをいくらにするかとか，重力加速度の大きさはいくらだったかなとかに煩わされるにすむように，この系に特徴的な長さと時間で無次元化しておきます。
まず，特徴的な長さとして $x_m$，特徴的な時間として $\displaystyle \tau \equiv \frac{v_0}{g}$ を採用し，無次元化された座標および時間を以下のように定義する。
$$ X \equiv \frac{x}{x_m}, \quad Y \equiv \frac{y}{x_m}, \quad T \equiv \frac{t}{\tau} $$無次元化されたこれらの量で解を書くと
\begin{eqnarray} X &=& T \cos\theta \\ Y &=& T \sin\theta - \frac{1}{2} T^2 \end{eqnarray}また，投射された物体がふたたび地上に戻る時間 $t_1$ も無次元化して
$$ T_1 \equiv \frac{t_1}{\tau} = 2\sin\theta$$plot_parametric()¶したがって SymPy 的には，投射角 $\theta$ を媒介変数として，以下のように plot_parametric() すればよい。
def X(T, th):
return T * cos(th*float(pi)/180)
def Y(T, th):
return T * sin(th*float(pi)/180) - 0.5 * T**2
def T1(th):
return 2 * sin(th*float(pi)/180)
p1 = plot_parametric(
(X(T, 45), Y(T, 45), (T, 0, T1(45))),
(X(T, 44), Y(T, 44), (T, 0, T1(44))),
(X(T, 46), Y(T, 46), (T, 0, T1(46))), legend=True,
title="空気抵抗がない場合"
)
SymPy の plot の残念なところは，複数のグラフを書く時，デフォルトでは同じ色の線で描くことです。これでは凡例をつけても区別できません。
そこで，以下のように3本の線に個別に凡例と線の色を設定します。
p1[0].label="θ = 45°"
p1[1].label="θ = 44°"
p1[2].label="θ = 46°"
p1[0].line_color="red"
p1[1].line_color="blue"
p1[2].line_color="green"
p1.show()
落下点付近を拡大表示します。確かに $\theta=45$° の時が最大水平到達距離（無次元化された座標では X = 1）であることがわかります。
p1.xlim=(0.998, 1.001)
p1.ylim=(0, 0.002)
p1.axis_center=(0.998, 0)
p1.show()
それぞれの線に凡例と線色を設定する別解です。
p2 = plot_parametric(
(X(T, 45), Y(T, 45), (T, 0, T1(45))), legend=True,
line_color="red", label="θ = 45°",
title="空気抵抗がない場合", show=False)
p2.extend(
plot_parametric(
(X(T, 44), Y(T, 44), (T, 0, T1(44))),
line_color="blue", label="θ = 44°", show=False)
)
p2.extend(
plot_parametric(
(X(T, 46), Y(T, 46), (T, 0, T1(46))),
line_color="green", label="θ = 46°", show=False)
)
p2.show()
SymPy の plot() や plot_parametric() には，それ自体ではグリッド線を表示させるオプションが見当たらないので，以下のように seaborn を import し，スタイル設定してみます。
seaborn を import 後，日本語フォント部分がトーフになるようなら，以下のようにplt.rcParams[] で日本語フォントを再設定します。以下の例では Noto Sans CJK JP を設定していますが，各自の環境に合わせて設定してください。
# seaborn のスタイルを利用してグリッドを表示
import seaborn as sns
# sns.set_style("whitegrid", {'grid.linestyle': '--'})
sns.set_style("whitegrid")
# seaborn を import したら日本語フォントを再設定
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Noto Sans CJK JP'
p2.show()
空気抵抗がない場合と同様に変数の無次元化を行います。また，空気抵抗の比例定数 $\beta$ は時間の逆数の次元をもつので，以下のように無次元化します。
$$ b \equiv \beta \, \tau = \beta \, \frac{v_0}{g} $$無次元化された変数で書いた解は，
\begin{eqnarray} X &=& \frac{1-e^{-b T}}{b} \, \cos\theta \\ Y &=& \frac{1-e^{-b T}}{b} \, \sin\theta - \frac{b T - (1-e^{-b T})}{b^2} \end{eqnarray}plot_parametric()¶def x(T, b, th):
return (1-exp(-b*T))/b*cos(th*float(pi)/180)
def y(T, b, th):
return (1-exp(-b*T))/b*sin(th*float(pi)/180)-(b*T-(1-exp(-b*T)))/(b**2)
b = 0.5;
p3 = plot_parametric(
(x(T, b, 45), y(T, b, 45)),
(x(T, b, 44), y(T, b, 44)),
(x(T, b, 46), y(T, b, 46)), (T, 0, 3), legend=True,
xlim = (0, 1.1), ylim = (0, 0.3),
title="空気抵抗がある場合 b = 0.5", show=False
)
p3[0].label="θ = 45°"
p3[1].label="θ = 44°"
p3[2].label="θ = 46°"
p3[0].line_color="red"
p3[1].line_color="blue"
p3[2].line_color="green"
p3.show()
落下点付近を拡大表示します。
p3.xlim=(0.662, 0.68)
p3.ylim=(0, 0.003)
p3.axis_center=(0.662, 0)
p3.show()
$44$° つまり $\theta < 45$° のほうが水平到達距離が大きくなっているようです。もう少し調べてみます。
b = 0.5;
p4 = plot_parametric(
(x(T, b, 45), y(T, b, 45)),
(x(T, b, 42), y(T, b, 42)),
(x(T, b, 41), y(T, b, 41)),
(x(T, b, 40), y(T, b, 40)),
(x(T, b, 39), y(T, b, 39)),
(x(T, b, 38), y(T, b, 38)), (T, 0, 3), legend=True,
xlim = (0, 1.1), ylim = (0, 0.3),
title="空気抵抗がある場合 b = 0.5", show=False
)
p4[0].label="θ = 45°"
p4[1].label="θ = 42°"
p4[2].label="θ = 41°"
p4[3].label="θ = 40°"
p4[4].label="θ = 39°"
p4[5].label="θ = 48°"
p4[0].line_color="red"
p4[1].line_color="blue"
p4[2].line_color="lightblue"
p4[3].line_color="darkblue"
p4[4].line_color="pink"
p4[5].line_color="lightgreen"
p4.xlim=(0.665, 0.682)
p4.ylim=(0, 0.003)
p4.axis_center=(0.665, 0)
p4.show()
$b = 0.5$ の場合，最大水平到達距離となるのは $\theta$ が 39° から 40° くらいの時であることがわかります。これをもう少し定量的に求めてみます。
nsolve()¶$T = 0$ で投射された物体がふたたび地上に戻る時間 $T_1$ は $ Y = 0$ を解いて求める必要がありますが，
$$ Y = \frac{1-e^{-b T}}{b} \, \sin\theta - \frac{b T - (1-e^{-b T})}{b^2} = 0 $$は解析的に解けないので，数値的に解く関数 nsolve() を使ってみます。
nsolve() によって $Y = 0$ の解である $T_1 (>0)$ がもとまったら，この時間を $X$ に代入して水平到達距離を求めます。
$\theta = 45$° のときの水平到達距離を求める例：
b = 0.5;
th = 40;
T1 = nsolve(y(T, b, th), T, 1, prec=16)
x1 = x(T1, b, th)
print("投射角 θ = ", th, "° のときの水平到達距離は ", x1)
$\theta$ の値を 0.01° きざみで変えて計算し，水平到達距離が最大になる $\theta$ の値とそのときの最大水平到達距離 $X_{\rm max}$ を出力します。
b = 0.5;
Xm = 0;
for i in range(1501):
th = 30 + 0.01*i
T1 = nsolve(y(T, b, th), T, 1, prec=16)
x1 = x(T1, b, th)
if x1 > Xm:
thm = th
Xm = x1
tm = T1
print("最大水平到達距離は θ = ", thm, "°のとき，Xmax = ", Xm)
空気抵抗がある場合の水平到達距離が最大となる投射角 $\theta_m$ を，$b$ の値をいろいろ変えて求めてみよ。またそのときの最大水平到達距離も求め，$b$ との間にどのような関係が成り立つか調べてみよ。