計算物理学２第 12 回：微分方程式の解法その２ver. 2018/7/12
定常状態にある 1 次元 Schrödinger 方程式−h̄2 d2 ψ+ V (x)ψ(x) = Eψ(x)2m dx2(1)の解法について考えます。二階の常微分方程式ですが、初期値が与えられていないこと、また、固有エネル
ギー E も未知であるところが前回扱った微分方程式と異なります。

k 2 (x) = 2m[E − V (x)]/h̄2 とおくと Schrödinger 方程式はd2 ψ= −k 2 (x)ψ(x)dx2(2)となります。
まずは簡単な場合として無限に深い井戸型ポテンシャル{0, (|x| ≤ a),
V (x) =∞ (|x| > a)

(3)を考えますが、式では V (x) を残しておきます。井戸型ポテンシャルでは境界条件は ψ(a) = ψ(−a) = 0 です。また、k(x) は |x| ≤ a で定数となります。

1 Runge-Kutta 法による解法
前回説明した Runge-Kutta 法 (以下では 4 次のものを使います) を Schrödinger 方程式に対して使うこともできます。二階微分方程式ですので ψ(x) の一階微分に対応する関数 ϕ(x) を導入し、連立一階微分方程式dψ=ϕ(x),dxdϕ= − k 2 (x)ψ(x)dx(4)(5)を解きます。k(x) の値には未知量の E が入っているので適当な値をセットします。
初期値としては、井戸型ポテンシャルの場合は例えば x = −a で ψ(−a) = 0, ϕ(−a) = 1 とできます。

ψ(−a) = 0 は井戸型ポテンシャルの境界条件ですが、ϕ(−a) は任意の値を用いて構いません。用いる値によっては規格化されていない波動関数が求まるだけですので、必要に応じてあとで波動関数を規格化します。
この初期条件で Runge-Kutta 法で方程式を x = −a から解き進め、x = a での波動関数の値を計算します。
固有値となるエネルギー E に対応する k が使われていれば、ψ(a) = 0 となります。k の値が解でない場合はψ(a) の値がゼロからずれます (図１左)。Runge-Kutta 法で x = −a から解き進めていったときの ψ(a) の値を用いた k の関数とし、二分法を用いることで ψ(a) = 0 となる解 k を求めます (図 1 右)。
この方法で、エネルギー E と境界条件を満たす波動関数 ψ(x) を求めることができますが、波動関数の初期値は ϕ(−a) = 1 として計算しており規格化されていません。規格化は Runge-Kutta 法で求まった波動関数を ψ ′ (x) とすると、

ψ(x) = √∫ψ ′ (x)aψ ′∗ (x′ )ψ ′ (x′ )dx′−a1(6)211.5
1ψ(a)

1.5
ψ(x)

2

k=0.0/ak=0.5/ak=1.0/ak=1.5/ak=2.0/a0.5

0.5
0-0.5-10-1.5-0.5-2-1-0.5

00.5

10x/a2468101214k

図 1 左:井戸型ポテンシャルの場合、様々な k を用いて x = −a から計算した波動関数。右:k の値を変えたときの ψ(a) の値。ψ(a) = 0 となる k が Schrödinger 方程式の解 E に対応する。

で与えられます。
対称性を使うと計算領域を半分に減らすこともできます。微分方程式 (2) および境界条件は x = 0 に対して対称ですので解は x = 0 での反転に対して正のものと負のものに分けることができます。パリティが正の解ψ(x) = ψ(−x) では x = 0 での一階微分がゼロとなるので ψ(0) = 1, ϕ(0) = 0 とできます。パリティが負の解 ψ(x) = −ψ(−x) 値がゼロになるので ψ(0) = 0, ϕ(0) = 1 とできます。x = 0 からこれらの初期条件を用いて Runge-Kutta 法で解き進め、x = a での波動関数の値がゼロとなるような k の値を決定します。

2 Numerov 法
1 次元 Schrödinger 方程式のように、微分方程式が波動関数の一階微分を含まない場合は Numerov 法という Runge-Kutta 法よりも精度のよい方法を用いることができます。x 周りので波動関数の Taylor 展開を考えます。(ψ (n) (x) ≡ (dn /dxn )ψ(x))

h2 (2)h3h4
ψ (x) + ψ (3) (x) + ψ (4) (x) + O(h5 ),2!
3!
4!
234hhhψ(x − h) =ψ(x) − hψ (1) (x) + ψ (2) (x) − ψ (3) (x) + ψ (4) (x) + O(h5 ),
2!
3!
4!

ψ(x + h) =ψ(x) + hψ (1) (x) +(7)(8)

これらを足し合わせると奇数次の項が消え

h4 (4)
ψ (x) + O(h6 )12(9)ψ(x + h) + ψ(x − h) − 2ψ(x) h2 (4)− ψ (x) + O(h4 )h2
12(10)

ψ(x + h) + ψ(x − h) = 2ψ(x) + h2 ψ (2) (x) +となります。Schrödinger 方程式の二階微分はこの式を用いると

ψ (2) (x) =とかけます。これは二階微分の3 点公式でした。これを代入する前に Schrödinger 方程式 (式 (2)) に)(1+h2 d2
12 dx2(1+h2 d212 dx2を作用させます。

)()
d2h2
h2 d2 22
ψ(x)+k[k (x)ψ(x)] = 0(x)ψ(x)= ψ (2) (x) + ψ (4) (x) + k 2 (x)ψ(x) +2dx
1212 dx2(11)2

この式の ψ (2) (x) に式 (10) を代入すると ψ (4) (x) の項が消えます。

ψ(x + h) + ψ(x − h) − 2ψ(x)h2 d2 22+k(x)ψ(x)+[k (x)ψ(x)] + O(h4 ) = 0h2
12 dx2(12)

最後の項に対しても二階微分の 3 点公式

d2 2k 2 (x + h)ψ(x + h) + k 2 (x − h)ψ(x − h) − 2k 2 (x)ψ(x)[k (x)ψ(x)] =+ O(h2 )2dxh2(13)

を用いると、][][][
5 2 2h2 2h2 2
1 + k (x + h) ψ(x + h) + −2 + h k (x) ψ(x) + 1 + k (x − h) ψ(x − h) + O(h6 ) = 0 (14)
12612
より

[]
[]
5h2
2 − h2 k 2 (x) ψ(x) − 1 + k 2 (x − h) ψ(x − h)612
ψ(x + h) =+ O(h6 )h2 2
1 + k (x + h)12(15)となり、ψ(x), ψ(x − h) と k(x), k(x − h), k(x + h) の値を用いて ψ(x + h) の値を計算することができます。

(井戸型ポテンシャルの場合は k は定数です) Runge-Kutta 法では誤差は O(h5 ) だったので、Runge-Kutta法よりも高い精度で計算することができます。

3 差分法
Schrödinger 方程式は ψ(x) に関して線形になっています。求めたい関数に対して方程式が線形の場合は差分法が有効です。
差分法では座標 x を離散化し、波動関数の値を離散化された格子点での値で表現します。求めたいのは

ψ(x) の −a < x < a での値なので、これらの領域を N + 1 等分した格子点 i = 0, 1, 2, · · · N, N + 1 を考えます。格子点上での x の値は

xi = −a +2ai
≡ −a + i∆xN +1(16)となります。(x0 = −a, xN +1 = a)。また格子点の間隔を ∆x = 2a/(N + 1) とおいています。ψi ≡ ψ(xi ) とおくと、微分方程式 (1) は i = 1, 2, · · · N の点で−h̄2 ψi+1 − 2ψi + ψi−1h̄2 d2+V(x)ψ(x)≃−ψ(x)+ V (xi )ψi = Eψiii2m dx22m(∆x)2x=xi(17)とかけます。ここで x = xi での二階微分は 3 点公式を使って xi−1 と xi+1 での値を用いて評価しています。
差分法では境界での振る舞いに注意する必要があります。i = 1 での式では ψ0 が出てきますが、この波動関数の値は井戸型ポテンシャルの場合は ψ0 = ψ(−a) = 0 です。また、i = N での ψN +1 = ψ(a) も同様にゼロです。以上の式をまとめると i = 1 のときは−h̄2 ψ2 − 2ψ1+ V (x1 )ψ1 = Eψ1 ,2m (∆x)23(18)2 ≤ i ≤ N − 1 のときは−h̄2 ψi+1 − 2ψi + ψi−1+ V (xi )ψi = Eψi2m(∆x)2(19)−h̄2 −2ψN + ψN −1+ V (xN )ψN = EψN ,2m(∆x)2(20)i = N のときはです。これらを行列の形でまとめると、t =h̄2を用いて2m(∆x)2
2t + V (x1 )−t0···
−t2t + V (x2 )−t0
0−t2t + V (x3 ) −t

.
..
..
..
.
.
0
.
..

···
0−t0······0






ψ1
ψ1
  ψ2  ψ2 ···  ψ3  ψ3   ..  . 0  .  = E  .. ..
 . 
 . .
0  ..  ..  ψN −1 ψN −1 2t + V (xN −1 )−tψN
ψN−t2t + V (xN )
(21)
···

0..
.
..
.となります。左辺の行列は Hamiltonian 演算子を座標表示で書いたものに対応し (Hamiltonian 演算子を

h̄2 d2+ V (x̂) とするとこの行列要素は Hij = ⟨xi |Ĥ|xj ⟩)、波動関数 ψi (= ⟨xi |ψ⟩) は固有ベクトル、2m dx2
エネルギー E は固有値に対応した固有値方程式

Ĥ = −Hψ = Eψ(22)になっています (ψ = (ψ1 , ψ2 , · · · , ψN )T )。井戸型ポテンシャルの場合は V (xi ) = 0 とおけるので、Hamil-tonian 行列はより簡単な形になります。
左辺の行列は三重対角の対称行列ですので LAPACK の三重対角行列、対称行列、一般行列の固有値問題のルーチン (DSTEV, DGESV, DGEEV など) を用いて固有値 E と対応する固有ベクトル ψi を求めることができます。ただし、波動関数の規格化は別途行う必要があります。ライブラリの出力結果の固有ベクトル ψi′は N 次元ベクトルとしての規格化N∑|ψi′ |2 = 1(23)i=1を満たしたものとなりますが量子力学での波動関数は∫ a−a|ψ(x)|2 dx =N∑|ψi |2 ∆x = 1(24)i=1で規格化しますので (積分を単純な和に置き換えています)、式 (23) と比較して

ψi = √とすれば波動関数の規格化条件 (24) を満たします。

41ψi′∆x(25)

Numerov 法と比較すると、差分法による解法では高い精度が得られません (演習問題)。これは二階微分を3 点公式で近似しているためで、9 点公式のような高次のものを使うことで、格子点の数 N を増やすことなく精度をあげることができます。Numerov 法は 1 次元の座標変数での方程式では威力を発揮しますが、差分法は 2 次元、3 次元の座標空間での Schrödinger 方程式も、格子点の次元を拡張することによって解くことができます。しかしながら、3 次元に拡張すると格子点の数は N から N 3 へと増えるため、行列の次元も増大し、対角化に時間がかかります。N = 100 のまま 3 次元にすると格子点の数は N 3 = 106 (100 万) となり、

Hamiltonian 行列の対角化に必要な計算時間は行列の次元の 3 乗に比例するため、行列の次元が N = 100 の1 次元問題のときと比べると、3 次元問題の N 3 = 106 の場合は N 6 = 1012 (1 兆) 倍の計算時間がかかり、容易ではないことがわかります。

4 演習問題
井戸型ポテンシャルでの Schrödinger 方程式を解き、エネルギーの低い順に 10 つの解 k を求めよ。ただし

a = 1.0、N = 100 とする。
(42) Runge-Kutta 法で解を求めよ。
(43) Numerov 法で解を求めよ(44) 差分法で求めよ。
(45) 調和振動子型ポテンシャル
V (x) =1 2kx
2(26)

での 1 次元 Schrödinger 方程式を差分法を用いて解け。十分に x が大きいところ (x = 5) で波動関数はゼロになると近似する境界条件を用いてよい。また、h̄ = 1, m = 1, k = 1 とせよ。

5