メモ：繰りこみヌメロフ法による 1 次元シュレーディンガー方程式の数値計算法
堀越宗一

2021 年 2 月 13 日

1 問題設定次の 1 次元シュレーディンガー方程式を考える。

[]
h̄2 d2−+ V (x) Ψ(x) = EΨ(x)2m dx2(1)

ここで、2m[E − V (x)]h̄2(2)]
d2+Q(x)Ψ(x) = 0dx2(3)

Q(x) =と定義すると、(1) 式は以下のように書ける。

[問題設定は、与えれた任意の Q(x) に対して (3) 式を解き、固有値 E と固有状態 Ψ(x) を得ることである。これに対する解法の一つが「繰りこみヌメロフ法」と呼ばれるものである。この手法は学生さんでも簡単に使えるのでこのメモで簡単に紹介します。

2 ヌメロフ法 (Numerov method)x 座標軸上で xi から xf まで h の等間隔で (N+1) 点に離散化されているとする。離散化された座標を
xn ≡ xi + nh, (n = 0, 1, . . . , N ) とし、離散化された関数をそれぞれ、Ψn ≡ Ψ(xn ), Qn ≡ Q(xn ) とする。
Ψn+1 と Ψn−1 は x = xn まわりでのテーラー展開で、
Ψn+1 =∞
∑hk
k=0k!

Ψ(k)
n , Ψn−1 =∞∑(−h)kk=0k!

Ψ(k)nと記述できる [1]。これらを用いると、

11
1 4 (4)1[Ψn+1 + Ψn−1 ] = Ψn + h2 Ψ(2)h Ψn + h6 Ψ(6)n +n + ...
224!
6!
この 2 階微分は、

1 4 (6)1 (2)1 2 (4)(2)[Ψn+1 + Ψn−1 ] = Ψ(2)h Ψn + . . .
n + h Ψn +224!
1(4)(k)となる。ここで Ψn は Ψ(x) の x = xn における k 階の微分である。これらと (3) 式を用いることにより、

Ψ(x) が従うべき微分方程式は以下の形式で書ける。
1 h2
5h6 (6)1(Ψn+1 + Ψn−1 ) −(−Qn+1 Ψn+1 − Qn−1 Ψn−1 ) = Ψn − h2 Qn Ψn −
Ψ + ...
22 1212480 n(5)ここで、Tn = −h2Qn12(6)と定義し、微分されていない波動関数を左辺に、6 階以上の微分項を右辺に移動すると以下のようになる。

(1 − Tn+1 )Ψn+1 − (2 + 10Tn )Ψn + (1 − Tn−1 )Ψn−1 = −h6 (6)
Ψ + ...
240 n(7)

ここで波動関数の 6 階微分以降をゼロと近似することにより、(3) 式の微分方程式の Ψ(x) が従うべき方程式は以下の近似式で与えられる。

(1 − Tn+1 )Ψn+1 − (2 + 10Tn )Ψn + (1 − Tn−1 )Ψn−1 ≈ 0(8)

この近似式を境界条件を満たすように解くのがヌメロフ法である。

3 繰りこみヌメロフ法 (Renormalized Numerov method)
参考文献 [2] に従って繰りこみヌメロフ法を説明する。

Fn = (1 − Tn )Ψn(9)と定義することにより、(8) 式は以下のようにすっきりとした形で書ける。

Fn+1 − Un Fn + Fn−1 = 0
ここで Un は、Un =2 + 10Tn1 − Tn

(10)(11)である。つまり (10) 式を解いて Fn さえ決まってしまえば、(9) 式によって波動関数 Ψn が直ちに求まる。

(10) 式を x0 から外側に積分していくため、x の大きい方向に一つずれた位置との比を次のように与える。
Ln =すると (10) 式は 2 点問題となり、

Fn+1
Fn

Ln = Un − L−1n−1(12)(13)となる。境界条件が Ψ0 = 0, Ψ1 ̸= 0 の場合、Ln の初期値は L0 = ∞ である。この初期値により、(13) 式に従って x0 から外側に Ln が計算できる。Ln は左側から計算しているため”’Left’ の L の文字を用いている。
同様に、(10) 式を xN から内側に積分していくため、x の小さい方向に一つずれた位置との比を次のように与える。

Rn =これを用いると (10) 式はFn−1
Fn−1
Rn = Un − Rn+12(14)(15)となる。境界条件が ΨN = 0, ΨN −1 ̸= 0 の場合、Rn の初期値は RN = ∞ である。この初期値により、(15)式に従って xN から内側に Rn が計算できる。Rn は右側から計算しているため”’Right’ の R の文字を用いている。

Rn の振る舞いを見ると、ある点で Rn ≤ 1 となる点が現れる。この点が波動関数の微分値が初めてゼロになる場所に相当し、この点を xm とする。波動関数の連続性の要請より、Ln と Rn は xm で滑らかに接続される必要がある。xm における両者のズレの評価は以下の評価関数によって行う [2]。

[]−1
D(E) = Am+1 (Rm+1− Lm ) − Am−1 (Rm − L−1m−1 ) (1 − Tm )(16)An = (1/2 − Tn )(1 − Tn )−1(17)ここで、である。D(E) = 0 となる試行エネルギー E が固有値に相当し、波動関数が xm において滑らかに接続される。
固有値で Ln と Rn を計算し、Fm = 1(18)とすることにより、(12), (14) 式に従い Fn が x0 ≤ x ≤ xN の領域で得られる。（ここでは波動関数の規格化∫は気にしない）これを (9) 式に代入し、 Ψ(x)2 dx = 1 となるように振幅を調整すれば規格化された波動関数

Ψ(x) が得られる。
与えたポテンシャル中に複数の束縛状態が存在する場合、(16) 式がゼロとなるエネルギーは束縛状態の数だけ存在する。ある固有値における波動関数の節の数が α 個の場合、その状態は基底状態から α 番目の状態に相当する。

4 等間隔でない離散化の手法
分子ポテンシャルを用いて解離極限近傍の固有値を計算する際、波動関数は近距離で激しく振動し遠距離で緩やかに振動する。そのため原子間距離に依存して離散化の間隔を調整すると、ヌメロフ法の計算精度の向上と計算コストの軽減が実現できる。ここでは参考文献 [3] の (5), (6) 式を用いた方法をまとめておく。
ここでは V (R) のポテンシャルの固有値を計算することを想定し、Rmin ≤ R ≤ Rmax の範囲で n 点離散データを取るとする。ここではデータ点の取り方を、

R = U (y) = bea(y−1)(19)とし、y = 1, 2, . . . , n に n 点サンプリングするとする。

a=1
lnn−1

(RmaxRmin), b = Rminとすると、

R(y = 1) = Rmin , R(y = n) = Rmax ,となり、R = Rmin 近傍で密に、R = Raxn 近傍で疎にサンプリングすることができる。U (R) の逆関数をu(y) とし、
1y = u(R) = lna3(Rb)+1(20)である。
等間隔な y 座標系でヌメロフ法を適応し、R 座標系の固有値と固有関数を求める。参考文献 [3] の (6) 式にU (y) = bea(y−1) を代入すると、(1) 式の 1 次元シュレーディンガー方程式は y 座標系で以下の様に書き換えることができる。

[ここで、Q(y) ≡]d2+ Q(y) Φ(y) = 0dy 2)
a22µ (a(y−1) 2E−Ṽ(y)(abe)−,4h̄2
Φ(y) ≡ϕ(y)ea(y−1)(21)(22)(23)

で定義しており、Ṽ (y) = V (U (y)) である。(21) 式の y 座標系は等間隔にサンプリングされているので従来のヌメロフ法を適応することができ、(3) 式と同様に固有値 E を計算することができる。
波動関数は、参考文献 [3] の (5) 式より、

Ψ(R) = √√ϕ(y)abea(y−1)=ea(y−1)Φ(y)ab(24)となり、規格化条件は (24) 式より以下の通りである。

∫|Ψ(R)|2dR
dy =dy∫(Φ(y)ea(y−1) )2 dy = 1(25)5 まとめ
以上の手順より、与えられたポテンシャル V (x) に対する固有値と固有状態を計算することができる。この手法はプログラムが単純なだけでなく、安定で精度が高い。

謝辞
ヌメロフ法と離散化手法を教えてくださった北海道大学の小林淳氏に感謝申し上げます。

[参考文献]
1. John M. Blatt ”Practical points concerning the solution of the Schrödinger equation.” Journal ofComputational Physics 1.3 (1967): 382-396.
2. Bernard R. Johnson ”New numerical methods applied to solving the one ‐ dimensional eigenvalueproblem.”, The Journal of Chemical Physics 67.9 (1977): 4086-4093.
3. Eite Tiesinga, Carl J. Williams, and Paul S. Julienne ”Photoassociative spectroscopy of highly
excited vibrational levels of alkali-metal dimers: Green-function approach for eigenvalue solvers.”
Phys. Rev. A 57, 4257 (1998).

4