二次元波動方程式のグラフ

Question

二次元の波動方程式を数値的に時に、時間による波の形の推移をアニメーションに表示するプログラムを作成しました。このプログラムにおいて、t=0.5,1など、指定した時間における波の状態のみを画像で出力して保存したいのですが、うまくいきません。現状アニメーションまではちゃんと動作しています。下のプログラムについて、三変数関数であるuに対して、変数であるtの部分に時刻を指定してやればいいと思ったのですが、tに具体的な数字を入れるだけでは作動してくれません。どうすればいいでしょうか？

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

#パラメータ設定
dx = 0.02
dy = dx
dt = 0.01
tmin = 0.0
tmax = 3.0

#正方形エリアの設定
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 1.0

v = 1.0
rsq = (v*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,nx)
Y = np.linspace(ymin,ymax,ny)
X,Y = np.meshgrid(Y,X)

u = np.zeros((nt,nx,ny))

#初期条件
u_0 = np.exp(-(X-0.5)**2*3)*np.exp(-(Y-0.5)**2*3)*2
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 x in range(1,nx-1):
        for y in range(1,ny-1):
            u[t+1,x,y] = 2*(1-2*rsq)*u[t,x,y]-u[t-1,x,y]+rsq*(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,x,0] = 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
fig = plt.figure()
fig.set_dpi(100)
ax = Axes3D(fig)

def animate(i):
    ax.clear()
    ax.plot_surface(X,Y,u[i],rstride=1,cstride=1,cmap=plt.cm.coolwarm)
    ax.set_zlim(-3,3)

anim = animation.FuncAnimation(fig,animate,frames=nt-1,interval=10,repeat=False)
anim.save("animation1.gif",writer="pillow")


plt.show()

Comments

[y_waiwai]

このままではコードが見づらいので、質門を編集し、<code>ボタンで、出てくる’’’の枠の中にコードを貼り付けてください

[jbpb0]

pythonのコードの一番最初の行のすぐ上に
```python
だけの行を追加してください
また、pythonのコードの一番最後の行のすぐ下に
```
だけの行を追加してください

現状、コードがとても読み辛いです
質問にコードを載せる際に上記をやってくれたら、他人がコードを読みやすくなり、コードの実行による現象確認もやりやすくなるので、回答されやすくなります

[lucy-0333]

申し訳ありません。修正いたしました。

[jbpb0]

> uに対して、変数であるtの部分に時刻を指定

u[int(時刻/dt)]
で、どうでしょうか？
時刻はtmin〜tmaxのどれか

[lucy-0333]

plt.plot(X,Y,u[int(0.5/dt),X,Y])
このような感じでつけてみましたが、「arrays used as indices must be of integer (or boolean) type」というエラーが発生しました。

[jbpb0]

plot()じゃなくてplot_surface()です

また、u[int(0.5/dt),X,Y]の「,X,Y」は要らないはず

anim =...
anim.save(...
の行を無効にして、その代わりに
plt.show()
のすぐ上に
animate(int(時刻/dt))
を追加したら、どうでしょうか？

[lucy-0333]

解決しました！本当にありがとうございました！

Answer 1

uに対して、変数であるtの部分に時刻を指定

anim = animation.FuncAnimation(fig,animate,frames=nt-1,interval=10,repeat=False)
anim.save("animation1.gif",writer="pillow")

↓

animate(int(時刻/dt))

時刻はtmin〜tmaxのどれか

Comments

[lucy-0333]

ありがとうございました！