万有引力による三体問題のシミュレーションと新しい周期解の探索
概要
万有引力によって引き合う三つの質点に関するニュートンの運動方程式を数値的に解くことによって三体問題をシミュレーションし, 複数の新しい周期解を発見した. この成果をまとめ, 第 ◦◦ 回日本物理学会 Jr. セッションにおいて口頭発表した. また, 運動の軌跡を示すシミュレーション動画を作成し, 20 ◦◦ 年の学園祭において解説した.
背景と目的
万有引力によって引き合う二つの物体に関するニュートンの運動方程式は一般に解くことができ, 運動の軌跡は重心系において円, 楕円, 放物線, 双曲線のいずれかになる. 一方で, 万有引力によって引き合う三つの物体に関するニュートンの運動方程式 (三体問題)は, 一般に解くことができない. 初期条件を適当に選ぶと, 運動の軌跡は非常に複雑になり, 規則性のないカオス的な振る舞いを示す. しかし, 運動が周期的になるような特殊解も存在することが古くから知られている. 本活動では, ニュートンの運動方程式を数値的に解くことによって万有引力による三体問題をシミュレーションし, 新しい周期解を探索することを目的とした.
三体問題
三つの質点が同じ質量 m を持ち, それらの位置ベクトルを r1 , r2 , r3 とすると, 解くべ
き運動方程式は

Gm2
Gm2(r−r)+(r3 − r1 ),21|r2 − r1 |3|r3 − r1 |3Gm2
Gm2(r−r)+(r1 − r2 ),mr̈2 =
32|r3 − r2 |3|r1 − r2 |3
Gm2
Gm2mr̈3 =(r−r)+(r2 − r3 )13|r1 − r3 |3|r2 − r3 |3mr̈1 =となる [1]. ここで G は万有引力定数であり, 簡単のため z = 0 の xy 平面内に限られた運動だけを考えると, r1 = (x1 , y1 , 0) などとなる. また, 以下では重心系 (r1 + r2 + r3 = 0)を考え, m = G = 1 とおく.

11.5

1.5

1.0

1.00.5

0.5

0.0

0.0

-0.5

-0.5

-1.0-1.0

-1.5
-1.5

-1.0-0.5

0.0

0.5

1.0

-1.5
-1.5

1.5

図 1: カオス的な軌跡-1.0-0.5

0.0

0.5

1.0

1.5

図 2:「８の字」的な軌跡

三体問題の一般解を解析的に得ることはできないので, 運動方程式を数値的に解くことにする. 常微分方程式の数値解法として○○○法を採用し, プログラムを Python で実装
して高校の情報教室や自宅のパソコン上で実行した [2, 3]. 例えば, 初期条件を適当に

r1 = (1, 0, 0),r2 = (0, 0, 0),ṙ1 = (0.3, 0.5, 0),ṙ2 = −2ṙ1 ,r3 = (−1, 0, 0),ṙ3 = ṙ1と選ぶと, 運動の軌跡は xy 平面において図 1 のようになる. r1 , r2 , r3 の軌跡はそれぞれ赤線, 緑線, 青線で示され, 非常に複雑で規則性のないカオス的な振る舞いとなる. 初期条
件を適当に変えても, 結果は同様である.
しかし, 運動が周期的になるような特殊解も存在することが古くから知られている. 例
えば, オイラーによって発見された(√ ) )(√ )5
r1 = cos 2 t , sin 25 t , 0 ,(r2 = (0, 0, 0),((√ )(√ ) )
r3 = − cos 25 t , − sin 25 t , 0
は運動方程式の解であり, 三つの質点が一直線上に並んだまま等速円運動するので「直線2解」と呼ばれる. また, ラグランジュによって発見された( ( t )( t ) )r1 = cos 31/4, sin 31/4,0 ,( ( t)() )t2π,sin+r2 = cos 31/4 + 2π
1/4
3
3 ,0 ,3( ( t)( t) )4π
4π
r3 = cos 31/4 + 3 , sin 31/4 + 3 , 0
も運動方程式の解であり, 三つの質点が正三角形の頂点に並んだまま等速円運動するので「正三角形解」と呼ばれる. さらに, 初期条件を

r1 = (1, 0, 0),r2 = (0, 0, 0),ṙ1 = (0.347111, 0.532728, 0),r3 = (−1, 0, 0),ṙ2 = −2ṙ1 ,ṙ3 = ṙ1と選ぶと, 三つの質点が図 2 のような同一の軌道上を周期的に運動する. この解は「８の字解」と呼ばれ, 1993 年にムーアによって数値的に発見された [4].
新しい周期解本活動の目的は, 新しい周期解を探索することである. そこで文献 [5] を参考にして, 以下の方針を立てた.

1. ○○○○○
2. ○○○○○
3. ○○○○○
この方針に基づいて運動が周期的になるような初期条件を数値的に探索した結果, ◦ 個もの新しい周期解を発見することができた. 例として, 初期条件

r1 = (1, 0, 0),r2 = (0, 0, 0),ṙ1 = (0.080584, 0.588836, 0),r3 = (−1, 0, 0),ṙ2 = −2ṙ1 ,ṙ3 = ṙ1と

r1 = (1, 0, 0),r2 = (0, 0, 0),ṙ1 = (0.417343, 0.313100, 0),r3 = (−1, 0, 0),ṙ2 = −2ṙ1 ,ṙ3 = ṙ1から得られた軌道を図 3 と図 4 に示す.
以上の成果を研究レポートとしてまとめ, 第 ◦◦ 回日本物理学会 Jr. セッションに応募し, 口頭発表に採択された. また, 運動の軌跡を示すシミュレーション動画を作成し,

3

1.5

1.5

1.0

1.0

0.5

0.5

0.0

0.0

-0.5

-0.5

-1.0

-1.0

-1.5
-1.5

-1.0

-0.5

0.0

0.5

1.0

-1.5
-1.5

1.5

図 3: 新しい周期的軌道 1-1.0-0.5

0.00.51.01.5

図 4: 新しい周期的軌道 220 ◦◦ 年の学園祭において展示して来場者に解説した. 本活動は ◦ 年次の科目である「課題研究」の一環として, 担当の○○先生の指導のもとで行ったものである. 研究テーマと研究に必要な文献は○○先生から与えられたが, プログラムの実装と周期解の探索は独力で行った. 特に, ◦ ヶ月間の授業期間が終わった後も自主的に探索を続けた結果, 新しい周期解を発見することができた.*1
参考文献
三体問題とその特殊解については文献 [1] を参照した. また, 常微分方程式の数値解法については文献 [2], Python については文献 [3] を参照した.

[1] ○○○○○[2] ○○○○○[3] ○○○○○[4] C. Moore, “Braids in classical dynamics,” Phys. Rev. Lett. 70, 3675-3679 (1993).
[5] M. Šuvakov and V. Dmitrašinović, “A guide to hunting periodic three-body orbits,”
Am. J. Phys. 82, 609-619 (2014).

*1 この文書の内容はM. Šuvakov and V. Dmitrašinović, “Three classes of Newtonian three-body planar periodicorbits,” Phys. Rev. Lett. 110, 114301 (2013); http://three-body.ipb.ac.rs/を参考にして作成された例であり, 研究上の新規性はない.

4