3つの質点が描く軌跡が面白い「ピタゴラス三体問題」のシミュレーションをgnuplotのみで再現したく、試行錯誤をしています。
軌跡を全て残した場合のアニメーション。t=50ぐらいから質量比3の青色質点が第2象限を通るけど先行研究にはなかった。 pic.twitter.com/YUlWzydarT
— Hiro (@hiroloquy) February 7, 2022
gnuplotのコマンドmultiplotを使い、3質点の軌跡と系のエネルギーを描くスクリプトは用意できました。そしてピタゴラス三体問題っぽい結果が得られました。
ただ、$t\geq 60$で赤色と黒色の2質点による連星、青色質点のエスケープを確認できなかったので、未完成です。
プログラムの記述、運動方程式の立式、数値計算手法の選定など、原因の候補はいくつもあるため、シミュレーションを作るまでの過程を振り返りつつ、原因を探そうと思います。
ピタゴラス三体問題の設定・初期条件
「ピタゴラス三体問題」がどんな問題なのかを確認します。
問題設定
まず、辺の比3：4：5の直角三角形の頂点上に質量比3：4：5の質点3つを配置します。質点の速度はすべてゼロとし、この状態を初期条件として時間発展を考える問題です。質点同士は万有引力が相互に作用しています。
質点の位置関係については、下図のように辺の比と質量比の値が対辺・対角の関係で一致するよう配置します。$m_i$、$l_i$は順に質量、距離を表します。
一般性を持たせるなら、距離も質量も「比」にした方が良いので、任意正数$k,\ n$を用いて
$$
\begin{align}
l_1=3k, l_2=4k, l_3=5k\\
m_1=3n, m_2=4n, m_3=5n
\end{align}
$$
とするのが望ましいですが、正直面倒であり、後の数値計算を考えると計算回数はできるだけ少なくしたいので、いずれも$k=1, \ n=1$としました。
初期条件
ここで、質量$m_i$（$i=1, 2, 3$）の質点を第$i$体と呼び、位置を$\boldsymbol{r}_i=(x_i, y_i)$とします。$m_i$は前述のとおりで、$r_i$は$l_i$の関係さえ守れればよいので、以下の位置を初期位置として考えます。
$$
\boldsymbol{r}_1=(1, 3), \boldsymbol{r}_2=(-2, -1), \boldsymbol{r}_3=(1, -1)
$$
そして速度$\boldsymbol{v}_i=d\boldsymbol{r}_i /dt=(\dot{x}_i, \dot{y}_i)$はいずれもゼロです（$t$は時刻）。
$$
\boldsymbol{v}_1=0, \boldsymbol{v}_2=0, \boldsymbol{v}_3=0
$$
（微分演算子$\mathrm{d}$は立体で記述すべきですが、texで毎回\mathrm{}と入力するのが大変なので、ここではあえて斜体で記述します。texで使うマクロ定義\newcommand{}のような機能をWordpressで使う方法をご存知の方はコメント等で教えてください。）
運動方程式および力学的エネルギーの導出
シミュレーションを作るためにも、各質点の運動方程式を求める必要があります。また、シミュレーション結果の誤差・精度等を評価するために、力学的エネルギーも考えてみます。
ニュートンの運動方程式
今回のシステムでは質点同士の万有引力が相互作用するので、第$i$体の運動方程式は次式により与えられます（$G$は万有引力定数）。
$$
m_i\frac{d^2\boldsymbol{r}_i}{dt^2}=-\sum_{j=1, j\neq i}^{3}G m_i m_j\frac{\boldsymbol{r}_i-\boldsymbol{r}_j}{\left|\boldsymbol{r}_i-\boldsymbol{r}_j\right|^3} \label{eq:eom1}
$$
三体問題なので$\Sigma$計算はありますが、立式するだけならカンタンです。
ただし、 式(\ref{eq:eom1})の左辺$\boldsymbol{r}_i$はベクトルなので、$x, y$座標それぞれで上式が成立します。また、$i=1, 2, 3$なので運動方程式の合計は、2成分×3質点=6式 になります。
数値計算に向けた運動方程式の修正
数値計算することを踏まえ、式($\ref{eq:eom1}$)に関して以下の4点の修正を加えます。
上記をもとに式($\ref{eq:eom1}$)を修正したのがこちらの2式。
$$
\begin{align}
\frac{d\boldsymbol{r}_i}{dt}&=\boldsymbol{v}_i \label{eq:eom21}\\
\frac{d\boldsymbol{v}_i}{dt}&=\sum_{j=1, j\neq i}^{3}m_j\frac{\boldsymbol{r}_j-\boldsymbol{r}_i}{\left|\boldsymbol{r}_j-\boldsymbol{r}_i\right|^3} \label{eq:eom22}
\end{align}
$$
これにより、1回の連立常微分方程式($\ref{eq:eom21}$)($\ref{eq:eom22}$)が得られ、より数値計算に適した形の式に修正することができました。
運動エネルギー / ポテンシャルエネルギー / 力学的エネルギー
数値計算で直接使うことはないですが、シミュレーション結果の妥当性の指標として、運動エネルギー、ポテンシャルエネルギー、力学的エネルギーも計算しようと思います。順に$K, U, H$とすると、それぞれ以下の式で表されます。
$$
\begin{align}
K &= \sum_{i=1}^{3}\frac{1}{2}m_i\left|\boldsymbol{v}_i\right|^2\\
U &= -\sum_{i\neq j}G\frac{m_i m_j}{\left|\boldsymbol{r}_i-\boldsymbol{r}_j\right|}\\
E &= K+U
\end{align}
$$
さいごに
続いては数値計算手法について検討します。
本文は日本語、数式はTexを用いて書いているため、読点（、）とカンマ（，）が混在していますが、ご容赦ください。
コメント