初めてYoutube埋め込んだ
はじめに
ロボティクス・バイオメカニクスを専攻する端くれとして、弾性・粘性を変化させたときの質点を挙動を可視化した。
徐々に自分の専攻に関わる分野で、遊べそうなものを可視化したい。
対象読者:
- Pythonの基本的な文法（変数、リスト、ループ、関数）を少し知っている方
- 物理シミュレーションやデータ可視化に興味がある方
- Matplotlibでアニメーションを作ってみたい方
目的
- アニメーション描画の振り返り
- ルンゲクッタ法の振り返り
コード全体
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FFMpegWriter
# パラメータ
m = 60.0 # 質量 [kg]
k = 30.0 # バネ定数 [N/m]
c = 2.0 # 減衰係数 [Ns/m]
g = 9.81 # 重力加速度 [m/s^2]
# 初期条件
x0 = 10.0 # 初期変位 [m]
v0 = 0.0 # 初期速度 [m/s]
# 時間設定
t_max = 15.0
dt = 0.01
steps = int(t_max / dt)
t = np.linspace(0, t_max, steps)
# 配列準備
x = np.zeros(steps)
v = np.zeros(steps)
x[0], v[0] = x0, v0
# 運動方程式
def dx_dt(x, v): return v
def dv_dt(x, v): return -(c/m)*v - (k/m)*x - g
# RK4 で解を事前計算
for i in range(steps - 1):
k1_x = dt * dx_dt(x[i], v[i])
k1_v = dt * dv_dt(x[i], v[i])
k2_x = dt * dx_dt(x[i] + 0.5*k1_x, v[i] + 0.5*k1_v)
k2_v = dt * dv_dt(x[i] + 0.5*k1_x, v[i] + 0.5*k1_v)
k3_x = dt * dx_dt(x[i] + 0.5*k2_x, v[i] + 0.5*k2_v)
k3_v = dt * dv_dt(x[i] + 0.5*k2_x, v[i] + 0.5*k2_v)
k4_x = dt * dx_dt(x[i] + k3_x, v[i] + k3_v)
k4_v = dt * dv_dt(x[i] + k3_x, v[i] + k3_v)
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
# 図の準備：左右2分割
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# 左：バネ―重りアニメーション
ax1.set_xlim(0.5, 1.5)
ax1.set_ylim(-100, 15)
ax1.set_title('Spring-weight motion')
ax1.set_ylabel('Displacement [m]')
line_mass, = ax1.plot([], [], 'o-', lw=2)
# 左上にパラメータ表示
param_text = (
r"$x_0 = %.2f\ \mathrm{m}$" "\n"
r"$v_0 = %.2f\ \mathrm{m/s}$" "\n"
r"$m = %.2f\ \mathrm{kg}$" "\n"
r"$k = %.2f\ \mathrm{N/m}$" "\n"
r"$c = %.2f\ \mathrm{Ns/m}$"
) % (x0, v0, m, k, c)
ax1.text(
0.70, 0.95, param_text,
transform=ax1.transAxes,
fontsize=12,
linespacing=1.8,
verticalalignment='top',
horizontalalignment='left',
bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.8)
)
# 右：時間 vs 変位
ax2.set_xlim(0, t_max)
ax2.set_ylim(np.min(x)*1.1, np.max(x)*1.1)
ax2.set_title('Time-Displacement graph')
ax2.set_xlabel('Time [s]')
ax2.set_ylabel('Displacement [m]')
ax2.grid()
ax2.axhline(0, color='black', lw=0.5, ls='--')
line_graph, = ax2.plot([], [], lw=1)
# 初期化関数
def init():
line_mass.set_data([], [])
line_graph.set_data([], [])
return line_mass, line_graph
# 更新関数
def update(frame):
# 左側：バネを固定点(1,0)─(1,x)で描画
pos = x[frame]
line_mass.set_data([1, 1], [0, pos])
# 右側：時刻 0→t[frame] の変位履歴
line_graph.set_data(t[:frame], x[:frame])
return line_mass, line_graph
# アニメーション生成
ani = animation.FuncAnimation(
fig, update,
frames=range(0, steps, 5),
init_func=init,
blit=True
)
# MP4 保存設定
metadata = dict(artist='Me', comment='Damped oscillator')
writer = FFMpegWriter(
fps=20,
codec='libx264',
bitrate=1800,
metadata=metadata,
extra_args=[
'-profile:v', 'baseline', # Baseline プロファイル指定
'-level', '3.0', # レベル指定（オプション）
'-pix_fmt', 'yuv420p' # 再生互換性の高いピクセルフォーマット
]
)
ani.save("damped_oscillator_dual2.mp4", writer=writer)
plt.close(fig)
print("アニメーションが damped_oscillator_dual.mp4 に保存されました。")
コードのメモ
1. ライブラリのインポート
最初のブロックでは、シミュレーションやグラフ描画に必要な「ライブラリ」を読み込む。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FFMpegWriter
from matplotlib.animation import FFMpegWriter:
作成したアニメーションを動画ファイル（MP4など）として保存するために使う FFMpegWriter クラスを直接読み込む。これを使うには、自身の環境に FFmpeg という動画処理ソフトウェアがインストールされている必要がある。
応用ポイント:
アニメーションをGIFで保存したい場合は PillowWriter や ImageMagickWriter を使うこともある。（from matplotlib.animation import PillowWriter など）
2. シミュレーションのパラメータ設定
ここでは、シミュレーション対象である「バネとおもり」の物理的な特性値を設定。これらの値を変えることで、振動の様子がどのように変化するかを後で試すことができる。
# パラメータ
m = 60.0 # 質量 [kg]
k = 30.0 # バネ定数 [N/m]
c = 2.0 # 減衰係数 [Ns/m]
g = 9.81 # 重力加速度 [m/s^2]
- これらのパラメータの値を色々変えて実行
-
cを0にするとどうなる？（減衰のない振動） -
kを大きく/小さくすると？（振動の速さ） -
mを大きく/小さくすると？（動きやすさ）
-
- 自分の興味のある物理現象に合わせて、パラメータの種類や値を変えることで、様々なシミュレーションに応用できる。
3. 初期条件の設定
シミュレーションを開始する時点（時刻 t=0）でのおもりの状態を設定。
# 初期条件
x0 = 10.0 # 初期変位 [m]
v0 = 0.0 # 初期速度 [m/s]
-
x0(初期変位): シミュレーション開始時のおもりの位置（つり合いの位置からのずれ）。- 単位はメートル[m]。
-
v0(初期速度): シミュレーション開始時のおもりの速度。- 単位はメートル毎秒[m/s]。
-
0なら静かに手を離した状態、正や負の値なら初めにある速度を与えた状態。
4. 運動方程式の定義
シミュレーションの核となる、おもりの運動を表す「運動方程式」をPythonの関数として定義。今回は、変位 x と速度 v の時間変化を表す2つの微分方程式を定義。
# 運動方程式
def dx_dt(x, v): return v
def dv_dt(x, v): return -(c/m)*v - (k/m)*x - g
-
def dx_dt(x, v): return v:- これは変位
xの時間微分 $ \frac{dx}{dt} $ を計算する関数 - 物理学的に、速度 $v$ は変位 $x$ の時間微分なので、この関数は単純に引数で受け取った速度
vをそのまま返す
- これは変位
-
def dv_dt(x, v): return -(c/m)*v - (k/m)*x - g:- これは速度
vの時間微分 $ \frac{dv}{dt} $ （つまり加速度）を計算する関数 - これは、バネダンパー系の運動方程式 $ m \ddot{x} + c \dot{x} + kx = -mg $を $ \ddot{x} = \frac{dv}{dt} = -\frac{c}{m}v - \frac{k}{m}x - g $ と変形したもの
- $\ddot{x}$ :加速度、$ \dot{x} $:速度
- 関数は現在の変位
xと速度vを引数として受け取り、パラメータc,m,k,gを使って次の瞬間の速度変化率（加速度）を計算して返す
- これは速度
応用ポイント:
- シミュレーションしたい物理現象に合わせて、この運動方程式の関数を書き換えることで、様々な現象をシミュレーションできる
- 例えば、振り子の運動、惑星の軌道計算など
- 重要なのは、状態の時間微分（変化率）を計算する関数を定義すること
5. 数値積分 (4次ルンゲ＝クッタ法)
- 運動方程式（微分方程式）をコンピューターで解くために、「数値積分」という手法を使う
- ここでは、比較的精度が高い「4次ルンゲ＝クッタ法 (RK4)」を用いて、時刻
dt後の変位xと速度vを計算 - それを
steps回繰り返して全体の動きをシミュレーションする
# RK4 で解を事前計算
for i in range(steps - 1): # steps-1回繰り返す (0からsteps-2まで)
# 現在時刻 t[i] での x[i], v[i] を使って計算
# k1: 現在地点での傾き
k1_x = dt * dx_dt(x[i], v[i])
k1_v = dt * dv_dt(x[i], v[i])
# k2: k1を使ってステップの中間地点での傾きを予測
k2_x = dt * dx_dt(x[i] + 0.5*k1_x, v[i] + 0.5*k1_v)
k2_v = dt * dv_dt(x[i] + 0.5*k1_x, v[i] + 0.5*k1_v)
# k3: k2を使ってステップの中間地点での傾きを再予測
k3_x = dt * dx_dt(x[i] + 0.5*k2_x, v[i] + 0.5*k2_v)
k3_v = dt * dv_dt(x[i] + 0.5*k2_x, v[i] + 0.5*k2_v)
# k4: k3を使ってステップの終点での傾きを予測
k4_x = dt * dx_dt(x[i] + k3_x, v[i] + k3_v)
k4_v = dt * dv_dt(x[i] + k3_x, v[i] + k3_v)
# 4つの傾きを重み付け平均して、次の時刻 t[i+1] の値を計算
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
6. グラフ描画の準備
- シミュレーション結果をアニメーションで表示するための準備
- Matplotlibを使って、図の全体構成や各グラフの初期設定を行う
# 図の準備：左右2分割の描画エリアを作成
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# fig: 図全体 / ax1: 左の描画エリア / ax2: 右の描画エリア
左側のグラフ (バネとおもりの動き)
# 左：バネ―重りアニメーション
ax1.set_xlim(0.5, 1.5) # x軸の表示範囲
ax1.set_ylim(-100, 15) # y軸の表示範囲
ax1.set_title('Spring-weight motion') # グラフタイトル
ax1.set_ylabel('Displacement [m]') # y軸ラベル
# アニメーションで動かす線オブジェクトを空で初期化
line_mass, = ax1.plot([], [], 'o-', lw=2, markersize=8)
# 'o-': マーカー(点)と実線 / lw: 線の太さ / markersize: マーカーの大きさ
左上のパラメータ表示
# 左上にパラメータ表示
# 表示する文字列を準備 (LaTeX形式使用)
param_text = (
r"$x_0 = %.2f\ \mathrm{m}$" "\n" # 初期変位
r"$v_0 = %.2f\ \mathrm{m/s}$" "\n" # 初期速度
r"$m = %.2f\ \mathrm{kg}$" "\n" # 質量
r"$k = %.2f\ \mathrm{N/m}$" "\n" # バネ定数
r"$c = %.2f\ \mathrm{Ns/m}$" # 減衰係数
) % (x0, v0, m, k, c) # %演算子で数値を文字列に埋め込む
# テキストを描画
ax1.text(
0.70, 0.95, param_text, # 表示位置(相対座標 0.7, 0.95)と表示文字列
transform=ax1.transAxes, # 座標系をAxes基準にする
fontsize=12, # フォントサイズ
linespacing=1.8, # 行間
verticalalignment='top', # 垂直方向の揃え（上揃え）
horizontalalignment='left', # 水平方向の揃え（左揃え）
# テキストを囲むボックスの設定
bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.8)
)
-
param_text = (...) % (x0, v0, m, k, c):- 表示したいパラメータの文字列を組み立てる
-
r"$...$": LaTeX形式の数式を使えるようにr(raw string) と$で囲む -
%.2f: これは文字列フォーマット指定子で、対応する変数の値を小数点以下2桁までの浮動小数点数として埋め込む -
% (x0, v0, m, k, c): 文字列中の%fに、括弧内の変数の値を順番に埋め込む
右側のグラフ (時間 vs 変位)
# 右：時間 vs 変位
ax2.set_xlim(0, t_max) # x軸(時間)の表示範囲
# y軸(変位)の表示範囲を計算結果に基づいて設定 (最小値・最大値の1.1倍)
ax2.set_ylim(np.min(x)*1.1, np.max(x)*1.1)
ax2.set_title('Time-Displacement graph') # グラフタイトル
ax2.set_xlabel('Time [s]') # x軸ラベル
ax2.set_ylabel('Displacement [m]') # y軸ラベル
ax2.grid(True) # グリッド線を表示
# y=0の位置に黒の破線を引く (つり合いの位置の目安)
ax2.axhline(0, color='black', lw=0.5, ls='--')
# アニメーションで動かす線オブジェクトを空で初期化
line_graph, = ax2.plot([], [], lw=1)
-
ax2.set_xlim(0, t_max): 右側のグラフ(ax2)のx軸（時間軸）の範囲を0からt_maxまでに設定 -
ax2.set_ylim(np.min(x)*1.1, np.max(x)*1.1): y軸（変位軸）の範囲を設定- シミュレーション結果
xの最小値と最大値を取得し (np.min(x),np.max(x)を使用)、グラフが見やすくなるように少し余裕を持たせて（1.1倍して）範囲を決定 - パラメータを変えて振動の幅が変わっても、グラフが適切に表示されやすくなる
- シミュレーション結果
-
line_graph, = ax2.plot([], [], lw=1): 左側のグラフと同様に、アニメーションで変位の履歴を描画するための空の線オブジェクトline_graphを作成
7. アニメーション用関数の定義
アニメーションを作成するために、FuncAnimation に渡す2つの重要な関数 init と update を定義。
# 初期化関数: アニメーション開始時に呼ばれる
def init():
# 各線オブジェクトのデータを空にする
line_mass.set_data([], [])
line_graph.set_data([], [])
# 更新対象のオブジェクトを返す (blit=True のため)
return line_mass, line_graph
# 更新関数: 各フレーム描画時に呼ばれる
def update(frame):
# frame: 現在のフレーム番号 (0, 5, 10, ... と変化)
# --- 左側のグラフ更新 ---
# 現在のフレームに対応する変位を取得
pos = x[frame]
# おもりの位置を表す線(点)のデータを更新
# x座標は常に1, y座標は天井(0)から現在の位置(pos)まで
line_mass.set_data([1, 1], [0, pos])
# --- 右側のグラフ更新 ---
# 時刻0から現在のフレームまでの時間と変位のデータを取得
current_t = t[:frame+1] # フレーム番号までの時間
current_x = x[:frame+1] # フレーム番号までの変位
# 変位履歴グラフのデータを更新
line_graph.set_data(current_t, current_x)
# 更新対象のオブジェクトを返す (blit=True のため)
return line_mass, line_graph
-
def init()::- アニメーションの開始時、またはアニメーションがループするときに最初に呼び出される関数
- グラフ上のアニメーション要素（ここでは
line_massとline_graph）を初期状態（何も描画されていない状態）に戻す役割を持つ -
line.set_data([], [])は、線オブジェクトlineのxデータとyデータを空のリストに設定することで、線を見えなくする -
return line_mass, line_graph:blit=Trueを使う場合（後述）、init関数とupdate関数は更新された描画オブジェクト（Artist）の イテラブル（タプルやリストなど） を返す必要がある
-
def update(frame)::- アニメーションの各フレーム（コマ）を描画するために呼び出される関数
- 引数
frameには、現在のフレーム番号（FuncAnimationのframes引数で指定した値が順番に入ってくる。- このコードでは
0, 5, 10, 15, ...）が渡されます。このframeを使って、シミュレーション結果の配列xや時間配列tから対応するデータを取り出す
- このコードでは
-
左側の更新:
-
pos = x[frame]: 現在のフレーム番号frameに対応する変位x[frame]をposに代入 -
line_mass.set_data([1, 1], [0, pos]):line_massオブジェクトのデータを更新 - x座標は常に
1（固定）、y座標は天井の0から現在の変位posまで、というデータを設定することで、おもりの位置を表現
-
-
右側の更新:
-
current_t = t[:frame+1]とcurrent_x = x[:frame+1]: 時間tと変位xの配列を、最初から現在のフレームframeまで スライス ([:frame+1]) して取得-
frameまで含めるために+1する
-
-
line_graph.set_data(current_t, current_x):line_graphオブジェクトのデータを更新- これにより、時間が進むにつれてグラフが過去の履歴を含めて右に伸びていくように見える
-
-
return line_mass, line_graph: 更新した描画オブジェクトを返す