テスト

1自由度ダフィング運動方程式

Python

import numpy as np import matplotlib.pyplot as plt def model(): ''' Define model ''' # 質量・減衰・剛性の集中定数を設定する M = 0.5 Ka = 1000 Kb = 100000000 C = 5 return M, Ka, Kb, C def f(x, v, M, Ka, Kb, C, force): ''' Equation ''' # 減衰強制振動モデル y = - (Ka / M) * x - (Kb / M) * x ** 3 - (C / M) * v + force return y def sim(M, Ka, Kb, C): ''' Solve using 4th-order Runge-Kutta method ''' # 解析条件 t0 = 0.0 t1 = 5.0 dt = 1e-5 # 初期条件 x0 = 0 v0 = 0 x, v = x0, v0 t_axis = np.arange(t0, t1, dt) x_sol = [] v_sol = [] # 外力を定義 A = 1000 freq = 10 force = A * np.sin(2 * np.pi * freq * t_axis) # 4次のRunge-Kutta法による数値解析 iteration = 0 for t in t_axis: # 1自由度目の変位と速度をappend x_sol.append(x) v_sol.append(v) # イタレーションをモニタ if iteration % 100 == 0: print('iteration=', iteration, 'time=', t) # Runge-Kuttaのメイン k11 = f(x, v, M, Ka, Kb, C, force[iteration]) * dt k12 = v * dt k21 = f(x + k11 / 2, v + k12 / 2, M, Ka, Kb, C, force[iteration]) * dt k22 = (v + k12 / 2) * dt k31 = f(x + k21 / 2, v + k22 / 2, M, Ka, Kb, C, force[iteration]) * dt k32 = (v + k22 / 2) * dt k41 = f(x + k31, v + k32, M, Ka, Kb, C, force[iteration]) * dt k42 = (v + k32) * dt v += (k11 + 2 * k21 + 2 * k31 + k41) / 6 x += (k12 + 2 * k22 + 2 * k32 + k42) / 6 iteration += 1 return t_axis, x_sol, v_sol def plot(t, x, v): ''' plot ''' # フォントの種類とサイズを設定する。 plt.rcParams['font.size'] = 14 plt.rcParams['font.family'] = 'Times New Roman' # 目盛を内側にする。 plt.rcParams['xtick.direction'] = 'in' plt.rcParams['ytick.direction'] = 'in' # グラフの上下左右に目盛線を付ける。 fig = plt.figure(figsize=(10, 6)) ax1 = fig.add_subplot(221) ax1.yaxis.set_ticks_position('both') ax1.xaxis.set_ticks_position('both') ax2 = fig.add_subplot(223) ax2.yaxis.set_ticks_position('both') ax2.xaxis.set_ticks_position('both') ax3 = fig.add_subplot(122) ax3.yaxis.set_ticks_position('both') ax3.xaxis.set_ticks_position('both') # 軸のラベルを設定する。 ax1.set_xlabel('Time[s]') ax1.set_ylabel('Displacement[m]') ax2.set_xlabel('Time[s]') ax2.set_ylabel('Velocity[m/s]') ax3.set_xlabel('Displacement[m]') ax3.set_ylabel('Velocity[m/s]') # スケールの設定をする。 #ax3.set_xlim(-0.05, 0.05) #ax3.set_ylim(-10, 10) # データプロットの準備とともに、ラベルと線の太さ、凡例の設置を行う。 ax1.plot(t, x, label='Displacement', lw=1, color='red') ax2.plot(t, v, label='Velocity', lw=1, color='red') ax3.plot(x, v, label='Phase plane', lw=1, color='red') # レイアウト設定 fig.tight_layout() # グラフを表示する。 plt.show() plt.close() if __name__ == '__main__': ''' メイン処理 ''' M, Ka, Kb, C = model() t, x, v = sim(M, Ka, Kb, C) plot(t, x, v)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127

import numpy as np import matplotlib.pyplot as plt   def model():     ''' Define model '''       # 質量・減衰・剛性の集中定数を設定する     M = 0.5     Ka = 1000     Kb = 100000000     C = 5       return M, Ka, Kb, C   def f(x, v, M, Ka, Kb, C, force):     ''' Equation '''       # 減衰強制振動モデル     y = - (Ka / M) * x - (Kb / M) * x ** 3 - (C / M) * v + force       return y   def sim(M, Ka, Kb, C):     ''' Solve using 4th-order Runge-Kutta method '''       # 解析条件     t0 = 0.0     t1 = 5.0     dt = 1e-5       # 初期条件     x0 = 0     v0 = 0     x, v = x0, v0       t_axis = np.arange(t0, t1, dt)     x_sol = []     v_sol = []       # 外力を定義     A = 1000     freq = 10     force = A * np.sin(2 * np.pi * freq * t_axis)       # 4次のRunge-Kutta法による数値解析     iteration = 0     for t in t_axis:           # 1自由度目の変位と速度をappend         x_sol.append(x)         v_sol.append(v)           # イタレーションをモニタ         if  iteration % 100 == 0:             print('iteration=', iteration, 'time=', t)           # Runge-Kuttaのメイン         k11 = f(x, v, M, Ka, Kb, C, force[iteration]) * dt         k12 = v * dt         k21 = f(x + k11 / 2, v + k12 / 2, M, Ka, Kb, C, force[iteration]) * dt         k22 = (v + k12 / 2) * dt         k31 = f(x + k21 / 2, v + k22 / 2, M, Ka, Kb, C, force[iteration]) * dt         k32 = (v + k22 / 2) * dt         k41 = f(x + k31, v + k32, M, Ka, Kb, C, force[iteration]) * dt         k42 = (v + k32) * dt         v += (k11 + 2 * k21 + 2 * k31 + k41) / 6         x += (k12 + 2 * k22 + 2 * k32 + k42) / 6         iteration += 1         return t_axis, x_sol, v_sol     def plot(t, x, v):     ''' plot '''       # フォントの種類とサイズを設定する。     plt.rcParams['font.size'] = 14     plt.rcParams['font.family'] = 'Times New Roman'       # 目盛を内側にする。     plt.rcParams['xtick.direction'] = 'in'     plt.rcParams['ytick.direction'] = 'in'       # グラフの上下左右に目盛線を付ける。     fig = plt.figure(figsize=(10, 6))     ax1 = fig.add_subplot(221)     ax1.yaxis.set_ticks_position('both')     ax1.xaxis.set_ticks_position('both')     ax2 = fig.add_subplot(223)     ax2.yaxis.set_ticks_position('both')     ax2.xaxis.set_ticks_position('both')     ax3 = fig.add_subplot(122)     ax3.yaxis.set_ticks_position('both')     ax3.xaxis.set_ticks_position('both')       # 軸のラベルを設定する。     ax1.set_xlabel('Time[s]')     ax1.set_ylabel('Displacement[m]')     ax2.set_xlabel('Time[s]')     ax2.set_ylabel('Velocity[m/s]')     ax3.set_xlabel('Displacement[m]')     ax3.set_ylabel('Velocity[m/s]')       # スケールの設定をする。     #ax3.set_xlim(-0.05, 0.05)     #ax3.set_ylim(-10, 10)       # データプロットの準備とともに、ラベルと線の太さ、凡例の設置を行う。     ax1.plot(t, x, label='Displacement', lw=1, color='red')     ax2.plot(t, v, label='Velocity', lw=1, color='red')     ax3.plot(x, v, label='Phase plane', lw=1, color='red')       # レイアウト設定     fig.tight_layout()       # グラフを表示する。     plt.show()     plt.close()   if __name__ == '__main__':     ''' メイン処理 '''       M, Ka, Kb, C = model()     t, x, v = sim(M, Ka, Kb, C)       plot(t, x, v)