技術の手帖

三体｜Netflixを数学的に見る

3月21日(木)からNetflixオリジナルシリーズとして長編SF小説「三体」の公開が始まりました。実はこの三体に触れるのは3度目で、1度目はテクノエッジ編集部の松尾氏からのおすすめ、2度目は同氏がパーソナリティーを務めるテック系Podcast「backspace.fm」のマネージャーからのおすすめ(正確には「三体０」の書籍)、そして今回というわけです。

先に申し上げておくと、本記事は「三体の〇〇が素晴らしくて」とか「△△は原作と違って」などを書くつもりはありません。なぜなら私は原作を読破しているわけでもなく、三体｜Netflixも1話目しか見ていないからです。書籍についてはこのところ活字離れがひどく、5ページほどで手が止まっています。

ではなぜ三体の記事を書いているのかというと、今回Netflixで公開中の1話目を見るにあたり、「そもそも題名である『三体問題』とはなんぞや？」という疑問が湧いたためです（物語よりも先に数式に興味が出てしまう短所は心得ております）。よって本記事では「三体問題を実際に解くこと」を最終目標として三体｜Netflixを数学的に予習していきます。

三体問題とは

三体問題に入る前に、先ずは万有引力について思い出してみます。

ニュートン（I. Newton）が発見した物体間に働く引力の法則。質量が\(m\)と\(M\)の2つの質点が距離\(r\)だけ離れて存在しているとき、この2つの質点には、お互いを引っ張り合う方向に

$$F=G\frac{mM}{r^2}$$

ベクトル表記では

$$\vec{F}=-G\frac{mM}{r^3}\vec{r}$$

の力が働く。

”万有引力（重力）の法則”天文学辞典：公益社団法人日本天文学会より抜粋

図として描くと以下のような状態にあるわけです。

なお、ここでは仮に質量\(m\)の質点の位置ベクトルを\(\vec{x}_1\)、質量\(M\)の質点の位置ベクトルを\(\vec{x}_2\)としています。この位置ベクトル\(\vec{x}_i\)は2次元であれば\((x_i, y_i)\)であり、3次元であれば\((x_i, y_i, z_i)\)となり同様に2つの質点間距離\(r\)もそれぞれ以下のような式となります。

$$r=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2} \tag{2D}$$
$$r=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2} \tag{3D}$$

さて、実は後述する三体問題とは「互いに重力が働く三質点系がどのような運動をするのか」を問う問題なのです。つまりここまでで述べた万有引力を用いれば少なくとも二質点の運動を考える問題、いわば「二体問題」を取り扱うことができるわけです。それぞれの質点の運動方程式は

（力）=（質量）✖️（加速度）

で求まるため、上述の万有引力、質点の質量、そして位置を時間\(t\)で2階微分した値（位置の時間変化量は速度、速度の時間変化量は加速度）を代入することで、それぞれの質点の運動方程式は以下のように書くことができます。

$$\frac{\mathrm{d}^2\vec{x}_1}{\mathrm{d}t^2}=-G\frac{M}{r^3}\vec{r} \tag{1}$$

$$\frac{\mathrm{d}^2\vec{x}_2}{\mathrm{d}t^2}=-G\frac{m}{r^3}\vec{r} \tag{2}$$

さて、同様に以下の図に従って三体問題を考える起点となる三質点の運動方程式を導出してみるとそれぞれは次のように立式できます。

$$\frac{\mathrm{d}^2\vec{x}_1}{\mathrm{d}t^2}=-G\frac{m_2}{r_{12}^3}\vec{r_{12}}-G\frac{m_3}{r_{31}^3}\vec{r_{31}} \tag{3}$$
$$\frac{\mathrm{d}^2\vec{x}_2}{\mathrm{d}t^2}=-G\frac{m_3}{r_{23}^3}\vec{r_{23}}-G\frac{m_1}{r_{12}^3}\vec{r_{12}} \tag{4}$$
$$\frac{\mathrm{d}^2\vec{x}_3}{\mathrm{d}t^2}=-G\frac{m_1}{r_{31}^3}\vec{r_{31}}-G\frac{m_2}{r_{23}^3}\vec{r_{23}} \tag{5}$$

この運動方程式を基とし、歴史に名高い数学者たちが三質点の運動の軌跡を解こうとするわけです。

三体問題は解けない？

三体問題をネットで検索するとしばしば「解ける？解けない？」という関連ワードが見られます。さてここで重要なのは「どのように解くか？」なのです。数学では「解く」ことについて3つの概念が存在します。それぞれは

  • 代数的に解く
  • 解析的に解く
  • 数値的に解く

です。もう少し具体的な問題を用いて「解く」ということを見ていきましょう。

代数的に解く

代数的解法は、方程式や不等式を解くために代数の原則や操作（加算、減算、乗算、除算、因数分解、方程式の変形など）を用いる方法です。このアプローチは、未知数を含む方程式に対して直接的な解を求めることが目的で、具体的な数値解を得ることができます。代数的方法は、特に線形方程式や二次方程式など、形式がよく定義された数学的構造を持つ問題に適しています。代数的解法を用いると、方程式が具体的な解を持つかどうか、持つとすればその解が何であるかを明確にすることができます。

例：二次方程式\(ax^2+bx+c=0\)の解の公式\(x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}\)の導出は正しく代数的に解くに該当します。

解析的に解く

解析的解法とは、方程式や系を微分方程式の形で表し、これを解いて関数の形で表される厳密な解を求める方法です。解析的解法は、物理現象や工学問題など、連続的な変化を表す問題に対してよく用いられます。例えば、物体の運動を記述する微分方程式や、流体の流れ、電磁場の分布などの物理的状況を解析する際に解析的解法が適用されます。解析的解は、問題の基本的な特性や挙動を深く理解する上で非常に有用であり、可能であれば解析的に問題を解くことが望ましいとされています。しかし、複雑な系や非線形系では、解析的な解を得ることが困難な場合が多いです。

例：ばね定数\(k\)（ばねに作用する力とばねの変形量の関係を表す定数）をもつばねの先についた質量\(m\)をもつ質点の運動方程式\(m\frac{\mathrm{d}^2x}{\mathrm{d}t^2}=-kx\)の一般解\(x(t)=Acos(\omega t+\phi)\)の導出は解析的に解くに該当します。

三体問題に対して数多の数学者たちが挑んだのはこの解析的に解くということです（なので正確には三体問題は現状、解析的に解けないと言っていただきたい）。三体問題が解析的に解けないということについての説明は少々議論が細かくなっていくため簡単に割愛するならば、「方程式が足りない」のです。運動の軌跡というのは3次元であれば質点のその瞬間の位置\((x,y,z)\)と速度\((v_x,v_y,v_z)\)です。三体問題では三質点の軌跡を求めたいので合計で18個の解を求めることとなります。そしてこの18個の解を求めるには関連する18個の式が必要になるわけです。ここまでに導出した三体問題の式は(3)~(5)の3式であり、これらに加えて例えば三質点のエネルギー保存則（1式）、運動量保存則（3式）、角運動量保存則（3式）、重心位置（1式）を式に加えてもあと7式足りません。結果現在では三体問題にいくらか制限を加えた「制限三体問題」の基に導出されたラグランジュ点などの特殊解を除いて三体問題を解析的に解き、一般解を得ることができないのです。

数値的に解く

数値的解法は、解析的に解くことが困難または不可能な問題に対して用いられるアプローチです。この方法では、方程式や系を直接解くのではなく、コンピュータを使用して問題を近似的に解きます。数値解析では、問題を離散化して小さなステップや要素に分割し、各ステップでの計算を行いながら全体の解を近似的に求めます。数値的方法には、微分方程式の数値積分、最適化問題の反復法、偏微分方程式の有限要素法など、様々な技術があります。数値解法は、複雑な工学問題、気象予測、天体物理学の問題など、実際の多くの科学技術問題に適用されています。

例：世に出回るシミュレーションという言葉はこの数値的に解くに該当します。解析的に解くとの対比として、数値的に解くとはある瞬間\(t\)における位置や速度がわかっている場合に経過時間\(\Delta t\)後に位置や速度が運動方程式等に従ってどのように変化するかを取り扱います。解析解の一般解を\(x(t)\)と書くならば数値解は常に\(x(t,\Delta t)\)のように瞬間\(t\)と経過時間\(\Delta t\)の要素が必要になります。

さて、三体問題が解析的に解けないことはわかったので最後は数値的に解いて三体問題への理解を深めていこうと思います。ここまでに導出した「二体問題」と「三体問題」をGoogle Colabolatory上でPythonにてコードに起こし、実行していきます。ただし3次元を取り扱うのは面倒なので全て2次元の制限問題として取り扱っていきます。対象は我々に最も身近な三体問題である「太陽、地球、月」とします（二体問題は「地球と月」です）。実際に取り扱ってみたい方は以下のコードをGoogle Colabolatoryにコピーして実行してみてください。

まずは二体問題から

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

# 万有引力定数
G = 6.67430e-11  # m^3 kg^-1 s^-2

# 質量
m1 = 5.972e24  # 地球の質量、kg
m2 = 7.348e22  # 月の質量、kg

# 初期位置 (地球と月の位置を仮定)
# 地球を原点とし、月をx軸上に配置
r1_initial = np.array([0, 0])  # m
r2_initial = np.array([384.4e6, 0])  # m, 地球から月までの平均距離

# 初期速度
# 地球は静止していると仮定し、月の軌道速度を設定
v1_initial = np.array([0, 0])  # m/s
v2_initial = np.array([0, 1022])  # m/s, 月の軌道速度

# 初期条件ベクトル
initial_conditions = np.concatenate((r1_initial, v1_initial, r2_initial, v2_initial))

# 時間範囲
t_span = (0, 2 * 10**6)  # 2 million seconds
t_eval = np.linspace(t_span[0], t_span[1], 500)  # 時間ステップ

def two_body_equations(t, y):
    r1 = y[:2]
    v1 = y[2:4]
    r2 = y[4:6]
    v2 = y[6:8]
    
    r = np.linalg.norm(r2 - r1)
    force = G * m1 * m2 / r**3
    
    dr1dt = v1
    dv1dt = force * (r2 - r1) / m1
    dr2dt = v2
    dv2dt = force * (r1 - r2) / m2
    
    return np.concatenate((dr1dt, dv1dt, dr2dt, dv2dt))

# 数値積分を実行
solution = solve_ivp(two_body_equations, t_span, initial_conditions, method='RK45', t_eval=t_eval)

# 解の抽出
r1_solution = solution.y[:2, :]
r2_solution = solution.y[4:6, :]

# アニメーションの作成
fig, ax = plt.subplots()
ax.set_xlim(-4.5e8, 4.5e8)
ax.set_ylim(-4.5e8, 4.5e8)
line1, = ax.plot([], [], 'o-', markersize=5, label="Earth")
line2, = ax.plot([], [], 'o-', markersize=3, label="Moon")
ax.legend()

def animate(i):
    line1.set_data([r1_solution[0, i]], [r1_solution[1, i]])
    line2.set_data([r2_solution[0, i]], [r2_solution[1, i]])
    return line1, line2,

ani = FuncAnimation(fig, animate, frames=len(t_eval), interval=20, blit=True)
gif_path = 'two-boy.gif'
ani.save(gif_path, writer='pillow', fps=30)

実行結果

月（オレンジ）が地球（青）の周りを回っている様子が描かれました。注意深く見ていただきたいのは地球の位置も若干変化している点で、単なる楕円軌道を描いているのではなく本当に数値的に解いているということがわかります。

では次にお待ちかねの三体問題です（実装される方は上述の二体問題に続いて記述してください）。

# 太陽の質量と初期条件を追加
m3 = 1.989e30  # 太陽の質量、kg
r3_initial = np.array([-1.496e11, 0])  # m, 太陽の位置 (地球からの距離を負の方向に設定)
v3_initial = np.array([0, -29.78e3])  # m/s, 太陽の軌道速度（地球と相対的に）

# 初期条件ベクトルに太陽の条件を追加
initial_conditions_3body = np.concatenate((r1_initial, v1_initial, r2_initial, v2_initial, r3_initial, v3_initial))

def three_body_equations(t, y):
    r1 = y[:2]
    v1 = y[2:4]
    r2 = y[4:6]
    v2 = y[6:8]
    r3 = y[8:10]
    v3 = y[10:12]
    
    r12 = np.linalg.norm(r2 - r1)
    r13 = np.linalg.norm(r3 - r1)
    r23 = np.linalg.norm(r3 - r2)
    
    force12 = G * m1 * m2 / r12**3
    force13 = G * m1 * m3 / r13**3
    force23 = G * m2 * m3 / r23**3
    
    dr1dt = v1
    dv1dt = force12 * (r2 - r1) / m1 + force13 * (r3 - r1) / m1
    dr2dt = v2
    dv2dt = force12 * (r1 - r2) / m2 + force23 * (r3 - r2) / m2
    dr3dt = v3
    dv3dt = force13 * (r1 - r3) / m3 + force23 * (r2 - r3) / m3
    
    return np.concatenate((dr1dt, dv1dt, dr2dt, dv2dt, dr3dt, dv3dt))

# 数値積分を実行（三体問題）
solution_3body = solve_ivp(three_body_equations, t_span, initial_conditions_3body, method='RK45', t_eval=t_eval)

# 解の抽出
r1_solution_3body = solution_3body.y[:2, :]
r2_solution_3body = solution_3body.y[4:6, :]
r3_solution_3body = solution_3body.y[8:10, :]

# アニメーションの作成（三体問題）
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))

# 全体の軌道用の設定
ax1.set_xlim(-1.5e11, 1.5e11)
ax1.set_ylim(-1.5e11, 1.5e11)
ax1.set_title("Overall View")
line1, = ax1.plot([], [], 'o-', markersize=5, label="Earth")
line2, = ax1.plot([], [], 'o-', markersize=3, label="Moon")
line3, = ax1.plot([], [], 'o-', markersize=15, label="Sun")
ax1.legend()

# 地球付近を拡大した軌道用の設定
# ax2.set_xlim(-5e10, 5e10)
# ax2.set_ylim(-5e10, 5e10)
ax2.set_title("Close-up View")
line4, = ax2.plot([], [], 'o-', markersize=15, label="Earth")
line5, = ax2.plot([], [], 'o-', markersize=10, label="Moon")
ax2.legend()

def animate_3body(i):
    # 地球と月の現在位置
    earth_x, earth_y = r1_solution_3body[0, i], r1_solution_3body[1, i]
    moon_x, moon_y = r2_solution_3body[0, i], r2_solution_3body[1, i]
    
    # 全体ビューの更新
    line1.set_data([earth_x], [earth_y])
    line2.set_data([moon_x], [moon_y])
    line3.set_data([r3_solution_3body[0, i]], [r3_solution_3body[1, i]])
    
    # Close-upビューの更新
    line4.set_data([earth_x], [earth_y])
    line5.set_data([moon_x], [moon_y])
    
    # 地球と月の最大距離に基づいて軸範囲を設定
    max_distance = np.sqrt((moon_x - earth_x)**2 + (moon_y - earth_y)**2) * 1.2  # 余裕を持たせる
    margin = max(max_distance, 1.0e7)  # 最小マージンを設定

    # 軸範囲を地球を中心に設定
    ax2.set_xlim(earth_x - margin, earth_x + margin)
    ax2.set_ylim(earth_y - margin, earth_y + margin)
    
    return line1, line2, line3, line4, line5,

ani_3body = FuncAnimation(fig, animate_3body, frames=len(t_eval), interval=20, blit=True)

# GIFとして保存
gif_path = 'three-body-close-up.gif'
ani_3body.save(gif_path, writer='pillow', fps=30)

実行結果（全体表示だと地球と月が見えにくいので拡大表示）

注意としては太陽の位置と速度は地球に対する相対的な位置と速度となっています。実際は太陽の質量が地球や月に対して圧倒的に大きいので地球や月が太陽から受ける引力の影響がとてつもなく大きいです。

ストーリーの数学的考察

ということで、ここまで長々とお付き合いいただきましたが結論としては「三体問題を解析的に解くには方程式が足りない」ことがわかりました（本結論に関しては数学者の皆様からクレームが来そうですが便宜的に取り扱わせていただきます）。このことからストーリーに何某かの影響を与えるならば、三体問題を解析的に解く上でまだ人類が発見できていない法則やこれを観測できる世界線が我々の世界、あるいは並行世界として存在し、三体の物語の中で未知の領域に触れることになるのではないか、と邪推しています。

この考察の内容を確かめるためこれから2話目を見ようと思います。

3月21日(木)からNetflixオリジナルシリーズとして長…

タグ: ,
Previous Post
Next Post

アーカイブ