有限差分法

有限差分法による数値解法 #

有限差分法は、微分方程式を差分方程式で近似して解く方法です。

Webツールを下記のリンク先に実装しました。

差分法で解く、1次元時間依存しないシュレーディンガー方程式

1. シュレーディンガー方程式 #

原子単位系 ($m=1, \hbar=1$) に置ける、量子系の波動関数に関するシュレーディンガー方程式\eqref{eq1}を考えます。

$$ \begin{align}\label{eq1} \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x) \end{align} $$

ここで、$V(x)$はポテンシャル、$E$はエネルギー固有値を表します。また、$x$は区間$[x_a, x_b]$で定義され、境界条件として固定端条件を考えます。つまり、

$$ \begin{align} \psi(x_a)=\psi(x_b)=0 \end{align} $$

が波動関数に課されている場合を考えます。

1.1. 有限差分による近似 #

式\eqref{eq1}を解くために、波動関数を離散的に扱います。まず、座標$x$を下記のように離散的に決めます1

$$ \begin{gather} x_n = n\Delta x + x_a,~~~(n=0,1,\cdots,N) \\ \Delta x \equiv \frac{x_b-x_a}{N} \end{gather} $$

また上記の$n=0,N$における両端では、境界条件のため $\psi(x_0)=\psi(x_{N+1})=0$ が満たされています。

これらの離散点を用いると、二階微分は差分として下記のように近似されます。

$$ \begin{align} \frac{d^2f}{dx^2}\approx\frac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{\Delta^2} + O(\Delta^2) \end{align} $$

これを式\eqref{eq1}に代入すると、

$$ \begin{align} -\frac{1}{2\Delta x^2}\psi(x_{n+1})+\left(\frac{1}{\Delta x^2}+V(x_n)\right)\psi(x_n)-\frac{1}{2\Delta x^2}\psi(x_{n-1})=E\psi(x_n) \end{align} $$

を得ます。それぞれの離散点$n=0,1,\cdots,N$における波動関数について上記が成り立つため、$\psi_n\equiv\psi(x_n)$と書けば、

$$ \begin{gather} \psi_0 = \psi_{N} = 0\\ -\frac{1}{2\Delta x^2}\psi_{n+1}+\left[\frac{1}{\Delta x^2}+V(x_n)\right]\psi_n-\frac{1}{2\Delta x^2}\psi_{n-1}=E\psi_n,~~~(n=1,2,\cdots,N-1) \label{eq2} \end{gather} $$

となります。式\eqref{eq2}を行列の形で書けば、

$$ \begin{align} \begin{pmatrix} d_1 & e_2 & 0 & \cdots & 0 \\ e_2 & d_2 & e_3 & \cdots & 0 \\ 0 & e_3 & d_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & e_{N-1} \\ 0 & 0 & 0 & e_{N-1} & d_{N-1} \end{pmatrix} \begin{pmatrix} \psi_1 \\ \psi_2 \\ \psi_3 \\ \vdots \\ \psi_{N-1} \end{pmatrix} =E \begin{pmatrix} \psi_1 \\ \psi_2 \\ \psi_3 \\ \vdots \\ \psi_{N-1} \end{pmatrix}\label{eq3} \end{align} $$

となります。ここで

$$ \begin{align} d_i &\equiv \frac{1}{\Delta x^2}+V(x_i), &&(i=1,2,\cdots, N-1) \\ e_i &\equiv -\frac{1}{2\Delta x^2}, &&(i=2,3,\cdots, N-1) \end{align} $$

と置きました。あとはこれを対角化すればエネルギーと固有関数を得ます。

1.2. 固有関数として解釈 #

差分を固有関数による展開として考えても、式\eqref{eq3}に辿り着きます。それを示しましょう。

波動関数を固有関数の線形結合として表現します。つまり、波動関数を未知の係数$c_i$を用いて

$$ \begin{align}\label{eq4} \psi(x)=\sum_i c_i \varphi_i(x) \end{align} $$

として展開します。また、固有関数$\varphi_i(x)$を$i=0,1,\cdots,N$について下記の通り定義します。

$$ \begin{eqnarray} \varphi_i(x)= \left\{ \begin{aligned} & \frac{1}{\sqrt{\Delta x}} &&(x_i-\Delta x/2 \lt x \lt x_i+\Delta x/2) \\ & \frac{1}{2\sqrt{\Delta x}} &&(x = x_i\pm \Delta x/2) \\ & 0 &&(\text{otherwise}) \end{aligned} \right. \end{eqnarray} $$

また、$\varphi_i(x)$を含む積分は下記の通りに計算できます。

$$ \begin{align} \int_{x_a}^{x_b} f(x)\varphi_i(x)dx = \frac{1}{\sqrt{\Delta x}}\int_{x_i-\Delta x/2}^{x_i+\Delta x/2} f(x)dx \end{align} $$

ここで$f(x)$は任意の関数であることを仮定しています。

1.2.1. 内積 #

行列要素を計算するために、任意の関数$f(x)$に対する内積を計算しておきます。すると

$$ \begin{align} \int_{x_a}^{x_b} \varphi_i(x)f(x)\varphi_j(x)dx &=\frac{1}{\sqrt{\Delta x}}\int_{x_i-\Delta x/2}^{x_i+\Delta x/2} f(x)\varphi_j(x)dx \nonumber\\ &=\frac{1}{\Delta x}\int_{x_i-\Delta x/2}^{x_i+\Delta x/2} f(x)dx \delta_{i,j}\nonumber\\ &\approx \frac{1}{\Delta x}\int_{x_i-\Delta x/2}^{x_i+\Delta x/2} f(x_i)dx \delta_{i,j}\nonumber\\ &= f_i\delta_{i,j} \label{eq5} \end{align} $$

となります。特に$f(x)=1$について考えると、

$$ \begin{align} \int_{x_a}^{x_b} \varphi_i(x)\varphi_j(x)dx = \delta_{i,j} \end{align} $$

となるため、基底関数$\varphi_i(x)$は正規直交基底であることが分かります。

微分演算子に関する内積は下記の通りになります。

$$ \begin{align} \int_{-\infty}^{\infty} \varphi_i(x)\frac{d}{dx}\varphi_j(x)dx &=\frac{1}{\sqrt{\Delta x}}\int_{x_i-\frac{\Delta x}{2}}^{x_i+\frac{\Delta x}{2}} \frac{d}{dx}\varphi_j(x)dx \nonumber\\ &=\frac{1}{\sqrt{\Delta x}}\int_{\varphi(x_i-\frac{\Delta x}{2})}^{\varphi(x_i+\frac{\Delta x}{2})} d\varphi_j(x) \nonumber\\ &=\frac{1}{\sqrt{\Delta x}}\left[\varphi_j(x_i+\frac{\Delta x}{2})-\varphi_j(x_i-\frac{\Delta x}{2})\right] \nonumber\\ &=\frac{1}{\sqrt{\Delta x}}\left[\frac{1}{2\sqrt{\Delta x}}\left(\delta_{j,i+1}+\delta_{j,i}-\delta_{j,i-1}-\delta_{j,i}\right)\right] \nonumber\\ &=\frac{1}{2\Delta x}\left(\delta_{j,i+1}-\delta_{j,i-1}\right) \label{eq6} \end{align} $$

二階の微分演算子に関する内積は下記の通りになります。

$$ \begin{align} \int_{-\infty}^{\infty} \varphi_i(x)\frac{d^2}{dx^2}\varphi_j(x)dx &=\frac{1}{\sqrt{\Delta x}}\biggl[\frac{d\varphi_j}{dx}\biggr|_{x=x_i+\Delta x/2}-\frac{d\varphi_j}{dx}\biggr|_{x=x_i-\Delta x/2}\biggr] \nonumber\\ &=\frac{1}{\Delta x}\Bigl[\delta(x_i-x_j+\Delta x)-2\delta(x_i-x_j)+\delta(x_i-x_j-\Delta x)\Bigr] \nonumber\\ &\to \frac{1}{\Delta x^2}\left(\delta_{j,i+1}-2\delta_{j,i}+\delta_{j,i-1}\right)\label{eq7} \end{align} $$

ここで、下記の結果を用いました。

$$ \begin{gather} \frac{d\varphi_j}{dx}=\frac{1}{\sqrt{\Delta x}}\Bigl[\delta\left(x-(x_j-\Delta x/2)\right)-\delta(x-(x_j+\Delta x/2))\Bigr] \\ \delta(x-x’)\xrightarrow{離散化}\frac{\delta_{i,j}}{\Delta x} \end{gather} $$

1.2.2. ハミルトニアンの行列要素 #

以上の計算結果を用いて、式\eqref{eq4}の係数$c_n$について変分原理に従ってハミルトニアン$\hat{H}=-\frac{1}{2}\frac{d^2}{dx^2}+V(x)$の行列要素を考えます。 式\eqref{eq5}, \eqref{eq7}を用いれば、

$$ \begin{align} H_{ij}\equiv\langle\varphi_i|\hat{H}|\varphi_j\rangle &= -\frac{1}{2}\langle\varphi_i|\frac{d^2}{dx^2}|\varphi_j\rangle + \langle\varphi_i|V(x)|\varphi_j\rangle \nonumber \\ &= -\frac{1}{2\Delta x^2}\left(\delta_{j,i+1}-2\delta_{j,i}+\delta_{j,i-1}\right) + V(x_i)\delta_{i,j} \end{align} $$

を得ます。よって、

$$ \begin{align} \begin{pmatrix} d_1 & e_2 & 0 & \cdots & 0 \\ e_2 & d_2 & e_3 & \cdots & 0 \\ 0 & e_3 & d_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & e_{N-1} \\ 0 & 0 & 0 & e_{N-1} & d_{N-1} \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_{N-1} \end{pmatrix} =E \begin{pmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_{N-1} \end{pmatrix}\label{eq8} \end{align} $$

を得ます。これは、式\eqref{eq3}と同じになります。

2. 数値計算 #

Javascriptとhtmlで差分法で解く方法を、下記のページに実装作成しました。ユーザーがポテンシャルを指定して計算を実行できます。

差分法で解く、1次元時間依存しないシュレーディンガー方程式

ページの様子は下記のような様子です。

tise webtool image

Javascriptがそんなに早くないか実装の仕方が悪いため、行列の要素数が600x600程度まで行くと計算時間が遅いと感じ始めると思います。


  1. 等間隔である必要はありませんが、簡単のために等間隔に決めます。 ↩︎