波の伝搬と干渉の３次元シミュレーションをしてみた（波動方程式、ダブルスリット、python）

こんにちにゃんです。
水色桜（みずいろさくら）です。
今回は波の伝搬と干渉の３次元シミュレーションについて書いていこうと思います。
実行結果は下のgifみたいな感じになります。

まず背景について説明いたします。
波の挙動は波動方程式によって記述されます。

波動方程式はスライドに示すような定数係数二階線形偏微分方程式です。そのままでは機械で分析ができないため、スライドに示すように数値的に微分を行います。
結果として一番下のような形になります。
この式を元にしてプログラムを書いていきます。
境界条件はディレクレ境界条件（端の値を常に０とする）とします。

プログラムは次の通りです。

wave.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as anim
from mpl_toolkits.mplot3d import Axes3D

dx=0.1 #空間、x軸の分解能
dy=dx #空間、y軸の分解能
dt=0.05 #時間の分解能
tmin=0.0 #最初の時刻
tmax=6.0 #終了の時刻

xmin=0.0 #x軸の始点
xmax=3 #ｘ軸の終点
ymin=0.0 #y軸の始点
ymax=3 #y軸の終点

c=1.0 #波の伝搬速度
k=(c*dt/dx)**2 #計算で使う定数

nx=int((xmax-xmin)/dx)+1 #ｘ軸の分割数
ny=int((ymax-ymin)/dy)+1 #ｙ軸の分割数
nt=int((tmax-tmin)/dt)+2 #時間の分割数

X=np.linspace(xmin,xmax*3,nx) #ｘ軸の定義
Y=np.linspace(ymin,ymax*3,ny) #ｙ軸の定義
X,Y=np.meshgrid(Y,X) #メッシュの生成

u=np.zeros((nt,nx,ny)) #波を記録する配列

u_0=np.exp(-((X-0.5)**2)*10)*np.exp(-((Y-4.5)**2)*10)*20#最初に波のもとを与える
u_1=np.zeros((nx,ny)) #初期条件（波のもととなる地点の隣）

u[0]=u_0
u[1]=u[0]+dt*u_1

for t in range(1,nt-1):
     for y in range(1,ny-1): 
        for x in range(1,nx-1):
            if (ny/2-1<y<ny/2+1) and (0<x<nx/4 or nx/3<x<2*nx/3 or 3*nx/4<x<nx):
                    u[t,x,y]=0 # 壁がある部分では波の高さは０
            else: 
                u[t+1,x,y]=(2*(1-2*k)*u[t,x,y]-u[t-1,x,y]+k*(u[t,x-1,y]+u[t,x+1,y]+u[t,x,y-1]+u[t,x,y+1]))
    # 境界条件
     for x in range(1,nx-1):
        u[t+1,x,0]=u[t+1,x,1]
        u[t+1,x,ny-1]=u[t+1,x,ny-2]
     for y in range(1,ny-1):
        u[t+1,0,y]=u[t+1,1,y]
        u[t+1,nx-1,y]=u[t+1,nx-2,y]
     u[t,0,0]=(u[t,1,0]+u[t,0,1])/2
     u[t,nx-1,0]=(u[t,nx-2,0]+u[t,nx-1,1])/2
     u[t,0,ny-1]=(u[t,1,ny-1]+u[t,0,ny-2])/2
     u[t,nx-1,ny-1]=(u[t,nx-2,ny-1]+u[t,nx-1,ny-2])/2
     
v=np.zeros((nt,nx,ny))
for t in range(1,nt-1):
     for x in range(nx): 
        for y in range(ny):
            if (ny/2-1<y<ny/2+1) and (0<x<nx/4 or nx/3<x<2*nx/3 or 3*nx/4<x<nx):
                v[t,x,y]=u[t,x,y]+5 # 波を遮る壁を設置
            else:
                v[t,x,y]=u[t,x,y]

# 三次元のグラフを定義
fig=plt.figure(figsize=(6,7.5))
fig.set_dpi(100)
ax=Axes3D(fig)

# アニメーションの定義
def animate(i):
    ax.clear()
    ax.plot_surface(X,Y,v[i],cmap='cool') # cmapは表面の色合いを定義
    ax.set_zlim(-3,20)
    ax.set_xlim(0,9)
    ax.set_ylim(0,9)
    k=int(i*10)
    ax.set_title('Time: '+str(k*dt/10)+' s',size=10) #タイトルの設定

ani=anim.FuncAnimation(fig,animate,frames=nt-1,interval=10,repeat=True)
plt.show()
ani.save('wave_rect.gif',writer='pillow')

まず、必要なライブラリのインポートと定数の定義を行います。

intro.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as anim
from mpl_toolkits.mplot3d import Axes3D

dx=0.1 #空間、x軸の分解能
dy=dx #空間、y軸の分解能
dt=0.05 #時間の分解能
tmin=0.0 #最初の時刻
tmax=6.0 #終了の時刻

xmin=0.0 #x軸の始点
xmax=3 #ｘ軸の終点
ymin=0.0 #y軸の始点
ymax=3 #y軸の終点

c=1.0 #波の伝搬速度
k=(c*dt/dx)**2 #計算で使う定数

nx=int((xmax-xmin)/dx)+1 #ｘ軸の分割数
ny=int((ymax-ymin)/dy)+1 #ｙ軸の分割数
nt=int((tmax-tmin)/dt)+2 #時間の分割数

X=np.linspace(xmin,xmax*3,nx) #ｘ軸の定義
Y=np.linspace(ymin,ymax*3,ny) #ｙ軸の定義
X,Y=np.meshgrid(Y,X) #メッシュの生成

u=np.zeros((nt,nx,ny)) #波を記録する配列

u_0=np.exp(-((X-0.5)**2)*10)*np.exp(-((Y-4.5)**2)*10)*20#最初に波のもとを与える
u_1=np.zeros((nx,ny)) #初期条件（波のもととなる地点の隣）

その後、波動方程式を数値的に解いていきます。

calc.py

for t in range(1,nt-1):
     for y in range(1,ny-1): 
        for x in range(1,nx-1):
            if (ny/2-1<y<ny/2+1) and (0<x<nx/4 or nx/3<x<2*nx/3 or 3*nx/4<x<nx):
                    u[t,x,y]=0 # 壁がある部分では波の高さは０
            else: 
                u[t+1,x,y]=(2*(1-2*k)*u[t,x,y]-u[t-1,x,y]+k*(u[t,x-1,y]+u[t,x+1,y]+u[t,x,y-1]+u[t,x,y+1]))
    # 境界条件
     for x in range(1,nx-1):
        u[t+1,x,0]=u[t+1,x,1]
        u[t+1,x,ny-1]=u[t+1,x,ny-2]
     for y in range(1,ny-1):
        u[t+1,0,y]=u[t+1,1,y]
        u[t+1,nx-1,y]=u[t+1,nx-2,y]
     u[t,0,0]=(u[t,1,0]+u[t,0,1])/2
     u[t,nx-1,0]=(u[t,nx-2,0]+u[t,nx-1,1])/2
     u[t,0,ny-1]=(u[t,1,ny-1]+u[t,0,ny-2])/2
     u[t,nx-1,ny-1]=(u[t,nx-2,ny-1]+u[t,nx-1,ny-2])/2
     
v=np.zeros((nt,nx,ny))
for t in range(1,nt-1):
     for x in range(nx): 
        for y in range(ny):
            if (ny/2-1<y<ny/2+1) and (0<x<nx/4 or nx/3<x<2*nx/3 or 3*nx/4<x<nx):
                v[t,x,y]=u[t,x,y]+5 # 波を遮る壁を設置
            else:
                v[t,x,y]=u[t,x,y]

最後にグラフを描画し、アニメーションを作製します。

gra.py

# 三次元のグラフを定義
fig=plt.figure(figsize=(6,7.5))
fig.set_dpi(100)
ax=Axes3D(fig)

# アニメーションの定義
def animate(i):
    ax.clear()
    ax.plot_surface(X,Y,v[i],cmap='cool') # cmapは表面の色合いを定義
    ax.set_zlim(-3,20)
    ax.set_xlim(0,9)
    ax.set_ylim(0,9)
    k=int(i*10)
    ax.set_title('Time: '+str(k*dt/10)+' s',size=10) #タイトルの設定

ani=anim.FuncAnimation(fig,animate,frames=nt-1,interval=10,repeat=True)
plt.show()
ani.save('wave_rect.gif',writer='pillow')

以上がpythonで波の伝搬と干渉の三次元シミュレーションを行う方法です。
是非皆さんもシミュレーションを行ってみてください。
では、ばいにゃん～

参考にさせていただいたもの
1次元波動方程式を差分法で解く（python)
https://qiita.com/t--shin/items/b109f960b0a109dc358e