Python: 並列処理実装例（応答スペクトル計算２）

はじめに

先般、以下の記事を投稿した。

Python: 並列処理実装例（応答スペクトル計算）

これは、線形加速度法により１質点系の運動方程式を解き、応答スペクトルを作図するものであった。
その後、興味と練習を兼ねて、Runge-Kutta法とFFTによる解法を用いてスペクトル図を書くプログラムを作ってみた。
その結果、線形加速度法では、１質点系の固有周期が大きくなると、他の方法に比べて精度が著しく低下することが判明。対応策を考えた。

現在の環境は以下の通り。

MacBook Pro 14" (Apple M1 Pro, Memory 16GB)
macOS Sequoia 15.1
Python 3.13.0

加速度応答スペクトルの比較

以下に各方法による加速度応答スペクトル図を示す。線形加速度法とRunge-Kutta法におけるmdは、0.01秒ピッチで計測された入力加速度の、ピッチ分割数（0.01秒間の内挿点数）である。
Runge-Kutta法では、数値積分の安定条件を満足するため、md=5としている。線形加速度法では、安定条件の問題はないためmd=1としている。
以下の図より、線形加速度法では、周期が１秒程度より大きくなると、他の手法と比べ、応答加速度が過大に計算されていることがわかる。

線形加速度法（md=1）

Rumge-Kutta法（md=5）

FFT

加速度応答時刻歴の比較

何が起こっているか確認するため、いくつかの周期での、加速度応答時刻歴を出力してみた。
以下、Runge-Kutta法では、md=5で固定、線形加速度法ではmdの値を変化させている。

周期0.3秒（線形加速度法md=1）

周期10秒（線形加速度法md=1）

周期10秒（線形加速度法md=50）

加速度高騰履歴より、以下のことがわかる。

• Runge-Kutta法とFFT法では、いずれの周期でも、同様の加速度時刻歴が得られている。
• 線形加速度法では、周期が小さければ、他の方法と同等の加速度応答時刻歴が得られるが、周期が大きくなると、加速度が過大評価される。
• 線形加速度法の加速度過大評価に対する対策として、入力加速の線形内挿点数を多くすると、他手法と同等の加速度が得られるようになる。周期10秒の場合、内挿点数を50点とすれば、概ね同等の加速度応答時刻歴が得られる。
• このことから、線形加速度違法においては、計算を行う系の固有周期に応じて入力加速度の内挿点数（md）を増加させることが、対策として考えられる。

以上を踏まえ、入力加速度の内挿点数をmd=int(tp*5)+1として、周期20秒における加速度維持国歴の比較を以下に示す。ここに、tpは系の固有周期である。

いずれの方法の応答加速度も、同等の履歴を示す事が確認できた。

以下に、修正した方法（md=int(tp*5)+1）を適用した線形加速度法による、加速度応答スペクトルを示す。

うまくいっているようである。

計算時間比較

加速度応答スペクトルの計算時間については以下の通りである。

計算法	md	計算時間	備考
線形加速度法	1	5.862 sec	不採用
線形加速度法	int(tp*5)+1	43.447 sec sec	適用可能
Runge-Kutta法	5	38.515 sec sec	適用可能
FFT法	-	2.404 sec	最速

上表より、for文で時間ステップ毎に計算を行う、線形加速度法やRunge-Kutta法に比べ、ベクトルをFFTに放り込んで計算を行うFFT法が時間的には圧倒的に有利であることがわかる。

コード

線形加速度法による応答スペクトル計算（修正版）

import numpy as np
import time
from multiprocessing import Pool
import pandas as pd


def inpdat(fnameR):
    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=np.zeros(ndata,dtype=np.float64)
    for i in range(0,ndata):
        text=fin.readline();text=text.strip();wave[i]=float(text)
    fin.close()
    return dt,ndata,wave


def iacc(dt,acc):
    #Integral of acceleration
    nn=len(acc)
    vel=np.zeros(nn,dtype=np.float64)
    dis=np.zeros(nn,dtype=np.float64)
    for i in range(1,nn):
        vel[i]=vel[i-1]+(acc[i-1]+acc[i])*dt/2.0
        dis[i]=dis[i-1]+vel[i-1]*dt+(acc[i-1]/3.0+acc[i]/6.0)*dt*dt
    return vel,dis


def crac(dt,acc,vel,dis):
    #Amendment of base line
    accmax=np.max(np.abs(acc))
    nn=len(acc)
    tt=float(nn-1)*dt
    t=0.0
    for i in range(0,nn):
        dis[i]=dis[i]*(3.0*tt-2.0*t)*t*t
        t=t+dt
    dsum=(dis[0]+dis[nn-1])/2.0
    for i in range(1,nn-1):
        dsum=dsum+dis[i]
    dsum=dsum*dt
    a1=28.0/13.0/tt/tt*(2.0*vel[nn-1]-15.0/tt**5.0*dsum)
    a0=vel[nn-1]/tt-a1/2.0*tt
    t=0.0
    for i in range(0,nn):
        acc[i]=acc[i]-a0-a1*t
        t=t+dt
    acmax=np.max(np.abs(acc))
    coef=accmax/acmax
    acc=acc*coef
    return acc


def newmark(t,dt,ainp,tp):
    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

    #Parameter setting
    hh=0.05
    #tp=0.3
    am=1.0
    ak=4*np.pi**2*am/tp**2
    ac=2*hh*np.sqrt(ak*am)
    f=-am*ainp
    beta=1.0/6.0

    #Numerical integration
    md=int(tp*5)+1
    ddt=dt/float(md)
    rac=np.zeros(len(t),dtype=np.float64)
    vel=np.zeros(len(t),dtype=np.float64)
    dis=np.zeros(len(t),dtype=np.float64)
    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
    arac=np.max(np.abs(rac))
    avel=np.max(np.abs(vel))
    adis=np.max(np.abs(dis))
    return arac,avel,adis


def wrapper(args):
    return newmark(*args)


def calc():
    for nnn in [1,2,3]:
        match nnn:
            case 1: fnameR='dat_acc_TCGH16_EW2.txt'
            case 2: fnameR='dat_acc_TCGH16_NS2.txt'
            case 3: fnameR='dat_acc_TCGH16_UD2.txt'
        dt,ndata,acc=inpdat(fnameR)
        vel,dis=iacc(dt,acc) # Integral of acceleration
        ainp=crac(dt,acc,vel,dis) #Modification of base line
        t=np.arange(0,dt*ndata,dt) # time step
        tw=np.linspace(np.log10(0.02),np.log10(10),200)
        ttp=10**(tw) # period

        val=[]
        for tp in ttp:
            val.append([t,dt,ainp,tp])
        result=Pool(8).map(wrapper,val)

        result=np.array(result)
        df=pd.DataFrame({
            't': ttp,
            'acc': result[:,0],
            'vel': result[:,1],
            'dis': result[:,2]
        })
        print(df)
        match nnn:
            case 1: fnameW='df_n_res1.csv'
            case 2: fnameW='df_n_res2.csv'
            case 3: fnameW='df_n_res3.csv'
        df.to_csv(fnameW)


def main():
    start=time.time()
    calc()
    end=time.time()-start
    print('{0:10.3f}(sec)'.format(end))
    

#---------------
# Execute
#---------------
if __name__ == '__main__': main()

Rungte-Kutta法による応答スペクトル計算

import numpy as np
import time
from multiprocessing import Pool
import pandas as pd


def inpdat(fnameR):
    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=np.zeros(ndata,dtype=np.float64)
    for i in range(0,ndata):
        text=fin.readline();text=text.strip();wave[i]=float(text)
    fin.close()
    return dt,ndata,wave


def iacc(dt,acc):
    #Integral of acceleration
    nn=len(acc)
    vel=np.zeros(nn,dtype=np.float64)
    dis=np.zeros(nn,dtype=np.float64)
    for i in range(1,nn):
        vel[i]=vel[i-1]+(acc[i-1]+acc[i])*dt/2.0
        dis[i]=dis[i-1]+vel[i-1]*dt+(acc[i-1]/3.0+acc[i]/6.0)*dt*dt
    return vel,dis


def crac(dt,acc,vel,dis):
    #Amendment of base line
    accmax=np.max(np.abs(acc))
    nn=len(acc)
    tt=float(nn-1)*dt
    t=0.0
    for i in range(0,nn):
        dis[i]=dis[i]*(3.0*tt-2.0*t)*t*t
        t=t+dt
    dsum=(dis[0]+dis[nn-1])/2.0
    for i in range(1,nn-1):
        dsum=dsum+dis[i]
    dsum=dsum*dt
    a1=28.0/13.0/tt/tt*(2.0*vel[nn-1]-15.0/tt**5.0*dsum)
    a0=vel[nn-1]/tt-a1/2.0*tt
    t=0.0
    for i in range(0,nn):
        acc[i]=acc[i]-a0-a1*t
        t=t+dt
    acmax=np.max(np.abs(acc))
    coef=accmax/acmax
    acc=acc*coef
    return acc


def runge(t,dt,ainp,tp):
    def func(am,ac,ak,ff,x,y):
        dxdt=y
        dydt=ff/am-ac/am*y-ak/am*x
        return dxdt,dydt

    #Parameter setting
    hh=0.05
    #tp=10
    am=1.0
    ak=4*np.pi**2*am/tp**2
    ac=2*hh*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
    arac=np.max(np.abs(rac))
    avel=np.max(np.abs(vel))
    adis=np.max(np.abs(dis))
    return arac,avel,adis


def wrapper(args):
    return runge(*args)


def calc():
    for nnn in [1,2,3]:
        match nnn:
            case 1: fnameR='dat_acc_TCGH16_EW2.txt'
            case 2: fnameR='dat_acc_TCGH16_NS2.txt'
            case 3: fnameR='dat_acc_TCGH16_UD2.txt'
        dt,ndata,acc=inpdat(fnameR)
        vel,dis=iacc(dt,acc) # Integral of acceleration
        ainp=crac(dt,acc,vel,dis) #Modification of base line
        t=np.arange(0,dt*ndata,dt) # time step
        tw=np.linspace(np.log10(0.02),np.log10(10),200)
        ttp=10**(tw) # period

        val=[]
        for tp in ttp:
            val.append([t,dt,ainp,tp])
        result=Pool(8).map(wrapper,val)

        result=np.array(result)
        df=pd.DataFrame({
            't': ttp,
            'acc': result[:,0],
            'vel': result[:,1],
            'dis': result[:,2]
        })
        print(df)
        match nnn:
            case 1: fnameW='df_r_res1.csv'
            case 2: fnameW='df_r_res2.csv'
            case 3: fnameW='df_r_res3.csv'
        df.to_csv(fnameW)


def main():
    start=time.time()
    calc()
    end=time.time()-start
    print('{0:10.3f}(sec)'.format(end))


#---------------
# Execute
#---------------
if __name__ == '__main__': main()

FFT法による応答スペクトル計算

import numpy as np
from matplotlib import pyplot as plt
import time
from multiprocessing import Pool
import pandas as pd


def inpdat(fnameR):
    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=np.zeros(ndata,dtype=np.float64)
    for i in range(0,ndata):
        text=fin.readline();text=text.strip();wave[i]=float(text)
    fin.close()
    return dt,ndata,wave


def iacc(dt,acc):
    #Integral of acceleration
    nn=len(acc)
    vel=np.zeros(nn,dtype=np.float64)
    dis=np.zeros(nn,dtype=np.float64)
    for i in range(1,nn):
        vel[i]=vel[i-1]+(acc[i-1]+acc[i])*dt/2.0
        dis[i]=dis[i-1]+vel[i-1]*dt+(acc[i-1]/3.0+acc[i]/6.0)*dt*dt
    return vel,dis


def crac(dt,acc,vel,dis):
    #Amendment of base line
    accmax=np.max(np.abs(acc))
    nn=len(acc)
    tt=float(nn-1)*dt
    t=0.0
    for i in range(0,nn):
        dis[i]=dis[i]*(3.0*tt-2.0*t)*t*t
        t=t+dt
    dsum=(dis[0]+dis[nn-1])/2.0
    for i in range(1,nn-1):
        dsum=dsum+dis[i]
    dsum=dsum*dt
    a1=28.0/13.0/tt/tt*(2.0*vel[nn-1]-15.0/tt**5.0*dsum)
    a0=vel[nn-1]/tt-a1/2.0*tt
    t=0.0
    for i in range(0,nn):
        acc[i]=acc[i]-a0-a1*t
        t=t+dt
    acmax=np.max(np.abs(acc))
    coef=accmax/acmax
    acc=acc*coef
    return acc


def calfft(nn,xx,sp,wk1,tp):
    h=0.05
    #tp=10
    w0=2*np.pi/tp
    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
    arac=np.max(np.abs(rac))
    avel=np.max(np.abs(vel))
    adis=np.max(np.abs(dis))
    return arac,avel,adis


def wrapper(args):
    return calfft(*args)


def calc():
    for nnn in [1,2,3]:
        match nnn:
            case 1: fnameR='dat_acc_TCGH16_EW2.txt'
            case 2: fnameR='dat_acc_TCGH16_NS2.txt'
            case 3: fnameR='dat_acc_TCGH16_UD2.txt'
        dt,ndata,acc=inpdat(fnameR)
        vel,dis=iacc(dt,acc) # Integral of acceleration
        ainp=crac(dt,acc,vel,dis) #Modification of base line
        tw=np.linspace(np.log10(0.02),np.log10(10),200)
        ttp=10**(tw) # period

        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.float64)
        fk=k/nn/dt
        wk1=2*np.pi*fk

        val=[]
        for tp in ttp:
            val.append([nn,xx,sp,wk1,tp])
        result=Pool(8).map(wrapper,val)

        result=np.array(result)
        df=pd.DataFrame({
            't': ttp,
            'acc': result[:,0],
            'vel': result[:,1],
            'dis': result[:,2]
        })
        print(df)
        match nnn:
            case 1: fnameW='df_f_res1.csv'
            case 2: fnameW='df_f_res2.csv'
            case 3: fnameW='df_f_res3.csv'
        df.to_csv(fnameW)


def main():
    start=time.time()
    calc()
    end=time.time()-start
    print('{0:10.3f}(sec)'.format(end))


#---------------
# Execute
#---------------
if __name__ == '__main__': main()

作図

import numpy as np
from matplotlib import pyplot as plt
import time
import pandas as pd


def add_ax_rt(xmin,xmax,ymin,ymax,fsz):
    # Right y-axis
    plt.twinx()
    plt.yscale('log')
    plt.ylabel('Response displacement (cm)')
    plt.ylim([ymin,ymax])
    sd2=np.array([2,5,10,20,50,100,200,500])
    yy=2*np.pi/xmax*sd2
    sy=[str(n) for n in sd2]    
    plt.yticks(yy,sy)
    labels = plt.gca().get_yticklabels()
    plt.setp(labels, rotation=-45, fontsize=fsz-1)
    plt.gca().tick_params(axis='y',which='both',length=0,pad=1)
    sd1=np.array([0.01,0.02,0.05,0.1,0.2,0.5,1])
    for ssd in sd1:
        ys=ymin
        xs=2*np.pi/ys*ssd
        ss=str(ssd)
        plt.text(xs,ys,ss,ha='right',va='bottom',fontsize=fsz-1,rotation=-45)
    sd=np.concatenate([sd1,sd2])
    for ssd in sd:
        x1,x2=xmin,xmax
        y1,y2=2*np.pi/x1*ssd,2*np.pi/x2*ssd
        plt.plot([x1,x2],[y1,y2],'-',color='#999999',lw=0.5)


def add_ax_tp(xmin,xmax,ymin,ymax,fsz):
    # Top x-axis
    plt.twiny()
    plt.xscale('log')
    plt.xlabel('Response acceleration (gal)')
    plt.xlim([xmin,xmax])
    sa2=np.array([500,1000,2000,4000,10000,20000,50000,100000])
    xx=ymax*2*np.pi/sa2
    sx=[str(n) for n in sa2]
    plt.xticks(xx,sx)
    labels = plt.gca().get_xticklabels()
    plt.setp(labels, rotation=45, fontsize=fsz-1)
    plt.gca().tick_params(axis='x',which='both',length=0,pad=1,zorder=10)
    sa1=np.array([1,2,5,10,20,50,100,200])
    for ssa in sa1:
        xs=xmax
        ys=xs/2/np.pi*ssa
        ss=str(ssa)
        plt.text(xs,ys,ss,ha='right',va='top',fontsize=fsz-1,rotation=45)
    sa=np.concatenate([sa1,sa2])
    for ssa in sa:
        x1,x2=xmin,xmax
        y1,y2=x1/2/np.pi*ssa,x2/2/np.pi*ssa
        plt.plot([x1,x2],[y1,y2],'-',color='#999999',lw=0.5)


def fig_v(ttp1,vel1,ttp2,vel2,ttp3,vel3,fnameF):
    fsz=8
    xmin,xmax=0.02,10
    ymin,ymax=1,500
    iw,ih=4,4
    plt.figure(figsize=(iw,ih),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.gca().tick_params(axis='both',which='both',length=0)
    plt.xscale('log')
    plt.yscale('log')
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    ax_x=[0.02,0.05,0.1,0.2,0.5,1,2,5,10]
    ax_y=[1,2,5,10,20,50,100,200,500]
    ax_sx=[str(n) for n in ax_x]
    ax_sy=[str(n) for n in ax_y]
    plt.xticks(ax_x,ax_sx)
    plt.yticks(ax_y,ax_sy)
    plt.gca().hlines(y=ax_y,xmin=xmin,xmax=xmax,color='#000000',lw=0.5)
    plt.gca().vlines(x=ax_x,ymin=ymin,ymax=ymax,color='#000000',lw=0.5)
    plt.xlabel('Period (sec)')
    plt.ylabel('Response velocity (kine)')
    plt.title('TCGH16 (h=0.05)',loc='left',fontsize=fsz)
    plt.plot(ttp1,vel1, '-', color='#ff0000',lw=1,label='NW2')
    plt.plot(ttp1,vel2, '--',color='#ff0000',lw=1,label='NS2')
    plt.plot(ttp3,vel3, '--',color='#0000ff',lw=1,label='UD2')
    plt.legend(loc='upper left')

    add_ax_rt(xmin,xmax,ymin,ymax,fsz) # add right axis
    add_ax_tp(xmin,xmax,ymin,ymax,fsz) # add top axis

    plt.tight_layout()
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    #plt.show()


def fig_a(ttp1,acc1,ttp2,acc2,ttp3,acc3,fnameF):
    fsz=8
    xmin,xmax=0.01,10
    ymin,ymax=10,10000
    iw,ih=4,4
    plt.figure(figsize=(iw,ih),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.gca().tick_params(axis='both',which='both',length=0)
    plt.xscale('log')
    plt.yscale('log')
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    ax_x=[0.01,0.1,1,10]
    ax_y=[10,100,1000,10000]
    ax_sx=[str(n) for n in ax_x]
    ax_sy=[str(n) for n in ax_y]
    plt.xticks(ax_x,ax_sx)
    plt.yticks(ax_y,ax_sy)
    plt.grid(which="both")
    plt.xlabel('Period (sec)')
    plt.ylabel('Response acceleration (gal)')
    plt.title('TCGH16 (h=0.05)',loc='left',fontsize=fsz)
    plt.plot(ttp1,acc1, '-', color='#ff0000',lw=1,label='NW2')
    plt.plot(ttp2,acc2, '--',color='#ff0000',lw=1,label='NS2')
    plt.plot(ttp3,acc3, '--',color='#0000ff',lw=1,label='UD2')
    plt.legend()
    plt.tight_layout()
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    #plt.show()


def main():
    for nnn in [1,2,3]:
        match nnn:
            case 1: ss='n'
            case 2: ss='r'
            case 3: ss='f'
        df1 = pd.read_csv('df_'+ss+'_res1.csv')
        df2 = pd.read_csv('df_'+ss+'_res2.csv')
        df3 = pd.read_csv('df_'+ss+'_res3.csv')
        ttp1=df1['t'].to_numpy()
        acc1=df1['acc'].to_numpy()
        vel1=df1['vel'].to_numpy()

        ttp2=df2['t'].to_numpy()
        acc2=df2['acc'].to_numpy()
        vel2=df2['vel'].to_numpy()

        ttp3=df3['t'].to_numpy()
        acc3=df3['acc'].to_numpy()
        vel3=df3['vel'].to_numpy()

        fnameFa='fig_'+ss+'_ars_a.jpg'
        fnameFv='fig_'+ss+'_ars_v.jpg'
        fig_a(ttp1,acc1,ttp2,acc2,ttp3,acc3,fnameFa) # Acceleration response spectrum
        fig_v(ttp1,vel1,ttp2,vel2,ttp3,vel3,fnameFv) # Velocity response spectrum


#---------------
# Execute
#---------------
if __name__ == '__main__': main()

最後にFFT法による速度応答スペクトル（３重応答スペクトル図）を示す。

以　上