計算物理学２第 13 回：微分方程式の解法その３ver. 2018/7/20
偏微分方程式 (partial diﬀerential equation) の数値解法について説明します。
一般的に二階二変数の偏微分方程式は

A(x, y)

∂2U
∂2U
∂U
∂2U
∂U+ C(x, y) 2 + D(x, y)+ E(x, y)+ F (x, y) = 0+ 2B(x, y)
2
∂x
∂x∂y
∂y∂x
∂y(1)と書けますが、これは B 2 (x, y) − A(x, y)C(x, y) の符号により 3 つの型に分類されます。
すべての x, y について B 2 (x, y) − A(x, y)C(x, y) < 0 のものは楕円型偏微分方程式 (elliptic partialdiﬀerential equation) と呼ばれます。例としては Poisson 方程式∂2U∂2U+= −4πρ(x, y)∂x2∂y 2(2)があります。Poisson 方程式は与えられた U (x, y) の境界条件を用いて問題を解きます。

B 2 (x, y) − A(x, y)C(x, y) > 0 のものは双曲型偏微分方程式 (hyperbolic partial diﬀerential equation) と呼ばれ、例えば 1 次元の波動方程式∂2U∂2U− c2 2 = 02∂t∂x(3)が該当します。

B 2 (x, y) − A(x, y)C(x, y) = 0 のものは放物型偏微分方程式 (parabolic partial diﬀerential equation) と呼ばれ、1 次元の熱方程式∂U
λ ∂2U=∂tCV ∂x2(4)が該当します。

1 楕円型偏微分方程式
Poisson 方程式は時間に依存しない関数の空間分布を求める問題です。通常境界値が与えられており、差分法によって解くことができます。
例えば三点公式を用いると

U (xi+1 , yj ) − 2U (xi , yj ) + U (xi−1 , yj )∂2U(xi , yj ) =,∂x2(∆x)2∂2U
U (xi , yj+1 ) − 2U (xi , yj ) + U (xi , yj−1 )(xi , yj ) =,∂y 2(∆y)2(5)(6)などのように U (xi , yj ), (i = 1, · · · , Nx , j = 1, · · · , Ny ) を用いてすべての点での二階微分が表現できます。
これら Nx Ny 個の成分を一次元に並べたものを U とすると Poisson 方程式は

AU = b(7)の形の連立一次方程式となります。b には ρ(xi , yj ) からの寄与の他、境界条件も入ります。連立一次方程式の問題は以前取り扱った Jacobi, Gauss-Seidel, SOR, 逆行列などの方法を用いて解くことができます。

12 双曲型偏微分方程式
次回時間があれば取り扱う予定です。

3 放物型偏微分方程式
放物型偏微分方程式として 1 次元時間依存 Schrödinger 方程式を扱います。波動関数は複素数であることに注意してください。方程式は[]∂h̄2 ∂ 2ih̄ ψ(x, t) = −+ V (x) ψ(x, t)∂t2m ∂x2(8)の形で、初期値としては時刻 t = 0 での波動関数 ψ(x, t = 0) が与えられており、その時間発展を方程式に従って数値計算で計算します。
例えば初期値としてはガウス波束

ψ(x, 0) = Ce−(x−x0 )2a2eik0 x(9)を考えます。係数 C は規格化因子です。a は波束の広がりのパラメター、x0 は波束の中心、p0 = h̄k0 は初期運動量です。
空間 0 < x < L を考え、波動関数はその内側に存在するとします。(外側ではポテンシャルが無限大と考えます) 空間を両端を除いた N 個の格子点で分割し、(xi = −L + i∆x,2L∆x =, i = 0, 1, 2, · · · , N + 1
N +1)(10)とします。時間についても適当なタイムステップ ∆t を設定し、tj = j∆t,

(j = 0, 1, 2, · · · )(11)とします。
波動関数の時間発展は形式的に[]iψ(x, t) = exp − Ĥt ψ(x, t = 0)h̄(12)とかけます。ここで Hamiltonian 演算子はĤ = −h̄2 ∂ 2+ V (x)2m ∂x2(13)で与えられます。

3.1 陽解法
式 (12) を用いて t = tj から tj+1 までタイムステップ ∆t だけ時刻を進めることを考えます。このとき、Taylor 展開の 1 次までをとると)(iψ(xi , tj+1 ) = 1 − Ĥ∆t ψ(xi , tj ) + O((∆t)2 )h̄2(14)

とかけます。これは、時間依存 Schrödinger 方程式の時間に関する偏微分を前進差分で置き換えることと同じで、さらに空間に関する二階微分を三点公式を用いて表すと、ih̄h̄2 ψ(xi+1 , tj ) − 2ψ(xi , tj ) + ψ(xi−1 , tj )ψ(xi , tj+1 ) − ψ(xi , tj )+ O(∆t) = −+ V (xi )ψ(xi , tj ) + O((∆x)2 ),∆t2m(∆x)2(15)となります。この方程式は t = tj での波動関数を用いて t = tj+1 の波動関数を計算することができる式となっており、原理的にこれで時間ステップを進めていくことができます。

[]∆th̄2 ψ(xi+1 , tj ) − 2ψ(xi , tj ) + ψ(xi−1 , tj )ψ(xi , tj+1 ) = ψ(xi , tj ) − i−+ V (xi )ψ(xi , tj ) + O((∆t)2 , ∆t(∆x)2 )h̄2m(∆x)2(16)これは陽解法 (explicit method) と呼ばれます。Euler 法や Runge-Kutta 法も陽解法です。(特にこの放物型偏微分方程式の場合は FTCS(Forward-Time Central-Space) 法と呼ばれます。時間について前進差分、空間について中心差分、の意味です) 陽解法は最も簡単に実装できますが、時間発展の間に解が安定しないのが問題です。実際に式 (16) を用いて時間発展を行うと誤差が積もって発散してしまいます。少し改良して時間微分を中心差分で置き換えると安定して解が求められます。

ih̄h̄2 ψ(xi+1 , tj ) − 2ψ(xi , tj ) + ψ(xi−1 , tj )ψ(xi , tj+1 ) − ψ(xi , tj−1 )+ O((∆t)2 ) = −+ V (xi )ψ(xi , tj ) + O((∆x)2 ),2∆t2m(∆x)2(17)ψ(xi , tj+1 ) =ψ(xi , tj−1 ) − i[]h̄2 ψ(xi+1 , tj ) − 2ψ(xi , tj ) + ψ(xi−1 , tj )2∆t−+V(x)ψ(x,t)+ O((∆t)3 , ∆t(∆x)2 )ii jh̄2m(∆x)2(18)この式では t = tj+1 での波動関数を計算するためには t = tj と tj−1 の２つの時刻での波動関数が必要になります。時刻 t = ∆t でのみ tj−1 の波動関数が準備できないため、初回だけ式 (16) を使い、その後は式 (18) を用いることで時間発展が計算できます。
時間依存 Schrödinger 方程式の場合に特に問題となるのは時間発展で波動関数の規格化が徐々に破れていくことです。t = 0 では

∫ Ldx|ψ(x, t)|2 = 1(19)−Lは満たされるように規格化定数 C を決めますが、時間ステップを進めていくとこの条件が徐々に破れていきます。また、∆x と ∆t の値の選び方によっても解が発散することがあります。

3.2 陰解法
陰解法 (implicit method) は演算子を含む式で書くと(ψ(xi , tj+1 ) =i
1 + Ĥ∆th̄)−1ψ(xi , tj )(20)となり、時間に関して後退差分を用います。

ih̄

ψ(xi , tj+1 ) − ψ(xi , tj )h̄2 ψ(xi+1 , tj+1 ) − 2ψ(xi , tj+1 ) + ψ(xi−1 , tj+1 )+ V (xi )ψ(xi , tj+1 ) (21)=−∆t2m(∆x)23

この式を解く際に既知の波動関数は ψ(xi , tj ) のみですので、これを右辺に移動した()h̄2h̄2h̄h̄2h̄ψ(x,t)+−+V(x)+iψ(x,t)+ψ(xi−1 , tj+1 ) = i ψ(xi , tj )i+1j+1iij+12222m(∆x)m(∆x)∆t2m(∆x)∆t(22)
を解きます。これは ψ(xi , tj+1 ) (i = 1, 2, · · · N ) に関する連立一次方程式となっており、時間 tj ごとに連立方程式を解く必要がありますが、陽解法に比べて解は安定します。

3.3 Crank-Nicolson 法
Crank-Nicolson 法では時間発展の演算子を(expi∆tĤh̄)1i
Ĥ∆t2 h̄=1i
Ĥ∆t1+2 h̄
1−(23)と近似します。これを用いると tj での波動関数と tj+1 での波動関数の関係は()()
1i1iĤ∆t ψ(xi , tj+1 ) = 1 −Ĥ∆t ψ(xi , tj )1+2 h̄
2 h̄(24)となります。Crank-Nicolson 法は時間微分に関しては tj と tj+1 の波動関数を用いたもの、Hamiltonian に関しては tj , tj+1 で評価した値の平均を用いたものに対応します。

ih̄[]h̄2 ψ(xi+1 , tj+1 ) − 2ψ(xi , tj+1 ) + ψ(xi−1 , tj+1 )ψ(xi , tj+1 ) − ψ(xi , tj ) 1= −+V(x)ψ(x,t)
ii j+1∆t22m(∆x)2][1h̄2 ψ(xi+1 , tj ) − 2ψ(xi , tj ) + ψ(xi−1 , tj )+ −+ V (xi )ψ(xi , tj )(25)
22m(∆x)2

この式も陰解法のときと同様に連立方程式の形にして時刻 tj ごとに解く必要があります。(Crank-Nicolson法も陰解法です) Crank-Nicolson 法では波動関数のユニタリー性が満たされます。

4 可視化
時間依存する波動関数 ψ(x, t) をアニメーションとして表示してみましょう。次の演習問題 (46) に対応したサンプルプログラムを講義資料ページに載せていますのでこれを使ってください。データファイルが大量に生成されますので、(46) を実行する前に mkdir で新しい専用のディレクトリを作り、そこで作業をすることを推奨します。

4.1 データファイルの処理
まず時刻 t ごとに波動関数のデータファイルを作成します。1 つのファイルにまとめても構いませんが、時刻ごとにデータファイルを分けるのが扱いやすいでしょう。(46) のプログラムを実行すると、

wavefunction 0000000.dat という名前で (0 のところには時刻ステップが入ります) x と Reψ(x, t), Imψ(x, t)が 3 列に並んだファイルが生成されます。
時刻ごとに出力するデータファイルの名前を変更するにはプログラムの中で少しテクニックが必要です。

Fortran では時間ステップ数を表す整数型変数 ti を 7 桁の整数で表示し (I7.7)、文字型変数 tindex に WRITE
4

文で変換しています。これを用いてファイル名を文字列の連結で作成しています。Fortran での文字列の結合//で行います。文字列は十分なサイズ (この例では 50 文字) で定義して、定義したサイズ分使わない場合は余った右側に空白が入ります。これは文字列を連結するときに邪魔ですので TRIM 関数で取り除きます。C
言語では sprintf 関数を用いると簡単に整数型変数の情報を含むファイル名が作成できます。

4.2 gnuplot を用いた可視化データファイルのセットが出来上がったら gnuplot で gif アニメーションを作成します。アニメーションは時刻を変えたデータファイルを読み込んで何度もグラフを描画することで作成されます。２つの設定ファイルが必要となります。plot.gpl はメインのファイルで、data.gpl は繰り返し実行されるデータファイルを描画する部分です。

plot.gpl の冒頭でset term gif animateとします。出力ファイル名は拡張子”.gif” を使います。時刻での描画ごとに x 軸や y 軸のスケールがずれるとアニメーションになりませんので set xrange や set yrange であらかじめ固定しておきます。ファイル名を指定する変数 tmin, tmax, dt(増え幅) を定義しておき、data.gpl を呼び出します。

data.gpl ではもし t が設定されていなければ tmin が設定され、それを用いてデータファイルの名前を設定します。データファイルの内容をプロットし、t を dt 分増やし、tmax を超えない間 reread コマンドでdata.gpl ファイルを最初から実行します。tmin から tmax までのループを実行していることに対応します。
ターミナルで% gnuplot plot.gplとすると gif アニメーションが作成されます。gif ファイルを開くには画像ビューアやブラウザが使えます。
ファイルマネージャーでディレクトリを開き、ダブルクリックしてください。

4.3 linux でのファイル操作
多くのデータファイルを作成しますのでディレクトリの中に大量のファイルが作成されます。容量等にも注意を払ってください。

% ls -lhと-h オプションをつけると各ファイルのサイズが表示されます。また、今いるディレクトリのデータ容量を調べるには% du -h .とします。”.” はカレントディレクトリの意味です。
複数のファイルをまとめて操作する場合はワイルドカード*を用います。*は該当するものすべてを表します。たとえば% ls -lh wavefunction_00003*.dat5とすると、時刻が 0000300 から 0000399 のファイルだけ表示されます。削除したいときは% rm wavefunction_*.datとすると、wavefunction からはじまって.dat で終わるすべてのファイルが削除されます。

% rm *とするとディレクトリにある全てのファイルが削除されますので気をつけてください。特に wavefunc-tion *.dat とするつもりがスペースが入って wavefunction * .dat などとすると、スペースで区切られた３つのファイルと解釈されるため、すべてのファイルが消えます。

5 演習問題(46) h̄2 /2m = 20 (MeV)、∆t = 10−4 (h̄/MeV), 幅 L = 100 (fm), 空間の分割数 N = 1000 での井戸型ポテンシャルで t = 0 での初期波動関数をガウス波束 (a = 5 fm, x0 = 20 (fm), k = 5 (fm−1 )) として時間依存

Schrödinger 方程式を陽解法で解いたサンプルプログラムを実行し、t = 5000(h̄/MeV) までの波動関数の時間発展を表示せよ。

(47) x = L/2 に幅 b (fm)、高さ V0 (MeV) のポテンシャル障壁を設置し、量子トンネル現象が見えるように b と V0 の値を調整せよ。ポテンシャルも図中に書き込むようにするとよい。

(48) 調和振動子ポテンシャルを設置した場合はどうか。
(49) 周期境界条件 ψ(x + L, t) = ψ(x, t) に境界条件を変更せよ。
(50) Crank-Nicolson 法を用いるように改良せよ。

6