# 数値計算

単振り子の運動について考える．

まずは数値計算に必要なライブラリのインポートから．[NumPy](http://www.numpy.org/)は科学計算を行うための基本的なライブラリである．[SciPy](https://www.scipy.org/)も科学計算を行う上で用いられるライブラリであるが，NumPyを用いて数学や工学などへの応用性のあるクラスが多い．

[math](https://docs.python.jp/3/library/math.html) はいつでも使えるデフォルトの数学関数系のライブラリである．


```python
import numpy as np
from scipy.integrate import odeint
from math import sin
```

## `math.sin` と `numpy.sin` について
`numpy` ライブラリの `sin` 関数を使わず， `math` ライブラリから `sin` 関数をインポートしたのは，計算を高速化させるためである．`math`ライブラリで定義される正弦関数と`numpy`ライブラリで定義される正弦関数には違いがある．次の例で示すように，スカラー値の正弦の計算は `numpy.sin` よりも`math.sin` の方が早い．


```python
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)) 
```

    Sine calculation by math: 0.193746 [sec]
    Sine calculation by numpy: 1.112725 [sec]


## `scipy.integrate.odeint` を用いた常微分方程式の計算

まず，運動方程式で用いる定数 ($m$，$l$，$g$) を定義しておく．


```python
''' 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秒を何分割するかを意味するパラメータである．


```python
''' 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}$ を与える．


```python
''' 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]
```




    array([0, 1, 0, 1])



そして，$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`](https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.integrate.odeint.html) は [SciPy](https://www.scipy.org/) ライブラリの `integrate` クラスで定義される，常微分方程式を解く関数である．


```python
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.')
```

    ODE calculation finished.


## 注意点：`odeint` 中では **行ベクトル** で解を扱う

振り子のモデリングにおいては，解$s = \begin{bmatrix} \theta & \dot{\theta} \end{bmatrix}^\mathrm{T}$は **列ベクトル** として扱っている．しかし，`odeint` 関数の引数の `s_init` や `odefunc`　の返り値は **行ベクトル** でなくてはならない．
