FFTを利用したTDSEの数値計算

FFTを利用した時間依存するシュレーディンガー方程式の数値計算 #

時間依存するシュレーディンガー方程式 (TDSE) を、

  • 空間に関する微分をFFTで計算
  • 時間に関する微分を適当な時間発展アルゴリズム（例えば適応刻み幅4次ルンゲクッタ法）

で解く方法を紹介します。計算に特別な工夫はありませんが、愚直で、ある程度一般的な計算が可能です。

1. 問題設定 #

次の時間依存するハミルトニアンを持つ1次元時間依存するシュレーディンガー方程式を解くことを考えます。

$$ \begin{align} \label{e1} i\hbar\frac{\partial}{\partial t}\psi=\hat{H}(t)\psi(t) \end{align} $$

ここで、$i$は虚数単位、$\hbar$はプランク定数であり、時間に依存する演算子$\hat{H}(t)$は次のように書けているとします。

$$ \begin{align}\label{e2} \hat{H} = \sum_{s=0} C_s(x,t)\frac{\partial^s}{\partial x^s} \end{align} $$

ここで、$\frac{\partial^0}{\partial x^0}\equiv 1$としています。

2. 計算方法 #

偏微分方程式を解くにあたって、式\eqref{e1}の右辺を計算することができれば、それは時刻$t$に対する1変数の微分方程式を解くことと同じになります。

つまり問題となるのは、関数$f(x)$が与えられた時に波動関数の空間に対する一階微分

$$ \begin{align}\label{e3} \frac{d}{dx} f(x) \end{align} $$

を計算する手立てがあれば良いわけです。高階微分は一階微分を繰り返し適用すれば得られるので、これが問題の本質となります。
また、本来の問題は偏微分ですが、時刻は今注目していないので通常の微分として書いておきます。

微分を計算する最も簡単な数値的な方法は、微分を適当な次数の有限差分に置き換えることでしょう。例えば前進差分、中心差分、後進差分等で置き換えて行列の形で書き表すことができます12。

ここでは微分を計算する方法としてフーリエ変換を介する方法を紹介します。

2.1. フーリエ変換を介した微分の計算 #

フーリエ変換後の空間では、微分演算子が単なる数値の掛け算になることを利用して関数の微分を計算します。

位置空間の計算区間$[0, L]$で周期境界条件が課されている場合を考えます。 この区間で任意の関数$f(x)$がある場合、周期境界条件を満たす離散的な波数$k=k_n$を持つ複素指数関数$e^{i2\pi k_n x}$の線形結合で表現することができます。 つまり、係数$f(k_n)$を用いて

$$ \begin{align}\label{e4} f(x)=\frac{1}{N}\sum_{n=0}^{N-1}f(k_n)e^{i2\pi k_n x} \end{align} $$

と書くことにします。ここで係数$\frac{1}{N}$は後の離散フーリエ変換を見越して導入しており、$k_n=n\Delta k=n/L, \hspace{0.2em}(n\text{: integer})$です。 $f(x)$のフーリエ変換は次のように書けます。

$$ \begin{align}\label{e4a} f(k_n)=\frac{N}{L}\int_{0}^{L} dx f(x)e^{-i2\pi k_n x} \end{align} $$

式\eqref{e4}の両辺にそれぞれ微分演算子を作用させると、

$$ \begin{align} \frac{d}{dx}f(x)&=\frac{1}{N}\sum_{n=0}^{N-1}f(k_n)\frac{d}{dx}e^{i2\pi k_n (x-x_a)} \\ &=\frac{1}{N}\sum_{n=0}^{N-1}(i2\pi k_n) f(k_n)e^{i2\pi k_n (x-x_a)}\label{e5} \end{align} $$

となります。式\eqref{e5}の$f(k_n)$は式\eqref{e4a}から得られます。

これを$\int \to \sum \Delta x$として離散化すれば

$$ \begin{align} f(k_n) &=\frac{N}{L}\int_{0}^{L} dx f(x)e^{-i2\pi k_n x} \\ &\rightarrow \sum_{m=0}^{N-1} f(x_m)e^{-i2\pi k_n x_m} \end{align} $$

となり、結局$f(k_n)$は$f(x)$の離散フーリエ変換で得られることが分かります。ここで、$x_m=m\Delta x=m\frac{L}{N}, \hspace{0.2em}(m\text{: integer})$として書いています。 より一般化して$s$階微分について考えれば

$$ \begin{align} \frac{d^s}{dx^s}f(x)&=\sum_{n=0}^{N-1}f(k_m)(i2\pi k_n)^s e^{i2\pi k_n (x-x_a)} \end{align} $$

と書けます。

2.2. まとめ #

これまでの話をまとめます。関数$f(x)$の$s$階微分$\frac{d^s f(x)}{dx^s}$を計算するためには、次の手順に従って計算することができます。

  1. $f(x_m)$を離散フーリエ変換し、$f(k_n)$を得る

    $$ \begin{align} f(k_n) = \sum_{m=0}^{N-1} f(x_m)e^{-i2\pi k_n x_m},\hspace{0.5em}(n=0,1,\cdots,N-1) \end{align} $$

  2. $f(k)$に$(i2\pi k_n)^s$を掛け合わせ、$\tilde{f}(k_n)$を計算する

    $$ \begin{align} \tilde{f}(k_n) = (i2\pi k_n)^s f(k_n),\hspace{0.5em}(n=0,1,\cdots,N-1) \end{align} $$

  3. $\tilde{f}(k_n)$を逆離散フーリエ変換する。この結果は$\frac{d^s f(x)}{dx^s}$となっている

    $$ \begin{align} \frac{d^s}{dx^s}f(x_m) = \frac{1}{N}\sum_{n=0}^{N-1}\tilde{f}(k_n)e^{i2\pi k_n x_m},\hspace{0.5em}(m=0,1,\cdots,N-1) \end{align} $$

フーリエ変換を介することで、微分演算子$d/dx$を作用させる代わりに$i2\pi k$を掛けるだけで計算することができます。 本方法は、離散フーリエ変換を使用しているため周期境界条件を満たす数値計算方法との相性が良いです。

特に、シュレーディンガー方程式を解く場合、周期境界条件を持つFourier Grid Hamiltonian法と非常に相性が良いです。

3. プログラム #

Pythonのプログラムを下記に置きます。

tdse_solver_v1.00.zip (31kB)

本プログラムは、1次元のシュレーディンガー方程式の数値解法を3段階で実施するシミュレーションです。

  1. TDSEの初期状態を決めるために、TISEを解いて固有関数を求めます。多くの場合でTDSEの初期状態として固有関数を用いますが、TISEを解かずに任意の関数を直接設定しても良いです。
    • もしTISEを解く場合、Fourier‐Grid Hamiltonian (FGH) 法を用いて時間依存しないシュレーディンガー方程式 (TISE) 数値的に解き、固有値および固有関数を求めます。
  2. scipyのsolve_ivpとFFTを利用して空間微分を計算することで、時間依存するシュレーディンガー方程式 (TDSE) を解いて波動関数の時間発展を求めます。
  3. 得られた波動関数を固有状態に射影する解析を行い、その結果や波動関数の実部・虚部の変化を時系列プロットおよびGIFアニメーションとして出力します。

TISEは下記の形を想定しています。

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

TDSEは下記の形を想定しています。

$$ \begin{align} i\frac{\partial}{\partial t}\psi(x,t)=\left[C_2(t,x,\psi)\frac{\partial^2}{\partial x^2} + C_1(t,x,\psi)\frac{\partial}{\partial x} + C_0(t,x,\psi)\right]\psi(x,t) \end{align} $$

デフォルトの計算パラメータは次のようになっています。

変数 デフォルト値 説明
$x$ $x=[-10,10]$ 位置の定義域
$N_x$ $2^6=64$ 位置の分割数
$t$ $t=[0,4\pi]$ 時間の計算区間
$N_t$ $33$ 時間の結果出力数（積分精度に影響しない）
rtol $10^{-6}$ solve_ivpの許容相対誤差
atol $10^{-6}$ solve_ivpの許容絶対誤差
method RK45 solve_ivpの手法
$V(x)$ $\frac{1}{2}x^2$ TISEのポテンシャル
$\psi(x,t=0)$ $\frac{1}{\sqrt{2}}[\phi_0(x)+\phi_1(x)]$ TDSEの初期状態 ($\phi_n(x)$はTISEの$n$番目の固有値に属する固有状態)
$C_2(t,x,\psi)$ $-\frac{1}{2}$ TDSEの係数
$C_1(t,x,\psi)$ $0$ TDSEの係数
$C_0(t,x,\psi)$ $\frac{1}{2}x^2$ TDSEの係数

実際に実行すると次の結果を得ます。

PS C:\work\python\tdse\tdse_solver_v1.00> python .\main.py
  Elapsed time TISE: 0.002004 seconds
  Elapsed time TDSE: 0.290015 seconds
  Elapsed time  GIF: 5.557767 seconds
PS C:\work\python\tdse\tdse_solver_v1.00>

これくらいの分点数と時間の計算区間であれば、0.3秒程度で計算できていることが分かります。

また、出力は次の結果が出力されます。

出力結果 出力されるファイル名 説明
wavefunction.gif 時間発展の様子
projection.png TISEの固有状態への射影の絶対値 $ | \langle \psi(x,t) | \phi_n\rangle | $
(ただし左図では0,1が重なっています)
eigenfunctions.png TISEの固有状態
テキスト tise_eigenvalues.csv TISEの固有値 (テキスト)

tise_eigenvalues.csvは下記のように書かれます。

0,0.5000000000000019
1,1.4999999999999434
2,2.49999999999999
3,3.500000000000009
4,4.500000000000011
5,5.500000000000007
...

4. 参考文献 #

  1. §3.4 Differentiation -DLMF ↩︎

  2. Finite difference -wikipedia ↩︎