gnuplotでピタゴラス三体問題 Part 2 – モデリング・運動方程式

3つの質点が描く軌跡が面白い「ピタゴラス三体問題」のシミュレーションをgnuplotのみで再現したく、試行錯誤をしています。

軌跡を全て残した場合のアニメーション。t=50ぐらいから質量比3の青色質点が第2象限を通るけど先行研究にはなかった。 pic.twitter.com/YUlWzydarT

— Hiro (@hiroloquy) February 7, 2022

gnuplotのコマンドmultiplotを使い、3質点の軌跡と系のエネルギーを描くスクリプトは用意できました。そしてピタゴラス三体問題っぽい・・結果が得られました。

ただ、$t\ge 60$で赤色と黒色の2質点による連星、青色質点のエスケープを確認できなかったので、未完成です。

プログラムの記述、運動方程式の立式、数値計算手法の選定など、原因の候補はいくつもあるため、シミュレーションを作るまでの過程を振り返りつつ、原因を探そうと思います。

gnuplotでピタゴラス三体問題 - 記事まとめ

解析的に解くことが難しいとされる「三体問題」の中でも、3つの質点が描く軌跡が面白い「ピタゴラス三体問題」のシミュレーショ...

hiroloquy.com

2022.02.09

ピタゴラス三体問題の設定・初期条件

「ピタゴラス三体問題」がどんな問題なのかを確認します。

問題設定

まず、辺の比3：4：5の直角三角形の頂点上に質量比3：4：5の質点3つを配置します。質点の速度はすべてゼロとし、この状態を初期条件として時間発展を考える問題です。質点同士は万有引力が相互に作用しています。

質点の位置関係については、下図のように辺の比と質量比の値が対辺・対角の関係で一致するよう配置します。$m_i$、$l_i$は順に質量、距離を表します。

一般性を持たせるなら、距離も質量も「比」にした方が良いので、任意正数$k,n$を用いて

$$\begin{array}{}\text{(1)}& l_1=3k,l_2=4k,l_3=5k\text{(2)}& m_1=3n,m_2=4n,m_3=5n\end{array}$$

とするのが望ましいですが、正直面倒であり、後の数値計算を考えると計算回数はできるだけ少なくしたいので、いずれも$k=1,n=1$としました。

初期条件

ここで、質量$m_i$（$i=1,2,3$）の質点を第$i$体と呼び、位置を${\mathit{r}}_i=(x_i,y_i)$とします。$m_i$は前述のとおりで、$r_i$は$l_i$の関係さえ守れればよいので、以下の位置を初期位置として考えます。

$$\begin{array}{}\text{(3)}& {\mathit{r}}_1=(1,3),{\mathit{r}}_2=(-2,-1),{\mathit{r}}_3=(1,-1)\end{array}$$

そして速度${\mathit{v}}_i=d{\mathit{r}}_i/dt=({\dot{x}}_i,{\dot{y}}_i)$はいずれもゼロです（$t$は時刻）。

$$\begin{array}{}\text{(4)}& {\mathit{v}}_1=0,{\mathit{v}}_2=0,{\mathit{v}}_3=0\end{array}$$

（微分演算子$\mathrm{d}$は立体で記述すべきですが、texで毎回\mathrm{}と入力するのが大変なので、ここではあえて斜体で記述します。texで使うマクロ定義\newcommand{}のような機能をWordpressで使う方法をご存知の方はコメント等で教えてください。）

運動方程式および力学的エネルギーの導出

シミュレーションを作るためにも、各質点の運動方程式を求める必要があります。また、シミュレーション結果の誤差・精度等を評価するために、力学的エネルギーも考えてみます。

ニュートンの運動方程式

今回のシステムでは質点同士の万有引力が相互作用するので、第$i$体の運動方程式は次式により与えられます（$G$は万有引力定数）。

$$\begin{array}{}\text{(5)}& m_i\frac{d^2{\mathit{r}}_i}{d{t}^2}=-\sum _{j=1,j\ne i}^{3}G{m}_i{m}_j\frac{{\mathit{r}}_i-{\mathit{r}}_j}{{|{\mathit{r}}_i-{\mathit{r}}_j|}^3}\end{array}$$

三体問題なので$\mathrm{\Sigma }$計算はありますが、立式するだけならカンタンです。

ただし、 式($\text{5}$)の左辺${\mathit{r}}_i$はベクトルなので、$x,y$座標それぞれで上式が成立します。また、$i=1,2,3$なので運動方程式の合計は、2成分×3質点=6式 になります。

数値計算に向けた運動方程式の修正

数値計算することを踏まえ、式($\text{5}$)に関して以下の4点の修正を加えます。

上記をもとに式($\text{5}$)を修正したのがこちらの2式。

$$\begin{array}{}\text{(6)}& \frac{d{\mathit{r}}_i}{dt}& ={\mathit{v}}_i\text{(7)}& \frac{d{\mathit{v}}_i}{dt}& =\sum _{j=1,j\ne i}^3{m}_j\frac{{\mathit{r}}_j-{\mathit{r}}_i}{{|{\mathit{r}}_j-{\mathit{r}}_i|}^3}\end{array}$$

これにより、1回の連立常微分方程式($\text{6}$)($\text{7}$)が得られ、より数値計算に適した形の式に修正することができました。

運動エネルギー / ポテンシャルエネルギー / 力学的エネルギー

数値計算で直接使うことはないですが、シミュレーション結果の妥当性の指標として、運動エネルギー、ポテンシャルエネルギー、力学的エネルギーも計算しようと思います。順に$K,U,H$とすると、それぞれ以下の式で表されます。

$$\begin{array}{}\text{(8)}& K& =\sum _{i=1}^3\frac{1}{2}{m}_i{\left|{\mathit{v}}_i\right|}^2\text{(9)}& U& =-\sum _{i\ne j}G\frac{m_i{m}_j}{|{\mathit{r}}_i-{\mathit{r}}_j|}\text{(10)}& E& =K+U\end{array}$$

さいごに

続いては数値計算手法について検討します。

本文は日本語、数式はTexを用いて書いているため、読点（、）とカンマ（，）が混在していますが、ご容赦ください。