【python】ルンゲクッタ法による2重振り子のカオス挙動シミュレーション

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

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

前回、多自由度系の運動方程式をルンゲクッタ法でシミュレーションする方法について記載しました。
mori-memo.hateblo.jp

今回は、それを応用してよくある2重振り子のカオス挙動をシミュレーションしてみましょう。
wikipediaにて、『2重振り子』と検索すると2重振り子のアニメーションが出てくるので、今回の目標はそれを作成することとします
ja.wikipedia.org

2重振り子の運動方程式

対象のシステムは下記の通りとします。

このシステムの運動方程式は、下記のような2階の常微分方程式で表すことができます。

$(m_{1} + m_{2})l_{1}\ddot{\theta_{1}} + m_{2}l_{2}\ddot{\theta_{2}}\cos(\theta_{1} - \theta_{2}) + m_{2}l_{2}\dot{\theta_{2}}^{2}\sin(\theta_{1} - \theta_{2}) + (m_{1} + m_{2})g\sin\theta_{1} = 0$
$l_{1}l_{2}\ddot{\theta_{1}}\cos(\theta_{1} - \theta_{2}) + l_{2}^{2}\ddot{\theta_{2}} - l_{1}l_{2}\dot{\theta_{1}}^{2}\sin(\theta_{1} - \theta_{2}) + gl_{2}\sin\theta_{2} = 0$

式変形

ここで、いろいろと式変形を行い、2階の常微分方程式を1階の常微分方程式の形で表します。
簡単のため、下記の変数を定義しておきます。

$C = \cos(\theta_{1} - \theta_{2})$
$S = \sin(\theta_{1} - \theta_{2})$
$M = m_{1} + m_{2}$

$\dot{\theta_{1}}, \dot{\theta_{2}}$ をそれぞれ $\omega_{1}, \omega_{2}$ と置くと、上記の運動方程式は、下記の1階の連立常微分方程式に変換できる。

$\dot{\theta_{1}} = \omega_{1}$
$\dot{\theta_{2}} = \omega_{2}$
$\dot{\omega_{1}} = \frac{g(C\sin\theta_{2} - \frac{M}{m_{2}}\sin\theta_{1}) - (l_{1}\omega_{1}^{2}C + l_{2}\omega_{2}^{2})S}{l_{1}(\frac{M}{m_{2}} - C^{2})}$
$\dot{\omega_{2}} = \frac{g\frac{M}{m_{2}}(C\sin\theta_{1} - \sin\theta_{2}) + (\frac{M}{m_{2}}l_{1}\omega_{1}^{2}C + l_{2}\omega_{2}^{2}C)S}{l_{2}(\frac{M}{m_{2}} - C^{2})}$

これにより $\omega_{1}, \omega_{2}, \theta_{1}, \theta_{2}$ を陽解法で解ける形になりました。
以前せっかく作ったので、ルンゲクッタ法を用いてシミュレーションしてみましょう。

Pythonで実装

コード全文は以下です。

import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
from solver import runge_kutta_method
import numpy as np

def func(omega, theta):
    '''運動方程式'''
    # 各種変数定義
    C = np.cos(theta[0] - theta[1])
    S = np.sin(theta[0] - theta[1])
    M = np.sum(m)
    # マトリクスの生成
    v = np.zeros((2))
    v[0] = (g * (np.sin(theta[1]) * C - M * np.sin(theta[0]) / m[1]) \
            - (l[0] * C * omega[0] ** 2 + l[1] * omega[1] ** 2) * S) \
            / (l[0] * (M / m[1] - C ** 2))
    v[1] = (g * (M / m[1]) * (C * np.sin(theta[0]) - np.sin(theta[1])) \
            + ((M * l[0] * omega[0] ** 2) / m[1] + (l[1] * omega[1] ** 2) * C) * S) \
            / (l[1] * (M / m[1] - C ** 2))
    return v


l = np.array((0.6, 0.5))                             # 振り子の長さ
m = np.array((0.5, 0.7))                             # 振り子の重さ
g = 9.8                                              # 重力加速度
omega = np.array((0, 0))                              # 初期の角速度(今回は0)
theta = np.array((np.pi / 2, np.pi / 2))              # 初期の角度(今回は90度)
t_max = 20                        # 終了時間
dt = 0.001                    # サンプリング間隔
t = np.arange(0, t_max, dt)       # 時間データ

# プロット(各自由度の初期間隔を2としてプロットする)
x_rk = runge_kutta_method(t, dt, theta, omega, func)         # theta
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
line2, = ax.plot([], [], color='k')         # 軌跡のトレース
im, = ax.plot([], [], color='b', marker='o', \
                            markersize=10, linestyle='None')         # 振り子の球の部分
line, = ax.plot([], [], color='b')         # 振り子の棒の部分
x_track = l[0] * np.sin(x_rk[:, 0]) + l[1] * np.sin(x_rk[:, 1])
y_track = - l[0] * np.cos(x_rk[:, 0]) - l[1] * np.cos(x_rk[:, 1])
ax.set_xlim(-1.2, 1.2)
ax.set_ylim(-2.2, 0.2)
ax.set_aspect('equal')

# 初期化関数
def init():  # only required for blitting to give a clean slate.
    im.set_xdata([np.nan, np.nan])
    im.set_ydata([np.nan, np.nan])
    line.set_xdata([np.nan, np.nan, np.nan])
    line.set_ydata([np.nan, np.nan, np.nan])
    line2.set_xdata([np.nan])
    line2.set_ydata([np.nan])
    return line, line2, im,

# グラフ更新関数
def update_anim(f):
    theta1 = x_rk[f, 0]
    theta2 = x_rk[f, 1]
    x = [l[0] * np.sin(theta1), l[0] * np.sin(theta1) + l[1] * np.sin(theta2)]
    y = [- l[0] * np.cos(theta1), - l[0] * np.cos(theta1) - l[1] * np.cos(theta2)]
    line2.set_data(x_track[:f], y_track[:f])
    im.set_data(x, y)
    line.set_data([0, x[0], x[1]], [0, y[0], y[1]])
    return line, line2, im,

ani = FuncAnimation(
      fig,  # Figureオブジェクト
      update_anim,  # グラフ更新関数
      init_func=init,
      frames = range(0, 20000, 100),  # フレームを設定
      interval = 100,  # 更新間隔(ms)
      repeat = True,  # 描画を繰り返す
      blit = True,
      )

# アニメーションのgif形式での保存
ani.save("anim_test.gif", writer='pillow')
plt.show()

実際に出力された結果はこちら。

プログラムの説明

以前の記事でアニメーションを作った際と題材がほぼ同じなので、違うところだけ説明します。
mori-memo.hateblo.jp

ルンゲクッタ法のモジュール化

まず、solver.pyというファイルを作って以下の記事で作ったコードを入れています。ルンゲクッタ法の関数です
mori-memo.hateblo.jp
solver.pyの中身は以下です。

import numpy as np

def runge_kutta_method(t, dt, x_0, v_0, func):
    '''ルンゲクッタ法'''
    n = len(x_0)
    x = np.zeros((len(t), n))
    v = np.zeros((len(t), n))
    for i in range(len(t) - 1):
        if i == 0:
            x[i, :] = x_0
            v[i, :] = v_0

        k1_x = v[i, :] * dt
        k1_v = func(v[i, :], x[i, :]) * dt

        k2_x = (v[i, :] + k1_v / 2) * dt
        k2_v = func(v[i, :] + k1_v / 2, x[i, :] + k1_x / 2) * dt

        k3_x = (v[i, :] + k2_v / 2) * dt
        k3_v = func(v[i, :] + k2_v / 2, x[i, :] + k2_x / 2) * dt

        k4_x = (v[i, :] + k3_v) * dt
        k4_v = func(v[i, :] + k3_v / 2, x[i, :] + k3_x) * dt

        x[i + 1, :] = x[i, :] + (k1_x + 2 * k2_x + 2 * k3_x + k4_x) / 6
        v[i + 1, :] = v[i, :] + (k1_v + 2 * k2_v + 2 * k3_v + k4_v) / 6
    return x

入力には終了時刻、タイムステップ、変位の初期値、速度の初期値、運動方程式の関数をいれます。
出力は各自由度の各タイムステップでの変位が2次元配列で出力されます。
今回は変位ではなく角度ですね。

球の軌跡のプロット

今回は球の軌跡のプロットを新しく導入しました。
ルンゲクッタ法で $\theta_{1}, \theta_{2}$ を求めた後、変位に変換して各タイムステップでの変位を計算しておきます。

x_track = l[0] * np.sin(x_rk[:, 0]) + l[1] * np.sin(x_rk[:, 1])
y_track = - l[0] * np.cos(x_rk[:, 0]) - l[1] * np.cos(x_rk[:, 1])

その後、グラフ更新関数(update_anim)内で現在のタイムステップまでの軌跡をスライスで指定してline2に.set_dataメソッドでデータを渡しています。

def update_anim(f):
    # 略
    line2.set_data(x_track[:f], y_track[:f])
    # 略

他は全部これまでやったコードを使いまわしているだけなのですぐ理解できるかと思います。

終わりに

初期位置や、棒の長さ、重さ等パラメータを変えていろいろなカオスを楽しみましょう！

リンク

morikomorou’s blog

自分が学んだことなどの備忘録的なやつ