空気抵抗がある場合の斜方投射の解
空気抵抗がある場合の無次元化された斜方投射の解は，
\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}
/* 度をラジアンに変換する関数を準備しておく */
radians(deg):= deg * %pi/180$
X(T, th, b):= (1 - exp(-b*T))/b * cos(radians(th))$
Y(T, th, H, b):= H + (1 - exp(-b*T))/b * sin(radians(th))
+ (1 - b*T - exp(-b*T))/(b**2)$
滞空時間 $T_1$
$T = 0$ で打ち出された粒子が地面 $Y=0$ に落ちるまでの滞空時間 $T_1$ は以下を満たす。
\begin{eqnarray}
Y(T_1, \theta, H, b) &=& H + \frac{1-e^{-b T_1}}{b} \, \sin\theta + \frac{1 – b T_1 -e^{-b T_1}}{b^2} = 0
\end{eqnarray}
この式は $T_1$ について超越方程式であり，解析的に解を求めることはできない。
そこで Maxima で方程式の数値解を求めるために，find_root() を使う。
/* どの範囲で find_root するか。ここでは 0.1 ~ 3.0 と広めに */
T1(th, H, b):= find_root(Y(T, th, H, b)=0, T, 0.1, 3.0)$
水平到達距離 $L$
$Y(T_1, \theta, H, b) =0$ の数値解 $T_1(\theta, H, b)$ が求まったら，この滞空時間の間の水平方向の到達距離 $L$ は，
$$L = X(T_1(\theta, H, b), \theta, b)$$
L(th, H, b):= X(T1(th, H, b), th, b)$
$H = 0$ の場合の軌道
空気抵抗の係数 $b=0.5$ の場合に，$\theta = 45^{\circ}$ 前後の角度での斜方投射の軌道をグラフにしてみる。
とりあえずは $\theta = 45^{\circ}$ の場合を1本，グラフにしてみる。
th1: 45$
H1: 0$
b1: 0.5$
draw2d(
parametric(X(T, th1, b1), Y(T, th1, H1, b1),
T, 0, T1(th1, H1, b1))
)$
オプションを設定してそれらしいグラフに。
th1: 45$
H1: 0$
b1: 0.5$
draw2d(
/* 縦横比。*/
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,
/* 線の太さ */
line_width = 2,
/* 数値を取り込んだタイトルの設定 */
title = printf(false, "b = ~,1f, H = ~,1f の場合の斜方投射", b1, H1),
/* 数値を取り込んだ凡例の設定 */
key = printf(false, "θ = ~2d°", th1),
/* 媒介変数表示で曲線を描く */
parametric(X(T, th1, b1), Y(T, th1, H1, b1),
T, 0, T1(th1, H1, b1))
)$
$\theta$ の値を $38^{\circ}$ から $45^{\circ}$ まで $1^{\circ}$ 刻みで変えて，複数の曲線を一度に描く例。
あらかじめ makelist() でグラフにするリストを作っておきます。。
なお，細かいことですが，グラフの最高点は $\theta$ の値が大きいときですので，先に $45^{\circ}$ から描きはじめて，$1^{\circ}$ ずつ減らしてみます。
H1: 0$
b1: 0.5$
/* 複数のグラフを一度に描く場合は... */
/* 最初に makelist で作っておくといいかも */
thmax: 45$
thmin: 38$
lines: makelist(
[key = printf(false, "θ = ~2d°", th),
color = mod(thmax, th),
parametric(X(T, th, b1), Y(T, th, H1, b1),
T, 0, T1(th, H1, b1))
],
th, thmax, thmin, -1)$
draw2d(
/* 縦横比。*/
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 = printf(false, "b = ~,1f, H = ~,1f の場合の斜方投射", b1, H1),
/* 線の数だけ羅列する必要なく，作成済みのリストを書くだけ */
lines
)$
落下点付近を拡大表示します。
/* 着地点付近を拡大表示 */
draw2d(
/* 落下点付近を拡大表示 */
xrange = [0.66, 0.682], yrange = [0, 0.01],
grid = true,
xlabel = "X", ylabel = "Y",
title = printf(false, "b = ~,1f, H = ~,1f の場合の斜方投射", b1, H1),
lines
)$
$b = 0.5$ の場合，最大水平到達距離となるのは $\theta$ が 39° から 40° くらいの時であることがわかります。これをもう少し定量的に求めてみます。
水平到達距離 $L$ のグラフ
定量的に調べるために，水平到達距離 L(b, th, H) をグラフにしてみる。
draw2d(
grid = true,
xlabel = "打ち出し角度 θ (°)",
ylabel = "規格化された水平到達距離 L",
title = printf(false,
"b = ~,1f, H = ~,1f の場合の投射角度と水平到達距離", b1, H1),
line_width = 2,
explicit(L(th, H1, b1), th, 38, 42)
)$
最大水平到達距離となる角度 $\theta_{\rm max}$ を数値微分で求める
$L(\theta, H, b)$ が最大となるのは $\displaystyle \frac{d L}{d \theta} = 0$ となる角度のときです。
$L(\theta, H, b)$ は解析的な関数として与えられていないので，ここでは以下のような数値微分で求めることにします。
$$\frac{d L}{d \theta} \simeq \frac{L(\theta+\delta, H, b)-L(\theta-\delta, H, b)}{2 \delta}$$
$\displaystyle \frac{d L}{d \theta} = 0$ となる角度 $\theta$ を方程式の解を数値的に解く find_root() を使って求めます。
dL(th, H, b, delta):=
(L(th+delta, H, b) - L(th-delta, H, b))/(2*delta)$
/* 方程式の数値解 */
/* h を変えてみる */
for delta in [1.e-3, 1.e-4, 1.e-5, 1.e-6, 1.e-7] do(
print(find_root(dL(th, H1, b1, delta) = 0, th, 38, 45))
)$
上の結果から，delta = 1e-4 にして以下のような関数を定義します。
thmax(H, b):= find_root(dL(th, H, b, 1e-4)=0, th, 10, 45)$
$H=0, b = 0.5$ の場合
最大水平到達距離と投射角度
H1: 0$
b1: 0.5$
printf(true,"b = ~,1f, H = ~,1f のとき，~%", b1, H1)$
printf(true,"L が最大となる角度 θmax は ~9,6f°~%", thmax(H1, b1))$
printf(true,"そのときの最大水平到達距離 Lmax は ~,6f", L(thmax(H1,b1),H1,b1))$
軌道のグラフ
H1: 0$
b1: 0.5$
thMax: thmax(H1, b1)$
draw2d(
/* 縦横比。*/
proportional_axes = xy,
/* 表示範囲 */
xrange = [0, 0.8], yrange = [0, 0.4], grid = true,
xtics = 0.1, ytics = 0.1,
xlabel = "X", ylabel = "Y",
title = printf(false, "H = ~1d の場合の斜方投射の最大水平到達距離",H1),
line_width = 2,
key = printf(false,
"b = ~,1f, θ_{max} = ~,4f°, L_{max} = ~,4f",
b1, thMax, L(thMax,H1, b1)),
parametric(X(T, thMax, b1), Y(T, thMax, H1, b1),
T, 0, T1(thMax, H1, b1))
)$
問題 1
空気抵抗がある場合の高さ $H=0$ からの斜方投射の水平到達距離が最大となる打ち出し角度と，そのときの最大水平到達距離がわかるようなグラフを次の規格化された抵抗係数 $b$ の場合についてまとめて描け。
$$ b = 0.5, \ 1.0, \ 1.5, \ 2.0$$
解答例
bl: makelist(float(5*i/10), i, 1, 4);
/* 一時的に6桁で thml を表示する。*/
fpprintprec: 6$
thml: makelist(
thmax(H1, bl[i]), i, 1, length(bl)
);
H1: 0$
/* 先に描く線をリストにしておく */
lines: makelist(
[color=i,
key = printf(false,
"b = ~,1f, θ_{max} = ~,4f°, L_{max} = ~,4f",
bl[i], thml[i], L(thml[i],H1,bl[i])),
parametric(X(T, thml[i], bl[i]), Y(T, thml[i], H1, bl[i]),
T, 0, T1(thml[i], H1, bl[i]))
],
i, 1, length(bl)
)$
draw2d(
/* 縦横比。*/
proportional_axes = xy,
/* 表示範囲 */
xrange = [0, 0.8], yrange = [0, 0.4], grid = true,
xtics = 0.1, ytics = 0.1,
xlabel = "X", ylabel = "Y",
title = printf(false, "H = ~1d の場合の斜方投射の最大水平到達距離",H1),
line_width = 2,
lines
)$
問題 2
問題 1 で，$H = 0$ の場合，空気抵抗がないとき，つまり $b=0$ のときには $\theta=45^{\circ}$ のときに最大水平到達距離になるが，このままでは b=0 の場合のグラフを追加しようとするとエラーになるはずである。
その原因をつきとめ，改善し，$b = 0, 0.5, 1.0, 1.5$ の場合のグラフをまとめて描け。
問題 3
空気抵抗 $b = 0.5$ の場合の高さ $H$ からの斜方投射の水平到達距離が最大となる打ち出し角度と，そのときの最大水平到達距離がわかるようなグラフを次の $H$ の場合についてまとめて描け。
$$ H = 0.0, 0.2, 0.4, 0.6, 0.8, 1.0$$
解答例
Hl: makelist(float(i/10), i, 0, 10, 2);
b1: 0.5$
/* 一時的に6桁で thml を表示する。*/
fpprintprec: 6$
thml: makelist(
thmax(Hl[i], b1), i, 1, length(Hl)
);
/* 先に描く線をリストにしておく */
lines: makelist(
[color=i,
key = printf(false,
"H = ~,1f, θ_{max} = ~,4f°, L_{max} = ~,4f",
Hl[i], thml[i], L(thml[i],Hl[i],b1)),
parametric(X(T, thml[i], b1), Y(T, thml[i], Hl[i], b1),
T, 0, T1(thml[i], Hl[i], b1))
],
i, length(Hl), 1, -1
)$
draw2d(
/* 縦横比。*/
proportional_axes = xy,
/* 表示範囲 */
xrange = [0, 2], yrange = [0, 1.2], grid = true,
xtics = 0.2, ytics = 0.2,
xlabel = "X", ylabel = "Y",
title = concat("空気抵抗 b = ", b1,
" の場合の斜方投射の最大水平到達距離"),
line_width = 2,
lines
,file_name = Mbtosha08,terminal = 'svg
)$
