夏休みの時期ですねー
夏休みの課題は決まりましたか?

今回は2重振り子のシミュレーションです


解説動画はこちら



2重振り子のシミュレーション


スクリーンショット 2025-07-26 16.28.53
こんな感じの2重振り子
よく有りますよね

支点があって
2個の振り子がついてて
ブルンブルンしちゃうやつです。

今回はこれをPythonで計算して
動画にするやつです。



ルンゲ・クッタ法

2重振り子をシミュレーションするには
2点の位置

2つのx,y軸の座標を求める必要があります。

この計算を行うために使用するのが
ルンゲ・クッタ法というものです。


4次のルンゲ・クッタ法(Runge-Kutta 4次法)

微分方程式の数値解法の一つで
特に4次までを考慮したテイラー展開を用いることで
より精度の高い近似解を求める方法だそうです。


この方法では、変化率の推定値を4通り計算し
それらを加重平均することで
次のステップの値を計算することになります。

計算手順は以下のようになります。

k1 → 現在の変化率
k2 → 半歩進んだ時点の変化率（k1からの推定）
k3 → さらに半歩進んだ時点の変化率（k2からの推定）
k4 → 1ステップ進んだ時点の変化率

これらを 1:2:2:1 の比率で足し合わせて
現在の状態を更新します。

𝑑𝑡 / 6 * (𝑘1 + 2*𝑘2 + 2*𝑘3 + 𝑘4) 

4次のルンゲ・クッタ法は、傾きを4回サンプリングし
テイラー展開の4次項まで考慮した近似になるため、
長時間シミュレーションでも精度が高いみたいです。

早速これをコードに落とし込んでいきましょう。


シミュレーションコード

Google colabで実行できるコードになっているので
コピペして実行することができると思います。

まずはライブラリのインポートです。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import warnings
warnings.filterwarnings('ignore')


次にパラメータ設定と
数値計算部分のコードです。
# パラメータ設定
g = 9.81   # 重力加速度
L1, L2 = 1.0, 1.0  # 振り子の長さ
m1, m2 = 1.0, 1.0  # 振り子の質量

# 初期条件（角度と角速度）
theta1 = np.pi/2
theta2 = np.pi/2 + 0.1
omega1 = 0.0
omega2 = 0.0
state = np.array([theta1, omega1, theta2, omega2])

# 時間設定
dt = 0.05
t_max = 25
steps = int(t_max/dt)

# 保存用リスト
x1_list, y1_list = [], []
x2_list, y2_list = [], []

# 二重振り子の運動方程式
def derivatives(state):
    theta1, omega1, theta2, omega2 = state
    delta = theta2 - theta1
    denom1 = (m1 + m2) * L1 - m2 * L1 * np.cos(delta)**2
    denom2 = (L2/L1) * denom1
    a1 = (m2*L1*omega1**2*np.sin(delta)*np.cos(delta) +
          m2*g*np.sin(theta2)*np.cos(delta) +
          m2*L2*omega2**2*np.sin(delta) -
          (m1+m2)*g*np.sin(theta1)) / denom1
    a2 = (-m2*L2*omega2**2*np.sin(delta)*np.cos(delta) +
          (m1+m2)*(g*np.sin(theta1)*np.cos(delta) -
                   L1*omega1**2*np.sin(delta) -
                   g*np.sin(theta2))) / denom2
    return np.array([omega1, a1, omega2, a2])

# Runge-Kutta 4次法で数値計算
for _ in range(steps):
    k1 = derivatives(state)
    k2 = derivatives(state + dt*k1/2)
    k3 = derivatives(state + dt*k2/2)
    k4 = derivatives(state + dt*k3)
    state += (dt/6)*(k1 + 2*k2 + 2*k3 + k4)
    theta1, omega1, theta2, omega2 = state

    # 座標変換
    x1 = L1 * np.sin(theta1)
    y1 = -L1 * np.cos(theta1)
    x2 = x1 + L2 * np.sin(theta2)
    y2 = y1 - L2 * np.cos(theta2)
    x1_list.append(x1)
    y1_list.append(y1)
    x2_list.append(x2)
    y2_list.append(y2)
for文のところがルンゲクッタ法の部分です。
ここでk1 - k4までを計算し
次のステップに状態を更新しています。

最後に描画用に 2点のx,y 軸としてデータ化します。

動画の生成部分です。
fig, ax = plt.subplots(figsize=(6, 6))
ax.set_xlim(-2.2, 2.2)
ax.set_ylim(-2.2, 2.2)
ax.set_aspect('equal')
ax.axis("off")
line, = ax.plot([], [], 'o-', lw=2)       # 振り子本体
trace, = ax.plot([], [], 'r-', alpha=0.5) # 軌跡

# 軌跡用リスト
trace_x, trace_y = [], []
def update(i):
    # 点の位置（始点 → 中点 → 先端）
    thisx = [0, x1_list[i], x2_list[i]]
    thisy = [0, y1_list[i], y2_list[i]]
    line.set_data(thisx, thisy)
    trace_x.append(x2_list[i])
    trace_y.append(y2_list[i])
    trace.set_data(trace_x, trace_y)
    return line, trace

ani = FuncAnimation(fig, update, frames=len(x1_list), interval=40, blit=True)
ani.save("double_pendulum.mp4", fps=30, dpi=150)
plt.close(fig)
matplotlibで描画を計算して
最終的にはmp4にします。


ファイル置き場にできた動画を見てみると
from IPython.display import Video
Video("double_pendulum.mp4", embed=True)
スクリーンショット 2025-07-26 16.50.32
こんな感じで軌跡が赤く残るような
動画が作成できます。

なかなか面白いので
長い時間で作ってみるのも
面白いかもしれません。

頑張れば振り子の数を増やしたり
色々なエフェクトを加えたりして
より面白くすることもできるかもしれません。

夏休みの課題が無くて困っているご家庭は
ぜひ面白い改造に挑戦してみてください。

今日はここまでです
それでは。

タグ ：
#プログラミング
#Python
#2重振り子
#ルンゲ・クッタ法
#シミュレーション