はじめに

自然科学は、膨大に積み上げられてきた観測データに対して、それらをシンプルに説明する一般則を発見する闘いであるといえる。例えば、同じ恒星系（今回の場合は地球を含んでいる太陽系）のすべての惑星が例外なくケプラーの第3法則を満たす。

自分が高校生の時、そのことを知ってとても不思議に感じた。何故なら、同じ恒星の周りを回っている惑星という条件だけで、それら全ての惑星は厳密にケプラーの第3法則を満たすからである。

さらに、ニュートンは、それら先人が見つけた法則や膨大な観測データから、万有引力についての仮説を提唱した。それは、2つの物体には双方の質量の積に比例し、距離の2乗に反比例する引力が働いているという説である。
これによって、ケプラーの第3法則は完璧に説明することができた。もちろん、ケプラーの第1,第2法則も同様である。

そこで、当時は手計算で行っていた惑星の運動のシミュレーションだが、今回は計算機を用いて当時のシミュレーションを再現する。

具体的には、恒星とその周りを回る惑星との間にはその距離の2乗に比例する万有引力が働くものと仮定する。その仮定のもとで、膨大な惑星を恒星の周りに回らせる仮想的な恒星系考える。そして、それらの公転周期と長半径（楕円軌道の長い方の半径）にはある関係性があるということを視覚的に理解すること目指す。

万有引力の法則

質量$M$の恒星と質量$m$の惑星には以下に定義される万有引力が作用する。

F(\vec{r})=G\frac{Mm}{r^2}\frac{\vec{r}}{r}

つまり、引き合う方向に、以下の大きさの引力が働く。

F(r)=G\frac{Mm}{r^2}

ただし、$G$を万有引力定数と呼ぶ。

このもとで、仮想的な恒星系を考え、その周りを回っている膨大な惑星の公転周期と長半径のデータを集めプロットしていくときに、ケプラーの第3法則が成立するかどうかを調べる。

惑星の楕円軌道について

上記の万有引力を認め、$m<<M$つまり、惑星よりも恒星の方が圧倒的に質量が大きいという条件ならば、恒星は固定されたものとみなせるので、惑星は、以下のような楕円軌道をする。

ケプラーの第3法則について

ある恒星系を考えたとき、その周りを回っている全ての惑星$(m_1,m_2,\cdot \cdot \cdot ,m_k,,\cdot \cdot \cdot ,m_n)$の公転周期$(T_1,T_2,\cdot \cdot \cdot ,T_k,\cdot \cdot \cdot ,T_n)$と長半径$(a_1,a_2,\cdot \cdot \cdot ,a_k,\cdot \cdot \cdot ,a_n)$には、以下の関係性が成立する。

T^2_{k}=ka^3_{k}

ただし、$1\le k\le n$である。

このことを以下の数値解散によって検証する。

数値計算

以下の記事のプログラムを改造する。

ただし、上記事と同様に恒星と惑星間に働く万有引力として、以下の関係性を使う。

F(r)\propto \frac{1}{r^2}

初期設定

まず、$xy$平面の$(0,0)$ 座標に恒星を配置する。そして、その周りに$n$個の惑星を配置する。ただし、初速度と初期位置は、惑星が楕円運動を行うことのできる範囲とする。

具体的には、初速度は$y$方向で$0.5\le v_k\le 0.8$の範囲でランダムな値とし、初期位置は$1\le x \le 2$で$y=0$とする。

また、楕円軌道をするある惑星の長半径は、原点と天体の最小距離と最大距離の平均である。
また、公転周期は、原点からの動径角$\theta=0$となるタイミングから求められる。

プログラム

上記の議論から、以下のようなプログラムを作成した。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#微小時間
h=1.0*10**-3
# 中心力の定数
alpha=-2.0

# 惑星の軌道計算
def planet(x,y,v_x,v_y):

    # 計算結果を入れるためのリスト
    t_ary=[]
    a_x_ary=[]
    a_y_ary=[]
    v_x_ary=[]
    v_y_ary=[]
    x_ary=[]
    y_ary=[]
    R_ary=[]
    Theta_ary=[]
    Theta_ary2=[]
    T_ary=[]
    ## 時間
    t=0
    while t<=50:
        r=math.sqrt(x**2+y**2)
        theta=np.angle(x+y*1j)
        a_x=-r**alpha*x/r
        a_y=-r**alpha*y/r
        v_x=v_x+a_x*h
        v_y=v_y+a_y*h
        x=x+v_x*h
        y=y+v_y*h
        
        t=t+h
        if theta>0 and Theta_ary[-1]<0:
            T_ary.append(t)
        t_ary.append(t)
        a_x_ary.append(a_x)
        a_y_ary.append(a_y)
        v_x_ary.append(v_x)
        v_y_ary.append(v_y)
        x_ary.append(x)
        y_ary.append(y)
        R_ary.append(r)
        Theta_ary.append(theta)
        Theta_ary2.append(theta**2)   
    XX=max(x_ary)-min(x_ary)
    YY=max(y_ary)-min(y_ary)
    RR=max(XX,YY)/2
    T_value=T_ary[0]
    
    RR3=RR**3
    T_value2=T_value**2
    return RR3, T_value2

RR3_ary=[]
T_value2_ary=[]
#Nは惑星の数
N=100
#初期位置と速度の配列を作成
x_ini_ary=np.linspace(1,2,N)
y_ini_ary=np.zeros(N)

v_x_ini_ary=np.zeros(N)
#初速度は、楕円運動になるような範囲内でランダムに生成する。
v_y_ini_ary=np.random.uniform(0.5,0.8,N)

for i in range(N):
    x_ini=x_ini_ary[i]
    y_ini=y_ini_ary[i]
    v_x_ini=v_x_ini_ary[i]
    v_y_ini=v_y_ini_ary[i]
    
    RR3, T_value2 = planet(x_ini,y_ini,v_x_ini,v_y_ini)
    RR3_ary.append(RR3)
    T_value2_ary.append(T_value2)
#グラフの描画
plt.scatter(RR3_ary,T_value2_ary,color='blue')

plt.xlabel('楕円軌道の長半径の3乗')  
plt.ylabel('公転周期の2乗')
plt.savefig('ケプラーの第3法則に関する数値計算.png')
plt.show()    

これを実行すると以下のようなグラフが出力される。

このように、惑星軌道の長半径の3乗と公転周期の2乗には厳密な比例関係が成立する。これは、非常に興味深いことであり、最初にこの関係をプロットして確かめた学者（おそらくケプラー？）は大変驚いたことであろう。これが、実測されたデータから算出されたものだとしたら、確かに当時の学者たちは科学や宗教は万能な学問であると『錯覚』してしまうだろう。

まとめ

今回は、計算機を用いて万有引力の法則を仮定することで、ケプラーの第3法則を視覚的に説明できるか調査した。結果、万有引力を仮定するだけで、見事に惑星の長半径と公転周期との関係性を説明することができた。もちろん、過去の物理学者が手計算で行っていたようなものを計算機を用いて『ズル』したのだが、そのおかげで高度な数学力や計算力が無くても上手くシミュレーションを行うことができた。なので、どのような分野でもそうだが、科学系の大学生は、専門科目以外にもプログラミングによる数値計算も学んでほしい。

参考文献