拡散方程式の数値計算
プログラムの説明
2015
年度数理工学演習
I
補足資料
この資料は、偏微分方程式の数値計算の中でも、特に「拡散方程式」や「自 由粒子の Schr¨odinger 方程式」などの、一次元の方物型 PDE を数値計算す る場合において、どのように数値計算を行っているのかのあらましを解説す るものである。拡散方程式を数値計算する自らの Python プログラムを元に 簡単に説明する。 熱方程式については以下の通りである。u(t, x) を位置 x、時刻 t での温度 とした場合、その時間発展は下記の方程式 ∂ ∂tu(x, t) = k ∂2 ∂x2u(t, x) で記述される。k は拡散係数である。詳しい話は教科書や資料等にまかせる として、この係数を、k = i¯h 2mとしたものが自由粒子の Schr¨odinger 方程式 に相当する。¯h は換算プランク定数で、m は自由空間中の粒子 1 個の質量で ある。正確に書くならば、ψ(x, t) を粒子の波動関数として、 i¯h∂ ∂tψ(x, t) = 1 2m ( −i¯h ∂ ∂x )2 ψ(x, t) となる。 このタイプの偏微分方程式については、いわゆる Cauchy 問題と呼ばれる 初期値問題、すなわち空間領域 x :−∞ ≤ x ≤ ∞ に対して、与えた初期状態 が方程式に従ってどのように時間発展していくか考える場合、Fourier 変換に よって一般解を記述することが可能であるので、少しここでまとめておく。 途中の計算についてはテキストや資料に詳細が記載されているので、ここ では述べないが、拡散方程式として上記の 一般的に、温度関数や波動関数、その他位置 x について|x| → ∞ のとき値 がゼロに向かう物理現象に Fourier 変換が適用される。 ∂ ∂tu(x, t) = k ∂2 ∂x2u(t, x) を、Fourier 変換対 F (ξ) =F[F (x)](ξ) ≡ √1 2π ∫ ∞ −∞ F (x)e−iξxdx F (x) =F−1[F (ξ)](x)≡ √1 2π ∫ ∞ −∞ F (ξ)eiξxdξ を用いてその一般解を求めた場合、その厳密解は、 u(x, t) =√ 1 4kπt ∫ ∞ −∞ e−(x−x′)24kt u(x′, 0)dx′
として畳み込み積分の形で解が記述される（積分変数は x′に置き換えた）。 ここで√1 4kπte −(x−x′)2 4kt ≡ G(x − x′, t) とまとめた場合、 u(x, t) = ∫ ∞ −∞ G(x− x′, t)u(x′, 0)dx′ となる。G はいわゆるインパルス応答関数、Green 関数にあたるものである。 同じタイプの方程式として、自由粒子の Schr¨odinger 方程式 i¯h∂ ∂tψ(x, t) = 1 2m ( −i¯h ∂ ∂x )2 ψ(x, t) を解く場合を考える。変換対として F (p) =F[F (x)](p) ≡ √1 2π¯h ∫ ∞ −∞ F (x)e−ipxh¯ dx F (x) =F−1[F (p)](x)≡ √1 2π¯h ∫ ∞ −∞ F (p)eipxh¯ dp を採用し（これは ¯h = 1 となるような単位系 (原子単位系) で計算するなら上 記変換対と同じ形である）、一般解を求めると下記を得る。 ψ(x, t) = √1 2π¯h ∫ ∞ −∞ √ m ite im(x−x′)22¯ht ψ(x′, 0)dx′ こちらについても、√1 2π¯h √m ite im(x−x′)22¯ht = U (x− x′, t) と置くと、上記の 一般解が ψ(x, t) = ∫ ∞ −∞ U (x− x′, t)ψ(x′, 0)dx′ と表される。U (x− x′, t) は先程の G(x− x′, t) に相当する。 このように、今回取り上げたタイプの偏微分方程式は空間領域として x : −∞ ≤ x ≤ ∞ を考える問題の場合、フーリエ変換を用いて厳密に一般解が 記述できる。但し具体的に初期状態として u|t=0や、ψ|t=0を与えた場合の時 間発展まで厳密解が記述できるかはそのときに上記の積分が実行可能か否か による。
上記のような無限領域での放物型 PDE の話はここまでにして、本題に戻 る。有限の空間領域での拡散問題を数値計算で解く場合を考えよう。実際に プログラムを見ながら、数値計算する場合どのように式が差分化されるのか を見ていく。この資料では、コード内で行っている計算のプリポスト処理（数 値計算を行う前後のデータの準備や可視化する際の処理）に関しては言及せ ず、数値計算に関わる箇所のみ簡単に解説する。まずは実際のコードを見て みよう。 下記のプログラムは以下のような拡散型の、初期値・境界地問題の数値解を FTCS 法 (後で解説する) で計算するコードであり、t = 0、t = 0.01、t = 0.02 で解をプロットしたものが図１である。
使い方であるが、Python と拡張モジュール「NumPy」「matplotlib」が PC にインストールされた状況で上記プログラムを実行すれば、自動的に計算結 果の図が表示される。
equation : ∂t∂u(x, t) = k∂x∂22u(t, x) (x : 0 < x < 1 , t : 0 < t≤ 0.02 , k = 0.5)
B.C : u(0, t) = u(1, t) = 1
I.C : u(x, 0) = sin(5πx) + sin(30πx) + sin(80πx) + 1
=======================================================
#-*-coding:utf-8-*-import numpy as np
from numpy import linalg as LA import matplotlib.pyplot as plt k=0.5 N=100 dt=0.0001 step=200 x_ini=0.0 x_end=1.0 dx=(x_end - x_ini)/N dx2=dx**2 r=k*(dt/dx2) def IC(): global u u=np.zeros(N+1) u[0]=1.0 u[N]=1.0 #IC:三角型 # for i in range(1,N): # x = x_ini + i*dx
# if 0<x and x<=0.5: # u[i]=x+1 # elif 0.5<x and x<1: # u[i]=-x+2 #IC:三角関数の線形和 x=np.linspace(x_ini,x_end,N+1)
y=np.sin(5*np.pi*x) + np.sin(30*np.pi*x) + np.sin(80*np.pi*x) +1 for i in range(1,N): u[i]=y[i] M=np.linspace(x_ini,x_end,N+1) plt.plot(M,u) def evolve_mat1(): #FTCS 法による計算 A=np.zeros((len(u),len(u)),dtype=float) A[0,0]=A[len(u)-1,len(u)-1]=1 for i in range(1,len(u)-1): A[i,i-1]=r A[i,i]= 1.0 - 2*r A[i,i+1]=r return A
def evolve_mat2(): #Crank-Nicolson 法による計算 myu=r/2 A=np.zeros((len(u),len(u)),dtype=float) B=np.zeros((len(u),len(u)),dtype=float) A[0,0]=A[len(u)-1,len(u)-1]=1 B[0,0]=B[len(u)-1,len(u)-1]=1 for i in range(1,len(u)-1): A[i,i-1]=1 B[i,i-1]=-1 A[i,i]= -2-1/myu B[i,i]= 2-1/myu A[i,i+1]=1 B[i,i+1]=-1 C=np.dot(LA.inv(A),B) return C def evolve(m,b): for i in range(1,step+1): t= i*dt u1=np.dot(m,b) if (i%100)==0 :
A=np.linspace(x_ini,x_end,N+1) plt.plot(A,u1) b=u1[:] else: plt.legend(["t=0.00","t=0.01","t=0.02"],loc=2) plt.xlabel("x") plt.ylabel("u") plt.axhline(y=0, ls=’-’, lw=1, color=’black’) plt.title("1-D diffusion-eq") plt.show() def main(): IC() Mat=evolve_mat1() # Mat=evolve_mat2() evolve(Mat,u) if __name__ == "__main__": main() ======================================================= 図 1: 上記の初期値・境界地問題の数値解
方程式を数値計算する上で、まずは微分を差分に離散化する。差分については 色々と方法があるが、一番簡単な例として FTCS（ForwardTimeCenterSpace） 法について説明する。位置 x を空間インデックス i 番目、時刻 t をステップ n 番目として、u(x, t) = uni と表現する場合、上記方程式の微分演算子の差 分は、 ∂ ∂tu n i = ( un+1i − uni)/∆t ∂2 ∂x2u n i = ( uni+1− 2uni + uni+1)/∆x2 となる。 このように方程式を離散化することで、 un+1i = uni + k∆t ∆x2 ( uni+1− 2u n i + u n i−1 ) として、1 ステップ先の各点各点を一点ずつこの差分式で計算することが 出来る（陽解法）。 空間全体を N 分割すると考えると、空間節点数は N+1 個であり、空間イ ンデックス i は i : 0≤ i ≤ N となる。空間右端の節点 i = 0 と左端の節点 i = N が境界に相当する。両境界上で u = C(定数) とディリクレ条件を課す とすれば、全てのステップ n で un 0,Nが不変であるから、全ての i について上 記の差分式を連立すれば、(k∆t ∆x2 ≡ r と置いて) un+10 un+11 .. . un+1i .. . un+1N−1 un+1N = 1 0 · · · 0 r 1− 2r r 0 · · · · · · 0 0 r 1− 2r r 0 · · · 0 .. . . .. ... 0 · · · · · · r 1− 2r r 0 0 · · · · · · · · · r 1− 2r r 0 · · · 1 un 0 un1 .. . un i .. . un N−1 un N 上記で行列を A とおけば、 − →un+1 = A−→un と表現される。すなわち A は時間発展演算子としての意味がある。 大抵のテキストでは境界値が不変の場合、−→u から u0,Nは除いて（−→un+1= A−→un+ −→b , 但し−→u や−→b は (i : 1≤ i ≤ N − 1) についての要素数 N − 1 個の一 次元ベクトルとして）、各ステップでの時間発展を計算している。そちらの方 が余計な演算箇所が減るため通常はそうすべきである。ただ、上記のように
書くと初めから各ステップのデータ列に境界値が並んでいるので、計算終了 後のポスト処理が容易になる為に敢えてこのまま計算している。
方程式を差分すると、この様な形を持っているということで敢えて書いた が、元々この離散化手法（時間 (time)：前進 (forward) 差分、空間 (space)： 中心 (center) 差分なので FTCS 法などと呼ばれる）だと、1 ステップ後の各 点の u を上記差分式から一点一点求めることが出来る（陽解法）為、1 ステッ プ後の各点の値を一斉に求める必要はない。 後はこの計算を何度も繰り返し行えば、各時刻における、各点上での u の 値が求まってゆく。計算終了後、横軸を空間座標、縦軸をある時刻 t での u として出力すれば、時刻 t における u の振る舞いが見られる。 ここでもう一度コードを見てほしい。 「def evolve_mat1():」から「return A」
までが、上記の行列 A を実際にコーディングしている部分である。 A=np.zeros((len(u),len(u)),dtype=float) という部分で、行と列いずれについ ても、初期の一次元要素列 u と同じ長さの全要素 0 の二次元要素列 A を用意 し、その下の行では、この要素列に上記で説明したような行列 A になるよう 値を入れている。 − →unから −→ un+1を求めるための方法として、他の離散化手法も存在する。 その 1 つに Crank-Nicolson 法がある。プログラムではこちら数値スキームで 時間発展を計算できるようにもしてある。これは時間精度が 2 になる解法で ある。3 年後期の計算科学 I の講義でも取り上げられるはずなので、ここで は解説はしない。 ところで、FTCS 法では与えるパラメータしだいでは演算を繰り返す毎に 計算される値が発散してゆく現象が見られる（数値不安定性）。これに関して は時間発展演算子としての行列 A の最大固有値 (スペクトル半径と呼ばれた りする) に注目すればよい。行列の固有値というのは、ベクトルに行列を作 用させる場合、その固有値に対応する固有ベクトルの方向に向かってそのベ クトルを何倍にして移すのかを表している。すなわち最大固有値として 1 よ りも大きい値を持つ場合、 −→un+1 >|−→un| となり、計算を進めるほど値が大 きくなってゆく。FTCS 法で数値不安定性を避けて計算を進める為には、パ ラメータが満たすべき条件として、 r≤ 1 2 が課される。対して、Crank-Nicolson 法に関してはどのようなパラメータに関 しても無条件安定となる。これらの値については差分方程式に「VonNeumann の方法」というのを用いれば得られる。この方法は解をいろいろな波数成分 の波に分解して、一つ一つの成分の増幅率を調べるといった意味合いがある。
おまけで最後の方に、今回の問題で各々の数値スキームの場合の時間発 展演算子 A の全固有値の実部、虚部をプロットした図を載せた。と、いっ ても拡散方程式には複素数が含まれていないため、取り得るのは実固有値 のみである。黄色いのは単位円の周および内部を色付けしたもの。図内の
「|max(EigV alue)|, |min(EigV alue)|」は、プロットされている全固有値の
ノルムをとった場合の最大値と最小値である。 図 2: FTCS 法を用いる場合の時間発展演算子としての行列 A の全固有値 図 3: Crank-Nicolson 法を用いる場合の時間発展演算子としての行列 A の全 固有値 図２，３から、FTCS 法に関しては負の固有値も取りえるが、Crank-Nicolson 法の場合は完全に固有値は正の値のみ取るという違いが見られる。しかしい ずれにしても、時間発展演算子 A の固有値の大きさは全て１以下となってい
るので、時間発展を計算し続けても値が発散することが無いことが視覚的に 分かる。
̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶-BOX (A)(B)(C) (A) (B)Python Numpy FTCS Crank-Nicolson (C) Schr ¨odinger 58Page Gauss dtype=complex ̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶-1996 I 1996-1998 2008 2011 ¡http://www.somedalab.c.u-tokyo.ac.jp/Classes/NeuRD_H20_1.pdf ¿ 2015-6-4
