"""
二重振り子のシミュレーションとアニメーションを描画するモジュール。
詳細説明:
このモジュールは、重力下で運動する二重振り子の運動方程式をSciPyを用いて数値的に解き、
その結果をMatplotlibを使用してアニメーションとして表示します。
関連リンク:
:doc:`double_pendulum_usage`
"""
import argparse
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import odeint
# 運動方程式の設定
[ドキュメント]
def derivs(state, t, g, L1, L2, m1, m2):
"""
二重振り子の運動方程式を計算する関数。
詳細説明:
現在の状態（各振り子の角度と角速度）から、ラグランジュ力学に基づく運動方程式を用いて
角度と角速度の時間微分を計算します。
:param state: list, 現在の状態 [theta1, theta2, theta1_dot, theta2_dot]
:param t: float, 現在の時間
:returns: list, 状態の時間微分値 [dtheta1_dt, dtheta2_dt, dtheta1_dot_dt, dtheta2_dot_dt]
"""
dtheta1_dt = state[2]
dtheta2_dt = state[3]
delta = state[1] - state[0]
denominator1 = (m1 + m2) * L1 - m2 * L1 * np.cos(delta) * np.cos(delta)
denominator2 = (L2 / L1) * denominator1
dtheta1_dot_dt = ((m2 * L1 * state[2] ** 2 * np.sin(delta) * np.cos(delta) +
m2 * g * np.sin(state[1]) * np.cos(delta) +
m2 * L2 * state[3] ** 2 * np.sin(delta) -
(m1 + m2) * g * np.sin(state[0])) / denominator1)
dtheta2_dot_dt = ((-m2 * L2 * state[3] ** 2 * np.sin(delta) * np.cos(delta) +
(m1 + m2) * g * np.sin(state[0]) * np.cos(delta) -
(m1 + m2) * L1 * state[2] ** 2 * np.sin(delta) -
(m1 + m2) * g * np.sin(state[1])) / denominator2)
return [dtheta1_dt, dtheta2_dt, dtheta1_dot_dt, dtheta2_dot_dt]
[ドキュメント]
def parse_args():
parser = argparse.ArgumentParser(
description="二重振り子のシミュレーションとアニメーションを描画します。"
)
# 定数の設定
parser.add_argument("--g", type=float, default=9.81,
help="重力加速度 (m/s^2)")
parser.add_argument("--L1", type=float, default=1.0,
help="上部の振り子の長さ (m)")
parser.add_argument("--L2", type=float, default=1.0,
help="下部の振り子の長さ (m)")
parser.add_argument("--m1", type=float, default=1.0,
help="上部の振り子の質量 (kg)")
parser.add_argument("--m2", type=float, default=1.0,
help="下部の振り子の質量 (kg)")
# 初期条件
parser.add_argument("--theta1", type=float, default=np.pi / 2,
help="上部の振り子の初期角度 (ラジアン)")
parser.add_argument("--theta2", type=float, default=np.pi / 2,
help="下部の振り子の初期角度 (ラジアン)")
parser.add_argument("--theta1_dot", type=float, default=0.0,
help="上部の振り子の初期角速度 (rad/s)")
parser.add_argument("--theta2_dot", type=float, default=0.0,
help="下部の振り子の初期角速度 (rad/s)")
parser.add_argument("--tmax", type=float, default=120.0,
help="シミュレーション時間の最大値")
parser.add_argument("--tstep", type=float, default=0.05,
help="時間刻み")
parser.add_argument("--interval", type=int, default=10,
help="アニメーションのフレーム間隔 (ms)")
return parser.parse_args()
[ドキュメント]
def main():
args = parse_args()
# 初期状態と時間配列の設定
state = [args.theta1, args.theta2, args.theta1_dot, args.theta2_dot]
t = np.arange(0.0, args.tmax, args.tstep)
# 運動方程式を数値的に解く
y = odeint(
derivs,
state,
t,
args=(args.g, args.L1, args.L2, args.m1, args.m2),
)
# シミュレーション結果の表示
x1 = args.L1 * np.sin(y[:, 0])
y1 = -args.L1 * np.cos(y[:, 0])
x2 = args.L2 * np.sin(y[:, 1]) + x1
y2 = -args.L2 * np.cos(y[:, 1]) + y1
xmax = args.L1 + args.L2
ymax = args.L1 + args.L2
fig = plt.figure()
ax = fig.add_subplot(
111,
autoscale_on=False,
xlim=(-xmax, xmax),
ylim=(-ymax, ymax),
)
ax.set_aspect('equal')
ax.grid()
line, = ax.plot([], [], 'o-', lw=2)
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
def init():
"""
アニメーションの初期化関数。
詳細説明:
プロットされる線オブジェクトと時間表示テキストをクリアし、アニメーションの描画準備を行います。
:returns: tuple, 初期化された線オブジェクトとテキストオブジェクト (line, time_text)
"""
line.set_data([], [])
time_text.set_text('')
return line, time_text
def animate(i):
"""
アニメーションの各フレームを描画する関数。
詳細説明:
インデックスに基づいて二重振り子の各質点の位置座標を更新し、現在の時間をテキストとして表示します。
:param i: int, フレームのインデックス
:returns: tuple, 更新された線オブジェクトとテキストオブジェクト (line, time_text)
"""
thisx = [0, x1[i], x2[i]]
thisy = [0, y1[i], y2[i]]
line.set_data(thisx, thisy)
time_text.set_text(time_template % t[i])
return line, time_text
ani = animation.FuncAnimation(
fig,
animate,
np.arange(1, len(y)),
interval=args.interval,
blit=True,
init_func=init,
)
plt.show()
if __name__ == "__main__":
main()