単振り子 の運動をシミュレーションする．
一言に力学シミュレーションといっても，その方法は多様であり，いくつもの方法があるのだが，そのうちの一つを例としてあげる．シミュレーションを行うために必要な作業を
の3つのプロセスに分けて考える．
まず，一つ目の「モデル化」について．これは，例えば考える物理モデルからニュートンの運動方程式を導くプロセスのことである．
次に，二つ目の「数値計算」について．これは，Python プログラムを用いてコンピュータ上で運動を計算させるプロセスのことである．このプロセスでコンピュータは振り子の運動を与えられたモデル情報に基づいて計算し，その結果を数値で示す．例えば
シミュレーション開始0秒後，振り子は10 [rad]の位置にいる
シミュレーション開始1秒後，振り子は0 [rad]の位置にいる
シミュレーション開始2秒後，振り子は-10 [rad]の位置にいる
シミュレーション開始3秒後，振り子は0 [rad]の位置にいる
...
という情報が，
| 時間 | 位置 |
|---|---|
| 0 | 10 |
| 1 | 0 |
| 2 | -10 |
| 3 | 0 |
のような行列となって生成される．
最後に，「視覚化」について．「数値計算」によって得られた数値のデータは，あくまで数字の羅列であり，これだけでは物体の運動を理解しようとするのは難しい．そこで，運動を直感的に理解できるようにするために，この数値データからグラフやアニメーションを生成する．
図のような，よくある 単振り子の力学モデル を考える．振り子の角度を$\theta$，振り子の長さを$l$，振り子の先端のおもりの質量を$m$，鉛直下向きにはたらく重力の重力加速度を$g$としている．
おもりの運動方向にはたらく力は \begin{align} F = -mg\sin\theta \end{align} である．おもりの運動方向の加速度は$l\ddot{\theta}$なので，ニュートンの運動方程式 ($ma=F$) は， \begin{align} m l\ddot{\theta}= -mg\sin\theta \end{align} したがって， \begin{align} \ddot{\theta} = -\frac{g}{l}\sin\theta \end{align} である．
さて，単振り子のモデルを考える上で，時間の経過とともに値が変化する 変数 と，時間の経過に関わらず値が変化しない 定数 が存在した．$\theta$ が前者であり，$m$，$l$，$g$が後者である． 1. モデル化 のプロセスで達成することは，運動方程式を変数 について整理し，「（変数の時間の微分項）$=$（それ以外の項）」という式を得ることである．
1. モデル化 では変数$\theta$についての運動方程式を導いたが， 2. 数値計算 では，これを用いて 常微分方程式 を構築する．
まず，$\theta$自身とその一階微分である$\dot{\theta}$を縦に並べた列ベクトルをつくり，これを$s$とする．（GiuHub上で行列が上手く表示されないのは仕方のないことなのかな？）
\begin{align}
s =
\begin{bmatrix}
\theta \\ \dot{\theta}
\end{bmatrix}
\end{align}
そして，$s$の時間の一階微分を式で表す．ここで，$\ddot{\theta}$は先ほど求めた運動方程式によって，$\ddot{\theta}$を含まない形で表すことができる．
\begin{align}
\frac{d}{dt} s =
\begin{bmatrix}
\dot{\theta} \\ \ddot{\theta}
\end{bmatrix} =
\begin{bmatrix}
\dot{\theta} \\ -\frac{g}{l}\sin\theta
\end{bmatrix}
\end{align}
こうして得られた常微分方程式は，数値計算に適した形となっている．具体的には，$\theta = $ s[0] ，$\dot{\theta} = $ s[1] とすることで，
def odefunc(s, t):
theta = s[0]
dtheta = s[1]
ddtheta = -g/l*sin(theta) # <- Equation of motion
return np.r_[dtheta, ddtheta]のようにして常微分方程式をPythonの関数として記述したうえで
s = odeint(odefunc, s_init, t)を実行することで常微分方程式を数値的に解くことができる．ただし odeint は SciPy ライブラリの integrate クラスで定義される，常微分方程式を解く関数である．
In [1]:
import numpy as np
from scipy.integrate import odeint
from math import sin
''' constants '''
m = 1 # mass of the pendulum [kg]
l = 1 # length of the pendulum [m]
g = 10 # Gravitational acceleration [m/s^2]
''' time setting '''
t_end = 10 # simulation time [s]
t_fps = 50 # frame per second. This value means smoothness of produced graph and animation
t_step = 1/t_fps
t = np.arange(0, t_end, t_step)
''' initial value '''
theta_init = 0 # initial value of theta [rad]
dtheta_init = 1 # initial value of dot theta [rad/s]
s_init = np.r_[theta_init, dtheta_init]
def odefunc(s, t):
theta = s[0]
dtheta = s[1]
ddtheta = -g/l*sin(theta) # <- Equation of motion
return np.r_[dtheta, ddtheta]
s = odeint(odefunc, s_init, t)
print('ODE calculation finished.')
以上が数値計算を行うプログラムである．
print(np.c_[t, s])を実行することで，
| 時間 | 振子の角度 | 振子の角速度 |
|---|---|---|
| 0. | 0. | 1. |
| ... | ... | ... |
| 9.98 | -0.018 | 0.998 |
のように，それぞれの時間で振子の角度がどのようになっていて，また振子の角速度がどのようになっているのかを確認することができる．
In [2]:
print(np.c_[t, s])
In [3]:
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from math import cos
def animfunc(s, t):
''' Create mp4 movie file of a pendulum '''
plt.close()
fig = plt.figure()
plt.axis('scaled')
plt.xlim(-1, 1)
plt.ylim(-1.5, .5)
plt.grid('on')
draw_ceiling, = plt.plot([-2, 2], [0, 0], c='k', lw=2)
draw_pendulum, = plt.plot([], [], lw=4, c='b')
draw_mass, = plt.plot([], [], lw=2, marker='o', ms=20, mew=4, mec='b', mfc='c')
indicate_time = plt.text(-0.3, 0.25, [], fontsize=12)
def update_figure(i):
''' Set data of each movie frame '''
mass_x = l*sin(s[i, 0])
mass_y = - l*cos(s[i, 0])
pendlum_x = [0, mass_x]
pendlum_y = [0, mass_y]
draw_pendulum.set_data(pendlum_x, pendlum_y)
draw_mass.set_data(mass_x, mass_y)
indicate_time.set_text('t = {0:4.2f} [s]'.format(t[i]))
''' Create a movie file '''
line_ani = animation.FuncAnimation(fig, update_figure, frames=len(t))
line_ani.save('./pendulum.mp4', fps=t_fps)
print('pendulum.mp4 created')
animfuncを実行することで，pendulum.mp4のようなアニメーション動画が保存される．
In [4]:
animfunc(s, t)
一方，振子の動きをグラフに示したい場合は，例えば次のようにしてプロットを行えば良い．次のプログラムを実行することで，pendulum_graph.pngのようなイメージファイルが保存される．
In [5]:
plt.figure()
plt.plot(t, s[:, 0])
plt.xlabel('t [s]')
plt.ylabel('theta [rad]')
plt.savefig('pendulum_graph.png')
In [ ]: