LoginSignup
27
34

More than 5 years have passed since last update.

@damyarou

Pythonで運動方程式を解く

Last updated at Posted at 2016-06-11

はじめに

私はプログラマではありませんが，設計関係の仕事をしている関係上，色々な計算を行うことが多く，これを簡単なプログラムを作って行っています．言語は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*}
$$

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*}
$$

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$ 個目のデータを中心に対称形に格納するが，虚数部は符号を逆転させることに注意する．
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*}
$$

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秒

以上

27
34
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
27
34