🎉

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

2023/01/14に公開

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

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

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

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

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 #x軸の終点
ymin=0.0 #y軸の始点
ymax=3 #y軸の終点

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

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

X=np.linspace(xmin,xmax*3,nx) #x軸の定義
Y=np.linspace(ymin,ymax*3,ny) #y軸の定義
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 # 壁がある部分では波の高さは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 #x軸の終点
ymin=0.0 #y軸の始点
ymax=3 #y軸の終点

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

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

X=np.linspace(xmin,xmax*3,nx) #x軸の定義
Y=np.linspace(ymin,ymax*3,ny) #y軸の定義
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 # 壁がある部分では波の高さは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

Discussion