Pythonで運動方程式を解く

# はじめに
私はプログラマではありませんが，設計関係の仕事をしている関係上，色々な計算を行うことが多く，これを簡単なプログラムを作って行っています．言語はFortran90，コンパイラはGNU gfortran主体ですが，１年半ほど前からPython3をいじり始めました．なかなかおもしろく，最近はだいぶ慣れてきて，Pythonで仕事用の計算を行うことも多くなってきました．
そこで，Python利用の練習の一環として，運動方程式を解くプログラムを作成してみました．とはいっても元ネタはFortranで作ったものをPythonに直したというのが本当のところですが，Python独自のScipyの利用にも挑戦してみました．

プログラミング環境は以下のとおり．

+ マシン： MacBook Pro (Retina, 13-inch, Mid 2014)
+ OS：　OS X El Capitan Version 10.11.5
+ 言語： Python 3.4.3 (含む Numpy, Scipy, matplotlib)

# 解くべき運動方程式

$$
m\cdot\ddot{u}(t)+c\cdot\dot{u}(t)+k\cdot u(t)=f(t)
$$

上式は，よく知られた強制振動を受ける質点の運動方程式です．
ここで紹介するプログラムでは，$f(t)$が観測された地震加速度による外力である場合を考えます．

この方程式の数値解法として有名なのは，

+ 線形加速度法（Newmark $\beta$法）
+ Runge-Kutta法
+ 周波数応答関数とFFTを用いた方法

などです．ここでは上記に加え，Scipyのodeintを用いた方法にも挑戦してみました．


# 入力データ
入力データとして，観測された地震加速度を用います．
データ書式は以下のとおりで，

+ 一行目：コメント
+ 二行目：サンプリングピッチ（秒）
+ 三行目：データ数
+ 四行目以降：加速度値（ここでは全３００００個）

となっています．

```dat_acc.csv
2011/03/11 14:46:00 TCGH16 EW2 Max.acc=1196.697
dt,0.010
ndata,30000
-1.890
-1.888
-1.888
......
```

入力データは加速度 $a(t)$ なので，外力 $f(t)$ は，

$$
f(t)=-m \cdot a(t)
$$

の形でプログラム内で扱われます．


# プログラムソース
地震動の解析を行う場合，入力加速度のサンプリングピッチは0.01秒，計算したい系の最小固有周期は0.02秒という場合が多いと思います．このケースを考える場合，もし数値積分の時間刻みを加速度のサンプリングピッチ0.01秒に合わせると，系の固有周期が小さい場合，数値積分の安定条件を満足できない場合が出てきます．実際，Runge−Kutta法では固有周期0.02秒の場合，積分時間ピッチ0.01秒では安定条件を満足できませんので，プログラム内では積分ピッチは０．００２秒（加速度のサンプリングピッチの1/5）としています．この時間分割は，プログラム内の変数md(=5)で実現しています．
線形加速度法では積分ピッチ0.01秒でも0.02秒の固有周期の場合の応答計算は可能ですが，Runge-Kutta法と計算時間を比較するため，積分ピッチはRunge−Kutta法の場合と同じに設定しています．

## 線形加速度法
以下に示す Newmark $\beta$法による次ステップ応答値計算式において，$\beta=1/6$ とおくことにより，線形加速度法の式が得られます．($\beta=1/4$ と置けば平均加速度法)
$$
\begin{align}
\ddot{u}(t+\Delta t)=&
\cfrac{f(t+\Delta t)-c\cdot\left(\dot{u}(t)+\cfrac{\Delta t}{2}\cdot\ddot{u}(t)\right)-k\cdot\left(u(t)+\Delta t\cdot\dot{u}(t)+\left(\cfrac{1}{2}-\beta\right)\cdot(\Delta t)^2\cdot\ddot{u}(t)\right)}
{m+c\cdot\cfrac{\Delta t}{2}+k\cdot\beta\cdot(\Delta t)^2} \\\
\dot{u}(t+\Delta t)=&
\dot{u}(t)+\frac{\Delta t}{2}\cdot\ddot{u}(t)+\frac{\Delta t}{2}\cdot\ddot{u}(t+\Delta t) \\\
u(t+\Delta t)=&
u(t)+\Delta t\cdot\dot{u}(t)+\left(\cfrac{1}{2}-\beta\right)\cdot(\Delta t)^2\cdot\ddot{u}(t)+\beta\cdot(\Delta t)^2\cdot\ddot{u}(t+\Delta t)
\end{align}
$$

応答値を求める場合，通常外力$f$，質点系の固有周期$T$，減衰定数$h$が与えられるので，プログラム内では，$m$，$k$，$c$ は以下のようにおきます．

$$
\begin{equation*}
m=1 \qquad k=\frac{4\pi^2 m}{T^2} \qquad c=2h\sqrt{k m}
\end{equation*}
$$

```py3:py_eng_acc_nmk.py
import numpy as np
from matplotlib import pyplot as plt
import time

def nmk(beta,am,ac,ak,dt,ff,acc1,vel1,dis1):
    acc2=(ff-ac*(vel1+0.5*dt*acc1)-ak*(dis1+dt*vel1+(0.5-beta)*dt*dt*acc1))/(am+ac*0.5*dt+ak*beta*dt*dt)
    vel2=vel1+0.5*dt*(acc1+acc2)
    dis2=dis1+dt*vel1+(0.5-beta)*dt*dt*acc1+beta*dt*dt*acc2
    return acc2,vel2,dis2

def INPDAT():
    fnameR='dat_acc.csv'
    fin=open(fnameR,'r')
    text=fin.readline()
    text=fin.readline();text=text.strip();text=text.split(',');dt=float(text[1])
    text=fin.readline();text=text.strip();text=text.split(',');ndata=int(text[1])
    wave=[]
    for i in range(0,ndata):
        text=fin.readline();text=text.strip();wave=wave+[float(text)]
    fin.close()
    return dt,ndata,wave

#Main routine
start=time.time()
#Definition of time and acceleration
dt,ndata,wave=INPDAT()
t=np.arange(0,dt*ndata,dt)
ainp = np.array(wave)

#Parameter setting
h=0.05
T=0.3
am=1.0
ak=4*np.pi**2*am/T**2
ac=2*h*np.sqrt(ak*am)
f=-am*ainp

#Numerical integration
md=5
ddt=dt/float(md)
beta=1.0/6.0
rac=np.zeros(len(t))
vel=np.zeros(len(t))
dis=np.zeros(len(t))
acc1=0.0
vel1=0.0
dis1=0.0
for i in range(0,len(t)):
    f1=0.0
    if 1<=i: f1=f[i-1]
    f2=f[i]
    for k in range(0,md):
        fm=f1+(f2-f1)/float(md)*float(k)
        acc2,vel2,dis2=nmk(beta,am,ac,ak,ddt,fm,acc1,vel1,dis1)
        acc1=acc2
        vel1=vel2
        dis1=dis2
    rac[i]=acc1+ainp[i]
    vel[i]=vel1
    dis[i]=dis1

#Show the results
#for i in range(0,len(t)):
#    print('{0:10.3f}{1:10.3f}{2:10.3f}{3:10.3f}{4:10.3f}'.format(t[i],ainp[i],acc[i],vel[i],dis[i]))
print('{0:10.3f}{1:10.3f}'.format(np.max(ainp),np.min(ainp)))
print('{0:10.3f}{1:10.3f}'.format(np.max(rac),np.min(rac)))
print('{0:10.3f}{1:10.3f}'.format(np.max(vel),np.min(vel)))
print('{0:10.3f}{1:10.3f}'.format(np.max(dis),np.min(dis)))
print('{0:10.3f}(sec)'.format(time.time()-start))

#Show drawings
plt.subplot(2, 1, 1)
plt.plot(t,ainp, 'k-')
plt.title('Acceleration')
plt.ylabel('Acceleration (gal)')

plt.subplot(2, 1, 2)
plt.plot(t,rac, 'r-')
plt.xlabel('Time (s)')
plt.ylabel('Responce acceleration (gal)')

plt.show()
```

## Runge-Kutta法
Runge-Kutta法では，以下に示すように，二階微分方程式を，一階連立微分方程式に置き換えます．

$$
\begin{equation*}
\begin{cases}
\quad f_a(t,x,y)=\dot{x}=\cfrac{dx}{dt}=y \\
\quad \\\
\quad f_b(t,x,y)=\ddot{x}=\cfrac{dy}{dt}=\cfrac{f(t)}{m}-\cfrac{c}{m}\cdot y-\cfrac{k}{m}\cdot x
\end{cases}
\end{equation*}
$$

Runge-Kutta 法による計算手順は，以下に示すとおり．ここに $\Delta t$ は時間増分です．
変数 $t$, $x$, $y$ のそれぞれに初期値を与えることが必要です．

$$
\begin{equation*}
\begin{cases}
&k_{a1}=f_a(t_n,x_n,y_n) \\\
&k_{a2}=f_a(t_n+\Delta t/2,x_n+\Delta t/2\cdot k_{a1},y_n+\Delta t/2\cdot k_{a1}) \\\
&k_{a3}=f_a(t_n+\Delta t/2,x_n+\Delta t/2\cdot k_{a2},y_n+\Delta t/2\cdot k_{a2}) \\\
&k_{a4}=f_a(t_n+\Delta t,x_n+\Delta t\cdot k_{a3},y_n+\Delta t\cdot k_{a3}) \\\
&x_{n+1}=x_n+\Delta t/6\cdot (k_{a1}+2k_{a2}+2k_{a3}+k_{a4})
\end{cases}
\end{equation*}
$$

$$
\begin{equation*}
\begin{cases}
&k_{b1}=f_b(t_n,x_n,y_n) \\\
&k_{b2}=f_b(t_n+\Delta t/2,x_n+\Delta t/2\cdot k_{b1},y_n+\Delta t/2\cdot k_{b1}) \\\
&k_{b3}=f_b(t_n+\Delta t/2,x_n+\Delta t/2\cdot k_{b2},y_n+\Delta t/2\cdot k_{b2}) \\\
&k_{b4}=f_b(t_n+\Delta t,x_n+\Delta t\cdot k_{b3},y_n+\Delta t\cdot k_{b3}) \\\
&y_{n+1}=y_n+\Delta t/6\cdot (k_{b1}+2k_{b2}+2k_{b3}+k_{b4})
\end{cases}
\end{equation*}
$$

$$
\begin{equation*}
\begin{cases}
\quad \ddot{x}=\cfrac{f(t)}{m}-\cfrac{c}{m}\cdot y_{n+1}-\cfrac{k}{m}\cdot x_{n+1} & \text{(acceleration)}\\\
\quad \dot{x}=y_{n+1} & \text{(velocity)}\\\
\quad x=x_{n+1}       & \text{(displacement)}
\end{cases}
\end{equation*}
$$

プログラム内での，$m$，$k$，$c$ の定義は，線形加速度法と同様に以下のようにします．

$$
\begin{equation*}
m=1 \qquad k=\frac{4\pi^2 m}{T^2} \qquad c=2h\sqrt{k m}
\end{equation*}
$$

```py3:py_eng_acc_rng.py
import numpy as np
from matplotlib import pyplot as plt
import time

def func(am,ac,ak,ff,x,y):
    dxdt=y
    dydt=ff/am-ac/am*y-ak/am*x
    return dxdt,dydt

def INPDAT():
    fnameR='dat_acc.csv'
    fin=open(fnameR,'r')
    text=fin.readline()
    text=fin.readline();text=text.strip();text=text.split(',');dt=float(text[1])
    text=fin.readline();text=text.strip();text=text.split(',');ndata=int(text[1])
    wave=[]
    for i in range(0,ndata):
        text=fin.readline();text=text.strip();wave=wave+[float(text)]
    fin.close()
    return dt,ndata,wave

#Main routine
start=time.time()
#Definition of time and acceleration
dt,ndata,wave=INPDAT()
ainp=np.array(wave)
t=np.arange(0,dt*ndata,dt)

#Parameter setting
h=0.05
T=0.3
am=1.0
ak=4*np.pi**2*am/T**2
ac=2*h*np.sqrt(ak*am)
f=-am*ainp

#Numerical integration
md=5
ddt=dt/float(md)
rac=np.zeros(len(t))
vel=np.zeros(len(t))
dis=np.zeros(len(t))
acc1=0.0
vel1=0.0
dis1=0.0
for i in range(0,len(t)):
    f1=0
    if 1<=i: f1=f[i-1]
    f2=f[i]
    x1=dis1
    y1=vel1
    for k in range(0,md):
        fm=f1+(f2-f1)/float(md)*float(k)
        kx1,ky1=func(am,ac,ak,fm,x1,y1)
        kx2,ky2=func(am,ac,ak,fm,x1+ddt*kx1/2.0,y1+ddt*ky1/2.0)
        kx3,ky3=func(am,ac,ak,fm,x1+ddt*kx2/2.0,y1+ddt*ky2/2.0)
        kx4,ky4=func(am,ac,ak,fm,x1+ddt*kx3,y1+ddt*ky3)
        x2=x1+ddt/6.0*(kx1+2.0*kx2+2.0*kx3+kx4)
        y2=y1+ddt/6.0*(ky1+2.0*ky2+2.0*ky3+ky4)
        x1=x2
        y1=y2
    acc1=f2/am-ac/am*y1-ak/am*x1
    vel1=y1
    dis1=x1
    rac[i]=acc1+ainp[i]
    vel[i]=vel1
    dis[i]=dis1

#Show the results
#for i in range(0,len(t)):
#    print('{0:10.3f}{1:10.3f}{2:10.3f}{3:10.3f}{4:10.3f}'.format(t[i],ainp[i],acc[i],vel[i],dis[i]))
print('{0:10.3f}{1:10.3f}'.format(np.max(ainp),np.min(ainp)))
print('{0:10.3f}{1:10.3f}'.format(np.max(rac),np.min(rac)))
print('{0:10.3f}{1:10.3f}'.format(np.max(vel),np.min(vel)))
print('{0:10.3f}{1:10.3f}'.format(np.max(dis),np.min(dis)))
print('{0:10.3f}(sec)'.format(time.time()-start))

#Show drawings
plt.subplot(2, 1, 1)
plt.plot(t,ainp, 'k-')
plt.title('Acceleration')
plt.ylabel('Acceleration (gal)')

plt.subplot(2, 1, 2)
plt.plot(t,rac, 'r-')
plt.xlabel('Time (s)')
plt.ylabel('Responce acceleration (gal)')

plt.show()
```

## 周波数応答関数とFFTを用いた方法
運動方程式は，以下の形式とします．

$$
\begin{equation}
\ddot{u}(t)+2\cdot h\cdot\omega_0\cdot\dot{u}(t)+\omega_0^2\cdot u(t)=-a(t)
\end{equation}
$$

ここに， $u$ ：質点の相対変位， $a$ ：地盤加速度， $h$ ：系の減衰定数， $\omega_0$ ：系の固有円振動数です．
$\omega_0$と系の固有周期$T$の関係は，以下のとおり．

$$
\omega_0=\cfrac{2 \pi}{T}
$$

地盤加速度および質点の変位・速度・加速度をフーリエ逆変換を用いて表すと，

$$
\begin{align}
&a(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}A(\omega) e^{i\omega t}d\omega \\\
&u(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}U(\omega) e^{i\omega t}d\omega \\\
&\dot{u}(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}i\omega\cdot U(\omega) e^{i\omega t}d\omega \\\
&\ddot{u}(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}-\omega^2\cdot U(\omega) e^{i\omega t}d\omega
\end{align}
$$

上式群を，運動方程式に代入することにより，
$$
\begin{align}
&U(\omega)=\frac{A(\omega)}{(\omega^2-\omega_0^2)-2 h \omega_0 \omega\cdot i} \\\
&u(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{A(\omega)}{(\omega^2-\omega_0^2)-2 h \omega_0 \omega\cdot i} e^{i\omega t}d\omega \\\
&\dot{u}(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{i\omega A(\omega)}{(\omega^2-\omega_0^2)-2 h \omega_0 \omega\cdot i} e^{i\omega t}d\omega \\\
&\ddot{u}(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{-\omega^2 A(\omega)}{(\omega^2-\omega_0^2)-2 h \omega_0 \omega\cdot i} e^{i\omega t}d\omega
\end{align}
$$

として，応答値である変位・速度・加速度が表現できます．
ここで，

$$
\begin{equation}
H(\omega)=\frac{1}{(\omega^2-\omega_0^2)-2 h \omega_0 \omega\cdot i}
\end{equation}
$$

は伝達関数あるいは周波数応答関数と呼ばれている．

応答値の計算手順は以下のとおり．

+ 入力加速度$a$の複素フーリエ係数 $A(\omega)$ の算出
+ 伝達関数 $H(\omega)$ の算出
+ $A(\omega)\cdot H(\omega)$ のフーリエ逆変換による変位時刻歴 $u(t)$ の算出
+ $i\omega\cdot A(\omega)\cdot H(\omega)$ のフーリエ逆変換による速度時刻歴 $\dot{u}(t)$ の算出
+ $-\omega^2\cdot A(\omega)\cdot H(\omega)$ のフーリエ逆変換による加速度時刻歴 $\ddot{u}(t)$ の算出

実際の計算では高速フーリエ変換を用いるため，以下に注意します．

+ データ数は有限で $N$ 個（ $N$ は 2 の累乗）とし，データの時間間隔を $\Delta t$ とする．
+ 計算される振動数 $f_k$ および円振動数 $\omega$ は以下のとおりとなる．

$$
\begin{equation}
f_k=\frac{k}{N\cdot\Delta t} \qquad \omega_k=2\pi f_k \qquad k=0\sim N/2
\end{equation}
$$

+ ここでデータは $N$ 個あるが実態としての振動数は $N/2+1$ 個しか定義できないため，
$N/2+1$ 個目のデータを中心に対称形に $N/2+2\sim N$ 個目の円振動数を定義しておく．
この際，係数は $N/2+1$ 個目のデータを中心に対称形に格納するが，虚数部は符号を逆転させることに注意する．

```py3:py_eng_acc_fft.py
import numpy as np
from matplotlib import pyplot as plt
import time

def INPDAT():
    fnameR='dat_acc.csv'
    fin=open(fnameR,'r')
    text=fin.readline()
    text=fin.readline();text=text.strip();text=text.split(',');dt=float(text[1])
    text=fin.readline();text=text.strip();text=text.split(',');ndata=int(text[1])
    wave=[]
    for i in range(0,ndata):
        text=fin.readline();text=text.strip();wave=wave+[float(text)]
    fin.close()
    return dt,ndata,wave

#Main routine
start=time.time()
#Definition of time and acceleration
dt,ndata,wave=INPDAT()
tt=np.arange(0,dt*ndata,dt)
ainp = np.array(wave)

h=0.05
T=0.3
w0=2*np.pi/T

nn=2
while nn<ndata:
    nn=nn*2
xx=np.zeros(nn)
xx[0:ndata]=xx[0:ndata]+ainp[0:ndata]
sp=np.fft.fft(xx,nn)/nn

k=np.arange(0,int(nn/2)+1,dtype=np.float)
fk=k/nn/dt
wk1=2*np.pi*fk
hh1=1.0/((wk1**2-w0**2)-(2*h*w0*wk1)*1j)
_w=wk1[::-1];wk2=_w[1:int(nn/2)];wk=np.r_[wk1,wk2]
_h=hh1[::-1];hh2=_h[1:int(nn/2)];hh=np.r_[hh1,hh2]
hh.imag[int(nn/2)+1:nn]=-hh2.imag[0:int(nn/2)]

jdis=sp*hh
jvel=1j*wk*sp*hh
jacc=-wk**2*sp*hh

dis= np.real(np.fft.ifft(jdis*nn,nn))
vel=-np.imag(np.fft.ifft(jvel*nn,nn))
acc= np.real(np.fft.ifft(jacc*nn,nn))
rac=acc+xx

#Show the results
#for i in range(0,len(t)):
#    print('{0:10.3f}{1:10.3f}{2:10.3f}{3:10.3f}{4:10.3f}'.format(t[i],ainp[i],acc[i],vel[i],dis[i]))
print('{0:10.3f}{1:10.3f}'.format(np.max(ainp),np.min(ainp)))
print('{0:10.3f}{1:10.3f}'.format(np.max(rac),np.min(rac)))
print('{0:10.3f}{1:10.3f}'.format(np.max(vel),np.min(vel)))
print('{0:10.3f}{1:10.3f}'.format(np.max(dis),np.min(dis)))
print('{0:10.3f}(sec)'.format(time.time()-start))

#Show drawings
plt.subplot(2, 1, 1)
plt.plot(tt,ainp,'k-')
plt.ylabel('Input acc. (gal)')

plt.subplot(2, 1, 2)
plt.plot(tt,rac[0:ndata],'r-')
plt.xlabel('Time (s)')
plt.ylabel('Res. acc. (gal)')

plt.show()
```

## Scipyのodeintを用いた方法
Scipy の odeint を利用するため，Runge-Kutta法と同様に，二階微分方程式を，一階連立微分方程式に置き換えます．

$$
\begin{equation*}
\begin{cases}
\quad f_a(t,x,y)=\dot{x}=\cfrac{dx}{dt}=y \\
\quad \\\
\quad f_b(t,x,y)=\ddot{x}=\cfrac{dy}{dt}=\cfrac{f(t)}{m}-\cfrac{c}{m}\cdot y-\cfrac{k}{m}\cdot x
\end{cases}
\end{equation*}
$$

```py3:py_eng_acc_ode.py
import numpy as np
from scipy import integrate
from scipy import interpolate
from matplotlib import pyplot as plt
import time

def func(x,t,am,ac,ak,dt,f):
    def ff(t):
        return f[int(t/dt)]
    dxdt=x[1]
    dydt=ff(t)/am-ac/am*x[1]-ak/am*x[0]
    return [dxdt,dydt]

def INPDAT():
    fnameR='dat_acc.csv'
    fin=open(fnameR,'r')
    text=fin.readline()
    text=fin.readline();text=text.strip();text=text.split(',');dt=float(text[1])
    text=fin.readline();text=text.strip();text=text.split(',');ndata=int(text[1])
    wave=[]
    for i in range(0,ndata):
        text=fin.readline();text=text.strip();wave=wave+[float(text)]
    fin.close()
    return dt,ndata,wave

#Main routine
start=time.time()
#Definition of time and acceleration
md=5
dt,ndata,wave=INPDAT()
_w=np.array(wave)
_t=np.arange(0,dt*ndata,dt)
f1 = interpolate.interp1d(_t, _w)
ddt=dt/float(md)
t=np.arange(0,dt*ndata-1,ddt)
ainp = f1(t)
dt=ddt

#Parameter setting
h=0.05
T=0.3
am=1.0
ak=4*np.pi**2*am/T**2
ac=2*h*np.sqrt(ak*am)
f=-am*ainp

#Numerical integration
x0=np.array([0,0])
result=integrate.odeint(func,x0,t,args=(am,ac,ak,dt,f))
dis=result[:,0]
vel=result[:,1]
acc=f/am-ac/am*vel-ak/am*dis
rac=acc+ainp

#Show the results
#for i in range(0,len(t)):
#    print('{0:10.3f}{1:10.3f}{2:10.3f}{3:10.3f}{4:10.3f}'.format(t[i],ainp[i],acc[i],vel[i],dis[i]))
print('{0:10.3f}{1:10.3f}'.format(np.max(ainp),np.min(ainp)))
print('{0:10.3f}{1:10.3f}'.format(np.max(rac),np.min(rac)))
print('{0:10.3f}{1:10.3f}'.format(np.max(vel),np.min(vel)))
print('{0:10.3f}{1:10.3f}'.format(np.max(dis),np.min(dis)))
print('{0:10.3f}(sec)'.format(time.time()-start))

#Show drawings
plt.subplot(2, 1, 1)
plt.plot(t,ainp, 'k-')
plt.title('Acceleration')
plt.ylabel('Acceleration (gal)')

plt.subplot(2, 1, 2)
plt.plot(t,rac, 'r-')
plt.xlabel('Time (s)')
plt.ylabel('Responce acceleration (gal)')

plt.show()
```

# 実行結果
以下に，ある条件で計算を実行した結果を示します．
加速度単位はgal，速度単位はkine，変位単位はcmです．

```
（出力書式）
$ python3 プログラム名
  入力加速度最大値  入力加速度最小値
  加速度応答最大値  加速度応答最小値
  速度応答最大値    速度応答最小値
  変位応答最大値    変位応答最小値
  計算時間

$ python3 py_eng_acc_nmk.py
  1196.697  -905.670
  4319.940 -4187.217
   183.217  -182.928
     9.514    -9.829
     2.072(sec)

$ python3 py_eng_acc_rng.py
  1196.697  -905.670
  4329.333 -4197.259
   183.099  -182.670
     9.522    -9.816
     2.533(sec)

$ python3 py_eng_acc_fft.py
  1196.697  -905.670
  4344.584 -4205.964
   185.899  -180.384
     9.558    -9.832
     1.613(sec)

$ python3 py_eng_acc_ode.py
  1196.697  -905.670
  4332.335 -4199.847
   183.121  -183.094
     9.535    -9.835
    30.134(sec)
```

# 感想
## Qiitaへは第一回目の投稿です．
数式入力は，通常のMarkdown記法でいいかと思ったら，そうでもありません．
alignやcaseを使う場合の数式末のバックスラッシュは，２つではなく，３つ連続で入れる必要があるようです．少なくともこの文書ではそうしています．

## Pythonプログラムを動かした感想は．．．
Scipyのodeintを除いては，まあまあの速度です．
Pythonでこの手の数値計算のプログラムを書く上では，

+ できるだけリストは使わずNumpyの配列を使う
+ できるだけforを使わない

というのが計算時間をいたずらに長くしないコツのようです．
その点，NumpyのFFTモジュールを用いた方法が，今回の比較の中では最もPythonに適した方法といえそうです．

Scipyのodeintでは，入力加速度の受け渡し方法がよく解っていませんが，時刻ベクトルの関数にしなければならないことは理解できるため，下のコード（def ff(t)の所）のような方法をとっています．ここで時間を食っているのかもしれません．

```
def func(x,t,am,ac,ak,dt,f):
    def ff(t):
        return f[int(t/dt)]
    dxdt=x[1]
    dydt=ff(t)/am-ac/am*x[1]-ak/am*x[0]
    return [dxdt,dydt]
```

なお，いずれの方法も，方程式を１回解くのに１秒以上かかっており，このプログラムを，系の固有振動数を変えながら100回とか200回解く，加速度応答スペクトル解析などに使うのは，気の長い人意外はやめたほうがいいと思われます．
そのような場合は，素直にFortranとかCのプログラムを使うべきでしょう．

## iPad Pythonistaでも使える
Scipyのodeintを用いたプログラム以外は，iOS用のPythonista（iOS用Python2.7開発環境）でもそのまま動きました．Python2.7といいながら，Python3のprint文がそのまま動くようです．ありがたや．また，Pythonistaには，Numpyやmatplotlibも含まれているので，色々遊べます．

手持ちのiPadに入れているPythonistaを用い，Macの場合と同じ入力での計算時間は以下のとおりでした．
（iPad mini 4 Wi-Fi 64GB, 言語： Pythonista version 2.0）

+ 線形加速度法:     9.7秒
+ Runge-Kutta法:  11.5秒
+ FFT法:          7.7秒

以上