魅力的な二重回転振子のシミュレーションを探る

はじめに

複振り子は、物理学と数学における古典的な問題です。互いに取り付けられた 2 つの振り子が含まれ、2 番目の振り子の運動は 1 番目の振り子の運動に影響を受けます。この実験では、Python と Matplotlib を使って複振り子の運動をシミュレートします。

VM のヒント

VM の起動が完了したら、左上隅をクリックしてノートブックタブに切り替え、Jupyter Notebook を使って練習しましょう。

時々、Jupyter Notebook が読み込み終わるまで数秒待つ必要があります。Jupyter Notebook の制限により、操作の検証は自動化できません。

学習中に問題に遭遇した場合は、Labby にお問い合わせください。セッション後にフィードバックを提供してください。すぐに問題を解決いたします。

必要なライブラリをインポートする

まず、シミュレーションに必要なライブラリをインポートします。

from collections import deque
import matplotlib.pyplot as plt
import numpy as np
from numpy import cos, sin
import matplotlib.animation as animation

パラメータを設定する

次に、シミュレーション用のパラメータを定義します。これらのパラメータには、重力加速度、各振り子の長さと質量、およびシミュレーションの時間間隔が含まれます。

G = 9.8  ## 重力加速度、m/s^2
L1 = 1.0  ## 振り子 1 の長さ、m
L2 = 1.0  ## 振り子 2 の長さ、m
L = L1 + L2  ## 合成振り子の最大長さ
M1 = 1.0  ## 振り子 1 の質量、kg
M2 = 1.0  ## 振り子 2 の質量、kg
t_stop = 2.5  ## シミュレートする秒数
history_len = 500  ## 表示する軌道点の数

微分関数を定義する

任意の時刻におけるシステムの微分を計算する関数を定義します。この関数は、システムの現在の状態（各振り子の角度と角速度）を入力として受け取り、それらの値の微分を返します。

def derivs(t, state):
    dydx = np.zeros_like(state)

    dydx[0] = state[1]

    delta = state[2] - state[0]
    den1 = (M1+M2) * L1 - M2 * L1 * cos(delta) * cos(delta)
    dydx[1] = ((M2 * L1 * state[1] * state[1] * sin(delta) * cos(delta)
                + M2 * G * sin(state[2]) * cos(delta)
                + M2 * L2 * state[3] * state[3] * sin(delta)
                - (M1+M2) * G * sin(state[0]))
               / den1)

    dydx[2] = state[3]

    den2 = (L2/L1) * den1
    dydx[3] = ((- M2 * L2 * state[3] * state[3] * sin(delta) * cos(delta)
                + (M1+M2) * G * sin(state[0]) * cos(delta)
                - (M1+M2) * L1 * state[1] * state[1] * sin(delta)
                - (M1+M2) * G * sin(state[2]))
               / den2)

    return dydx

初期条件を設定して常微分方程式を積分する

次に、シミュレーションの初期条件を設定します。各振り子の初期角度と角速度、およびシミュレーションの時間間隔を設定します。その後、オイラー法を使って常微分方程式を積分し、各時刻での各振り子の位置と速度を求めます。

## create a time array from 0..t_stop sampled at 0.02 second steps
dt = 0.01
t = np.arange(0, t_stop, dt)

## th1 and th2 are the initial angles (degrees)
## w10 and w20 are the initial angular velocities (degrees per second)
th1 = 120.0
w1 = 0.0
th2 = -10.0
w2 = 0.0

## initial state
state = np.radians([th1, w1, th2, w2])

## integrate the ODE using Euler's method
y = np.empty((len(t), 4))
y[0] = state
for i in range(1, len(t)):
    y[i] = y[i - 1] + derivs(t[i - 1], y[i - 1]) * dt

各振り子の位置を計算する

次に、各時刻での各振り子の位置と速度を使って、各振り子の x 座標と y 座標を計算します。

x1 = L1*sin(y[:, 0])
y1 = -L1*cos(y[:, 0])

x2 = L2*sin(y[:, 2]) + x1
y2 = -L2*cos(y[:, 2]) + y1

グラフを設定する

次に、シミュレーション用のグラフを設定します。振り子の最大長に等しい x と y の範囲を持つグラフを作成し、アスペクト比を等しく設定し、グリッドを追加します。

fig = plt.figure(figsize=(5, 4))
ax = fig.add_subplot(autoscale_on=False, xlim=(-L, L), ylim=(-L, 1.))
ax.set_aspect('equal')
ax.grid()

アニメーション関数を定義する

二重回転振子の動きをアニメーション化する関数を定義します。この関数は、現在のタイムステップを受け取り、それを使って各振子の位置を更新します。また、2 番目の振子の軌跡も更新します。

line, = ax.plot([], [], 'o-', lw=2)
trace, = ax.plot([], [], '.-', lw=1, ms=2)
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
history_x, history_y = deque(maxlen=history_len), deque(maxlen=history_len)

def animate(i):
    thisx = [0, x1[i], x2[i]]
    thisy = [0, y1[i], y2[i]]

    if i == 0:
        history_x.clear()
        history_y.clear()

    history_x.appendleft(thisx[2])
    history_y.appendleft(thisy[2])

    line.set_data(thisx, thisy)
    trace.set_data(history_x, history_y)
    time_text.set_text(time_template % (i*dt))
    return line, trace, time_text

アニメーションを作成する

次に、Matplotlib の FuncAnimation 関数を使ってアニメーションを作成します。

ani = animation.FuncAnimation(
    fig, animate, len(y), interval=dt*1000, blit=True)
plt.show()

まとめ

この実験では、Python と Matplotlib を使って二重回転振子の動きをシミュレートしました。シミュレーション用のパラメータを設定し、システムの微分を計算する関数を定義し、オイラー法を使って常微分方程式を積分し、各時刻での各振子の位置を計算し、二重回転振子の動きのアニメーションを作成しました。