無次元化された解

空気抵抗がある場合の無次元化された解は，

\begin{eqnarray}
X &=& \frac{1-e^{-b T}}{b} \, \cos\theta \\
Y &=& H + \frac{1-e^{-b T}}{b} \, \sin\theta + \frac{1 – b T -e^{-b T}}{b^2}
\end{eqnarray}

詳細は以下のページ

空気抵抗がある場合の斜方投射を調べる準備

In [1]:

/* 度をラジアンに変換する関数を準備しておく */
radians(deg):= deg * %pi/180$

X(T, b, th):= (1 - exp(-b*T))/b * cos(radians(th))$
Y(T, b, th, H):= H + (1 - exp(-b*T))/b * sin(radians(th))
                   + (1 - b*T - exp(-b*T))/(b**2)$

空気抵抗がゼロの極限

本題とは直接関係ないですが，$b \rightarrow 0$ の極限で空気抵抗なしの斜方投射の解になることを確認しておきます。

$b$ が分母にあるので，単に $b=0$ とすると大変なことが起こりそうですが，きちんと極限をとると…

In [2]:

thet: theta*180/%pi$
limit(X(T, b, thet), b, 0);
limit(Y(T, b, thet, H), b, 0), expand;

Out[2]:

\[\tag{${\it \%o}_{5}$}T\,\cos \vartheta\]

Out[2]:

\[\tag{${\it \%o}_{6}$}T\,\sin \vartheta-\frac{T^2}{2}+H\]

滞空時間と水平到達距離

$T = 0$ で打ち出された粒子が地面 $Y=0$ に落ちるまでの滞空時間 $T_1$ は $Y(T, b, \theta, H) = 0$ となる $T$ を数値的に解いて求める。

この時間での水平方向の到達距離 $L$ は，

$$L = X(T_1, b, \theta)$$

In [3]:

/* どの範囲で find_root するか。ここでは 0.1 ~ 3.0 と広めに */
T1(b, th, H):= find_root(Y(T, b, th, H)=0, T, 0.1, 3.0)$
L(b, th, H):= X(T1(b, th, H), b, th)$

$H = 0$ の場合の軌道

$\theta = 45^{\circ}$ 前後の角度での斜方投射の軌道をグラフにしてみる。

In [4]:

/* θ = 45° のグラフ */
/* 必要なパラメータの数値設定 */
th1: 45$
b1: 0.5$
H1: 0$

draw2d(
  /* 滑らかにするためにサンプリングを多めに */
  nticks = 100, 
  /* 縦横比。*/
  proportional_axes = xy,
  /* 表示範囲とグリッド線 */
  xrange = [0, 0.8], yrange = [0, 0.4], grid = true,
  /* 軸のラベル */
  xlabel = "X", ylabel = "Y", 
  /* x軸 y軸の目盛 */
  xtics = 0.1, ytics = 0.1,
  /* 数値を取り込んだタイトルの設定 */
  title = concat("b = ", b1, ", H = ", H1, " の場合の斜方投射"),

  /* 数値を取り込んだ凡例の設定 */
  key = concat("θ_0 = ", th1, "°"), 
  /* 媒介変数表示で曲線を描く */
  parametric(X(T, b1, th1), 
             Y(T, b1, th1, H1), 
             T, 0, T1(b1, th1, H1))
)$

In [5]:

/* 複数のグラフを一度に描く例 */
b1: 0.5$
H1: 0$

draw2d(
  /* オプションの設定部分は同じ */
  /* 滑らかにするためにサンプリングを多めに */
  nticks = 100, 
  /* 縦横比。*/
  proportional_axes = xy,
  /* 表示範囲とグリッド線 */
  xrange = [0, 0.8], yrange = [0, 0.4], grid = true,
  /* 軸のラベル */
  xlabel = "X", ylabel = "Y", 
  /* x軸 y軸の目盛 */
  xtics = 0.1, ytics = 0.1,
  /* 数値を取り込んだタイトルの設定 */
  title = concat("b = ", b1, ", H = ", H1, " の場合の斜方投射"),

  /* 角度 th1 を変えて以下のように羅列すれば描けるが... */
  /* th = 45 */
  key = concat("θ_0 = ", 45, "°"), color = 1, 
  parametric(X(T, b1, 45), 
             Y(T, b1, 45, H1), 
             T, 0, T1(b1, 45, H1)),
  /* th = 44 */
  key = concat("θ_0 = ", 44, "°"), color = 2,  
  parametric(X(T, b1, 44), 
             Y(T, b1, 44, H1), 
             T, 0, T1(b1, 44, H1))
)$

In [6]:

/* 複数のグラフを一度に描く場合は... */
/* 最初に makelist で作っておくといいかも */
thmax: 45$
thmin: 38$

lines: makelist(
  [key = concat("θ_0 = ", th),
   color = mod(thmax, th), 
   parametric(X(T, b1, th), 
              Y(T, b1, th, H1), 
              T, 0, T1(b1, th, H1))
  ], 
  th, thmax, thmin, -1)$

draw2d(
  /* オプションの設定部分は同じ */
  /* 滑らかにするためにサンプリングを多めに */
  nticks = 100, 
  /* 縦横比。*/
  proportional_axes = xy,
  /* 表示範囲とグリッド線 */
  xrange = [0, 0.8], yrange = [0, 0.4], grid = true,
  /* 軸のラベル */
  xlabel = "X", ylabel = "Y", 
  /* x軸 y軸の目盛 */
  xtics = 0.1, ytics = 0.1,
  /* 数値を取り込んだタイトルの設定 */
  title = concat("b = ", b1, ", H = ", H1, " の場合の斜方投射"),

  /* 線の数だけ羅列する必要なく，作成済みのリストを書くだけ */
  lines
)$

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

In [7]:

/* 着地点付近を拡大表示 */
draw2d(
  nticks = 100, 
  /* 表示範囲を狭くして拡大表示する */
  xrange = [0.66, 0.682], yrange = [0, 0.01], grid = true,
  xlabel = "X", ylabel = "Y", 
  title = concat("b = ", b1, ", H = ", H1, " の場合の斜方投射"),

  lines
)$

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

水平到達距離 L

定量的に調べるために，水平到達距離 L(b, th, H) をグラフにしてみる。

In [8]:

draw2d(
  grid = true, 
  xlabel = "打ち出し角度 θ_0 (°)", 
  ylabel = "規格化された水平到達距離 L",
  title = concat("b = ", b1, ", H = ", H1, 
                 " の場合の打ち出し角度と水平到達距離"),
  
  explicit(L(b1, th, H1), th, 38, 42)
)$

最大水平到達距離を数値微分で求める

$L(b, \theta, H)$ が最大となるのは $\displaystyle \frac{\partial L}{\partial \theta} = 0$ となる角度のときです。

$L(b, \theta, H)$ は解析的な関数として与えられていないので，ここでは以下のような数値微分で求めることにします。

$$\frac{\partial L}{\partial \theta} \simeq \frac{L(b, \theta+h, H)-L(b, \theta-h, H)}{2 h}$$

In [9]:

/* 数値微分: h は小さいほどよい，というわけではない */
dL(b, th, H, h):= (L(b, th+h, H)-L(b, th-h, H))/(2*h)$

In [10]:

/* 方程式の数値解 */
/* h を変えてみる */
for h in [1.e-3, 1.e-4, 1.e-5, 1.e-6, 1.e-7] do(
  print(find_root(dL(b1, th, H1, h) = 0, th, 38, 45))
)$

$39.48760970738729$

$39.48760970505247$

$39.48760969486192$

$39.48760979510863$

$39.48760887562081$

In [11]:

/* 以後は h = 1e-5 として計算することにしよう */
thmax: 
find_root(dL(b1, th, H1, 1e-5) = 0, th, 10, 45)$

printf(true,"b = ~3,1f, H = ~3,1f のとき，", b1, H1)$
printf(true,"L が最大となる角度は ~9,6f°~%", thmax)$
printf(true,"    最大到達距離は ~9,6f", L(b1, thmax, H1))$

b = 0.5, H = 0.0 のとき，L が最大となる角度は 39.487610°
    最大到達距離は  0.679421

In [12]:

/* あらためて使う数値をまとめておく。*/
H1 = 0$
b1 = 0.5$
/* 到達距離が最大となる打ち出し角度 */
thmax: 
find_root(dL(b1, th, H1, 1e-5) = 0, th, 10, 45)$

draw2d(
  nticks = 100, 
  /* 縦横比。*/
  proportional_axes = xy,
  /* 表示範囲 */
  xrange = [0, 0.8], yrange = [0, 0.4], grid = true,
  xtics = 0.1, ytics = 0.1,
  xlabel = "X", ylabel = "Y",
  title = concat("H = ", H1, 
                 " の場合の斜方投射の最大水平到達距離"),

  key = printf(false, 
               "b =~,1f, θ_0 =~,4f°, L_{max} =~,4f", 
                b1,       thmax,       L(b1, thmax,H1)),
  parametric(X(T, b1, thmax), 
             Y(T, b1, thmax, H1), 
             T, 0, T1(b1, thmax, H1))
)$

問題 1

空気抵抗がある場合の高さ $H=0$ からの斜方投射の水平到達距離が最大となる打ち出し角度と，そのときの最大水平到達距離がわかるようなグラフを次の規格化された抵抗係数 $b$ の場合についてまとめて描け。

$$ b = 0.5, \ 1.0, \ 1.5, \ 2.0$$

解答例

In [13]:

bl: makelist(float(5*i/10), i, 1, 4);

Out[13]:

\[\tag{${\it \%o}_{32}$}\left[ 0.5 , 1.0 , 1.5 , 2.0 \right] \]

In [14]:

H1 = 0$

/* 到達距離が最大となる角度を追加するリスト */
thml:[]$

for b1 in bl do(
  thmax: find_root(dL(b1, th, H1, 1e-5) = 0, th, 20, 45),
  thml: append(thml, [thmax])
)$

/* 一時的に6桁で thml を表示する。*/
fpprintprec: 6$
thml;

Out[14]:

\[\tag{${\it \%o}_{37}$}\left[ 39.4876 , 35.5897 , 32.6181 , 30.245 \right] \]

In [15]:

/* 先に描く線をリストにしておく */
/* b1 = bl[i], thmax = thml[i] として... */
lines: makelist(
  [color=i, 
   key = printf(false, 
                "b = ~,1f,  θ_0 = ~,4f°,  L_{max} = ~,4f", 
                 bl[i],     thml[i],      L(bl[i], thml[i], H1)),
   parametric(X(T, bl[i], thml[i]), 
              Y(T, bl[i], thml[i], H1), 
              T, 0, T1(bl[i], thml[i], H1))
  ], 
  i, 1, length(bl)
)$

draw2d(
  /* 上下左右に主目盛・副目盛をつける例 */
  user_preamble = "set key samplen 1;
                   set xtics mirror; set ytics mirror;
                   set mxtics 10; set mytics 10; set grid;",
  /* 滑らかにするためにサンプリングを多めに */
  nticks = 100, 
  /* 縦横比。*/
  proportional_axes = xy,
  /* 表示範囲 */
  xrange = [0, 0.8], yrange = [0, 0.4], 
  xtics = 0.1, ytics = 0.1, 
  xlabel = "X", ylabel = "Y",
  title = "空気抵抗がある場合の斜方投射の最大水平到達距離",
  
  lines
)$

問題 2

空気抵抗 $b = 0.5$ の場合の高さ $H$ からの斜方投射の水平到達距離が最大となる打ち出し角度と，そのときの最大水平到達距離がわかるようなグラフを次の $H$ の場合についてまとめて描け。

$$ H = 0.0, 0.2, 0.4, 0.6, 0.8, 1.0$$

解答例

In [16]:

Hl: makelist(float(i/10), i, 0, 10, 2);

Out[16]:

\[\tag{${\it \%o}_{40}$}\left[ 0.0 , 0.2 , 0.4 , 0.6 , 0.8 , 1.0 \right] \]

In [17]:

b1: 0.5$

/* 到達距離が最大となる角度を追加するリスト */
thml:[]$

for H1 in Hl do(
  thmax: find_root(dL(b1, th, H1, 1e-5) = 0, th, 10, 45),
  thml: append(thml, [thmax])
)$

/* 一時的に6桁で thml を表示する。*/
fpprintprec: 6$
thml;

Out[17]:

\[\tag{${\it \%o}_{45}$}\left[ 39.4876 , 32.0487 , 26.956 , 23.1965 , 20.2851 , 17.9536 \right] \]

In [18]:

/* 先に描く線をリストにしておく */
lines: makelist(
  [color=i, 
   key = printf(false, 
                "{/=12 H = ~,1f,  θ_0 = ~,4f°,  L_{max} = ~,4f}", 
                 Hl[i],      thml[i],       L(b1, thml[i], Hl[i])),
   parametric(X(T, b1, thml[i]), 
              Y(T, b1, thml[i], Hl[i]), 
              T, 0, T1(b1, thml[i], Hl[i]))
  ], 
  i, length(Hl), 1, -1
)$

draw2d(
  /* 上下左右に主目盛・副目盛をつける例 */
  user_preamble = "set key samplen 1;
                   set xtics mirror; set ytics mirror;
                   set mxtics 4; set mytics 4; set grid;",
  /* 滑らかにするためにサンプリングを多めに */
  nticks = 100, 
  /* 縦横比。*/
  proportional_axes = xy,
  /* 表示範囲 */
  xrange = [0, 2], yrange = [0, 1.2], 
  xtics = 0.2, ytics = 0.2,
  xlabel = "X", ylabel = "Y",
  title = concat("空気抵抗 b = ", b1,
                 " の場合の斜方投射の最大水平到達距離"),
  
  line_width = 2,
  lines,file_name = mbtosha08a,terminal = 'svg
)$

Maxima で空気抵抗がある場合の斜方投射を調べる