　時間発展の方程式を数値計算するときに、物体の運動と物理量の時間変化を対応させてアニメーションで表示したいときがあります。例えば、単振り子のシュミレーションでは、振り子が運動している様子とその時の位相を並べてアニメーションで表示できれば、物理現象が視覚的に理解できます。

そこで、今回はmatplotlibのanimationを使ってそれを実現させるコードを紹介します。

matplotlib.animationの使い方

まずは、matplotlib.animationでsinカーブ上の点を動かすアニメーションの作り方を紹介します。

animate_sin.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

#グラフの体裁を整える
plt.rcParams['figure.figsize'] = (6,5)
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 16

#枠線の太さを指定する関数
def axes_set_linewidth(axes, t=1, b=1, r=1, l=1):
    axes.spines['top'].set_linewidth(t)
    axes.spines['bottom'].set_linewidth(b)
    axes.spines['right'].set_linewidth(r)
    axes.spines['left'].set_linewidth(l)

#変数設定
N = 100
x = np.linspace(-2 * np.pi, 2 * np.pi, N)
y = np.sin(x)

#アニメーション作成
fig, ax = plt.subplots()

ax.plot(x, y, color = 'blue', label = 'sin x')  #sinカーブをplot

axes_set_linewidth(ax, t=0, r=0, b=2, l=2)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_ylim(-1.2, 1.4)
ax.legend(frameon = False)
plt.subplots_adjust(top=0.9, bottom=0.2, right=0.9, left=0.2)

ims = []  #ここに1ステップごとのグラフを格納
for i in range(N):
    p = ax.plot(x[i], y[i], color = 'darkblue', marker = 'o', markersize = 10)
    ims.append(p)
    
ani = animation.ArtistAnimation(fig, ims, interval=100)  #ArtistAnimationでアニメーションを作成する。
ani.save('animate.gif', writer='imagemagick', dpi = 300)  #gifで保存

以下に保存したgifを載せておきます。

単振り子の運動と位相の同時アニメーション

半径Lの円弧上を質量mの質点が往復運動する振り子について考える。このときの質点の運動方程式は、
$$
mL \frac{d^2 \theta}{dt^2} = -mg \sin{\theta}
$$

となる。$\theta$が十分小さいとして$\sin\theta\simeq\theta$の近似を用いてこの微分方程式を解くと、一般解は次のようになる。
$$
\theta (t) = \theta_{max} \cos(\omega t + \alpha) ; ,\quad \omega= \sqrt{\frac{g}{L}}
$$
ここで、$\theta_{max}$は振幅、$\omega$は角振動数、$\alpha$は初期位相である。

この解をもとに以下のように、単振り子の簡単な数値計算を行います。

pendulum.py

#条件
alpha = 0
g = 9.8  #m/s2
L = 0.1  #m
omega = np.sqrt(g/L)
theta_max = np.pi / 6
dt = 0.01
t = np.arange(0, 2 + dt, dt)

#位相の計算
def theta(ts):
    return theta_max*np.cos(omega*ts + alpha)

#位置の計算
def theta_position(ts):
    x = L * np.sin(theta(ts))
    z = -L * np.cos(theta(ts))
    return x, z

それでは以下のコードを実行して、振り子の運動とその時の位相のアニメーションを並べて表示してみます。

pendulum_animate.py

fig, ax = plt.subplots(1,2, figsize = (14,6))
axes_set_linewidth(ax[0], t=0, b=0, r=0, l=0)
axes_set_linewidth(ax[1], t=0, b=2, r=0, l=2)

ax[0].set_xticks([])
ax[0].set_yticks([])
ax[0].set_xlim(-0.06, 0.06)
ax[0].set_ylim(-0.12, 0.02)
ax[0].vlines(0, -0.1, 0, linestyle = '--')
ax[0].hlines(0, -1, 1)

ax[1].plot(t, theta(t))
ax[1].set_xlim(0,2.2)
ax[1].set_ylim(-0.6,0.6)
ax[1].hlines(0,0,2.2, color='black', linestyle='--')
ax[1].set_xlabel('Time [s]')
ax[1].set_ylabel(r'$\theta$ [rad]')

plt.subplots_adjust(top=0.9, bottom=0.2, right=0.9, left=0.2, wspace = 0.3)

ims = []
for i in t:
    x, y = theta_position(i)
    p = ax[0].plot([0,x], [0,y], color = 'black') + ax[0].plot(x, y, marker = 'o', color = 'red', markersize = 15)\
        + ax[1].plot(i, theta(i), marker='o', color='darkblue', markersize=10)
    ims.append(p)

ani = animation.ArtistAnimation(fig, ims, interval=50, repeat=True)
ani.save('pendulum.gif', writer='imagemagick')

このようにax[0]とax[1]で連動したアニメーションを表示するときは、

p = ax[0].plot([0,x], [0,y], color = 'black') + ax[0].plot(x, y, marker = 'o', color = 'red', markersize = 15)\
        + ax[1].plot(i, theta(i), marker='o', color='darkblue', markersize=10)

のようにplotをプラスしていけばよいです。

最後に出力結果を示して終わりとしたいと思います。