Pythonで拡散方程式の数値計算
数値計算 (理工系の基礎数学 8) に拡散方程式の数値計算法について書かれていましたので、実際に手を動かして数値計算してみました。
拡散方程式の例：棒の温度分布の時間変化
拡散方程式をモデルとした物理現象として、棒の温度分布の時間変化が挙げられます。上式は、細長い棒の両端の温度を0に固定した時（上式③）、初期温度分布（上式②）が時間と共にどのような温度分布になるかを示しています。
xは棒の各点の位置、u(x,t)は位置xにおける温度です。
初期条件②、境界条件③のもと拡散方程式①を数値計算してみます。Pythonでコードを書いて数値計算した結果、下図のように時間が経つにつれ、棒全体の温度が両端の温度（0℃）に収束していく様子を確認することができました。本に書かれていた図と同じような結果が得られ満足です。
コード
上図を導出するためのコードは以下の通りです。Jupyter notebook上で表示することを想定しています。
初期条件や境界条件を変えて挙動がどう変わるか色々試せます。
%matplotlib inline import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation # 各種定数の設定 t0 = 0 tM = 1 M = 500 dt = (tM-t0)/M # 時間間隔 x0 = 0 xN = 1 N = 10 dx = (xN-x0)/N # 空間間隔 a = dt/(dx**2) # 差分方程式に出てくる係数 ux0 = 0 # 境界条件 uxN = 0 # 境界条件 A = 2 # 初期条件の式の振幅 # 変数の初期化 u = np.zeros([M+1, N+1]) # 行：時間、列：ｘ # 初期条件 (t=0でのuの分布) u(x)=A*x(1-x) (0<=x<=1) for xi in range(N+1): u[0,xi] = A*(xi*dx)*(1-xi*dx) # 差分解 for ti in range(M): # 境界条件（両端の温度を0に固定) u[ti,0] = ux0 u[ti,N] = uxN # 両端以外を差分方程式により更新 for xj in range(1,N): u[ti+1, xj] = a*u[ti,xj+1] + (1-2*a)*u[ti,xj] + a*u[ti,xj-1] # matplotlibで表示 x = np.linspace(x0, xN, N+1) # maplotlibで表示するためのx軸の設定 fig = plt.figure() ax = fig.add_subplot(1,1,1) # t=0 ax.plot(x, u[0,:], color='red') ax.scatter(x, u[0,:], color='red', label='t=0 ms') # t=20*dt (dt=1/500) ax.plot(x, u[20,:], color='orange') ax.scatter(x, u[20,:], color='orange', label='t=40 ms') # t=40*dt ax.plot(x, u[40,:], color='green') ax.scatter(x, u[40,:], color='green', label='t=80 ms') # t=100*dt ax.plot(x, u[100,:], color='cyan') ax.scatter(x, u[100,:],color='cyan', label='t=200 ms') # t=500*dt ax.plot(x, u[500,:], color='blue') ax.scatter(x, u[500,:], color='blue', label='t=1000 ms') # 軸ラベルの設定 ax.set_xlabel('x', fontsize=12) ax.set_ylabel('u', fontsize=12) # サブプロットタイトルの設定 ax.set_title('Temperature Profile') # 凡例の設定 ax.legend(loc='upper left', bbox_to_anchor=(1,1)) # 表示 plt.show()
