シュレディンガー方程式の数値解法

前方差分法（陽解法）

　シュレディンガー方程式

\[ \small i\hbar\frac{\partial \psi(x,t)}{\partial t} = -\frac{\hbar^2}{2m}\frac{\partial^2 \psi(x,t)}{\partial x^2} +V(x)\psi(x,t) \]

について、数値計算で波動関数の時間発展を求める方法を考えよう。いちいち方程式を解くのが面倒くさいし、解析解で解けない場合も少なくないからである。関数\(\small \psi(x,t)\)について、時間と空間を等間隔に分割し

\[ \small \begin{align*} &x_j = 0,\Delta x, 2\Delta x,\cdots, (M-1)\Delta x \\ &t_k = 0,\Delta t, 2\Delta t,\cdots, (N-1)\Delta t \end{align*} \]

と表す。\(\small \psi_{j,k} = \psi(j\Delta x, k\Delta t)\)と表すことにする。

　常微分方程式と同様に偏微分を差分近似して時間発展を計算する方法を考える。実際に差分近似を計算すると

\[ \small \begin{align*} &\frac{\partial \psi_{j,k}}{\partial t} = \frac{\psi_{j,k+1}-\psi_{j,k}}{\Delta t} \\ &\frac{\partial^2 \psi_{j,k}}{\partial x^2} = \frac{\psi_{j+1,k}+\psi_{j-1,k}-2\psi_{j,k}}{\Delta x^2} \end{align*} \]

と置くことができる。これを用いて偏微分方程式を差分方程式の形に書き直そう。

\[ \small i\hbar\frac{\psi_{j,k+1}-\psi_{j,k}}{\Delta t} = -\frac{\hbar^2}{2m}\frac{\psi_{j+1,k}+\psi_{j-1,k}-2\psi_{j,k}}{\Delta x^2} +V(x_{j})\psi_{j,k} \]

であるから、式を整理すると

\[ \small \psi_{j,k+1} = \psi_{j,k}+i\Delta t\left(\frac{\hbar}{2m}\frac{\psi_{j+1,k}+\psi_{j-1,k}-2\psi_{j,k}}{\Delta x^2}-\frac{V(x_{j})}{\hbar}\psi_{j,k}\right) \]

と表すことができる。手前の時点から順番に計算していけば、偏微分方程式の解を求めることができる。これは常微分方程式におけるオイラー法に相当する手法と考えることができるだろう。ルンゲ・クッタ法のように多段階で計算することも可能であると考えられる。偏微分方程式の数値解法では、このような計算方法は陽解法（Explicit Method）といわれる。次節で説明するクランク-ニコルソン法と比較して収束性が悪いため、現実にはあまり使われないが、非常に複雑な偏微分方程式についても容易に計算することが可能であり、手っ取り早く解のイメージをつかむために利用しやすい手法である。

クランク-ニコルソン法

　前節では、時間について前方差分で偏微分方程式を解いたが、これは時間についてバイアスが生じることは容易に推測できるだろう。もちろん、グリッドを細かくしていけば、そういったバイアスは減少していくものの、バイアスがかからない形式で定式化できるのであれば、そちらを利用した方が良いだろう。

　前方差分の近似と、後方差分の近似の２つを考えて平均をとることを考える。時間の変化幅を\(\small \Delta t/2\)として、手前のグリッドから前方差分で近似した式と後ろのグリッドから後方差分で近似した式の2つを考える。具体的には

\[ \small \begin{align*} &i\hbar\frac{\psi_{j,k+1/2}-\psi_{j,k}}{\Delta t/2} = -\frac{\hbar^2}{2m}\left(\frac{\psi_{j+1,k}+\psi_{j-1,k}-2\psi_{j,k}}{\Delta x^2} \right) +V(x_{j})\psi_{j,k} \\ &i\hbar\frac{\psi_{j,k+1}-\psi_{j,k+1/2}}{\Delta t/2} = -\frac{\hbar^2}{2m}\left(\frac{\psi_{j+1,k+1}+\psi_{j-1,k+1}-2\psi_{j,k+1}}{\Delta x^2} \right) +V(x_{j})\psi_{j,k+1} \ \end{align*} \]

となる。二つの式を足すと

\[ \small \begin{align*} &i\hbar\frac{\psi_{j,k+1}-\psi_{j,k}}{\Delta t} = -\frac{\hbar^2}{4m}\Biggl(\frac{\psi_{j+1,k+1}+\psi_{j-1,k+1}-2\psi_{j,k+1}}{\Delta x^2}+\frac{\psi_{j+1,k}+\psi_{j-1,k}-2\psi_{j,k}}{\Delta x^2} \Biggl) \\ &\qquad\qquad\qquad\qquad+V(x_{j})\frac{\psi_{j,k}+\psi_{j,k+1}}{2} \end{align*} \]

を得る。このように座標に関する微分を計算する方法はクランク-ニコルソン法といわれる。後方差分を用いる解法は陽解法との対照から陰解法（Implicit Method）といわれ、クランク-ニコルソン法は陰解法の一種である。式を整理すると

\[ \small \begin{align*} &a_j\psi_{j-1,k+1}+b_j\psi_{j,k+1}+c_j\psi_{j+1,k+1} = -a_j\psi_{j-1,k}-d_j\psi_{j,k}-c_j\psi_{j+1,k}\\ &a_j = -\frac{\hbar\Delta t}{4m\Delta x^2}i\\ &b_j = 1+\frac{\hbar\Delta t}{2m\Delta x^2}i+i\frac{V(x_j)}{2\hbar}\Delta t \\ &c_j = -\frac{\hbar\Delta t}{4m\Delta x^2}i\\ &d_j = -1+\frac{\hbar\Delta t}{2m\Delta x^2}i+i\frac{V(x_j)}{2\hbar}\Delta t \end{align*} \]

と表すことができる。行列の形式で表せば、以下のとおりである。

\[ \small \left[ \begin{array}{cccccc} b_{1} & c_{1} & 0 & 0 & \cdots & 0 & 0 \\ a_{2} & b_{2} & c_{2} & 0 & \cdots & 0 & 0 \\ 0 & a_{3} & b_{3} & c_{3} & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \ddots & b_{M-3} & c_{M-3} \\ 0 & 0 & 0 & 0 & \cdots & a_{M-2} & b_{M-2} \end{array} \right] \left[ \begin{array}{c} \psi_{1,k+1} \\ \psi_{2,k+1} \\ \psi_{3,k+1} \\ \vdots \\ \psi_{M-3,k+1} \\ \psi_{M-2,k+1} \end{array} \right]=\qquad\qquad\qquad\qquad \\
\small \left[ \begin{array}{l} -a_{1} \psi_{0,k} – d_{1} \psi_{1,k} – c_{1} \psi_{2,k} – a_{1} \psi_{0, k+1} \\ -a_{2} \psi_{1,k} – d_{2} \psi_{2,k} – c_{2} \psi_{3,k} \\ -a_{3} \psi_{2,k} – d_{3} \psi_{3,k} – c_{3} \psi_{4,k} \\ \qquad\qquad\qquad\vdots \\ -a_{M-3} \psi_{M-4,k} – d_{M-3} \psi_{M-3,k} – c_{M-3} \psi_{M-2,k} \\ -a_{M-2} \psi_{M-3,k} – d_{M-2} \psi_{M-2,k} – c_{M-2} \psi_{M-1,k} – c_{M-2} \psi_{M-1, k+1} \end{array} \right] \]

境界条件として、\(\small \psi_{j,0},\psi_{0,k},\psi_{M-1,k},\forall j,k\)が与えられれば、時間グリッドの手前の方からこの連立方程式を繰り返し解くことで波動関数の時間発展を計算することができる。以上が有限差分法と言われる偏微分方程式の数値解法である。シュレディンガー方程式を数値解法で解く場合は境界条件や波動関数のグリッド上の値は複素数として計算しなければならないことに注意する。

自由状態の場合

　実際の例として、自由状態のシュレディンガー方程式について有限差分法で解いてみたものを掲載しておく。こちらは解析式でもExactな解なので、まあ有限差分法で解いてもイメージ通りであろうと思われる。お試しで動画の形式にしてみたので、イメージがつかみやすくなっているのではないだろうか。

１次元のクーロンポテンシャル

　前回の記事

で怪しげな近似解を提出した１次元クーロンポテンシャルのシュレディンガー方程式についても、有限差分法で解いてみた。波動関数が\(\small x>0\)の領域と\(\small x<0\)の領域が混じるため、干渉縞が生じてカオスのような運動になった。両極を行ったり来たりする運動をするという意味では、ある程度予想は合っていたと言いたいけど、それにしても想定外な結果である。時間発展を動画にしてみたが、なんかおバカなネタっぽくなってしまった・・・