Pythonで描くカオスの世界：二重振り子シミュレーターを作ってみよう

「単純な仕組みなのに、予測不能な動きをする」。そんな物理現象の代表格が二重振り子です。

一重の振り子は規則正しく揺れますが、振り子の先にもう一つの振り子をつなげるだけで、その軌道は驚くほど複雑で、少しの初期値の違いが将来の挙動を劇的に変えてしまいます。これがいわゆる**「カオス理論」**の入り口です。

今回は、Pythonを使ってこの美しいカオスをシミュレーションし、ブラウザ上で自由にパラメータを操作できるツールを作成しました。

1. シミュレーターの仕組み

このプログラムは、大きく分けて3つの要素で構成されています。

① 物理演算（微分方程式の数値解法）

二重振り子の運動は、ニュートン力学（またはラグランジュ力学）から導かれる複雑な微分方程式で記述されます。コード内の deriv 関数がその心臓部です。
ここでは、scipy.integrate.odeint を使用して、刻一刻と変化する振り子の角度（）と角速度（）を計算しています。

② インタラクティブなUI

ipywidgets を活用することで、コードを書き換えなくても「棒の長さ」「重りの重さ」「初期角度」を画面上で調整できるようにしています。

• 初期角1・2: これを数度変えるだけで、数秒後の軌跡が全く別物になるのがカオスの面白いところです。
• 時間(秒): シミュレーションの長さを指定できます。

③ アニメーション描画

matplotlib.animation.FuncAnimation を使い、計算結果をパラパラ漫画のように動かします。赤い線で描かれる「軌跡（Trace）」が、カオスな動きを視覚的に強調してくれます。

---

2. 実際に遊んでみよう：カオスを体験するポイント

シミュレーターを実行したら、ぜひ以下の設定を試してみてください。

1. 高エネルギー状態: 初期角度（th1, th2）を共に 180 度（真上）に近い状態に設定してみてください。激しく回転し、予測不可能な動きが加速します。
2. 微小な差の観察: 一度シミュレーションを行い、次に初期角度を「0.1度」だけ変えて再度実行してみてください。最初は同じ動きに見えますが、途中から全く異なる軌道を描き始めます。これを「初期値鋭敏性」と呼びます。
3. 無重力実験: G（重力加速度）を 0 にするとどうなるでしょうか？慣性だけで動く不思議な挙動が観察できます。

---

ご提示いただいたコードをGoogle ColabやJupyter Notebookなどの環境で実行した際の**「出力結果（イメージ）」**と、コードのポイントを整理して出力します。

このコードは、実行すると対話的なダッシュボードが表示され、ボタンを押すことで物理シミュレーションが走る仕組みになっています。

---

1. 実行時の表示（UIレイアウト）

実行直後、以下のような入力フォームとボタンが表示されます。

カテゴリ	入力項目（デフォルト値）
長さ・質量	棒1: 1.2m / 棒2: 1.0m / 重り1: 1.0kg / 重り2: 1.0kg
環境・初期値	重力: 9.8 / 初期角1: 90.0度 / 初期角2: 45.0度 / 時間: 10秒
操作	[Simulation] (青ボタン) / [Reset] (黄ボタン)

2. シミュレーション結果（アニメーション）

[Simulation] ボタンをクリックすると、計算が走り、以下のようなアニメーションが出力されます。

• 青い線: 2本の棒と重りの現在の位置を示します。
• 赤い薄い線: 2番目の重りが通った「軌跡（Trace）」です。時間が経つにつれて、この線が複雑に絡み合い、カオスな模様を描き出します。

3. 数学的背景（微分方程式）

このプログラムが内部で解いているのは、以下のラグランジュ形式から導出される二階常微分方程式です（簡略化のため についてのみ例示）。

この複雑な計算を scipy.integrate.odeint が高速に処理しています。

---

実行用コード

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML, display
import ipywidgets as widgets

# --- 1. 物理計算ロジック ---
def deriv(y, t, L1, L2, M1, M2, G):
    theta1, z1, theta2, z2 = y
    delta = theta2 - theta1
    
    # 分母の計算
    den1 = (M1 + M2) * L1 - M2 * L1 * np.cos(delta)**2
    
    d_theta1 = z1
    d_z1 = (M2 * L1 * z1**2 * np.sin(delta) * np.cos(delta)
            + M2 * G * np.sin(theta2) * np.cos(delta)
            + M2 * L2 * z2**2 * np.sin(delta)
            - (M1 + M2) * G * np.sin(theta1)) / den1
            
    den2 = (L2 / L1) * den1
    d_theta2 = z2
    d_z2 = (-M2 * L2 * z2**2 * np.sin(delta) * np.cos(delta)
            + (M1 + M2) * G * np.sin(theta1) * np.cos(delta)
            - (M1 + M2) * L1 * z1**2 * np.sin(delta)
            - (M1 + M2) * G * np.sin(theta2)) / den2
            
    return d_theta1, d_z1, d_theta2, d_z2

# --- 2. UIの作成（単位を追加） ---
style = {'description_width': '100px'} # 幅を少し広げました
layout = widgets.Layout(width='220px')

# デフォルト値とラベルの設定
param_info = {
    'L1': {'val': 1.2, 'desc': '棒1長さ (m)'},
    'L2': {'val': 1.0, 'desc': '棒2長さ (m)'},
    'M1': {'val': 1.0, 'desc': '重り1質量 (kg)'},
    'M2': {'val': 1.0, 'desc': '重り2質量 (kg)'},
    'G':  {'val': 9.8, 'desc': '重力加速度 (m/s²)'},
    'th1':{'val': 90.0,'desc': '初期角1 (deg)'},
    'th2':{'val': 45.0,'desc': '初期角2 (deg)'},
    'dur':{'val': 10,  'desc': 'シミュ時間 (s)'}
}

inputs = {}
for k, v in param_info.items():
    if k == 'dur':
        inputs[k] = widgets.IntText(value=v['val'], description=v['desc'], style=style, layout=layout)
    else:
        inputs[k] = widgets.FloatText(value=v['val'], description=v['desc'], style=style, layout=layout)

btn_sim = widgets.Button(description="Simulation", button_style='primary', icon='play')
btn_reset = widgets.Button(description="Reset", button_style='warning', icon='undo')
output = widgets.Output()

# --- 3. ボタン動作 ---
def run_simulation(b):
    with output:
        output.clear_output(wait=True)
        print("計算中... 単位に基づいた物理演算を実行しています。")
        
        l1, l2, m1, m2 = inputs['L1'].value, inputs['L2'].value, inputs['M1'].value, inputs['M2'].value
        g, dur = inputs['G'].value, inputs['dur'].value
        t1, t2 = np.radians(inputs['th1'].value), np.radians(inputs['th2'].value)
        
        fps = 30
        t = np.linspace(0, dur, int(dur * fps))
        sol = odeint(deriv, [t1, 0, t2, 0], t, args=(l1, l2, m1, m2, g))
        
        x1 = l1 * np.sin(sol[:, 0])
        y1 = -l1 * np.cos(sol[:, 0])
        x2 = x1 + l2 * np.sin(sol[:, 2])
        y2 = y1 - l2 * np.cos(sol[:, 2])
        
        fig, ax = plt.subplots(figsize=(6, 6))
        limit = (l1 + l2) * 1.1
        ax.set_xlim(-limit, limit); ax.set_ylim(-limit, limit)
        ax.set_aspect('equal')
        ax.grid(True, linestyle='--', alpha=0.6)
        
        # グラフ軸に単位を追加
        ax.set_xlabel("Position x [m]")
        ax.set_ylabel("Position y [m]")
        ax.set_title("Double Pendulum Simulation")
        
        line, = ax.plot([], [], 'o-', lw=2, color='#1f77b4', markersize=8)
        trace, = ax.plot([], [], '-', lw=1, color='red', alpha=0.4)
        
        def update(i):
            line.set_data([0, x1[i], x2[i]], [0, y1[i], y2[i]])
            trace.set_data(x2[:i], y2[:i])
            return line, trace

        ani = FuncAnimation(fig, update, frames=len(t), interval=1000/fps, blit=True)
        plt.close()
        display(HTML(ani.to_jshtml()))

def reset_fields(b):
    for k, v in param_info.items():
        inputs[k].value = v['val']

btn_sim.on_click(run_simulation)
btn_reset.on_click(reset_fields)

# --- 4. 表示 ---
display(widgets.VBox([
    widgets.HBox([inputs['L1'], inputs['L2'], inputs['M1'], inputs['M2']]),
    widgets.HBox([inputs['G'], inputs['th1'], inputs['th2'], inputs['dur']]),
    widgets.HBox([btn_sim, btn_reset])
]), output)