gnuplotでピタゴラス三体問題 Part 1 – はじめに

gnuplot

2022.02.062022.02.09

記事内に広告、プロモーション、アフィリエイトリンクが含まれます。これらを通して購入すると売上の一部が当ブログに還元され、今後の活動の励みになります！

Hiro (@hiroloquy) です。

大学3年生のときに実験でグラフ描画ソフト”gnuplot“に触れたのをキッカケに、物理学や工学系のシミュレーションを趣味で作ってきました。これまで倒立振子や二重振り子のような複雑なシステムのプログラムを考えてきました。

今回は解析的に解くことが難しいとされる「三体問題」のシミュレーションを作ってみようと思います。三体問題の中でも3つの質点が描く軌跡が面白い「ピタゴラス三体問題」をgnuplotで考えます。

「ピタゴラス三体問題」とは？

ピタゴラス三体問題とは辺の比3：4：5の直角三角形の頂点上に質量比3：4：5の質点3つを配置して時間発展を考える問題です（質点の初速度はすべて0）。

このとき、下図のように辺の比と質量比の値が対辺・対角の関係で一致するよう配置します。

以下の2点を目標にシミュレーション作成を目指します。

数値計算と結果の描画をgnuplotのみで完結させる。
Szebehely, Peters(1967) の論文のような数値的解を再現する（精度は問わない）。具体的には、t=60 以降で質点3つのうち2つが連星となり、残り1つがエスケープするような挙動。

ピタゴラス三体問題を数値計算で解いた先行研究の論文がこちらです。

過去に一度作ろうと挑戦してみたことがあります。

当時は4次のルンゲ＝クッタ法で解きましたが、挙動は論文にあったものとは異なり、また系のエネルギーも急接近後に増えておかしな結果となりました。

gnuplotのみで3体問題を4次のルンゲクッタ使って解こうとしたけどエネルギーの総和が急接近したら増えてしまう。数値計算難しい。 pic.twitter.com/SKTCk4hTDI
— Hiro (@hiroloquy) May 31, 2020

このときは、

「微分方程式の数値解は4次のルンゲ＝クッタ法をとりあえず使えば求められる」

と安易に考えていましたが、実装ミスなのか時間刻みによる誤差の影響なのか、急接近すると結果はダメになりました。

この結果に対し、エネルギー保存性を意識した「シンプレクティック数値積分」というものをオススメしてもらったので、今回はこの手法を含め数値計算の勉強もしつつ、シミュレーションを作ろうと思います。

以下のことを書く予定にしています。記事を追加すればその都度更新していきます。