計算物理学２第 14 回：微分方程式の解法その４ver. 2018/7/23
前回は放物型の偏微分方程式の解法について解説を行いました。陽解法では時間について前進差分を行うと解が発散してしまうという問題がありました。今回はまず微分方程式の数値解の安定性を調べる方法を見ていきます。

1 von Neumann の安定性解析
時間依存の 1 次元 Schrödinger 方程式の陽解法による解として[]h̄2 ψ(xi+1 , tj ) − 2ψ(xi , tj ) + ψ(xi−1 , tj )∆t−+V(x)ψ(x,t),ii jh̄2m(∆x)2[]2∆th̄2 ψ(xi+1 , tj ) − 2ψ(xi , tj ) + ψ(xi−1 , tj )ψ(xi , tj+1 ) =ψ(xi , tj−1 ) − i−+ V (xi )ψ(xi , tj )h̄2m(∆x)2

ψ(xi , tj+1 ) =ψ(xi , tj ) − i(1)(2)の２つについて前回説明しました。式 (1) は時間微分を前進差分で近似したもの、式 (2) は中心差分を用いたものです。両者はほとんど同じ式に見えますが、式 (1) を使うとすぐに発散し、式 (2) だと安定して解を求めることができました。どの方程式 (アルゴリズム) を用いて時間発展を計算すれば安定して数値解が求められるかを分析するには von Neumann の安定性解析を用いて議論ができます。
ポテンシャルが一定であるとすると、方程式の固有解は平面波 (波数 k でラベルされる) でψ(xi , tj ) = ξ j eikxi(3)となります。任意の波動関数はこの固有関数の重ね合わせでかけ、時間依存 Schrödinger 方程式は ψ に関して線形ですのである k についての振る舞いを調べれば十分です。ここで、ξ は k に依存する複素関数で、時間ステップを進めるごとに、ξ の整数乗で変化する増幅係数 (amplification factor) と呼ばれるものです。誤差も同じ式で増えていきますので、もし |ξ| > 1 となる k が (固有モードが) 存在すれば、そのモードについては時間発展に対して不安定であることがわかります。この式 (3) を式 (1) に代入することによって ξ を求めることができます。

√ξ =1 − ir, |ξ| = 1 + r2([)]
k∆xh̄2∆t2
sin4+Vr≡ih̄
2m(∆x)22(4)(5)r は実数、ポテンシャルがなければ正の実数です。ポテンシャルがない場合でも、sin2 (k∆x/2) がゼロにならないような k については |ξ| > 1 となるため、式 (1) を用いて時間発展させると不安定であることがわかります。
式 (2) を用いて時間発展させた場合は、式 (3) を代入すると ξ に関する２次方程式となります。

ξ 2 + 2irξ − 1 = 0(6)この解は、|r| ≤ 1 のときは

ξ = −ir ±1√

1 − r2(7)この解は２つとも |ξ| = 1 となり、安定であることがわかります。
一方で、|r| > 1 のときは

ξ = i(−r ±√r2 − 1),|ξ|2 = 2r2 − 1 ± 2r√r2 − 1(8)となり、r > 1 の領域では + の符号の解は常に |ξ| > 1 であり不安定です。
これより式 (2) で安定に解を求めるにはポテンシャルがない場合は

r=2h̄ ∆t<12m (∆x)2(9)ととればよいことがわかります。これは座標の格子点 ∆x を小さくすれば、対応して ∆t も小さくする必要があることを示しています。演習問題 (46) では ∆t/h̄ = 10−4 (1/MeV), ∆x ≃ 10−1 (fm), h̄2 /2m = 20(MeVfm−2 ) としているため、この条件が満たされています。
また、陰解法の安定性についても見ておきます。陰解法ではih̄h̄2 ψ(xi+1 , tj+1 ) − 2ψ(xi , tj+1 ) + ψ(xi−1 , tj+1 )ψ(xi , tj+1 ) − ψ(xi , tj )=−+ V (xi )ψ(xi , tj+1 ) (10)∆t2m(∆x)2を用いて連立一次方程式を解くことで時間発展を計算します。これに式 (3) を代入すると、r を用いて

1= 1 + ir,
ξ|ξ| = √11 + r2(11)となり、常に |ξ| < 1 を満たして安定であることがわかります。陰解法では解は安定しますが、|ξ| < 1 である
ため、時間発展させると波動関数の絶対値が小さくなる問題があります。長時間時間発展を行うといずれ波動関数は消えてなくなってしまいますが、本来時間依存 Schrödinger 方程式は波動関数のノルムを保存します。
このような場合は |ξ| = 1 となるような解法が望まれます。

Crank-Nicolson 法では以下の式を用いて時間発展します。
[]
ψ(xi , tj+1 ) − ψ(xi , tj ) 1h̄2 ψ(xi+1 , tj+1 ) − 2ψ(xi , tj+1 ) + ψ(xi−1 , tj+1 )= −ih̄
+V(x)ψ(x,t)
iij+1∆t22m(∆x)2[]
1h̄2 ψ(xi+1 , tj ) − 2ψ(xi , tj ) + ψ(xi−1 , tj )+ −+V(x)ψ(x,t)(12)ii j22m(∆x)2
これに式 (3) を代入すると

ξ=r + 2i,−r + 2i|ξ| = 1(13)となり、任意の r に対して安定かつ |ξ| = 1 となりノルムも保存します。以上より、時間依存 Schrödinger 方程式を解く場合は陽解法で |r| < 1 の場合に式 (2) を用いるか、また、任意の r の値で Crank-Nicolson 法を用いるのがノルムが保存した安定解を求める方法であると言えます。von Neumann の解析の結果はもちろん解く微分方程式の形によって変わりますので注意してください。

2 双極型偏微分方程式
今回は 2 次元空間での波動方程式を扱います。

∂2U− c2∂t2(∂2U∂2U+∂x2∂y 22)=0(14)

の形をしています。これに t = 0 での初期関数 U (x, y, t = 0) を与え、その時間発展を見てみます。格子点

xi = i∆x, (i = 1, 2, · · · Nx ), yj = j∆y, (j = 1, 2, · · · Ny ) としこの点の上で関数を表現します。時間ステップ
も tk = k∆t, (k = 0, 1, 2, · · · ) とします。
波動方程式の二階の偏微分は三点公式を使うと差分に置き換えることができます。

∂2U
U (xi , yj , tk+1 ) − 2U (xi , yj , tk ) + U (xi , yj , tk−1 )(xi , yj , tk ) =,∂t2(∆t)2∂2U
U (xi+1 , yj , tk ) − 2U (xi , yj , tk ) + U (xi−1 , yj , tk )(xi , yj , tk ) =,∂x2(∆x)2U (xi , yj+1 , tk ) − 2U (xi , yj , tk ) + U (xi , yj−1 , tk )∂2U
(xi , yj , tk ) =,2∂y(∆y)2(15)(16)(17)偏微分ですので微分しない他の変数の格子点の値は変わりません。これらを波動方程式 (14) に代入すると

U (xi , yj , tk+1 ) =2U (xi , yj , tk ) − U (xi , yj , tk−1 )[]22 U (xi+1 , yj , tk ) − 2U (xi , yj , tk ) + U (xi−1 , yj , tk )+ c (∆t)(∆x)2][22 U (xi , yj+1 , tk ) − 2U (xi , yj , tk ) + U (xi , yj−1 , tk )+ c (∆t)(∆y)2(18)

この式は t = tk+1 での関数 U (x, y, tk+1 ) が t = tk と tk−1 での関数 U (x, y, t) を用いて計算できることを示しているので、この式をすべての xi , yj について解くことで波動方程式が解けます。
境界条件は固定端反射であれば、計算する格子点の外では関数 U はすべての時刻でゼロとしておきます。

U (x0 , yj , tk ) =U (xNx +1 , yj , tk ) = 0(for all j, k),(19)U (xi , y0 , tk ) =U (xi , yNy +1 , tk ) = 0(for all i, k),(20)

自由端反射の場合は、境界で一階偏微分がゼロとなります。

∂U∂U(x0 , yj , tk ) =(xNx +1 , yj , tk ) = 0,
∂x
∂x∂U
∂U
(xi , y0 , tk ) =(xi , yNy +1 , tk ) = 0,∂y∂y(for all j, k),(21)(for all i, k),(22)

計算する格子点の一つ外でも同じ値としておくことで、一階微分がゼロの条件が満たされます。(境界となる座標が固定端の時と少しずれますが)U (x0 , yj , tk ) = U (x1 , yj , tk ),U (xNx +1 , yj , tk ) = U (xNx , yj , tk ),(for all j, k),(23)

U (xi , y0 , tk ) = U (xi , y1 , tk ),U (xi , yNy +1 , tk ) = U (xi , yNy , tk ),(for all i, k),(24)

実際に波動方程式を解く前に von Neumann の安定性解析を行ってみます。固有関数

U (xi , yj , tk ) = ξ k eikx xi +iky yj(25)を代入すると、(())kx ∆xky ∆y22sinsin22r = c2 (∆t)2 +22(∆x)(∆y)
3(26)

を用いて ξ は 2 次方程式

ξ 2 − 2(1 − 2r)ξ + 1 = 0(27)の解となります。これが安定になるのは r < 1 のときで、この場合√ξ = 1 − 2r ± 2i r − r2 ,|ξ| = 1(28)]1
1+<1(∆x)2(∆y)2(29)となるため、[2c (∆t)2または ∆x = ∆y のときは (∆t/∆x)2 < 1/2c2 を満たすように時間ステップと格子間隔を決める必要があります。

3 可視化
gnuplot で 3 次元プロットを行うには splot を使います。gnuplot を起動後
splot "ファイル名" using 1:2:3 w lとすると 1 列目を x 座標、2 列目を y 座標、3 列目を z 座標とした 3 次元プロットが表示されます。with lines
オプションはすべてのデータ点を読み込んだ順番に線で結びますので計算結果をファイルに出力するときに連続していない点を出力するときは改行を入れておきます。

4 演習問題(51) 初期関数がガウス型[U (x, y, t = 0) = exp](y − y0 )2(x − x0 )2+,a2xa2y

U (x, y, t =∆t) = U (x, y, t = 0)

(30)
(31)

ax = ay = 3 (m), x0 = y0 = 50 (m), Lx = Ly = 100 (m), c = 10 (m/s), 分割数 Nx = Ny = 100,∆t = 10−3 (s) の場合の 2 次元波動方程式を自由端反射の境界条件で解くサンプルプログラムを実行し、結果を gif ファイルで図示せよ。

(52) 初期関数としてガウス型を２つおいた場合 (たとえば (x0 , y0 ) = (30, 30) と (30,70) においた場合) はどうなるか。

(53) 初期関数を x 方向のみに進む波に変更せよ。
(54) 境界条件を固定端反射に変更せよ (境界ですべての時刻で U (x, y, t) = 0)4