落ち続ける物体の運動

前項の補足。カテゴリー的にはmathsというよりphysicsだとは思うが。

\documentclass{jsarticle}
\begin{document}

\subsection*{万有引力下の物体の運動方程式}

質量$M$の物体との間に働く万有引力を受けて運動する質量$m$の物体について考える。
$m \ll M$であれば質量$M$の物体は動かないとみなしてよいだろう。

時刻$t$における質量$m$の物体の位置ベクトルを$\mathbf{r}_m(t)$、
質量$M$の物体の位置ベクトルを$\mathbf{r}_M(t)$、
万有引力定数を$G$とすると、
質量$m$の物体に働く万有引力の大きさは二つの物体の質量の積に比例し、
物体間の距離の二乗に反比例し、
その方向は質量$m$の物体から質量$M$の物体への向きとなるので、
\begin{eqnarray*}
m\ddot{\mathbf{r}}_m &=& \frac{GMm}{|\mathbf{r}_M-\mathbf{r}_m|^2} \frac{\mathbf{r}_M-\mathbf{r}_m}{|\mathbf{r}_M-\mathbf{r}_m|} \\
\ddot{\mathbf{r}}_m &=& \frac{GM}{|\mathbf{r}_M-\mathbf{r}_m|^3} (\mathbf{r}_M-\mathbf{r}_m)
\end{eqnarray*}
となる。
質量$M$の物体の位置ベクトルを原点とすると、
$\mathbf{r}_M=\mathbf{0}$であるから、
\[
\ddot{\mathbf{r}}_m = - \frac{GM}{|\mathbf{r}_m|^3} \mathbf{r}_m
\]
が質量$m$の物体の運動方程式になる。

$\mathbf{r}_m$を軌道面上の直交する二つの単位ベクトル$\mathbf{i}_x$、$\mathbf{i}_y$の線形結合で表現すると、
\begin{eqnarray*}
\mathbf{r}_m &=& x\mathbf{i}_x+y\mathbf{i}_y \\
|\mathbf{r}_m| &=& (x^2+y^2)^{1/2} \\
\ddot{\mathbf{r}}_m &=& \ddot{x}\mathbf{i}_x+\ddot{y}\mathbf{i}_y
\end{eqnarray*}
であるので、運動方程式は、
\begin{eqnarray*}
\ddot{x} &=& - \frac{GM}{((x^2+y^2)^{1/2})^3} x = - GMx(x^2+y^2)^{-3/2} \\
\ddot{y} &=& - \frac{GM}{((x^2+y^2)^{1/2})^3} y = - GMy(x^2+y^2)^{-3/2}
\end{eqnarray*}
のように、
$\mathbf{i}_x$と$\mathbf{i}_y$で定まる座標$x(t)$と$y(t)$の連立二階微分方程式として書くことができる。

したがって、これは4変数$x$、$y$、$\dot{x}$、$\dot{y}$の連立一階微分方程式として、
\begin{eqnarray*}
\frac{\mathrm{d}x}{\mathrm{d}t} &=& \dot{x} \\
\frac{\mathrm{d}y}{\mathrm{d}t} &=& \dot{y} \\
\frac{\mathrm{d}\dot{x}}{\mathrm{d}t} &=& - GMx(x^2+y^2)^{-3/2} \\
\frac{\mathrm{d}\dot{y}}{\mathrm{d}t} &=& - GMy(x^2+y^2)^{-3/2}
\end{eqnarray*}
と表せる。

\end{document}

最後の連立一階微分方程式の右辺それぞれがプログラム中の関数f1、f2、f3、f4になっている。
初期値は、原点からの距離が平均軌道半径の位置にあり、軌道の接線方向へ平均軌道速度の大きさの速度をもっているとした。