粘弾性をいじって質点挙動をPythonアニメーションで可視化したい

Last updated at 2025-05-02Posted at 2025-05-01

初めて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) と $ で囲んでいます。これにより、添え字 ($x_0$) や単位 ($\mathrm{m}$) をきれいに表示できます。（LaTeXの知識がなくても、見よう見まねで使えます）
- %.2f: これは文字列フォーマット指定子で、対応する変数の値を小数点以下2桁までの浮動小数点数として埋め込みます。
- \n: 改行文字です。
- % (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: 更新した描画オブジェクトを返す

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up