[Pythonによる科学・技術計算] FTCS法(陽解法)による1次元・2次元波動方程式の数値解法，双曲型偏微分方程式

#はじめに

空間に広がった物質・空間そのもの(媒質)中を伝わる振動を一般に波または波動という[1]。力学的な波動では弾性波や流体運動が，電磁気学的な波動は電磁波がそれに対応する。媒質中で振動する物理量$u(x,y,z,t)$が満たす方程式が波動方程式である。

本稿ではPythonを利用し2次元波動方程式(双曲型偏微分方程式)
$\frac{1}{c^2}\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} $

を陽解法の差分スキームであるFTCS法(Forward in Time and Centered Difference method)[補遺]により解く(追記:1次元の場合も入れた。内容(2)を参照)

なお， この方程式の力学における導出にあたって微少振動の仮定と波を伝える媒質の一様・等方性が仮定されている(位相速度cが定数)点に注意すること。

#内容

##(1) 2次元波動方程式
陽的差分スキームであるFTCS法により2次元波動方程式を解く。
計算条件は以下の通りである。

• 計算領域: Lx = 1, Ly=1の正方形
  x方向 (Nx)，y方向(Ny)共に41点刻み (Nx=Ny=41)

• 計算時間の最大値(t_max)を0.5としている。

• 時間刻み(delta_t):数値安定性を考慮した値をさらに半分にしたものを使っている。

• 波の位相速度: c=1

• 境界条件: 4辺とも固定端とする (振動せず動かない)

• 初期条件1: 時刻0で正方形の中心位置にガウス分布関数を設置。これで波を"引っ張る" (図1参照)。

• 初期条件2:$\frac{\partial u}{\partial t} = 0 \ (t=0)$ (初速は0 )

また，計算ではt=1ステップ目の波動$u$の情報を必要とする。本シミュレーションでは初期条件2を中心差分して得られる結果を用いて評価した[4].

図1. 初期条件として正方形中心にガウシアンを設置。

---

##(2) 1次元波動方程式
$\frac{1}{c^2}\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} $
FTCS法で解く。

境界条件は，端点で$u=0$ (固定端)とする。
初期条件は，$\frac{\partial u}{\partial t} = 0 \ (t=0)$ (初速は0 )である。

波形に対する初期条件は以下の図のとおり。

図2. 初期波形

---

#計算コード (1)2次元波動方程式

#%matplotlib nbagg
"""
2D Wave Equation
一様媒質: 位相速度が時間・空間的に一定
FTCS法
"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import ArtistAnimation # アニメーション作成のためのメソッドをインポート


c=1
print ("c=",c)
Lx = 1 # 長さ[m]
Ly = 1
t_max = 0.5

Nx = 41
Ny =41

delta_x = Lx/Nx
delta_y = Ly/Ny
delta_t=(1/(1/delta_x+1/delta_y))/c*0.5 #時間刻み幅。あるていど小さめにとる。

print("delta_t=",delta_t)

Nt = int(t_max/delta_t)
 

print("Nt=",Nt)

stx=delta_x/delta_t
sty=delta_y/delta_t



u =  np.empty((Nx,Ny,Nt),dtype=float)

#境界条件
u[0,:,:], u[-1,:,:] , u[:,0,:] , u[:,-1,:] = 0, 0, 0, 0

#初期条件
for i in range (Nx) :
    for ii in range(Ny):
        u[i,ii,0]  = (np.exp(-((i*delta_x-Lx/2)**2+(ii*delta_y-Ly/2)**2)/0.01))  #初期条件。ガウス分布関数の設置。

#波動 u[i,ii,1]を作るための計算。初期条件2を中心差分して得られる。
for i in range(1,Nx-1): 
    for ii in range(1,Ny-1):
        u[i,ii,1]= u[i,ii,0]+0.5*((delta_t/delta_x)**2)*(c**2)*(u[i+1,ii, 0]+u[i-1, ii,0]-2*u[i,ii,0])\
+0.5*((delta_t/delta_y)**2)*(c**2)*(u[i,ii+1, 0]+u[i,ii-1, 0]-2*u[i,ii,0])


# #波動の時間発展のためのループ
for j in range(1,Nt-1):
    for i in range(1,Nx-1):
        for ii in range(1,Ny-1):
            u[i,ii, j+1] = 2*u[i,ii,j]-u[i,ii,j-1]+((c/stx)**2)* (u[i+1,ii,j]+u[i-1,ii,j]\
-2*u[i,ii,j])+((c/sty)**2)* (u[i,ii+1,j]+u[i,ii-1,j]-2*u[i,ii,j]) 
       

得られた波動関数 u[x,y,t]を可視化するため，以下のコードを利用する。
ワイヤーフレーム表示[3]でプロットし，結果をアニメーションにしている。

#可視化のためのプログラム

%matplotlib nbagg 
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation # アニメーション作成のためのメソッドをインポート

fig = plt.figure(figsize = (8, 6))
ax = fig.gca(projection='3d')

x=list(range(Nx))
y=list(range(Ny))
X, Y = np.meshgrid(x,y)

def update(i,u,fig_title):
    if i != 0:
        ax.clear()
    ax.view_init(elev=60., azim=60.) # アングル設定
    #ax.plot_surface(X,Y,u[X,Y,i],rstride=1,cstride=1,cmap='Greys',linewidth=0.3) # サーフェスプロット
    ax.plot_wireframe(X,Y,u[X,Y,i],rstride=1,cstride=1,color='blue',linewidth=0.3) # ワイヤーフレーム表示
    
    ax.set_title(fig_title + 'i=' + str(i))
    ax.set_zlim(-0.4,0.4) # zの描画範囲の指定
    ax.set_xlabel('X') # ラベル
    ax.set_ylabel('Y')
    ax.set_zlabel('U')


    
    return ax


ani = animation.FuncAnimation(fig,update,fargs = (u,'Wave motion: time step='), interval =1, frames = Nt,blit=False)
fig.show()

ani.save("output.gif", writer="imagemagick")

#結果(1)
アニメーションにした結果を示す。中心位置に生じた波が周りへ伝播していく様子が見られる。

---

#計算コード(2): 1次元波動方程式

%matplotlib nbagg
"""
1D　Wave Equation
FTCS法

"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import ArtistAnimation # アニメーション作成のためのメソッドをインポート


T=40 # 張力 [N]
rho=0.01 # 密度 [Kg/m]
c = np.sqrt(T/rho) # 位相速度

print ("c=",c)
Lx = 1 # 長さ[m]
t_max = 0.1

Nx = 100

delta_x = Lx/Nx

delta_t=delta_x/c   # 安定条件を満たす時間刻みの最大値
#print("delta_t=",delta_t)

Nt = int(t_max/delta_t)
st=delta_x/delta_t 
print("c-st=",c-st) # 安定性条件のチェック (負数なら不安定)

u =  np.zeros([Nx,Nt])

#境界条件
u[0,:] = 0
u[-1,:] = 0

#初期条件
for i in range (Nx) :
    if i <= 0.8*Nx :
        u[i,0]=1.25*i/Nx
    
    else :  
        u[i,0] = 5.0*(1-i/Nx)

for i in range(1,Nx-1): # u[:,1]の設定
    u[i,1]  = u[i,0]+0.5*((delta_t/delta_x)**2)*(c**2)*(u[i+1,0]+u[i-1,0]-2*u[i,0])

for j in range(Nt-1):
    for i in range(1,Nx-1):
        u[i,j+1] = 2*u[i,j]-u[i,j-1]+(c**2/st**2)* (u[i+1,j]+u[i-1,j]-2*u[i,j])
        

x=list(range(Nx))
y=list(range(Nt))

X, Y = np.meshgrid(x,y)

def functz(u):
    z=u[X,Y]
    return z

Z = functz(u)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z, color='r')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('U')

plt.show()

#結果(2) 1次元波動方程式

図3. 一次元波動方程式の解

次に，アニメーションで結果を図示する。

%matplotlib nbagg 
from matplotlib.animation import ArtistAnimation # アニメーション作成のためのメソッドをインポート

fig = plt.figure()

anim = [] #アニメーション用に描くパラパラ図のデータを格納するためのリスト

for i in range(Nt):
    U=list(u[:,i])
    x=list(range(Nx))
    if i % int(Nt*0.01) ==0: 
        im=plt.plot(x,U, '-', color='red',markersize=10, linewidth = 2, aa=True)
        anim.append(im)

anim = ArtistAnimation(fig, anim) # アニメーション作成    
plt.xlabel('x')
plt.ylabel('U')

fig.show() 
anim.save("t.gif", writer='imagemagick')   #アニメーションをt.gifという名前で保存し，gifアニメーションファイルを作成する。

1次元波動方程式の解のアニメーション

---

#補遺

##(1) FTCS法
偏微分方程式における**FTCS法(Forward in Time and Centered Difference method)[2]**は，時間微分については前進差分，空間微分については中心差分をとる単純なスキームである。実装が容易であるが数値不安定性がでないように時間刻みを細かくとる必要がある

(2) ここで扱ったシミュレーションでは減衰効果を考慮していない。また，位相速度の時間・空間依存性を考慮していない。これらについて，参考文献[4]に初歩的な説明がある。

---

#参考文献

[1] 長岡 洋介,「振動と波」
[2] William H. Press,「ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ」
[3] Qiita記事: [Pythonによる科学・技術計算] 2次元(カラー)等高線等の描画，可視化，matplotlib
[4] Rubin H. Landau、「計算物理学 応用編」 (訳書)