葛西真寿　弘前大学大学院理工学研究科

空気抵抗がある場合の斜方投射¶

鉛直下向きの一様重力場中で，ある初速度をもって物体を空中に投げ出す。これを業界用語で「斜方投射」（Wikipedia 英語版タイトルでは Projectile motion）という。

空気抵抗が無視できる場合には，斜方投射された物体は放物線を描いて運動し，初速度一定の条件で水平な地上面から投射角 $\theta$ で投射された物体の水平到達距離が最大となるのは，$\theta = 45$° であることが知られている。

では，空気抵抗がある場合はどうなるか？

空気抵抗がある場合，初速度一定の条件で投射された物体の水平到達距離が最大となる投射角 $\theta$ を，ここでは Maxima を使って求めてみよう。

問題（3択）¶

まずは予想をたててみよう。

空気抵抗がある場合の斜方投射の問題。初速度一定の条件で投射された物体の水平到達距離が最大となる投射角 $\theta$ は...

45° より小さくなる。
やっぱり 45° である。
45° より大きくなる。

斜方投射の運動方程式¶

（以下では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$ とすると分母がゼロになって発散してしまいそうだけど...

ヒント：

Maxima では $\displaystyle \lim_{\beta\rightarrow 0} $ が使えます。以下では $y$ 成分について $\beta\rightarrow 0$ の極限をとっています。

ybeta: y_0 + 1/beta * v_0 * (1-exp(-beta*t))
           - 1/beta**2 * g * (beta*t - (1 - exp(-beta*t)))$

limit(ybeta, beta, 0)$
expand(%)$

運動が平面内に限られる理由¶

斜方投射の問題は，しばしば運動が $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$ の値もすぐわかる。

最大水平到達距離と投射角¶

空気抵抗がない場合¶

解は

\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$$

`parametric` で `plot2d`¶

したがって Maxima 的には，投射角 $\theta$ に対して，以下のように parametric として，媒介変数表示で plot2d すればよい。

/* Jupyter Notebook にインラインでグラフを表示させる場合，*/
/* 最初に以下のようにプロットオプションをセット。 */
set_plot_option([svg_file, "~/.maxplot.svg"])$

pi: float(%pi)$
X(T, th) := T * cos(th*pi/180)$
Y(T, th) := T * sin(th*pi/180) - 0.5 * T**2$
T1(th):= 2*sin(th*pi/180)$

plot2d([[parametric, X(T, 45), Y(T, 45), [T, 0, T1(45)]], 
        [parametric, X(T, 44), Y(T, 44), [T, 0, T1(44)]],
        [parametric, X(T, 46), Y(T, 46), [T, 0, T1(46)]]], 
       [x, 0, 1.1], [y, 0, 0.3], 
       [legend, "θ = 45°", "θ = 44°", "θ = 46°"], grid2d,
       [title, "空気抵抗がない場合"])$

落下点付近を拡大表示します。確かに $\theta=45$° の時が最大水平到達距離（無次元化された座標では X = 1）であることがわかります。

plot2d([[parametric, X(T, 45), Y(T, 45), [T, T1(45)*0.995, T1(45)]],
        [parametric, X(T, 44), Y(T, 44), [T, T1(44)*0.995, T1(44)]], 
        [parametric, X(T, 46), Y(T, 46), [T, T1(46)*0.995, T1(46)]]], 
       [x, 0.998, 1.001], [y, 0, 0.002], [nticks, 100],
       [legend, "θ = 45°", "θ = 44°", "θ = 46°"], grid2d,
       [title, "空気抵抗がない場合"])$

plot2d: some values were clipped.
plot2d: some values were clipped.
plot2d: some values were clipped.

空気抵抗がある場合¶

無次元化¶

空気抵抗がない場合と同様に変数の無次元化を行います。また，空気抵抗の比例定数 $\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}

`parametric` で `plot2d`¶

/* 無次元化された解ではあるが，
   空気抵抗なしの解と区別するため，小文字 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)$

b: 0.5$
plot2d([[parametric, x(T, b, 45), y(T, b, 45), [T, 0, 3]], 
        [parametric, x(T, b, 44), y(T, b, 44), [T, 0, 3]],
        [parametric, x(T, b, 46), y(T, b, 46), [T, 0, 3]]], 
       [x, 0, 1.1], [y, 0, 0.3], [nticks, 100],
       [legend, "θ = 45°", "θ = 44°", "θ = 46°"], grid2d,
       [title, "空気抵抗がある場合 b=0.5"])$

plot2d: some values were clipped.
plot2d: some values were clipped.
plot2d: some values were clipped.

落下点付近を拡大表示します。

b: 0.5$
plot2d([[parametric, x(T, b, 45), y(T, b, 45), [T, 1, 2]], 
        [parametric, x(T, b, 44), y(T, b, 44), [T, 1, 2]],
        [parametric, x(T, b, 46), y(T, b, 46), [T, 1, 2]]], 
       [x, 0.662, 0.68], [y, 0, 0.003], [nticks, 1000],
       [legend, "θ = 45°", "θ = 44°", "θ = 46°"], grid2d,
       [title, "空気抵抗がある場合 b=0.5"])$

plot2d: some values were clipped.
plot2d: some values were clipped.
plot2d: some values were clipped.

$44$° つまり $\theta < 45$° のほうが水平到達距離が大きくなっているようです。もう少し調べてみます。

b: 0.5$
plot2d([[parametric, x(T, b, 45), y(T, b, 45), [T, 1, 2]], 
        [parametric, x(T, b, 42), y(T, b, 42), [T, 1, 2]],
        [parametric, x(T, b, 41), y(T, b, 41), [T, 1, 2]],
        [parametric, x(T, b, 40), y(T, b, 40), [T, 1, 2]],
        [parametric, x(T, b, 39), y(T, b, 39), [T, 1, 2]],
        [parametric, x(T, b, 38), y(T, b, 38), [T, 1, 2]]], 
    [x, 0.666, 0.683], [y, 0, 0.003], [nticks, 2000],
    [legend, "θ = 45°", "θ = 42°", "θ = 41°", "θ = 40°", "θ = 39°", "θ = 38°"], 
    grid2d,
    [title, "空気抵抗がある場合 b=0.5"])$

plot2d: some values were clipped.
plot2d: some values were clipped.
plot2d: some values were clipped.
plot2d: some values were clipped.
plot2d: some values were clipped.
plot2d: some values were clipped.

$b = 0.5$ の場合，最大水平到達距離となるのは $\theta$ が 39° から 40° くらいの時であることがわかります。これをもう少し定量的に求めてみます。

`find_root()`¶

$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 $$

は解析的に解けないので，数値的に解く関数 find_root() を使ってみます。以下を参照。

Maxima で数値解析 - 弘前大学 Home Sweet Home

find_root() によって $Y = 0$ の解である $T_1 (>0)$ がもとまったら，この時間を $X$ に代入して水平到達距離を求めます。

$\theta = 45$° のときの水平到達距離を求める例：

b: 0.5$
th: 40$
T1: find_root(y(T, b, th)=0, T, 0.1, 2)$
x1: float(x(T1, b, th))$
print("投射角 θ = ", th, "° のときの水平到達距離は ", x1)$

$\theta$ の値を 0.01° きざみで変えて計算し，水平到達距離が最大になる $\theta$ の値とそのときの最大水平到達距離 $X_{\rm max}$ を出力します。

Maxima における条件分岐 (if) や繰り返し処理 (for) については以下も参照。

if 条件式 then のあとに複数の文を書く書き方がわからなくて苦労しましたが，以下の例のように，かっこでくくり，カンマで区切って書くようです。

if 条件式 then 
    (実行文1, 実行文2, 実行文3)

b: 0.5$
Xm: 0$

for i: 0 thru 1500 
    do(
        th: 30 + 0.01*i,
        T1: find_root(y(T, b, th)=0, T, 0.1, 2),
        x1: float(x(T1, b, th)),
        if x1 > Xm then 
            (thm: th, Xm: x1, tm: T1)
    )$
print("最大水平到達距離は θ = ", thm, "°のとき，Xmax = ", Xm)$

問題¶

空気抵抗がある場合の水平到達距離が最大となる投射角 $\theta_m$ を，$b$ の値をいろいろ変えて求めてみよ。またそのときの最大水平到達距離も求め，$b$ との間にどのような関係が成り立つか調べてみよ。

空気抵抗がある場合の斜方投射¶

問題（3択）¶

斜方投射の運動方程式¶

空気抵抗がない場合¶

空気抵抗がある場合¶

問題¶

運動が平面内に限られる理由¶

問題¶

最大水平到達距離と投射角¶

空気抵抗がない場合¶

無次元化¶

parametric で plot2d¶

空気抵抗がある場合¶

無次元化¶

parametric で plot2d¶

find_root()¶

問題¶

`parametric` で `plot2d`¶

`parametric` で `plot2d`¶

`find_root()`¶