空気抵抗がある場合の斜方投射の最大水平到達距離を調べる

「空気抵抗がある場合の斜方投射を調べる準備」を参照。

SymPy で空気抵抗がある場合の斜方投射の最大水平到達距離を調べる

空気抵抗がある場合の斜方投射の軌道を，Python の SymPy を使って（かつ SciPy や NumPy を使わずに）調べる。

詳細は以下のページ

SymPy の import

In [1]:
from sympy import *
# 1文字変数の Symbol の宣言が省略できる
from sympy.abc import *
# 円周率
from sympy import pi

init_printing()

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

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

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

In [2]:
def X(T, theta, b):
    return (1-exp(-b*T))/b * cos(theta)

def Y(T, theta, H, b):
    return H + (1-exp(-b*T))/b * sin(theta) + (1-b*T-exp(-b*T))/b**2
In [3]:
X(T, theta, b)
Out[3]:
$\displaystyle \frac{\left(1 – e^{- T b}\right) \cos{\left(\theta \right)}}{b}$
In [4]:
Y(T, theta, H, b)
Out[4]:
$\displaystyle H + \frac{\left(1 – e^{- T b}\right) \sin{\left(\theta \right)}}{b} + \frac{- T b + 1 – e^{- T b}}{b^{2}}$

滞空時間 $T_1$

滞空時間 $T_1\ \ (>0)$ は以下の式を満たす。

$$Y(T_1, \theta, H, b) = 0 \quad\Rightarrow\quad T_1 = T_1(\theta)$$

陰関数定理を使う

簡単な場合には，$T_1$ は $\theta$ の陽関数として求めることができるが，一般にはそうはいかないかもしれない。そこで，$T_1$ は $\theta$ の陰関数であるとして話を進める。

$Y(T_1(\theta), \theta) = 0$ であるから，

$$\frac{d}{d\theta} Y(T_1(\theta),\theta)=
\frac{\partial Y}{\partial T_1} \frac{d T_1}{d\theta} +
\frac{\partial Y}{\partial \theta}
= 0$$

より，$\displaystyle \frac{d T_1}{d\theta}$ を求めることができる。

In [5]:
T1 = Function('T1')(theta)

dY = diff(Y(T1, theta, H, b), theta)
sol1 = solve(dY, diff(T1, theta))
sol1
Out[5]:
$\displaystyle \left[ \frac{\left(1 – e^{b T_{1}{\left(\theta \right)}}\right) \cos{\left(\theta \right)}}{b \sin{\left(\theta \right)} – e^{b T_{1}{\left(\theta \right)}} + 1}\right]$
In [6]:
dT1 = sol1[0]
Eq(Derivative(T1, theta), dT1)
Out[6]:
$\displaystyle \frac{d}{d \theta} T_{1}{\left(\theta \right)} = \frac{\left(1 – e^{b T_{1}{\left(\theta \right)}}\right) \cos{\left(\theta \right)}}{b \sin{\left(\theta \right)} – e^{b T_{1}{\left(\theta \right)}} + 1}$

水平到達距離 $L$

滞空時間 $T_1(\theta)$ の間に水平方向に進む距離は

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

In [7]:
def L(theta):
    return X(T1, theta, b)

L(theta)
Out[7]:
$\displaystyle \frac{\left(1 – e^{- b T_{1}{\left(\theta \right)}}\right) \cos{\left(\theta \right)}}{b}$

$L$ が最大となる角度 $\theta_{\rm max}$

$\displaystyle \frac{dL}{d\theta} = 0$ となる角度 $\theta \equiv \theta_{\rm max}$ を求める。陰関数定理の結果を .subs() で代入して…

In [8]:
dL = diff(L(theta), theta).subs(diff(T1, theta), dT1)
dL = simplify(dL)
Eq(Derivative(L(theta), theta), dL)
Out[8]:
$\displaystyle \frac{\partial}{\partial \theta} \frac{\left(1 – e^{- b T_{1}{\left(\theta \right)}}\right) \cos{\left(\theta \right)}}{b} = – \frac{\left(e^{b T_{1}{\left(\theta \right)}} – 1\right) \left(b – e^{b T_{1}{\left(\theta \right)}} \sin{\left(\theta \right)} + \sin{\left(\theta \right)}\right) e^{- b T_{1}{\left(\theta \right)}}}{b \left(b \sin{\left(\theta \right)} – e^{b T_{1}{\left(\theta \right)}} + 1\right)}$

$L$ が最大となる角度は，連立方程式

\begin{eqnarray}
Y(T_1(\theta), \theta) &=& 0 \\
\frac{d}{d\theta} L(\theta) &=& 0
\end{eqnarray}

を，$\sin\theta$ と $T_1(\theta)$ について解いて求めることができる。

In [9]:
Eq(Y(T1, theta, H, b), 0)
Out[9]:
$\displaystyle H + \frac{\left(1 – e^{- b T_{1}{\left(\theta \right)}}\right) \sin{\left(\theta \right)}}{b} + \frac{- b T_{1}{\left(\theta \right)} + 1 – e^{- b T_{1}{\left(\theta \right)}}}{b^{2}} = 0$
In [10]:
Eq(dL, 0)
Out[10]:
$\displaystyle – \frac{\left(e^{b T_{1}{\left(\theta \right)}} – 1\right) \left(b – e^{b T_{1}{\left(\theta \right)}} \sin{\left(\theta \right)} + \sin{\left(\theta \right)}\right) e^{- b T_{1}{\left(\theta \right)}}}{b \left(b \sin{\left(\theta \right)} – e^{b T_{1}{\left(\theta \right)}} + 1\right)} = 0$

この連立方程式は，$T_1$ について超越方程式であり，solve() で解くことはできない。

そこでまず，$Y(T_1(\theta), \theta) = 0$ の式から $\sin\theta$ を $T_1$ で表し，$\sin\theta$ を消去して1変数 $T_1$ に関する方程式にする。

In [11]:
eqsin = solve(Y(T1, theta, H, b), sin(theta))
Eq(sin(theta), eqsin[0])
Out[11]:
$\displaystyle \sin{\left(\theta \right)} = – \frac{H b^{2} e^{b T_{1}{\left(\theta \right)}} – b T_{1}{\left(\theta \right)} e^{b T_{1}{\left(\theta \right)}} + e^{b T_{1}{\left(\theta \right)}} – 1}{b \left(e^{b T_{1}{\left(\theta \right)}} – 1\right)}$

これを $\displaystyle \frac{d}{d\theta} L(\theta) = 0$ の式に代入して，$T_1(\theta)$ に関する1本の方程式にする。

In [12]:
eq = factor(dL.subs(sin(theta), eqsin[0]))
Eq(eq, 0)
Out[12]:
$\displaystyle \frac{\left(e^{b T_{1}{\left(\theta \right)}} – 1\right)^{2} \left(H b^{2} e^{b T_{1}{\left(\theta \right)}} + b^{2} – b T_{1}{\left(\theta \right)} e^{b T_{1}{\left(\theta \right)}} + e^{b T_{1}{\left(\theta \right)}} – 1\right) e^{- 2 b T_{1}{\left(\theta \right)}}}{b^{2} \left(H b^{2} – b T_{1}{\left(\theta \right)} + e^{b T_{1}{\left(\theta \right)}} – 1\right)} = 0$
In [13]:
def eqT1(Ti, Hi, bi):
    tmp = eq.subs(T1, Ti).subs(H, Hi).subs(b, bi)
    return tmp

def T1max(H, b):
    sol = nsolve(eqT1(x, H, b), x, [0.1, 3], prec=14)
    return sol

$T_{1 \rm max}$ が数値的に求まったら，水平到達距離 $L$ を最大にする角度 $\theta_{\rm max}$ は…

In [14]:
var('T1_max theta_max')
Eq(sin(theta_max), eqsin[0].subs(T1, T1_max))
Out[14]:
$\displaystyle \sin{\left(\theta_{max} \right)} = – \frac{H b^{2} e^{T_{1 max} b} – T_{1 max} b e^{T_{1 max} b} + e^{T_{1 max} b} – 1}{b \left(e^{T_{1 max} b} – 1\right)}$

から…

In [15]:
Eq(theta_max, asin(eqsin[0]).subs(T1, T1_max))
Out[15]:
$\displaystyle \theta_{max} = – \operatorname{asin}{\left(\frac{H b^{2} e^{T_{1 max} b} – T_{1 max} b e^{T_{1 max} b} + e^{T_{1 max} b} – 1}{b \left(e^{T_{1 max} b} – 1\right)} \right)}$
In [16]:
def thetamax(Hi, bi):
    tmp = eqsin[0].subs(T1, T1max(Hi, bi)).subs(H, Hi).subs(b, bi)
    return asin(tmp)

def degrees(theta):
    from sympy import pi
    return theta*180/float(pi)
def thmax(H, b):
    ans = thetamax(H, b)
    return degrees(ans)

最大水平到達距離 $L_{\rm max}$

最大水平到達距離
$$L_{\rm max} = X(T_1(\theta_{\rm max}), \theta_{\rm max}, b)$$

は，$\theta_{\rm max}(H, b)$ であることから最終的に，$L_{\rm max}(H, b)$ となる。

In [17]:
def Lmax(H, b):
    return X(T1max(H, b), thetamax(H, b), b)

$H=0, b = 0.5$ の場合の最大水平到達距離と軌道

例として，$H=0, b=0.5$ の場合に水平到達距離 $L$ が最大となる角度 $\theta_{\rm max}$ と，そのときの最大水平到達距離 $L_{\rm max}$ を表示させてみる。

In [18]:
Htmp = 0
btmp = 0.5

print("H = %.1f, b = %.1fのとき，" % (Htmp, btmp))
print("L が最大となる角度は θmax=%.14f°" % thmax(Htmp, btmp))
print("そのときの最大水平到達距離は Lmax=%.14f" % Lmax(Htmp, btmp))

H = 0.0, b = 0.5のとき，
L が最大となる角度は θmax=39.48760970478555°
そのときの最大水平到達距離は Lmax=0.67942108743974

Maxima で空気抵抗がある場合の斜方投射の最大水平到達距離を調べる

空気抵抗がある場合の斜方投射の軌道を，Maxima を使って調べる。

詳細は以下のページ

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

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

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

In [1]:
X(T, theta, b):= (1-exp(-b*T))/b * cos(theta);
Y(T, theta, H, b):= H + (1-exp(-b*T))/b * sin(theta) + (1-b*T - exp(-b*T))/b**2;
Out[1]:
\[\tag{${\it \%o}_{1}$}X\left(T , \vartheta , b\right):=\frac{1-\exp \left(\left(-b\right)\,T\right)}{b}\,\cos \vartheta\]
Out[1]:
\[\tag{${\it \%o}_{2}$}Y\left(T , \vartheta , H , b\right):=H+\frac{1-\exp \left(\left(-b\right)\,T\right)}{b}\,\sin \vartheta+\frac{1-b\,T-\exp \left(\left(-b\right)\,T\right)}{b^2}\]

滞空時間 $T_1$

滞空時間 $T_1\ \ (>0)$ は以下の式を満たす。

$$Y(T_1, \theta, H, b) = 0 \quad\Rightarrow\quad T_1 = T_1(\theta)$$

陰関数定理を使う

簡単な場合には，$T_1$ は $\theta$ の陽関数として求めることができるが，一般にはそうはいかないかもしれない。そこで，$T_1$ は $\theta$ の陰関数であるとして話を進める。

$Y(T_1(\theta), \theta) = 0$ であるから，

$$\frac{d }{d\theta} Y(T_1(\theta), \theta)=
\frac{\partial Y}{\partial T_1} \frac{d T_1}{d\theta} +
\frac{\partial Y}{\partial \theta}
= 0$$

より，$\displaystyle \frac{d T_1}{d\theta}$ を求めることができる。

In [2]:
depends(T1, theta);

dY: diff(Y(T1, theta, H, b)=0, theta);
eqdT1: solve(dY, diff(T1, theta));
Out[2]:
\[\tag{${\it \%o}_{3}$}\left[ T_{1}\left(\vartheta\right) \right] \]
Out[2]:
\[\tag{${\it \%o}_{4}$}\frac{d}{d\,\vartheta}\,T_{1}\,e^ {- T_{1}\,b }\,\sin \vartheta+\frac{\left(1-e^ {- T_{1}\,b }\right)\,\cos \vartheta}{b}+\frac{\frac{d}{d\,\vartheta}\,T_{1}\,b\,e^ {- T_{1}\,b }-\frac{d}{d\,\vartheta}\,T_{1}\,b}{b^2}=0\]
Out[2]:
\[\tag{${\it \%o}_{5}$}\left[ \frac{d}{d\,\vartheta}\,T_{1}=-\frac{\left(e^{T_{1}\,b}-1\right)\,\cos \vartheta}{b\,\sin \vartheta-e^{T_{1}\,b}+1} \right] \]

水平到達距離 $L$

滞空時間 $T_1(\theta)$ の間に水平方向に進む距離は

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

In [3]:
L(theta, b):= X(T1, theta, b);
Out[3]:
\[\tag{${\it \%o}_{6}$}L\left(\vartheta , b\right):=X\left(T_{1} , \vartheta , b\right)\]

$L$ が最大となる角度 $\theta_{\rm max}$

$\displaystyle \frac{dL}{d\theta} = 0$ となる角度 $\theta \equiv \theta_{\rm max}$ を求める。陰関数定理の結果を ev() で代入して…

In [4]:
dL: diff(L(theta, b), theta);
dL: trigsimp(ev(dL, eqdT1));
Out[4]:
\[\tag{${\it \%o}_{7}$}\frac{d}{d\,\vartheta}\,T_{1}\,e^ {- T_{1}\,b }\,\cos \vartheta-\frac{\left(1-e^ {- T_{1}\,b }\right)\,\sin \vartheta}{b}\]
Out[4]:
\[\tag{${\it \%o}_{8}$}\frac{\left(e^{2\,T_{1}\,b}-2\,e^{T_{1}\,b}+1\right)\,\sin \vartheta-b\,e^{T_{1}\,b}+b}{b^2\,e^{T_{1}\,b}\,\sin \vartheta-b\,e^{2\,T_{1}\,b}+b\,e^{T_{1}\,b}}\]

$L$ が最大となる角度は，連立方程式

\begin{eqnarray}
Y(T_1(\theta), \theta) &=& 0 \\
\frac{d}{d\theta} L(\theta) &=& 0
\end{eqnarray}

を，$\sin\theta$ と $T_1(\theta)$ について解いて求めることができる。

In [5]:
Y(T1, theta, H, b) = 0;
dL = 0;
Out[5]:
\[\tag{${\it \%o}_{9}$}\frac{\left(1-e^ {- T_{1}\,b }\right)\,\sin \vartheta}{b}+\frac{-e^ {- T_{1}\,b }-T_{1}\,b+1}{b^2}+H=0\]
Out[5]:
\[\tag{${\it \%o}_{10}$}\frac{\left(e^{2\,T_{1}\,b}-2\,e^{T_{1}\,b}+1\right)\,\sin \vartheta-b\,e^{T_{1}\,b}+b}{b^2\,e^{T_{1}\,b}\,\sin \vartheta-b\,e^{2\,T_{1}\,b}+b\,e^{T_{1}\,b}}=0\]

この連立方程式は，$T_1$ について超越方程式であり，solve() で解くことはできない。

そこでまず，$Y(T_1(\theta), \theta) = 0$ の式から $\sin\theta$ を $T_1$ で表し，$\sin\theta$ を消去して1変数 $T_1$ に関する方程式にする。

In [6]:
eqsin: solve(Y(T1, theta, H, b) = 0, sin(theta));
Out[6]:
\[\tag{${\it \%o}_{11}$}\left[ \sin \vartheta=-\frac{\left(H\,b^2-T_{1}\,b+1\right)\,e^{T_{1}\,b}-1}{b\,e^{T_{1}\,b}-b} \right] \]

これを $\displaystyle \frac{d}{d\theta} L(\theta) = 0$ の式に代入して，$T_1(\theta)$ に関する1本の方程式にする。

In [7]:
ev(dL=0, eqsin[1]), factor;
Out[7]:
\[\tag{${\it \%o}_{12}$}\frac{e^ {- 2\,T_{1}\,b }\,\left(e^{T_{1}\,b}-1\right)^2\,\left(H\,b^2\,e^{T_{1}\,b}-T_{1}\,b\,e^{T_{1}\,b}+e^{T_{1}\,b}+b^2-1\right)}{b^2\,\left(e^{T_{1}\,b}+H\,b^2-T_{1}\,b-1\right)}=0\]

この式もまた $T_1$ に関する超越方程式であり，solve() では解けない。左辺の分子 num() の必要な部分のみを取り出し…

In [8]:
define(eqT1(H, b), 
  num(factor(ev(dL, eqsin[1])*exp(2*b*T1)/(1-exp(b*T1))**2)));
Out[8]:
\[\tag{${\it \%o}_{13}$}{\it eqT}_{1}\left(H , b\right):=H\,b^2\,e^{T_{1}\,b}-T_{1}\,b\,e^{T_{1}\,b}+e^{T_{1}\,b}+b^2-1\]

$T_1$ について find_root() で数値的に解くことにする。

$$H b^2 e^{b T_1}-T_1 b e^{b T_1}+e^{b T_1}+b^2-1 = 0$$を $T_1$ について解いた解を $T_{1 \rm{max}}$ とすると…

In [9]:
T1max(H, b):= find_root(eqT1(H, b)=0, T1, 0.1, 3);
Out[9]:
\[\tag{${\it \%o}_{14}$}{\it T1max}\left(H , b\right):={\it find\_root}\left({\it eqT}_{1}\left(H , b\right)=0 , T_{1} , 0.1 , 3\right)\]

$T_{1 \rm max}$ が数値的に求まったら，水平到達距離 $L$ を最大にする角度 $\theta_{\rm max}$ は…

In [10]:
eqsin[1];
Out[10]:
\[\tag{${\it \%o}_{15}$}\sin \vartheta=-\frac{\left(H\,b^2-T_{1}\,b+1\right)\,e^{T_{1}\,b}-1}{b\,e^{T_{1}\,b}-b}\]

上記の式から

$$\theta_{\rm max} = \arcsin\left(-\frac{\left(H\,b^2-T_{1 \rm{max}}\,b+1\right)\,e^{
T_{1 \rm{max}}\,b}-1}{b\,e^{T_{1 \rm{max}}\,b}-b}\right)$$

より求めることができる。

In [11]:
thetamax(H, b):= asin(ev(rhs(eqsin[1]), T1=T1max(H, b)));
/* 度とラジアンの変換関数を作っておく */
degrees(theta):= float(theta*180/%pi);
thmax(H, b):= degrees(thetamax(H, b));
Out[11]:
\[\tag{${\it \%o}_{16}$}{\it thetamax}\left(H , b\right):=\arcsin {\it ev}\left({\it rhs}\left({\it eqsin}_{1}\right) , T_{1}={\it T1max}\left(H , b\right)\right)\]
Out[11]:
\[\tag{${\it \%o}_{17}$}{\it degrees}\left(\vartheta\right):={\it float}\left(\frac{\vartheta\,180}{\pi}\right)\]
Out[11]:
\[\tag{${\it \%o}_{18}$}{\it thmax}\left(H , b\right):={\it degrees}\left({\it thetamax}\left(H , b\right)\right)\]

最大水平到達距離 $L_{\rm max}$

最大水平到達距離
$$L_{\rm max} = X(T_1(\theta_{\rm max}), \theta_{\rm max}, b)$$

は，$\theta_{\rm max}(H, b)$ であることから最終的に，$L_{\rm max}(H, b)$ となる。

In [12]:
Lmax(H, b):= X(T1max(H, b), thetamax(H, b), b);
Out[12]:
\[\tag{${\it \%o}_{19}$}{\it Lmax}\left(H , b\right):=X\left({\it T1max}\left(H , b\right) , {\it thetamax}\left(H , b\right) , b\right)\]

$H=0, b = 0.5$ の場合の最大水平到達距離と軌道

例として，$H=0, b=0.5$ の場合に水平到達距離 $L$ が最大となる角度 $\theta_{\rm max}$ と，そのときの最大水平到達距離 $L_{\rm max}$ を表示させてみる。

In [13]:
H: 0$
b: 0.5$

printf(true, "H = ~,f, b = ~,1fのとき，~%", H, b)$
printf(true, "L が最大となる角度は θmax=~,14f°~%", thmax(H, b))$
printf(true, "そのときの最大水平到達距離は Lmax=~,14f~%", Lmax(H, b))$

H = 0.0, b = 0.5のとき，
L が最大となる角度は θmax=39.48760970478557°
そのときの最大水平到達距離は Lmax=0.67942108743974