AAMの速度変化（解説）

読む人はあまりいないと思いますが書いてみます。SLGとは全く関係ありません。

質量mの物体の位置が時間tの関数として\(x(t)\)に沿って運動するとき、ニュートンの運動方程式に従います。ニュートンの運動方程式とは、加速度をaとし、物体にかかる力をFとして、\(ma(t)=F(t)\)です。aもFも時間の関数としておきます。

加速度は速度の時間微分で、速度は位置の時間微分です。つまり
\[a=\frac{d^2x}{dt^2}(t)\]
から
\[m\frac{d^2x}{dt^2}(t)=F(t)\]
です。

ミサイルの場合、F(t)は推力と空気抵抗ということにします。推力が一定値 F であるとして、空気抵抗は速度の二乗に比例するものとすれば
\[m\frac{d^2x}{dt^2}(t)=F-r v(t)^2\]
です。rは比例定数です。rとFとmはミサイルの種類によってそれぞれ定まるものです。

さて、\( v(t)=dx/dt \)ですから左辺を vで書き直すと
\[m\frac{dv}{dt}(t)=F-r v(t)^2\]
vだけの式になります。

右辺はvと定数だけの式だけですから、変数分離形の微分方程式です。

ニュートンの凄みは、物体の運動を位置の関数の微分方程式で書けることに気がついて、微分積分学を作り上げてしまったところにあります。

さて、簡単のためにm=F=r=1としてv(t)=vと書いて
\[\frac{dv}{dt}=1- v^2\]
この微分方程式を解きます。両辺を\(1-v^2\)で割りましょう
\[ \frac{1}{1-v^2}\frac{dv}{dt}=1\]

両辺をtで積分すると、Cを積分定数として
\[ \int \frac{1}{1-v^2}\frac{dv}{dt}dt=t+C\]

置換積分の公式から、左辺はvでの積分になります。
\[ \int \frac{1}{1-v^2}dv=t+C\]

左辺の積分には、次の関係式を使います。
\[\frac{2}{1-v^2}=\frac{1}{1-v}+\frac{1}{1+v}\]

これから、左辺の積分は次の通りに計算できます。
\[ \frac{1}{2}(-\log|1-v| + \log|1+v| )=\frac{1}{2}\log|\frac{1+v}{1-v}|\]

結局、ニュートン方程式から得られた微分方程式から次の関係式を得られました。
\[ \frac{1}{2}\log|\frac{1+v}{1-v}|=t+C\]

vについて整理すると、
\[ \log|\frac{1+v}{1-v}|=2(t+C)\]
\[ \frac{1+v}{1-v}=e^{2(t+C)}\]

\(D=e^{2C}\)とおいて、
 \[ \frac{1+v}{1-v}=De^{2t}\]

結局、次のようにv=v(t)を得られました。
\[ v(t)=\frac{De^{2t}-1}{De^{2t}+1}\]

AAMの速度変化

いりやっくさんのMAの日記（MA会員のみ）から。AAMの速度変化はポアソン分布の密度関数のような感じになるようです。じゃあということでニュートンの運動方程式を立ててみます。(符号等が間違っていたので諸処修正)

初速というか母機の速度を\(v_0\)とします。高校の物理で習ったように、高速で運動する物体の空気抵抗は速度の二乗に比例するとします。質量も空気抵抗の比例定数も正規化して1としましょう。発射時を時刻0, 推力\(a\gt 0\)は時刻\(t_0\)まで一定で、その後は0になるものとします。

速度を\(v(t)\)で表すと、時刻0から\(t_0\)までは次のニュートンの運動方程式が成立します。

\[
\frac{dv}{dt}(t)=a-v(t)^2
\]

変数分離形なので解くことにします。まず
\[
\frac{1}{2}\log\frac{\sqrt{a}+v(t)}{\sqrt{a}-v(t)}+C=t
\]
（絶対値は省略。最後に符号を調整。）

整理して、
 \[
v(t)=\sqrt{a}\frac{De^{e^{2t}}+1}{De^{2t}-1}
\]
ここで、\(D\)は\(v_0\)から次のように定まります。
\[
D=\frac{\sqrt{a}+v_0}{\sqrt{a}-v_0}
\]
\(t\)が正負の無限大で\(\pm\sqrt{a}\)に収束します。これが推力\(a\)で飛翔し続けた時の限界速度です。が、今は\(t_0\)で推力が止まるのです。

\(t_0\)以降は次のようになります。
\[
\frac{dv}{dt}(t)=-v(t)^2
\]
\[
\frac{1}{v(t)}-\frac{1}{v(t_0)} =t-t_0
\]
これより、
\[
v(t)=\frac{v(t_0)}{v(t_0)(t-t_0)+1}
\]

\(t_0\)では滑らかな形になりませんが、推力を細かく与えると解けるかどうかわからなくなりますのでここまで。

a=1, t0=2, v0=1/2 でグラフを描くとこんな感じです。t=t0で推力が不連続に変わりますので、グラフは滑らかになりません。x軸は0から10までに1/100します。これを描くのに使ったMATLABクローンのフリーソフトoctaveのソースをつけます。

a=1;
v0=1/2;
t0=2;
function y=v(x,a,v0,t0)
d=(sqrt(a)+v0)/(sqrt(a)-v0);
if( x < t0 )
y=(d*exp(2*x)-sqrt(a))/(d*exp(2*x)+sqrt(a));
else
y=1/(x-t0+1/((sqrt(a)*(d*exp(2*t0)-1))/(d*exp(2*t0)+1)));
endif
endfunction

z=[1:1000];

for x=[1:1000];
z(x)=v(x/100,1,1/2,2);
endfor

plot(z)