In [1]:
import numpy as np
from scipy.integrate import odeint
from math import sin
In [2]:
import time
import math
import numpy
''' Sine function of math library '''
start = time.time()
for i in range(1000000):
math.sin(i)
elapsed_time = time.time() - start
print ('Sine calculation by math: {0:f} [sec]'.format(elapsed_time))
''' Sine function of numpy library '''
start = time.time()
for i in range(1000000):
numpy.sin(i)
elapsed_time = time.time() - start
print ('Sine calculation by numpy: {0:f} [sec]'.format(elapsed_time))
In [3]:
''' constants '''
m = 1 # mass of the pendulum [kg]
l = 1 # length of the pendulum [m]
g = 10 # Gravitational acceleration [m/s^2]
次に，シミュレーションを行う時間に関する定数を定義しておく． t_fps はグラフの描画やアニメーションの生成のために1秒を何分割するかを意味するパラメータである．
In [4]:
''' 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)
$\theta$自身とその一階微分である$\dot{\theta}$を縦に並べた列ベクトルをつくり，これを$s$とする．（GiuHub上で行列が上手く表示されないのは仕方のないことなのかな？） \begin{align} s = \begin{bmatrix} \theta \\ \dot{\theta} \end{bmatrix} \end{align}
常微分方程式$\dot{s} = f(s)$を解くために，初期値 $s_\mathrm{init}$ を与える．
In [5]:
''' 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]
np.r_[s_init, s_init]
Out[5]:
そして，$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) # <- This is the equation of motion
return np.r_[dtheta, ddtheta]のようにして常微分方程式をPythonの関数として記述したうえで
s = odeint(odefunc, s_init, t)を実行することで常微分方程式を数値的に解くことができる．ただし odeint は SciPy ライブラリの integrate クラスで定義される，常微分方程式を解く関数である．
In [6]:
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.')