はじめに

ラグランジュ点について調べていると、youtube でシミュレーションの動画を見つけました。

これを実際にやってみようと、いろいろなサイトをめぐっても、なかなか実際のシミュレーションの方法が思いつかないという方もいたのではないでしょうか(私もそうでした)。そこで、この記事でこのシミュレーションをどうやってやったのかを書き残しておこうと思います。

運動方程式を立てる

公転系の主星と伴星(例えば地球と月)のまわりに、小惑星がいる状況を考えてみましょう。そして、小惑星の運動方程式を書き下します。これが作れれば、原理的にはそれを使ってシミュレーションするだけです。(すぐあとで書きますが、実際には少し工夫してからシミュレーションしています。)

まず、座標系を設定します。重心を原点とし、 $M_2$ の天体が公転する角速度 $\omega$ で同時に回転する回転座標系を考えます。この座標系を図1に示します。つまり、この座標系では、上記の動画のように、主星と伴星は静止して見えることになります。また、主星、伴星、及び小惑星の質量を、それぞれ $M_1, M_2, m$ 、位置を $\boldsymbol{r_1}, \boldsymbol{r_2}, \boldsymbol{r}$ と置きました。また、小惑星の質量が主星、伴星の質量と比べて無視できるほど小さいとしています。また、ラグランジュ点は主星と伴星からの引力と遠心力が釣り合う点で、コリオリ力はラグランジュ点の安定性に関与しているようです[1]。

f:id:shijimi29431:20200808222548p:plain — 図1 座標系

この座標系で、運動方程式は、

$\begin{eqnarray} \label{1} \begin{split} m \ddot{\boldsymbol{r}}&=&\left(\boldsymbol{r_1}-\boldsymbol{r}\right) G \frac{M_{1} m}{\left|\boldsymbol{r_1}-\boldsymbol{r}\right|^{3}}+\left(\boldsymbol{r_2}-\boldsymbol{r}\right) G \frac{M_{2} m}{\left|\boldsymbol{r_2}-\boldsymbol{r}\right|^{3}} \\ &&+m \boldsymbol{r}|\boldsymbol{\omega}|^{2}+2 m(\boldsymbol{r} \times \boldsymbol{\omega}) \end{split} \end{eqnarray}$

となります。ここで、 $G$ は万有引力定数です。右辺1, 2項は、小惑星が主星と伴星から受ける引力、3項目は遠心力、4項目はコリオリ力です。

座標系を、扱いやすい形に変える

運動方程式はできましたが、そのままでは値が大きすぎて、私の環境の Python などではうまく計算されないようです。そこで、扱いやすいように座標を取り直すことを考えます。

座標のスケールを、主星と伴星の距離 $L$ が 1、公転の角速度の大きさが 1 であるように設定します。この座標系での万有引力定数 $G$ を計算しましょう。

まず、主星（質量 $M_1$ の天体）に関して、伴星（質量 $M_2$ の天体）による引力と、主星が受ける遠心力は釣り合っていることから、

$\begin{eqnarray} \label{2} M_{1}\left|\boldsymbol{r}_{1}\right| \underbrace{|\boldsymbol{\omega}|^{2}}_{=1}=G \frac{M_{1} M_{2}}{\underbrace{L^{2}}_{=1}}. \end{eqnarray}$

また、 $\boldsymbol{r_1}$ （ $M_1$ の天体から重心までの距離）は、 $r_1$ を原点としたときの重心の方程式より

$\begin{eqnarray*} |r_{1}|=\frac{\overbrace{L}^{=1} M_{2}}{M_{1}+M_{2}} \end{eqnarray*}$

であり、この式を式(2)に代入すると、

$\begin{eqnarray*} \cancel{M_{1}} \frac{\cancel{M_{2}}}{M_{1}+M_{2}}=G \cancel{M_{1} M_{2}} \quad \therefore G=\frac{1}{M_{1}+M_{2}} \end{eqnarray*}$

となります。

この $G$ を式(1)に代入すると、

$\begin{eqnarray*} m \ddot{\boldsymbol{r}}&=&\left(\boldsymbol{r}_{1}-\boldsymbol{r}\right) \frac{M_{1} m}{M_{1}+M_{2}}\frac{1} {\left|\boldsymbol{r}_{1}-\boldsymbol{r}\right|^{3}} +\left(\boldsymbol{r}_{2}-\boldsymbol{r}\right) \frac{M_{2} m}{M_{1}+M_{2}} \frac{1}{\left|\boldsymbol{r}_{2}-\boldsymbol{r}\right|^{3}} \\ &&+m \boldsymbol{r} \underbrace{|\boldsymbol{\omega}|^{2}}_{=1}+2 m(\dot{\boldsymbol{r}} \times \underbrace{\boldsymbol{\omega}}_{=(0,0,1)}) \\ \therefore \ddot{r}&=& \frac{r_{1}-r}{\left|r_{1}-r\right|^{3}} \frac{M_{1}}{M_{1}+M_{2}}+\frac{r_{2}-r}{\left|r_{2}-r\right|^{3}} \frac{M_{2}}{M_{1}+M_{2}}+r+2(\dot{r} \times(0,0,1)) \end{eqnarray*}$

となります。この微分方程式を用いて、シミュレーションすることができます。また、 $\boldsymbol{\omega}=(0,0,1)$ としているのは、角速度ベクトルは位置ベクトルと速度ベクトルに直交し、今は大きさ 1 の角速度を考えているためです。

ソースコードとシミュレーション結果

ソースコードは、[2] のサイトでシミュレーションに使われているものを、少し書き換えて使わせていただきました。元のサイトで紹介されているものは、上の式を使ってシミュレーションしているようです。運動方程式は、4次のルンゲ・クッタ法を用いて解かれています。この記事のコードは、小惑星の初期位置を、公転軌道上にまんべんなく配置するように与えています。

"""
Created on Tue Aug  9 23:00:14 2016

@author: Yuki Miyake

Edited on Sat Aug 8 23:11:05 2020
@shijimi29431
"""

import math
import numpy as np
import matplotlib.animation as animation
import matplotlib.pyplot as plt

M1 = 5.972*10**24# 地球の質量 kg
M2 = 7.347673*10**22# 月の質量 kg

# 質点の個数
N =500

# 座標系の選択フラグ （True: 静止座標系, False: 回転座標系）
rotate = False

mu = M2/ (M1 + M2)

x_unit = np.array([1.,0.])
r1 = - mu * x_unit # M1 の天体の座標
r2 = r1 + np.array([1.,0.]) # M2 の天体の座標

def f_gravity(z):
    """質点の位置と速度を並べた行列zの時間微分dz/dt。
      z = [
          [x1,x2,...,xn],
          [y1,y2,...,yn],
          [vx1,vx2,...,vxn],
          [vy1,vy2,...,vyn]
          ]
    """
    f = np.empty_like(z)
    f_x = f[:2]
    f_v = f[2:]
    x = z[:2]
    v = z[2:]
    rho1 = np.linalg.norm(r1[:,np.newaxis]-x,axis=0)
    rho2 = np.linalg.norm(r2[:,np.newaxis]-x,axis=0)
    f_vx = 2 * v[1] - (1-mu) * (x[0]-r1[0]) / rho1**3 - mu * (x[0]-r2[0]) / rho2**3 + x[0]
    f_vy = -2 * v[0] - (1-mu) * (x[1]-r1[1]) / rho1**3 - mu * (x[1]-r2[1]) / rho2**3 + x[1]
    f_v[0] = f_vx
    f_v[1] = f_vy
    f_x[:,:] = v
    return f              
                
def RK4(dt, x, f):
    k1 = dt * f(x)
    k2 = dt * f(x+k1/2)
    k3 = dt * f(x+k2/2)
    k4 = dt * f(x+k3)
    return (k1+2*k2+2*k3+k4)/6


def gen():
    theta = np.random.rand(N)*2*math.pi
    x = np.zeros((2,N)) +  np.vstack((np.cos(theta),np.sin(theta))) + + np.random.randn(2,N) * 0.01
    v = np.zeros_like(x)
    z = np.vstack((x,v))
    dt = 0.005
    t = 0.
    for ti in range(100000):
        t = ti * dt
        z += RK4(dt, z, f_gravity)
        if ti % 5 == 0:
            yield t, z
            
def func(data):
    t, z = data

    nmax = 2
    x = z[:2]
    R = np.array([[np.cos(t),-np.sin(t)],
                   [np.sin(t),np.cos(t)]])
    
    if rotate:
        # 回転座標にする処理
        _x = R @ x
        _r1 = R @ r1
        _r2 = R @ r2
    else:
        _x, _r1, _r2 =  x, r1, r2
    
    ax.cla()
    ax.set_aspect("equal")
    
    ax.set_title("Number of revolutions = {:.3f}".format(t/(2*math.pi)))
    ax.set_xlim(-nmax,nmax)
    ax.set_ylim(-nmax,nmax)
    
    ax.set_yticks(np.arange(-2,3))
    ax.scatter(_x[0],_x[1], s=2, c="cyan", marker=".")
    ax.scatter(_r1[0],_r1[1], c="mediumblue", s=30) 
    ax.scatter(_r2[0],_r2[1], c="orange", s=6) 
    
fig = plt.figure()
ax = fig.add_subplot(111)
ani = animation.FuncAnimation(fig, func, gen, interval=10)
plt.show()

このコードを実行して得られる結果の例を、次に示します。

f:id:shijimi29431:20200808232629g:plain — ラグランジュ点付近の小惑星の軌道

はてなブログで、なぜか長い時間の gif が投稿できなかったので、 3 回公転するまでの過程を載せました。もっと長い時間やりたい方は、ぜひ試してみてください。

まとめ

ラグランジュ点付近の小惑星の運動をシミュレーションするため、運動方程式を立ててそれを解きました。その際、座標系を扱いやすいものに取り直してから、4次のルンゲ・クッタ法で運動方程式を解きました。

参考文献

[1] LagrangePoint詳解-3

[2] ラグランジュ点の安定性を調べてみた

しじみノート

ラグランジュ点付近の物体のシミュレーション

はじめに

運動方程式を立てる

座標系を、扱いやすい形に変える

ソースコードとシミュレーション結果

まとめ

参考文献

関連リンク