Python
1import numpy as np 2import pandas as pd 3import mpl_toolkits.axes_grid1 4from math import * 5from mpl_toolkits.axes_grid1 import make_axes_locatable 6from matplotlib import pyplot as plt 7from matplotlib import animation as animation 8# model N x M:PBC 9N, M=30,60 10#parameters 11pi=4*atan(1.0) #定数の定義 12dt, h=0.01,0.1 13h2=h*h 14DX,DY=0.0001,0.0001 15A,B=1.0,3.0 16maxX,maxY=0.0,0.0 17minX,minY=6.0,6.0 18T=20 19diagX, offX=1-4*DX*dt/h2, DX*dt/h2#差分法の定数部分の定義 20diagY, offY=1-4*DY*dt/h2, DY*dt/h2 21#preparation 22m=np.zeros([N+2,M+2],float) 23X=np.zeros([N+2,M+2],float) 24Y=np.zeros([N+2,M+2],float) 25XX=np.zeros([N+2,M+2],float) 26YY=np.zeros([N+2,M+2],float) 27rhsX=np.zeros([N+2,M+2],float) 28rhsY=np.zeros([N+2,M+2],float) 29fig = plt.figure() 30ax1 = fig.add_subplot(1, 2, 1) 31ax2 = fig.add_subplot(1, 2, 2) 32ax1.set_xlabel("x") 33ax1.set_ylabel("y") 34ax1.invert_yaxis() 35ax2.set_xlabel("x") 36ax2.set_ylabel("y") 37ax2.invert_yaxis() 38for j in range(N+2): 39 for k in range(M+2): 40 s,t=pi*j/N,pi*k/M 41 X[j][k]=2.0+0.4*sin(s)*sin(3*t) 42 Y[j][k]=1.0+0.4*sin(s)*sin(t) 43def SX(u,v): 44 return (A-(B+1)*u+u*u*v) 45def SY(u,v): 46 return (B*u-u*u*v) 47def update(frames): 48 for j in range (N+1): 49 for k in range (M+1): 50 rhsX[j][k]=dt*SX(X[j][k],Y[j][k]) 51 rhsY[j][k]=dt*SY(X[j][k],Y[j][k]) 52 rhsX[j][k]+=offX*(X[j+1][k]+X[j-1][k]+X[j][k-1]+X[j][k+1])+diagX*X[j][k] 53 rhsY[j][k]+=offY*(Y[j+1][k]+Y[j-1][k]+Y[j][k-1]+Y[j][k+1])+diagY*Y[j][k] 54 XX[j][k],YY[j][k]=rhsX[j][k],rhsY[j][k] 55 X, Y=XX, YY 56 X_T=np.array(X).T 57 Y_T=np.array(Y).T 58 im1=ax1.imshow(X_T,interpolation="nearest",cmap="jet") 59 im2=ax2.imshow(Y_T,interpolation="nearest",cmap="jet") 60 fig.tight_layout() 61 divider1 = mpl_toolkits.axes_grid1.make_axes_locatable(ax1) 62 divider2 = mpl_toolkits.axes_grid1.make_axes_locatable(ax2) 63 cax1 = divider1.append_axes('right', '5%', pad='3%') 64 cax2 = divider2.append_axes('right', '5%', pad='3%') 65 fig.colorbar(im1,cax=cax1) 66 fig.colorbar(im2,cax=cax2) 67ani=animation.FuncAnimation(fig,update,frames=T+1,interval=10) 68ani.save("BZ_reaction2.gif",writer="ffmpeg")
後述のArtistAnimationを用いての出力はうまく行ったのでFuncAnimationで計算結果をアニメーションで出力したいのですがうまく行かず困っています。
発生している問題・エラーメッセージ
line 50, in update rhsX[j][k]=dt*SX(X[j][k],Y[j][k]) UnboundLocalError: local variable 'X' referenced before assignment
試したこと
def update(frames)内におそらくXの定義がないといけないと思われるので定数の定義からすべてをdef update(frames)の中に入れてみてもうまくグラフは出力されませんでした。。
補足情報(ArtistAnimationではうまくいったのでそのコードを以下に示します。)
Python
1import numpy as np 2import pandas as pd 3import mpl_toolkits.axes_grid1 4from math import * 5from mpl_toolkits.axes_grid1 import make_axes_locatable 6from matplotlib import pyplot as plt 7from matplotlib import animation as animation 8 9 10# model N x M:PBC 11N, M=30,60 12#parameters 13pi=4*atan(1.0) #定数の定義 14dt, h=0.01,0.1 15h2=h*h 16DX,DY=0.0001,0.0001 17A,B=1.0,3.0 18maxX,maxY=0.0,0.0 19minX,minY=6.0,6.0 20T=10 21#preparation 22m=np.zeros([N+2,M+2],float) 23X=np.zeros([N+2,M+2],float) 24Y=np.zeros([N+2,M+2],float) 25XX=np.zeros([N+2,M+2],float) 26YY=np.zeros([N+2,M+2],float) 27rhsX=np.zeros([N+2,M+2],float) 28rhsY=np.zeros([N+2,M+2],float) 29for j in range(N+2): 30 for k in range(M+2): 31 s,t=pi*j/N,pi*k/M 32 X[j][k]=2.0+0.4*sin(s)*sin(3*t) 33 Y[j][k]=1.0+0.4*sin(s)*sin(t) 34diagX, offX=1-4*DX*dt/h2, DX*dt/h2#差分法の定数部分の定義 35diagY, offY=1-4*DY*dt/h2, DY*dt/h2 36def SX(u,v): 37 return (A-(B+1)*u+u*u*v) 38def SY(u,v): 39 return (B*u-u*u*v) 40fig = plt.figure() 41ax1 = fig.add_subplot(1, 2, 1) 42ax2 = fig.add_subplot(1, 2, 2) 43ims = [] 44for n in range(T+1): 45 for j in range (N+1): 46 for k in range (M+1): 47 rhsX[j][k]=dt*SX(X[j][k],Y[j][k]) 48 rhsY[j][k]=dt*SY(X[j][k],Y[j][k]) 49 rhsX[j][k]+=offX*(X[j+1][k]+X[j-1][k]+X[j][k-1]+X[j][k+1])+diagX*X[j][k] 50 rhsY[j][k]+=offY*(Y[j+1][k]+Y[j-1][k]+Y[j][k-1]+Y[j][k+1])+diagY*Y[j][k] 51 XX[j][k],YY[j][k]=rhsX[j][k],rhsY[j][k] 52 for j in range(0,N+2): 53 XX[j][0],XX[j][M+1]=XX[j][M],XX[j][1] 54 YY[j][0],YY[j][M+1]=YY[j][M],YY[j][1] 55 for k in range(0,M+2): 56 XX[0][k],XX[N+1][k]=XX[N][k],XX[1][k] 57 YY[0][k],YY[N+1][k]=YY[N][k],YY[1][k] 58 # renew 59 X, Y=XX, YY 60 #プロット 61 X_T=np.array(X).T 62 Y_T=np.array(Y).T 63 im1=ax1.imshow(X_T,interpolation='none',cmap='jet') 64 im2=ax2.imshow(Y_T,interpolation='none',cmap='jet') 65 ims.append([im1, im2]) 66ax1.set_xlabel("x") 67ax1.set_ylabel("y") 68ax1.invert_yaxis() 69ax2.set_xlabel("x") 70ax2.set_ylabel("y") 71ax2.invert_yaxis() 72divider1 = mpl_toolkits.axes_grid1.make_axes_locatable(ax1) 73divider2 = mpl_toolkits.axes_grid1.make_axes_locatable(ax2) 74cax1 = divider1.append_axes('right', '5%', pad='3%') 75cax2 = divider2.append_axes('right', '5%', pad='3%') 76fig.colorbar(im1,cax=cax1) 77fig.colorbar(im2,cax=cax2) 78fig.tight_layout() 79ani = animation.ArtistAnimation(fig, ims, interval=50) 80#ani.save("BZ_reaction.gif", writer="ffmpeg") 81plt.show()
陽解法を用いて反応拡散方程式を解いてそれをj×kの行列としてX,Yを出してそれをヒートマップとして出力しています。