振り子とシミュレーション

はじめに

単振り子や2重振り子の運動方程式は、ラグランジアンを用いることで求めることができる。しかし、単振り子は周期的な単純な運動をする一方で、2重振り子は複雑な運動を行う。今回は、Pythonを用いて運動方程式を解析的に解くことで単振り子や2重振り子の運動の様子をシミュレーションすることを目的とする。ただし、シミュレーションに用いた運動方程式は以下のサイトを参考にしたので、導出等を知りたい方は参考にして欲しい。

一番最後に以下の内容を追記した。

運動方程式を数値計算で解く場合、自分がよく用いる簡易なオイラー法では力学的なエネルギー保存されない場合があるという致命的な弱点がある。そこで、scipyを用いたルンゲクッタ法を用いることでこの問題を解決した。
ちなみに、プログラムは生成AIであるClaudeに出力してもらった。

単振り子

運動方程式

xy平面に質点$m_1$が存在していて、原点を支店として糸で繋げられている単振り子を考える。
重力は-y方向にかかっているとする。
質点の座標は糸の長さ$l_1$,重力方向との角度$\theta_1$を用いて以下のように考えることができる。

\begin{equation}
\left\{ \,
    \begin{aligned}
    & x_1 = l_1 sin \theta_1 \\
    & y_1 = -l_1 cos \theta_1 \\
    \end{aligned}
\right.
\end{equation}

一方で動径方向の運動方程式は以下のように表すことができる。

\frac{d^2 \theta_1}{dt^2}= -g sin \theta_1 

プログラム

これを元にプログラムを作成すると以下のようになる。

python simple_pendulum.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import japanize_matplotlib
import math

# アニメ
fig = plt.figure()

ims = []

# パラメータと定数
g=9.8


n=100




# 振り子の長さ
l_1=2


#時刻の初期値
t=0
#角度の初期値と角速度の初期値
theta_1 = math.pi/2
d_theta_dt_1=0
dd_theta_dd_t_1=0

#微小時間
dt =0.1

while(t<10):
 
  dd_theta_dd_t_1 = -g*np.sin(theta_1)

 
  d_theta_dt_1 = d_theta_dt_1 + dd_theta_dd_t_1*dt
  theta_1 = theta_1 + d_theta_dt_1*dt


  t = t+dt
  
  

  x_1=l_1*np.sin(theta_1)
  y_1=-l_1*np.cos(theta_1)
  


  X_1= np.linspace(0,x_1,n)
  Y_1= np.linspace(0,y_1,n)


  im1=plt.plot(X_1,Y_1,'-')
  
  ims.append(im1)


# 10枚のプロットを 100ms ごとに表示
ani = animation.ArtistAnimation(fig, ims, interval=100)
#保存
ani.save("furiko_simple.gif", writer="pillow")
plt.show()

このプログラムを実行すると以下のような動画が出力される。

考察

上記の動画のように単振り子は往復運動を示す。

2重振り子

運動方程式

運動方程式の導出は以下のサイトを参照して欲しい。本記事では、導出された結果を用いるのみとする。

\begin{equation}
\left\{ \,
    \begin{aligned}
    & \ddot{\theta_1} = -\frac{(m_1+m_2)gl_1 sin\theta_1+m_2 l_1 l_2 (\ddot{\theta_2} cos(\theta_1-\theta_2)+\ddot {\theta_2^2} sin(\theta_1-\theta_2))}{(m_1+m_2)l_1^2}\\
    & \ddot{\theta_2} = -\frac{m_2 gl_1 sin\theta_2+m_2 l_1 l_2 (\ddot{\theta_1} cos(\theta_1-\theta_2)+\ddot \theta_1^2 sin(\theta_1-\theta_2))}{m_2 l_2^2} \\
    \end{aligned}
\right.
\end{equation}

一方で、質点$m_1$と質点$m_2$の位置は以下のようになる。

\begin{equation}
\left\{ \,
    \begin{aligned}
    & x_1 = l_1 sin \theta_1 \\
    & y_1 = -l_1 cos \theta_1 \\
    \end{aligned}
\right.
\end{equation}

\begin{equation}
\left\{ \,
    \begin{aligned}
    & x_2 = l_1 sin \theta_1+ l_2 sin \theta_2 \\
    & y_2 = -l_1 cos \theta_1- l_2 cos \theta_2 \\
    \end{aligned}
\right.
\end{equation}

プログラム

以上のことに留意してプログラムを書くと以下のようになる。

python double_pendulum.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import japanize_matplotlib
import math

# アニメ
fig = plt.figure()

ims = []

# パラメータと定数
g=9.8


n=100

m_1=1.0
m_2=2.0

# 振り子の長さ
l_1=2
l_2=1

#時刻の初期値
t=0
#角度の初期値と角速度の初期値
theta_1 = math.pi/2
d_theta_dt_1=0
dd_theta_dd_t_1=0

theta_2 = 0
d_theta_dt_2=0
dd_theta_dd_t_2=0

#微小時間
dt =0.01

while(t<10):
 
  dd_theta_dd_t_1 = -((m_1+m_2)*g*l_1*np.sin(theta_1)+m_2*l_1*l_2*(dd_theta_dd_t_2*np.cos(theta_1-theta_2)+((d_theta_dt_2)**2)*np.sin(theta_1-theta_2)))/((m_1+m_2)*l_1**2)

  dd_theta_dd_t_2 = -(m_2*g*l_2*np.sin(theta_2)+m_2*l_1*l_2*(dd_theta_dd_t_1*np.cos(theta_1-theta_2)-((d_theta_dt_1)**2)*np.sin(theta_1-theta_2)))/(m_2*l_2**2)


  d_theta_dt_1 = d_theta_dt_1 + dd_theta_dd_t_1*dt
  theta_1 = theta_1 + d_theta_dt_1*dt

  d_theta_dt_2 = d_theta_dt_2 + dd_theta_dd_t_2*dt
  theta_2 = theta_2 + d_theta_dt_2*dt

  t = t+dt
  
  

  x_1=l_1*np.sin(theta_1)
  y_1=-l_1*np.cos(theta_1)
  
  x_2= l_1*np.sin(theta_1)+l_2*np.sin(theta_2)
  y_2= -l_1*np.cos(theta_1)-l_2*np.cos(theta_2)

  X_1= np.linspace(0,x_1,n)
  Y_1= np.linspace(0,y_1,n)

  X_2= np.linspace(x_1,x_2,n)
  Y_2= np.linspace(y_1,y_2,n)

  im1=plt.plot(X_1,Y_1,'-')
  im2=plt.plot(X_2,Y_2,'-')
  ims.append(im1+im2)


# 10枚のプロットを 100ms ごとに表示
ani = animation.ArtistAnimation(fig, ims, interval=100)
#保存
ani.save("furiko_2tai_simple.gif", writer="pillow")
plt.show()

考察

このプログラムを実行すると以下のようになる。

このように、2重振り子はカオスな運動となる。しかし、やはり運動の様子を正確に記述するためには十分な計算量が必要となる。

補足

ただ、この手法だと力学的エネルギーが保存されていない場所が存在するので、以下のようなプログラムを生成AIによって出力した。

python Pendukum_scipy.py

"""
単振り子・2重振り子 — scipy.integrate.solve_ivp による高精度シミュレーション
元記事: https://qiita.com/arairuca/items/a2cd950d7dc2a8d7fb74

【元のコードの問題点】
  単振り子:
    - 差分法（オイラー法）のため dt=0.1 では誤差蓄積が大きく、
      エネルギーが保存されず振幅が増大/減衰してしまう。

  2重振り子:
    - 運動方程式に根本的なバグあり:
        元コードの式は θ1_ddot, θ2_ddot が互いに依存する連立方程式なのに、
        単純に逐次代入しているだけ(θ2_ddot=0 の前の値を θ1_ddot に使っている)。
        正しくは連立方程式を解いて θ1_ddot, θ2_ddot を同時に求める必要がある。
    - その結果、θ2_ddot の式にも誤り: m_2*g*l_2 となっているが正しくは m_2*g*l_1。

【改善内容】
  1. scipy.integrate.solve_ivp + DOP853 (8次RK、自動ステップ幅制御) を使用
  2. 2重振り子: 連立方程式を解析的に解いて正しい加速度を導出
  3. エネルギー保存の確認プロットを追加
  4. 軌跡のトレイル表示でアニメーションを改善
"""

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


# ════════════════════════════════════════════════════════════════
# 1.  単 振 り 子
# ════════════════════════════════════════════════════════════════

def simple_pendulum_ode(t, y, g, l):
    """
    状態ベクトル y = [theta, omega]
    dydt = [omega, -g/l * sin(theta)]
    """
    theta, omega = y
    return [omega, -(g / l) * np.sin(theta)]


def run_simple_pendulum():
    g = 9.8
    l_1 = 2.0
    theta0 = np.pi / 2      # 初期角度 (元記事と同じ)
    omega0 = 0.0            # 初期角速度
    y0 = [theta0, omega0]
    t_end = 10.0

    sol = solve_ivp(
        simple_pendulum_ode,
        [0, t_end],
        y0,
        args=(g, l_1),
        method='DOP853',
        rtol=1e-10,
        atol=1e-12,
        dense_output=True
    )

    n_frames = 200
    t_eval = np.linspace(0, t_end, n_frames)
    Y = sol.sol(t_eval)    # shape (2, n_frames)
    theta = Y[0]

    # 座標変換
    x1 = l_1 * np.sin(theta)
    y1 = -l_1 * np.cos(theta)

    # ── アニメーション ──────────────────────────────────────
    fig, ax = plt.subplots(figsize=(5, 5))

    def update(i):
        ax.clear()
        ax.plot([0, x1[i]], [0, y1[i]], '-o', color='steelblue',
                lw=2, ms=10, markerfacecolor='tomato')
        ax.set_xlim(-l_1 * 1.3, l_1 * 1.3)
        ax.set_ylim(-l_1 * 1.3, l_1 * 0.3)
        ax.set_aspect('equal')
        ax.set_title(f'単振り子  t = {t_eval[i]:.2f} s  (DOP853)', fontsize=11)
        ax.set_xlabel('x [m]'); ax.set_ylabel('y [m]')
        ax.axhline(0, color='gray', lw=0.5)
        ax.axvline(0, color='gray', lw=0.5)

    ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=50)
    ani.save('simple_pendulum_scipy.gif', writer='pillow', fps=20)
    plt.close(fig)

    # ── エネルギー保存の確認 ──────────────────────────────
    omega = Y[1]
    KE = 0.5 * l_1**2 * omega**2          # 運動エネルギー (m=1)
    PE = g * l_1 * (1 - np.cos(theta))    # ポテンシャルエネルギー (m=1, 最下点を基準)
    E_total = KE + PE

    fig2, ax2 = plt.subplots(figsize=(7, 3))
    ax2.plot(t_eval, KE,      label='運動エネルギー', color='tomato')
    ax2.plot(t_eval, PE,      label='位置エネルギー', color='royalblue')
    ax2.plot(t_eval, E_total, label='全エネルギー',   color='seagreen', lw=2)
    ax2.set_xlabel('時間 [s]'); ax2.set_ylabel('エネルギー [J/kg]')
    ax2.set_title('単振り子 — エネルギー保存の確認 (DOP853)')
    ax2.legend(); ax2.grid(True, alpha=0.4)
    fig2.tight_layout()
    fig2.savefig('simple_pendulum_energy.png', dpi=120)
    plt.close(fig2)

    print('単振り子: simple_pendulum_scipy.gif / simple_pendulum_energy.png を保存しました')
    print(f'  全エネルギー変動: {E_total.max() - E_total.min():.2e} J/kg  (完全保存なら 0)')


# ════════════════════════════════════════════════════════════════
# 2.  2 重 振 り 子
# ════════════════════════════════════════════════════════════════
#
# 【正しい運動方程式の導出】
# ラグランジアンから導かれる連立方程式 (標準形):
#
#   (m1+m2) l1 θ1'' + m2 l2 θ2'' cos(θ1-θ2) + m2 l2 θ2'^2 sin(θ1-θ2) + (m1+m2) g sinθ1 = 0
#   m2 l2 θ2'' + m2 l1 θ1'' cos(θ1-θ2) - m2 l1 θ1'^2 sin(θ1-θ2) + m2 g sinθ2 = 0
#
# これを θ1'', θ2'' について連立して解くと (行列式を展開):
#
# 分母: D = (m1+m2)*l1^2 * m2*l2^2 - (m2*l1*l2*cos(θ1-θ2))^2
#        = m2*l1*l2 * [(m1+m2)*l1*l2 - m2*l1*l2*cos^2(Δ)]
#
# この解析解を ODE 関数内で毎ステップ計算する。

def double_pendulum_ode(t, y, m1, m2, l1, l2, g):
    """
    状態ベクトル y = [θ1, θ2, ω1, ω2]
    連立方程式を解析的に解いた正しい加速度を返す。
    """
    theta1, theta2, omega1, omega2 = y
    delta = theta1 - theta2
    cos_d = np.cos(delta)
    sin_d = np.sin(delta)

    # 分母
    denom = (m1 + m2) * l1 * l2 - m2 * l1 * l2 * cos_d**2
    # = l1*l2 * [(m1+m2) - m2*cos^2(delta)]

    # 分子 θ1''
    num1 = (
        m2 * l2 * omega2**2 * sin_d
        - (m1 + m2) * g * np.sin(theta1)
        + m2 * g * np.cos(delta) * np.sin(theta2)
        - m2 * l1 * omega1**2 * cos_d * sin_d  # クロスタームの寄与
    )
    # より正確な表現 (行列の逆行列で導出):
    # α = (m1+m2)*l1, β = m2*l2*cos_d, γ = m2*l2, δ_ = m1*l1*cos_d
    # [[α, β],[γ_,δ_]] [[θ1''],[θ2'']] = [f1, f2]
    alpha = (m1 + m2) * l1
    beta  = m2 * l2 * cos_d
    gamma = l1 * cos_d
    delta_ = l2

    f1 = -m2 * l2 * omega2**2 * sin_d - (m1 + m2) * g * np.sin(theta1)
    f2 =  l1 * omega1**2 * sin_d      - g * np.sin(theta2)

    det = alpha * delta_ - beta * gamma
    alpha2 = delta_ * f1 - beta * f2
    alpha3 = alpha * f2 - gamma * f1

    theta1_ddot = alpha2 / det
    theta2_ddot = alpha3 / det

    return [omega1, omega2, theta1_ddot, theta2_ddot]


def run_double_pendulum():
    g = 9.8
    m1, m2 = 1.0, 2.0
    l1, l2 = 2.0, 1.0
    theta1_0 = np.pi / 2    # 元記事と同じ
    theta2_0 = 0.0
    omega1_0 = 0.0
    omega2_0 = 0.0
    y0 = [theta1_0, theta2_0, omega1_0, omega2_0]
    t_end = 20.0            # 元記事より長時間でも安定

    sol = solve_ivp(
        double_pendulum_ode,
        [0, t_end],
        y0,
        args=(m1, m2, l1, l2, g),
        method='DOP853',
        rtol=1e-11,
        atol=1e-13,
        dense_output=True
    )

    n_frames = 400
    t_eval = np.linspace(0, t_end, n_frames)
    Y = sol.sol(t_eval)   # shape (4, n_frames)

    theta1 = Y[0]; theta2 = Y[1]
    omega1 = Y[2]; omega2 = Y[3]

    # 座標変換
    x1 = l1 * np.sin(theta1)
    y1 = -l1 * np.cos(theta1)
    x2 = x1 + l2 * np.sin(theta2)
    y2 = y1 - l2 * np.cos(theta2)

    # ── アニメーション ──────────────────────────────────────
    fig, ax = plt.subplots(figsize=(6, 6))
    trail_len = 80   # 質点2の軌跡フレーム数

    def update(i):
        ax.clear()
        lo = max(0, i - trail_len)
        # 質点2の軌跡
        ax.plot(x2[lo:i+1], y2[lo:i+1], '-', color='seagreen', alpha=0.4, lw=1)
        # ロッド
        ax.plot([0, x1[i], x2[i]], [0, y1[i], y2[i]],
                '-o', color='steelblue', lw=2, ms=8,
                markerfacecolor='tomato', markeredgecolor='steelblue')
        ax.plot(x2[i], y2[i], 'o', color='seagreen', ms=10)
        # 支点
        ax.plot(0, 0, 'k^', ms=10)

        span = (l1 + l2) * 1.2
        ax.set_xlim(-span, span)
        ax.set_ylim(-span, span * 0.3)
        ax.set_aspect('equal')
        ax.set_title(f'2重振り子  t = {t_eval[i]:.2f} s  (DOP853)', fontsize=11)
        ax.set_xlabel('x [m]'); ax.set_ylabel('y [m]')
        ax.axhline(0, color='gray', lw=0.5)
        ax.axvline(0, color='gray', lw=0.5)

    ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=50)
    ani.save('double_pendulum_scipy.gif', writer='pillow', fps=20)
    plt.close(fig)

    # ── エネルギー保存の確認 ──────────────────────────────
    # 運動エネルギー
    KE = (0.5 * m1 * l1**2 * omega1**2
          + 0.5 * m2 * (l1**2 * omega1**2
                        + l2**2 * omega2**2
                        + 2 * l1 * l2 * omega1 * omega2 * np.cos(theta1 - theta2)))
    # ポテンシャルエネルギー (最下点を基準)
    PE = (m1 * g * l1 * (1 - np.cos(theta1))
          + m2 * g * (l1 * (1 - np.cos(theta1)) + l2 * (1 - np.cos(theta2))))
    E_total = KE + PE

    fig2, ax2 = plt.subplots(figsize=(8, 3))
    ax2.plot(t_eval, KE,      label='運動エネルギー', color='tomato', alpha=0.8)
    ax2.plot(t_eval, PE,      label='位置エネルギー', color='royalblue', alpha=0.8)
    ax2.plot(t_eval, E_total, label='全エネルギー',   color='seagreen', lw=2)
    ax2.set_xlabel('時間 [s]'); ax2.set_ylabel('エネルギー [J]')
    ax2.set_title('2重振り子 — エネルギー保存の確認 (DOP853)')
    ax2.legend(); ax2.grid(True, alpha=0.4)
    fig2.tight_layout()
    fig2.savefig('double_pendulum_energy.png', dpi=120)
    plt.close(fig2)

    print('2重振り子: double_pendulum_scipy.gif / double_pendulum_energy.png を保存しました')
    print(f'  全エネルギー変動: {E_total.max() - E_total.min():.2e} J  (完全保存なら 0)')


# ════════════════════════════════════════════════════════════════
if __name__ == '__main__':
    print('=== 単振り子を計算中 ===')
    run_simple_pendulum()
    print()
    print('=== 2重振り子を計算中 ===')
    run_double_pendulum()
    print()
    print('完了')

以下の動画のように、かなり精度よく2重振り子の様子を描写することができた。

また、力学的エネルギーについても保存されていることが分かる。

まとめ

今回は、単振り子や2重振り子の運動方程式からPythonを用いて振り子のシミュレーションを行うことを目的としてプログラムを作成した。結果、微分方程式を解くために用いた差分法によるアルゴリズムは上手く動作しているようで、シミュレーションをすることができた。