両端が断熱的な導線の熱伝導を記述する, 次の斉次 Neumann 境界付きの熱方程式の初期値境界値問題
を考えます*1.
ここで, , を表します.
前進差分法とスペクトル法を用いる方法の2通りでの数値計算例のデモを作成しました. *2
コードがメインで, 数値解析の手法については多くを説明しません.
意味が通らない箇所があれば気軽にコメントをよろしくお願いします.
いずれの場合も初期値は
により与え, 時刻 を の範囲で解きます.
import matplotlib.pyplot as plt #line domain xmin = 0.0 xmax = np.pi nx = 64 dx = (xmax-xmin)/nx #time interval tmin = 0.0 tmax = 1.0 nt = 1024 dt = (tmax-tmin)/nt c = 1.0 #spreading coefficient M = 51 #maximam fourier number #initialize F = np.zeros(nx) T = np.linspace(tmin, tmax, nt, endpoint = False) X = np.linspace(0, np.pi, nx, endpoint = False)+dx/2 F = np.cos(4*X) for i in range(nx//2): F[i] = 1
1. 前進差分法で解こう
ここでは空間メッシュ数を 刻み幅を とし,
時間メッシュ数を , 時間刻み幅を とします.
熱方程式は解 がある程度滑らかだという仮定のもと, 時間微分を前進差分, 空間微分を中心差分に置き換えることで,
と離散化されます.
rsq = c*dt/(dx**2) #constant #solution v = np.zeros((nt,nx)) #initial data v[0] = F print(len(F)) #simulation for t in range(nt-1): for x in range(1,nx-1): v[t+1,x] = (1-2*rsq)*v[t,x]+rsq*(v[t,x-1]+v[t,x+1]) #Neumann condition v[t+1,0] = v[t+1,1] v[t+1,nx-1] = v[t+1,nx-2]
これを可視化すると次のような動画が得られます.
fig1 = plt.figure() fig1.set_dpi(100) ax = fig1.add_subplot(1,1,1) def animate(t): ax.clear() plt.ylim([-1.2,1.5]) p, = plt.plot(X,v[t]) time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes) time_text.set_text('t = %.2f/%.2f,' % ((t+1)/nt,tmax)) anim = animation.FuncAnimation(fig1,animate,frames=nt,interval=60,repeat=False) plt.show()
2.スペクトル法で解こう
記事
trtmfile.com
で導入したフーリエ変換を用いて数値解を求めます. まずは初期値をフーリエ変換しましょう.
F_M = [np.zeros(nx) for m in range(M)] #Fourier coefficients A = [0]*M for m in range(M): A[m] = (np.dot(np.cos(m*X),F))/nx F_M[0] = np.cos(0*X)*A[0] for m in range(1,M): F_M[m]=2*np.cos(m*X)*A[m] fig1 = plt.figure() fig1.set_dpi(100) ax = fig1.add_subplot(1,1,1) def animate(t): ax.clear() plt.ylim([-1.2,1.5]) p, = plt.plot(X,F_M[t]) time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes) time_text.set_text('t = %.2d/%.2d,' % (t,M-1)) anim = animation.FuncAnimation(fig1,animate,frames=M,interval=600,repeat=False) plt.show()
端点を含めて滑らかな関数なので, ギブス現象も起こらずによく近似されています.
スペクトル法を用いて解を求めましょう.
u = [[[0 for x in range(nx)] for m in range(M)] for t in range(nt)] #u[t][m][x] u[0]= F_M for m in range(M): for t in range(nt): u[t][m] = np.dot(np.exp(-t*dt*m*m),u[0][m]) for m in range(1,M): for t in range(nt): u[t][m] += u[t][m-1]
これを可視化すると次のような動画が得られます.
fig1 = plt.figure() fig1.set_dpi(100) ax = fig1.add_subplot(1,1,1) def animate(t): ax.clear() plt.ylim([-1.2,1.5]) p, = plt.plot(X,u[10*t][M-1]) time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes) time_text.set_text('t = %.2f/%.2f,' % ((t+1)/nt,tmax)) anim = animation.FuncAnimation(fig1,animate,frames=nt,interval=60,repeat=False) plt.show()
3. それぞれの解を比較しよう.
fig1 = plt.figure() fig1.set_dpi(100) ax = fig1.add_subplot(1,1,1) def animate(t): ax.clear() plt.ylim([-1.2,1.5]) p, = plt.plot(X,v[10*t]) p, = plt.plot(X,u[10*t][M-1]) time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes) time_text.set_text('t = %.2f/%.2f,' % ((t+1)/nt,tmax)) anim = animation.FuncAnimation(fig1,animate,frames=nt,interval=60,repeat=False) plt.show()
2つの数値解はそこそこ近いものになっていますね!
教えてくれたけますさん, ありがとうございます.
コードはこちらにございます.
github.com