万有引力の数値計算

地球上からロケットを発射すると、ロケットは地球の引力の影響を受けて運動する。 宇宙空間は3次元空間だが、地球の引力による運動では2次元平面上で考えることができる。

プログラム rocket_orbit.rb は、地球上から発射した推進力のないロケット（SpinLaunch）がどのように運動するかをシミュレーションする。

プログラムは rbCanvas で動作する。 rbCanvas editor に入力して[▶]をクリックしてみよう。

rocket_orbit.rb

# coding: utf-8
# ロケットの運動を描画する
# 動作環境：rbCanvas
# 作者：山本裕樹
 
# 画像サイズ
WIDTH = 800
HEIGHT = 800
 
# 画像上の中心座標
IOX = WIDTH / 2
IOY = HEIGHT / 2
 
# 空間のスケール [pix/m]
SCALE = 5e-6
 
# 物理定数
G = 6.67430e-11                 # 万有引力定数 [m^3/kg/s^2]
ME = 5.9724e24                  # 地球の質量 [kg]
RE = 6.3781e6                   # 地球の半径 [m]
 
# 経過時間 [s]
t = 0
 
# $dt：時間の刻み [s]
$dt = 60
 
# ロケットの初期値
# (x0, y0)：ロケットの発射位置 [m]
x0 = RE + 50
y0 = 0
 
# (x, y)：ロケットの空間座標 [m]
x = x0
y = y0
 
# (u, v)：ロケットの速度 [m/s]
u = 0
v = 0
 
# ロケットの軌跡
tr = []
trmax = 1000
 
# 原点から (x, y) までの距離
def r(x, y) return Math.sqrt(x ** 2 + y ** 2) end
 
# 運動方程式の右辺1
def f(x, y) return - G * ME * x / r(x, y) ** 3 end
 
# 運動方程式の右辺2
def g(x, y) return - G * ME * y / r(x, y) ** 3 end
 
# ku, kv, kx, ky の計算
def calc_k(u, v, x, y, du, dv, dx, dy)
  return f(x + dx, y + dy) * $dt,
         g(x + dx, y + dy) * $dt,
         (u + du) * $dt,
         (v + dv) * $dt
end
 
# ルンゲ=クッタ法で $dt 後の (x,y) と (u,v) を求める
def next_step(u, v, x, y)
  ku1, kv1, kx1, ky1 = calc_k(u, v, x, y, 0, 0, 0, 0)
  ku2, kv2, kx2, ky2 = calc_k(u, v, x, y,
                              ku1 * 0.5, kv1 * 0.5, kx1 * 0.5, ky1 * 0.5)
  ku3, kv3, kx3, ky3 = calc_k(u, v, x, y,
                              ku2 * 0.5, kv2 * 0.5, kx2 * 0.5, ky2 * 0.5)
  ku4, kv4, kx4, ky4 = calc_k(u, v, x, y, ku3, kv3, kx3, ky3)
  return u + (ku1 + 2 * ku2 + 2 * ku3 + ku4) / 6.0,
         v + (kv1 + 2 * kv2 + 2 * kv3 + kv4) / 6.0,
         x + (kx1 + 2 * kx2 + 2 * kx3 + kx4) / 6.0,
         y + (ky1 + 2 * ky2 + 2 * ky3 + ky4) / 6.0
end
 
# 空間座標 (x,y) から画像座標 (ix,iy) へ変換
def x2ix(x) return (IOX + x * SCALE).round end
def y2iy(y) return (IOY - y * SCALE).round end
 
# 背景の設定
Window.width   = WIDTH
Window.height  = HEIGHT
Window.bgcolor = 'black'
 
# 発射角度のスライダーとラベル
label1_angle = UI.label(WIDTH + 20, 120, 100, 30, '発射角度')
label2_angle = UI.label(WIDTH + 370, 120, 180, 30, '')
slider_angle = UI.slider(WIDTH + 150, 120, 200, [-90, 90, 30, 1])
 
# 発射速度のスライダーとラベル
label1_speed = UI.label(WIDTH + 20, 170, 100, 30, '発射速度')
label2_speed = UI.label(WIDTH + 370, 170, 180, 30, '')
slider_speed = UI.slider(WIDTH + 150, 170, 200, [5, 15, 10, 0.1])
 
# 発射ボタンとラベル
$launch = false
button_launch = UI.button(WIDTH + 20, 220, 100, 30, '発射')
button_launch.on_action do
  $launch = true
end
 
flight = false
 
Window.loop do
 
  # 発射時
  if $launch
    # 発射角度 [rad]
    alpha = slider_angle.value.to_f * Math::PI / 180
    # 発射速度 [m/s]
    speed = slider_speed.value.to_f * 1e3
 
    u = speed * Math.cos(alpha)
    v = speed * Math.sin(alpha)
    x = x0
    y = y0
    t = 0
    tr = []
    $launch = false
    flight = true
  end
 
  # ロケットの空間座標から画像座標への変換
  ix = x2ix(x)
  iy = y2iy(y)
 
  # 飛行中
  if flight
    # 軌跡の記録
    if tr.length > trmax
      tr.shift
    end
    tr.push([ix, iy])
  end
 
  # ロケットの描画
  Window.draw_circle_fill(ix, iy, 3, 'green')
 
  # 軌跡の描画
  (1..(tr.length - 1)).each do |i|
    Window.draw_dot(tr[i - 1][0], tr[i - 1][1], 'cyan', weight: 1)
  end
 
  # 地球の描画
  Window.draw_circle_fill(x2ix(0), y2iy(0), RE * SCALE, 'blue')
 
  # ラベルの更新
  label2_angle.value = "#{slider_angle.value}度"
  label2_speed.value = "#{slider_speed.value}km/s"
 
  # 経過時間の表示
  Window.draw_font(0, 25, sprintf("経過時間:%d秒", t),
                   Font.default, color: 'yellow')
 
  if flight
    # 次のステップ
    u, v, x, y = next_step(u, v, x, y)
    t += $dt
 
    # 地球に衝突ならストップ
    if r(x, y) < RE
      flight =false
    end
  end
end

ロケットの発射速度と発射角度を色々変えてみよう。

地球の中心を原点としたロケットの座標を $(x,y)$ とする。 ロケットは地球の引力を受けて運動するため、$(x,y)$ は時間 $t$ で変化する。

ロケットの速度を $(v_{x},v_{y})$ とすると、$(x,y)$ と以下の関係にある。 \begin{eqnarray*} \frac{dx}{dt}&=&v_{x}\\ \frac{dy}{dt}&=&v_{y} \end{eqnarray*} ここで $\frac{dx}{dt}$ と $\frac{dy}{dt}$ は $x$ と $y$ をそれぞれ時間 $t$ で微分することを意味する。

ロケットは地球の引力を受けて運動するため、万有引力の法則から次の運動方程式が成り立つ。 \begin{eqnarray*} m\frac{dv_{x}}{dt}&=&-\frac{GM_{\mathrm{E}}m}{r^{3}}x\\ m\frac{dv_{y}}{dt}&=&-\frac{GM_{\mathrm{E}}m}{r^{3}}y \end{eqnarray*} ここで $G$ は万有引力定数、$M_{\mathrm{E}}$ は地球の質量、$m$ はロケットの質量である。 \[ G=6.67430\times 10^{-11}\mathrm{[m^{3}kg^{-1}s^{-2}]} \] \[ M_{\mathrm{E}}=5.9724\times 10^{24}\mathrm{[kg]} \] 両辺を $m$ で割ると $m$ は消えるため、この計算にロケットの質量は不要である。 \begin{eqnarray*} \frac{dv_{x}}{dt}&=&-\frac{GM_{\mathrm{E}}}{r^{3}}x\\ \frac{dv_{y}}{dt}&=&-\frac{GM_{\mathrm{E}}}{r^{3}}y \end{eqnarray*} $r$ は地球の中心からロケットまでの距離で、三平方の定理（ピタゴラスの定理）から次の式で求められる。 \[ r=\sqrt{x^{2}+y^{2}} \]

以上の式を解くとロケットの座標 $(x,y)$ が時間 $t$ とともにどのように変化するか、ロケットの軌道を求めることができる。

以下の4つの方程式（微分方程式という）を解くとロケットの運動の軌道が求められる。 \begin{eqnarray*} \frac{dx}{dt}&=&v_{x}\\ \frac{dy}{dt}&=&v_{y}\\ \frac{dv_{x}}{dt}&=&f(x,y)\\ \frac{dv_{y}}{dt}&=&g(x,y) \end{eqnarray*} ここで \begin{eqnarray*} f(x,y)&=&-\frac{GM_{\mathrm{E}}}{r^{3}}x\\ g(x,y)&=&-\frac{GM_{\mathrm{E}}}{r^{3}}y \end{eqnarray*} とした。

微分方程式はオイラー法やルンゲ=クッタ法などを用いて数値的に計算できる。