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)