Pythonで波動方程式の数値計算②
以前作成した一次元波動方程式の数値計算の初期条件を変えてどのような挙動になるか見てみました。
例えば上のアニメーションのように弦の真ん中でインパルス上の波を起こすと、時間が経つにつれ波が左右に移動し（振幅は半分になる）、弦の両端で位相がひっくりかえり逆方向に進みます。二つの波がぶつかると合成し振幅が大きくなります。波の挙動がちゃんと表現できていました。面白いですね。
コード
# matplotlibのアニメーション表示用 %matplotlib nbagg import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation import math # 各種定数の設定 t0 = 0 tM = 2 M = 40 dt = (tM-t0)/M # 時間間隔 x0 = 0 xN = 1 N = 20 dx = (xN-x0)/N # 空間間隔 a = (dt/dx)**2 # 差分方程式に出てくる係数 ux0 = 0 # 境界条件 uxN = 0 # 境界条件 A = 2 # 初期条件の式の振幅 # 変数の初期化 u = np.zeros([M+1, N+1]) # 行：時間、列：ｘ # 初期条件1 (t=0でのuの分布) u(x)=A*exp(-5*(x-0.5)^2) for xi in range(N+1): u[0,xi] = A*math.exp(-((xi-N/2)**2)*5) # 初期条件2 (t=dtでのuの分布) u[1,0] = 0 u[1,N] = 0 for xi in range(1,N): u[1,xi] = u[0,xi] + dt*0 + a/2*(u[0,xi+1] -2*u[0,xi] + u[0,xi-1]) # 差分解 (t=2*dt~M*dt) for ti in range(1,M): # 両端の弦の変異は常に0 u[ti+1,0] = ux0 u[ti+1,N] = uxN # 両端以外を差分方程式により更新 for xj in range(1,N): u[ti+1, xj] = 2*u[ti,xj] -u[ti-1,xj] + a*(u[ti,xj+1] -2*u[ti,xj] + u[ti,xj-1]) # matplotlibで表示 x = np.linspace(x0, xN, N+1) # maplotlibで表示するためのx軸の設定 fig = plt.figure(figsize=(6,4)) ax = fig.add_subplot(1,1,1) # アニメ更新用の関数 def update_func(i): # 前のフレームで描画されたグラフを消去 ax.clear() ax.plot(x, u[i,:], color='blue') ax.scatter(x, u[i,:], color='blue') # 軸の設定 ax.set_ylim(-2.22, 2.22) # 軸ラベルの設定 ax.set_xlabel('x', fontsize=12) ax.set_ylabel('u', fontsize=12) # サブプロットタイトルの設定 ax.set_title('Time: ' + '{:.2f}'.format(tM*i/M)) ani = animation.FuncAnimation(fig, update_func, frames=M, interval=100,repeat=True) # アニメーションの保存 #ani.save('wave_eq3.gif', writer='imagemagick') # 表示 plt.show()
