2,3体問題と数値解析

はじめに

複数の惑星が存在した場合、その惑星が万有引力によりどのような軌跡を描くのかを調査する問題を多体問題という。一般的に3体問題以上の系は厳密に数式を用いて運動方程式を解くことができない。そこで、今回は2体問題や3体問題についてpythonを用いた差分法で運動方程式を解くことにより、適切な計算量以内の範囲でシミュレーションをすることを目的とする。最後に3体問題について本記事で扱った差分法よりも近似精度が上がるルンゲ・クッタ法を用いることによっても上手く描写することができるかも調査する。

最後に、以下のgif動画のように、生成AIの力を借りて、さらに効率よくルンゲクッタ法を用いて2,3体問題をシミュレーションする。

アルゴリズム

加速度と変位の関係式である運動方程式から変位を求めるために、以下のようなアルゴリズムを用いる。

速度と変位

微小区間($t=t$から$t=t+\Delta t$)に粒子の加速度$\textbf{a}$が与えられたとき
速度と変位は以下のように表すことができる。これは、微分の定義から来ている。

\textbf{v}(t+\Delta t)=\textbf{v}(t)+\textbf{a}(t)\Delta t

\textbf{x}(t+\Delta t)=\textbf{x}(t)+\textbf{v}(t)\Delta t

これを$t=0$から$t=t$まで加算することで$t=t$での位置$\textbf{x}(t)$を求めることができる。

説明

このアルゴリズムは差分法と呼ばれ、微分方程式を数値解析的に解くのに最も簡単な方法である。したがって、精度は粗くなってしまうので、$\delta t$を如何に小さくするかに懸かっている。しかし、それを行うと計算量が膨大になってしまうことから、上手くバランスを取ることがシミュレーションを適切に行う上で求められる。

2体問題

導入

以下のように、A,Bという惑星のみが存在する系を考える。

運動方程式

万有引力定数を$k$とすると運動方程式は以下のように表すことができる。

\begin{equation}
\left\{ \,
    \begin{aligned}
    & m_a \overrightarrow{a_a} = k\frac{m_a m_b}{|\overrightarrow{AB}|^2} \frac{\overrightarrow{AB}}{|\overrightarrow{AB}|} \\
    & m_b \overrightarrow{a_b} = k\frac{m_b m_a}{|\overrightarrow{BA}|^2} \frac{\overrightarrow{BA}}{|\overrightarrow{BA}|} \\
    \end{aligned}
\right.
\end{equation}

つまり、加速度と変位の関係は以下のようになる。

\begin{equation}
\left\{ \,
    \begin{aligned}
    &  \overrightarrow{a_a} = k\frac{ m_b}{|\overrightarrow{AB}|^2} \frac{\overrightarrow{AB}}{|\overrightarrow{AB}|} \\
    &  \overrightarrow{a_b} = k\frac{ m_a}{|\overrightarrow{BA}|^2} \frac{\overrightarrow{BA}}{|\overrightarrow{BA}|} \\
    \end{aligned}
\right.
\end{equation}

プログラム

上記の運動方程式に留意してプログラムを書くと以下のようになる。

2body-problem.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import japanize_matplotlib
#パラメータ
m_a=3
m_b=4
# 万有引力定数
k=1
# 初期条件
## 変位
x_a=0
x_b=1
## 速度
v_a=1j
v_b=-1j
## 加速度
a_a=0
a_b=0


# アニメを作る初期設定
fig = plt.figure()

ims = []


#初期時間
t=0
#微小時間の時間間隔(小さくするほど精度は上がるが計算時間が比例して増大する)
delta_t=0.1
while(t<30):
  im=[]
  #速度の微小変化
  v_a=v_a+a_a*delta_t
  #位置の微小変化
  x_a=x_a+v_a*delta_t

  #速度の微小変化
  v_b=v_b+a_b*delta_t
  #位置の微小変化
  x_b=x_b+v_b*delta_t

  a_a= (k*m_b/(abs(x_b-x_a))**3)*(x_b-x_a)
  a_b= (k*m_a/(abs(x_a-x_b))**3)*(x_a-x_b)
  
  im=plt.plot(x_a.real,x_a.imag,marker='o',color="red")
  im1=plt.plot(x_b.real,x_b.imag,marker='o',color="blue")
  
  ims.append(im+im1)
  #時間の更新
  t=t+delta_t
 
# 複数枚のプロットを 20ms ごとに表示
ani = animation.ArtistAnimation(fig, ims, interval=20)
#保存
ani.save("sample_2tai.gif", writer="pillow")

これを実行すると以下のようなgifファイルが生成される。

考察

上記のシミュレーション結果により、2つの天体は楕円運動をすることが分かる。これは、ケプラーの法則から上手く適切な計算精度でシミュレーションを行えているということを示唆している。ただし、天体があまりにも近づき過ぎると、位置エネルギーが発散してしまうため、力が無限大になり、加速度が無限大になり計算精度が悪化してしまうと考えられる。これは、後に述べる3体問題の計算精度がこの方法では悪くなってしまう原因であると考えられる。

3体問題

導入

以下のように、A,B,Cという惑星のみが存在する系を考える。

運動方程式

万有引力定数を$k$とすると運動方程式は以下のように表すことができる。

\begin{equation}
\left\{ \,
    \begin{aligned}
    & m_a \overrightarrow{a_a} = k\frac{m_a m_b}{|\overrightarrow{AB}|^2} \frac{\overrightarrow{AB}}{|\overrightarrow{AB}|} +k\frac{m_a m_c}{|\overrightarrow{AC}|^2} \frac{\overrightarrow{AC}}{|\overrightarrow{AC}|}\\
    & m_b \overrightarrow{a_b} = k\frac{m_b m_a}{|\overrightarrow{BA}|^2} \frac{\overrightarrow{BA}}{|\overrightarrow{BA}|} +k\frac{m_b m_c}{|\overrightarrow{BC}|^2} \frac{\overrightarrow{BC}}{|\overrightarrow{BC}|}\\
    & m_c \overrightarrow{a_c} = k\frac{m_c m_a}{|\overrightarrow{CA}|^2} \frac{\overrightarrow{CA}}{|\overrightarrow{CA}|} +k\frac{m_c m_b}{|\overrightarrow{CB}|^2} \frac{\overrightarrow{CB}}{|\overrightarrow{CB}|}\\
    \end{aligned}
\right.
\end{equation}

つまり、加速度と変位の関係は以下のようになる。

\begin{equation}
\left\{ \,
    \begin{aligned}
    & \overrightarrow{a_a} = k\frac{ m_b}{|\overrightarrow{AB}|^2} \frac{\overrightarrow{AB}}{|\overrightarrow{AB}|} +k\frac{ m_c}{|\overrightarrow{AC}|^2} \frac{\overrightarrow{AC}}{|\overrightarrow{AC}|}\\
    & \overrightarrow{a_b} = k\frac{ m_a}{|\overrightarrow{BA}|^2} \frac{\overrightarrow{BA}}{|\overrightarrow{BA}|} +k\frac{ m_c}{|\overrightarrow{BC}|^2} \frac{\overrightarrow{BC}}{|\overrightarrow{BC}|}\\
    &  \overrightarrow{a_c} = k\frac{ m_a}{|\overrightarrow{CA}|^2} \frac{\overrightarrow{CA}}{|\overrightarrow{CA}|} +k\frac{ m_b}{|\overrightarrow{CB}|^2} \frac{\overrightarrow{CB}}{|\overrightarrow{CB}|}\\
    \end{aligned}
\right.
\end{equation}

プログラム

上記の運動方程式に留意してプログラムを書くと以下のようになる。

3body-problem.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import japanize_matplotlib
#パラメータ
## 質量
m_a=1
m_b=1
m_c=1
# 万有引力定数
k=1
# 初期条件
## 変位
x_a=0
x_b=1
x_c=0.5+(0.5*3**0.5*1j)
## 速度
v_a=1j
v_b=-1j
v_c=-1
## 加速度
a_a=0
a_b=0
a_c=0

# アニメを作る初期設定
fig = plt.figure()

ims = []


#初期時間
t=0
#微小時間の時間間隔(小さくするほど精度は上がるが計算時間が比例して増大する)
delta_t=0.001
while(t<0.8):
  im=[]
  #速度の微小変化
  v_a=v_a+a_a*delta_t
  #位置の微小変化
  x_a=x_a+v_a*delta_t

  #速度の微小変化
  v_b=v_b+a_b*delta_t
  #位置の微小変化
  x_b=x_b+v_b*delta_t

  #速度の微小変化
  v_c=v_c+a_c*delta_t
  #位置の微小変化
  x_c=x_c+v_c*delta_t

  a_a= (k*m_b/(abs(x_b-x_a))**3)*(x_b-x_a)+(k*m_c/(abs(x_c-x_a))**3)*(x_c-x_a)
  a_b= (k*m_c/(abs(x_c-x_b))**3)*(x_c-x_b)+(k*m_a/(abs(x_a-x_b))**3)*(x_a-x_b)
  a_c= (k*m_a/(abs(x_a-x_c))**3)*(x_a-x_c)+(k*m_b/(abs(x_b-x_c))**3)*(x_b-x_c)

  im=plt.plot(x_a.real,x_a.imag,marker='o',color="red")
  im1=plt.plot(x_b.real,x_b.imag,marker='o',color="blue")
  im2=plt.plot(x_c.real,x_c.imag,marker='o',color="black")

  ims.append(im+im1+im2)
  #時間の更新
  t=t+delta_t
 
# 複数枚のプロットを 100ms ごとに表示
ani = animation.ArtistAnimation(fig, ims, interval=20)
#保存
ani.save("sample_3tai3.gif", writer="pillow")

これを実行すると以下のようなgifファイルが生成される。

考察

3体問題の場合、特に惑星同士が接近する場合において高い計算精度がもとめられる。したがって、今回のシミュレーションでは、あえて短時間の範囲で行った。なので、長時間のシミュレーションを行う場合は、より精度が高い微分方程式を解くための数値解析手法を用いるべきであろう。

ルンゲ・クッタ法

ルンゲ・クッタ法とは、微分方程式を積算する数値解析法の一種で計算量と計算精度をバランス良く追求できる方法である。

以下の記事を参考に、本記事では4次のルンゲ・クッタ法を用いて数値解析を行うものとする。

python 3body-problem2.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import japanize_matplotlib


#パラメータ
## 質量
m_a=1
m_b=1
m_c=1
# 万有引力定数
k=1
# 初期条件
## 変位
x_a=0
x_b=1
x_c=0.5+(0.5*3**0.5*1j)
## 速度
v_a=1j
v_b=-1j
v_c=-1
## 加速度
a_a=0
a_b=0
a_c=0

# アニメを作る初期設定
fig = plt.figure()

ims = []


#微小時間の時間間隔(小さくするほど精度は上がるが計算時間が比例して増大する)
delta_t=1.0*10**-3


def F1_a(x_a,x_b,x_c,v_a,v_b,v_c, t):
    return v_a
def F2_a(x_a,x_b,x_c,v_a,v_b,v_c, t):
    return (k*m_b/(abs(x_b-x_a))**3)*(x_b-x_a)+(k*m_c/(abs(x_c-x_a))**3)*(x_c-x_a)
def F1_b(x_a,x_b,x_c,v_a,v_b,v_c, t):
    return v_b
def F2_b(x_a,x_b,x_c,v_a,v_b,v_c, t):
    return (k*m_c/(abs(x_c-x_b))**3)*(x_c-x_b)+(k*m_a/(abs(x_a-x_b))**3)*(x_a-x_b)
def F1_c(x_a,x_b,x_c,v_a,v_b,v_c, t):
    return v_c
def F2_c(x_a,x_b,x_c,v_a,v_b,v_c, t):
    return (k*m_a/(abs(x_a-x_c))**3)*(x_a-x_c)+(k*m_b/(abs(x_b-x_c))**3)*(x_b-x_c)



t=0
while(t<1):
  im=[]

  #print(x_b,x_c)
  #Runge-Kutta法(a)
  k1_a = delta_t*F1_a(x_a,x_b,x_c,v_a,v_b,v_c, t)
  m1_a = delta_t*F2_a(x_a,x_b,x_c,v_a,v_b,v_c, t)
  #k2 = delta_t*F1_a(x+0.5*k1, v+0.5*m1, t+0.5*delta_t)
  k2_a = delta_t*F1_a(x_a+0.5*k1_a,x_b,x_c,v_a+0.5*m1_a,v_b,v_c, t+0.5*delta_t)
  m2_a = delta_t*F2_a(x_a+0.5*k1_a,x_b,x_c,v_a+0.5*m1_a,v_b,v_c, t+0.5*delta_t)
  k3_a = delta_t*F1_a(x_a+0.5*k2_a,x_b,x_c,v_a+0.5*m2_a,v_b,v_c, t+0.5*delta_t)
  m3_a = delta_t*F2_a(x_a+0.5*k2_a,x_b,x_c,v_a+0.5*m2_a,v_b,v_c, t+0.5*delta_t)
  k4_a = delta_t*F1_a(x_a+k3_a,x_b,x_c,v_a+m3_a,v_b,v_c, t+delta_t)
  m4_a = delta_t*F2_a(x_a+k3_a,x_b,x_c,v_a+m3_a,v_b,v_c, t+delta_t)
  x_a = x_a + (k1_a + 2*k2_a + 2*k3_a + k4_a)/6
  v_a = v_a + (m1_a + 2*m2_a + 2*m3_a + m4_a)/6

  #Runge-Kutta法(b)
  k1_b = delta_t*F1_b(x_a,x_b,x_c,v_a,v_b,v_c, t)
  m1_b = delta_t*F2_b(x_a,x_b,x_c,v_a,v_b,v_c, t)
  #k2 = delta_t*F1_a(x+0.5*k1, v+0.5*m1, t+0.5*delta_t)
  k2_b = delta_t*F1_b(x_a,x_b+0.5*k1_b,x_c,v_a,v_b+0.5*m1_b,v_c, t+0.5*delta_t)
  m2_b = delta_t*F2_b(x_a,x_b+0.5*k1_b,x_c,v_a,v_b+0.5*m1_b,v_c, t+0.5*delta_t)
  k3_b = delta_t*F1_b(x_a,x_b+0.5*k2_b,x_c,v_a,v_b+0.5*m2_b,v_c, t+0.5*delta_t)
  m3_b = delta_t*F2_b(x_a+0.5*k2_b,x_b,x_c,v_a+0.5*m2_b,v_b,v_c, t+0.5*delta_t)
  k4_b = delta_t*F1_b(x_a,x_b+k3_b,x_c,v_a,v_b+m3_b,v_c, t+0.5*delta_t)
  m4_b = delta_t*F2_b(x_a,x_b+k3_b,x_c,v_a,v_b+m3_b,v_c, t+0.5*delta_t)
  x_b = x_b + (k1_b + 2*k2_b + 2*k3_b + k4_b)/6
  v_b = v_b + (m1_b + 2*m2_b + 2*m3_b + m4_b)/6


  #Runge-Kutta法(c)
  k1_c = delta_t*F1_c(x_a,x_b,x_c,v_a,v_b,v_c, t)
  m1_c = delta_t*F2_c(x_a,x_b,x_c,v_a,v_b,v_c, t)
  #k2 = delta_t*F1_a(x+0.5*k1, v+0.5*m1, t+0.5*delta_t)
  k2_c = delta_t*F1_c(x_a,x_b,x_c+0.5*k1_c,v_a,v_b,v_c+0.5*m1_c, t+0.5*delta_t)
  m2_c = delta_t*F2_c(x_a,x_b,x_c+0.5*k1_c,v_a,v_b,v_c+0.5*m1_c, t+0.5*delta_t)
  k3_c = delta_t*F1_c(x_a,x_b,x_c+0.5*k2_c,v_a,v_b,v_c+0.5*m2_c, t+0.5*delta_t)
  m3_c = delta_t*F2_c(x_a,x_b,x_c+0.5*k2_c,v_a,v_b,v_c+0.5*m2_c, t+0.5*delta_t)
  k4_c = delta_t*F1_c(x_a,x_b,x_c+k3_c,v_a,v_b,v_c+m3_c, t+0.5*delta_t)
  m4_c = delta_t*F2_c(x_a,x_b,x_c+k3_c,v_a,v_b,v_c+m3_c, t+0.5*delta_t)
  x_c = x_c + (k1_c + 2*k2_c + 2*k3_c + k4_c)/6
  v_c = v_c + (m1_c + 2*m2_c + 2*m3_c + m4_c)/6


  im=plt.plot(x_a.real,x_a.imag,marker='o',color="red")
  im1=plt.plot(x_b.real,x_b.imag,marker='o',color="blue")
  im2=plt.plot(x_c.real,x_c.imag,marker='o',color="black")

  ims.append(im+im1+im2)
  #時間の更新
  t=t+delta_t
 
# 複数枚のプロットを 2ms ごとに表示
ani = animation.ArtistAnimation(fig, ims, interval=2)
#保存
ani.save("sample_3tai3_ran3.gif", writer="pillow")

これを実行すると以下のようになる。

ルンゲ・クッタ法を用いても簡易的な差分法とあまり差が見られなかった。従って、3体問題を扱うときは、物体が近くなったときのみ時間の粒度を細かくするといった手法を持ちた方がよいようである。

参考

上記の考察をもとにして、生成AI(claude)で2,3体問題の数値計算を行うプログラムを作成した。予想の斜め上の成果物で感動している。

python 2_3body_scipy.py

"""
2体・3体問題 — scipy.integrate.solve_ivp による高精度シミュレーション
元記事: https://qiita.com/arairuca/items/1393c2cfcf8e9a1b552e

改善点:
  1. scipy の DOP853 (8次 Runge-Kutta) + 自動ステップ幅制御で積分
  2. 全天体の状態ベクトルを同時に積分 (連成ODE を正しく扱う)
  3. 軌跡の trailing trail 表示でアニメーションを視認しやすく改善
  4. japanize_matplotlib 非依存 (フォント設定で日本語表示)
"""

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import matplotlib.animation as animation
from scipy.integrate import solve_ivp

# ─── 共通パラメータ ───────────────────────────────────────────
G = 1.0   # 万有引力定数（無次元）

# ─── ヘルパー: 複素数 ↔ 実数ベクトル 変換 ────────────────────
def c2r(z):
    """複素スカラー → [Re, Im]"""
    return np.array([z.real, z.imag])

def r2c(v):
    """[Re, Im] → 複素スカラー"""
    return v[0] + 1j * v[1]


# ════════════════════════════════════════════════════════════════
# 1.  2 体 問 題
# ════════════════════════════════════════════════════════════════

def make_2body_ode(m_a, m_b, G=G):
    """
    状態ベクトル y = [xa, ya, xb, yb, vxa, vya, vxb, vyb]
    dydt を返す関数を生成する
    """
    def ode(t, y):
        xa, ya, xb, yb, vxa, vya, vxb, vyb = y
        # 相対ベクトル a→b
        dx = xb - xa;  dy = yb - ya
        r3 = (dx**2 + dy**2)**1.5
        # 加速度
        axa =  G * m_b * dx / r3;  aya =  G * m_b * dy / r3
        axb = -G * m_a * dx / r3;  ayb = -G * m_a * dy / r3
        return [vxa, vya, vxb, vyb, axa, aya, axb, ayb]
    return ode

def run_2body():
    m_a, m_b = 3.0, 4.0
    # 初期条件 (元記事と同じ、複素数表現 → 実数ベクトルに変換)
    x_a0 = c2r(0+0j);     v_a0 = c2r(0+1j)
    x_b0 = c2r(1+0j);     v_b0 = c2r(0-1j)
    y0 = np.concatenate([x_a0, x_b0, v_a0, v_b0])

    t_end = 30.0
    # solve_ivp: method='DOP853' は8次RK + 自動ステップ制御
    sol = solve_ivp(
        make_2body_ode(m_a, m_b),
        [0, t_end],
        y0,
        method='DOP853',
        rtol=1e-10,   # 相対許容誤差
        atol=1e-12,   # 絶対許容誤差
        dense_output=True   # 補間解を保持
    )
    # アニメーション用に等間隔サンプリング
    n_frames = 300
    t_eval = np.linspace(0, t_end, n_frames)
    Y = sol.sol(t_eval)    # shape (8, n_frames)

    # ── アニメーション描画 ──────────────────────────────────
    fig, ax = plt.subplots(figsize=(6, 6))
    trail_len = 60   # 軌跡の長さ (フレーム数)

    def update(i):
        ax.clear()
        lo = max(0, i - trail_len)
        ax.plot(Y[0, lo:i+1], Y[1, lo:i+1], '-', color='tomato',   alpha=0.5, lw=1)
        ax.plot(Y[2, lo:i+1], Y[3, lo:i+1], '-', color='royalblue', alpha=0.5, lw=1)
        ax.plot(Y[0, i], Y[1, i], 'o', color='tomato',   ms=8, label=f'A (m={m_a})')
        ax.plot(Y[2, i], Y[3, i], 'o', color='royalblue', ms=8, label=f'B (m={m_b})')
        ax.set_xlim(-3, 3); ax.set_ylim(-3, 3)
        ax.set_aspect('equal')
        ax.legend(loc='upper right', fontsize=8)
        ax.set_title(f'2体問題  t = {t_eval[i]:.2f}  (DOP853)', fontsize=11)
        ax.set_xlabel('x'); ax.set_ylabel('y')

    ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=40)
    ani.save('2body_scipy.gif', writer='pillow', fps=25)
    plt.close(fig)
    print('2体問題: 2body_scipy.gif を保存しました')


# ════════════════════════════════════════════════════════════════
# 2.  3 体 問 題
# ════════════════════════════════════════════════════════════════

def make_3body_ode(m_a, m_b, m_c, G=G):
    """
    状態ベクトル y = [xa, ya, xb, yb, xc, yc,
                      vxa, vya, vxb, vyb, vxc, vyc]
    ※ 全天体を「同時に」評価するので連成系が正しく解ける
    """
    def ode(t, y):
        xa, ya, xb, yb, xc, yc = y[0:6]
        vxa, vya, vxb, vyb, vxc, vyc = y[6:12]

        def grav(xi, yi, xj, yj, mj):
            dx = xj - xi;  dy = yj - yi
            r3 = (dx**2 + dy**2)**1.5
            return G * mj * dx / r3, G * mj * dy / r3

        axa, aya = np.add(grav(xa, ya, xb, yb, m_b), grav(xa, ya, xc, yc, m_c))
        axb, ayb = np.add(grav(xb, yb, xa, ya, m_a), grav(xb, yb, xc, yc, m_c))
        axc, ayc = np.add(grav(xc, yc, xa, ya, m_a), grav(xc, yc, xb, yb, m_b))

        return [vxa, vya, vxb, vyb, vxc, vyc,
                axa, aya, axb, ayb, axc, ayc]
    return ode

def run_3body():
    m_a = m_b = m_c = 1.0
    # 初期条件 (元記事と同じ)
    x_a0 = c2r(0+0j);            v_a0 = c2r(0+1j)
    x_b0 = c2r(1+0j);            v_b0 = c2r(0-1j)
    x_c0 = c2r(0.5 + 0.5*3**0.5*1j);  v_c0 = c2r(-1+0j)
    y0 = np.concatenate([x_a0, x_b0, x_c0, v_a0, v_b0, v_c0])

    t_end = 3.0   # 元記事より長い時間でも安定して計算できる
    sol = solve_ivp(
        make_3body_ode(m_a, m_b, m_c),
        [0, t_end],
        y0,
        method='DOP853',
        rtol=1e-11,
        atol=1e-13,
        dense_output=True
    )

    n_frames = 400
    t_eval = np.linspace(0, t_end, n_frames)
    Y = sol.sol(t_eval)   # shape (12, n_frames)

    fig, ax = plt.subplots(figsize=(6, 6))
    trail_len = 80
    colors = ['tomato', 'royalblue', 'seagreen']
    labels = ['A', 'B', 'C']
    # 位置インデックス: (0,1), (2,3), (4,5)
    pos_idx = [(0, 1), (2, 3), (4, 5)]

    def update(i):
        ax.clear()
        lo = max(0, i - trail_len)
        for (xi, yi), c, lb in zip(pos_idx, colors, labels):
            ax.plot(Y[xi, lo:i+1], Y[yi, lo:i+1], '-', color=c, alpha=0.4, lw=1)
            ax.plot(Y[xi, i],      Y[yi, i],      'o', color=c, ms=8, label=lb)
        ax.set_xlim(-2.5, 2.5); ax.set_ylim(-2.5, 2.5)
        ax.set_aspect('equal')
        ax.legend(loc='upper right', fontsize=8)
        ax.set_title(f'3体問題  t = {t_eval[i]:.3f}  (DOP853)', fontsize=11)
        ax.set_xlabel('x'); ax.set_ylabel('y')

    ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=40)
    ani.save('3body_scipy.gif', writer='pillow', fps=25)
    plt.close(fig)
    print('3体問題: 3body_scipy.gif を保存しました')


# ════════════════════════════════════════════════════════════════
if __name__ == '__main__':
    print('=== 2体問題を計算中 ===')
    run_2body()
    print('=== 3体問題を計算中 ===')
    run_3body()
    print('完了')

やはり、生成AIの技術は人知を超えた速さである。

まとめ

今回は、差分法を用いた数値解析手法を用いて、2体問題や3体問題をシミュレーションすることを目的とした。結果、2体問題では、惑星は楕円運動をした。一方で3体問題では複雑な運動をした。ただし、今回のような問題は高い精度と膨大な計算量が必要になるので、アルゴリズムや性能の高い計算機を使うことが、より良い精度のシミュレーションを行うには必要であると考えられる。