葛西 真寿 弘前大学大学院理工学研究科
鉛直下向きの一様重力場中で，ある初速度をもって物体を空中に投げ出す。これを業界用語で「斜方投射」（Wikipedia 英語版タイトルでは Projectile motion）という。
空気抵抗が無視できる場合には，斜方投射された物体は放物線を描いて運動し，初速度一定の条件で水平な地上面から投射角 $\theta$ で投射された物体の水平到達距離が最大となるのは，$\theta = 45$° であることが知られている。
では，空気抵抗がある場合はどうなるか？
空気抵抗がある場合，初速度一定の条件で投射された物体の水平到達距離が最大となる投射角 $\theta$ を，ここでは gnuplot を使って求めてみよう。
まずは予想をたててみよう。
空気抵抗がある場合の斜方投射の問題。初速度一定の条件で投射された物体の水平到達距離が最大となる投射角 $\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$ とすると分母がゼロになって発散してしまいそうだけど...
斜方投射の問題は，しばしば運動が $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次関数として書けることを示せ。
解は
\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}set parametric で plot¶したがって gnuplot 的には，投射角 $\theta$ に対して，以下のように set parametric として，媒介変数表示で plot すればよい。
# th, つまりθを ° 単位で与える。
X(T, th) = T * cos(th*pi/180)
Y(T, th) = T * sin(th*pi/180) - 0.5 * T**2
set samples 5000 # 拡大表示するので多めに
set xrange [0:1.1] # 最大到達距離は無次元化された X で 1
set yrange [0:0.3] # y 軸の表示は地面にめり込まないようにゼロ以上に。
set xlabel "X"
set ylabel "Y"
set parametric # 媒介変数表示にする。
set grid # グリッドの表示
set title "空気抵抗がない場合"
plot [T=0:2] \
X(T, 44), Y(T, 44) title "θ = 44°", \
X(T, 45), Y(T, 45) title "θ = 45°" lc -1, \
X(T, 46), Y(T, 46) title "θ = 46°" ;
落下点付近を拡大表示します。確かに $\theta=45$° の時が最大水平到達距離（無次元化された座標では X = 1）であることがわかります。
set xrange [0.998:1.001]
set yrange [0:0.002]
replot
空気抵抗がない場合と同様に変数の無次元化を行います。また，空気抵抗の比例定数 $\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}set parametric で plot¶# 無次元化された解ではあるが，
# 空気抵抗なしの解と区別するため，小文字 x, y に。
x(T, b, th) = (1.-exp(-b*T))/b*cos(th*pi/180)
y(T, b, th) = (1.-exp(-b*T))/b*sin(th*pi/180)-(b*T-(1.-exp(-b*T)))/(b**2)
set samples 5000 # 拡大表示するので多めに
set xrange [0:1.1]
set yrange [0:0.3]
set xlabel "X"
set ylabel "Y"
set parametric
set grid
set title "空気抵抗がある場合 b = 0.5"
b = 0.5
plot [T=0:2] \
x(T, b, 44), y(T, b, 44) title "θ = 44°", \
x(T, b, 45), y(T, b, 45) title "θ = 45°" lc -1, \
x(T, b, 46), y(T, b, 46) title "θ = 46°" ;
落下点付近を拡大表示します。
set xrange [0.662:0.68]
set yrange [0:0.003]
replot
$44$° つまり $\theta < 45$° のほうが水平到達距離が大きくなっているようです。もう少し調べてみます。
set xrange [0.666:0.683]
set yrange [0:0.003]
set title "空気抵抗がある場合 b = 0.5"
b = 0.5
plot [T=0:2] \
x(T, b, 45), y(T, b, 45) title "θ = 45°" lc -1, \
x(T, b, 42), y(T, b, 42) title "θ = 42°", \
x(T, b, 41), y(T, b, 41) title "θ = 41°", \
x(T, b, 40), y(T, b, 40) title "θ = 40°", \
x(T, b, 39), y(T, b, 39) title "θ = 39°", \
x(T, b, 38), y(T, b, 38) title "θ = 38°";
$b = 0.5$ の場合，最大水平到達距離となるのは $\theta$ が 39° から 40° くらいの時であることがわかります。これをもう少し定量的に求めてみます。
$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 $$は解析的に解けないので，ニュートン法で数値的に解く関数を以下のように定義します。 ニュートン法については，例えば以下を参照。
# 再帰定義による解法。 y(t, b, th) = 0 となる t1 を求める。
# 最初，t には探索の初期値を入れる
# h は数値微分をするときに使う微小量，h*0.1 が解の収束条件である微小量。
newton_y(t, b, th, h) = \
(abs(y(t, b, th)) > h*0.1 ? \
newton_y(t-y(t,b,th)*2.*h/(y(t+h,b,th)-y(t-h,b,th)), b, th, h) : t)
ニュートン法によって $Y = 0$ の解である $t_1 (>0)$ がもとまったら，この時間を $X$ に代入して水平到達距離を求めます。
ニュートン法で求めた数値解の精度を確認するため，h の値をかえて計算します。
b = 0.5
th = 40
t1 = newton_y(1, b, th, 1.E-5); print x(t1, b, th)
t1 = newton_y(1, b, th, 1.E-6); print x(t1, b, th)
t1 = newton_y(1, b, th, 1.E-7); print x(t1, b, th)
h = 1.E-6 で小数点以下15桁の精度が出ているようですので，以下の計算ではこの値で続けます。
$\theta$ の値を 0.01° きざみで変えて計算し，水平到達距離が最大になる $\theta$ の値とそのときの最大水平到達距離 $X_{\rm max}$ を出力します。
gnuplot における条件分岐 (if) や繰り返し処理 (do for) については以下も参照。
Xm = 0
b = 0.5
do for [i = 0: 1500] {
th = 30 + 0.01*i
t1 = newton_y(1, b, th, 1.E-6)
x1 = x(t1, b, th)
if (x1 > Xm) {
Xm = x1
thm = th}
}
print "水平最大到達距離は θ = ", thm, "°のとき，Xmax = ", Xm
空気抵抗がある場合の水平到達距離が最大となる投射角 $\theta_m$ を，$b$ の値をいろいろ変えて求めてみよ。またそのときの最大水平到達距離も求め，$b$ との間にどのような関係が成り立つか調べてみよ。