Python (numpy)を用いて1次元ポテンシャル関数のフーリエ変換を行う際にスケーリングがわからなくなったのでメモ。
結論
はじめに結論だけまとめておく。
以下フーリエ変換の定義は
\tilde V(k) = \int_{-\infty}^\infty V(x)e^{-ikx} dx
である。
$\Delta x$を離散化幅として、ポテンシャル関数$V(x)$が$x_m=m\Delta x\ (m=0,1,\ldots,N-1)$での値$V(x_m)$として得られているとする。
この$V(x_m)$に対してPythonのnumpy.fft.fft関数で離散フーリエ変換を行う。その結果が$\tilde P(k_j)\ (k_j=j\Delta k)$として得られるとき、スケール補正した正しい値は
\tilde V(k_j) = \Delta x \tilde P(\frac{\Delta x}{2\pi}k_j)
で得られる。
コード
過去に同様の問題を取り扱った記事を参考に、ガウス関数のフーリエ変換を例に取り上げる。ポテンシャルとそのフーリエ変換は以下のように表される。
V(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{x^2}{2\sigma^2}\right)\\
\tilde V(k)=\exp\left(-\frac{\sigma^2k^2}{2}\right)
Pythonでフーリエ変換するコードは以下。
import numpy as np
import matplotlib.pyplot as plt
N = 4096
xlen = 30.0
dlt = 2.0*xlen/N
xr = np.linspace(-xlen,xlen,N)
sig = 0.1
pot = np.exp(-xr**2/(2.0*sig**2))/(np.sqrt(2.0*np.pi)*sig)
potk = np.fft.fft(pot)
k = np.fft.fftfreq(N)
ks = 2.0*np.pi/dlt*np.roll(k,int(N/2))
theory = np.exp(-0.5*ks**2*sig**2)
fig, ax = plt.subplots(figsize=(12,8))
ax.plot(ks,dlt*abs(np.roll(potk,int(N/2))),label="Python",color="red")
ax.plot(ks,theory,label="Theory",ls="dashed",color="blue")
ax.set_xlabel("k")
ax.set_ylabel(r"$\tilde{V}(k)$")
ax.legend()
plt.show()
上のコードでksが適切にスケーリングされた波数である。また、Pythonの結果をプロットする際にdlt（$\Delta x$）をかけていることに注目してほしい。
結果として、フーリエ変換は解析解とよく一致することが見て取れる。
矩形ポテンシャルについても同様の計算を行った。
V(x)=\begin{cases}1 & (-d\le x\le d) \\
0 & ({\rm otherwise})
\end{cases}\\
\tilde V(k)=2d\frac{\sin{kd}}{kd}=2d\ {\rm sinc}(kd)
pot = np.zeros(N)
d = 2.0
dst = int((xlen-d)/(2.0*xlen)*N)
dln = int((xlen+d)/(2.0*xlen)*N)
for i in range(dst,dln+1):
pot[i] = 1.0
potk = np.fft.fft(pot)
k = np.fft.fftfreq(N)
ks = 2.0*np.pi/dlt*np.roll(k,int(N/2))
theory = 2.0*d*np.sinc(d*ks/np.pi)
numpyでのsinc関数の定義に注意！（$\pi$で割るのを忘れて結果が合わず、しばらく悩んだ）
理論的な補足
フーリエ積分を数値積分の形式に書き換えると以下のようになる。
\begin{align}
\tilde V(k_j) &= \int_{-\infty}^\infty V(x)e^{-ik_jx} dx\\
&\approx \sum_{m=0}^{N-1} V(x_m)e^{-ik_jx_m} \Delta x\\
&=\Delta x \sum_{m=0}^{N-1} V(x_m) e^{-i \frac{jm}{N}}
\end{align}
ここで、$x_m=m\Delta x,\ k_j=j\Delta k$を用いた。$\Delta k$はk空間での離散化幅であり、$\Delta k=1/(N\Delta x)$で表される。
注意しなければならないのは、numpyのfft関数で計算される$\tilde P(k_j)$は$\boldsymbol\Delta \boldsymbol {x=1}$のもとでの
\tilde P(k_j)=\sum_{m=0}^{N-1} V(x_m) e^{-2\pi i \frac{jm}{N}}
であるということである(exp内の$2\pi$のファクターにも注意)。よって、ポテンシャル関数の表現に用いている実際の離散化幅$\Delta x'$に対応したスケーリングを施す必要がある。
そこで、波数$2\pi k_j'\ (k_j'=j\Delta k')$でのフーリエ積分を離散化幅$\Delta x'$の数値積分の形に直す。
\begin{align}
\tilde V(2\pi k_j') &= \int_{-\infty}^\infty V(x)e^{-2\pi ik_j'x} dx\\
&\approx \sum_{m=0}^{N-1} V(x_m')e^{-2\pi ik_j'x_m'} \Delta x'\\
&=\Delta x' \sum_{m=0}^{N-1} V(x_m') e^{-2\pi i \frac{jm}{N}}\\
&= \Delta x' \tilde P(k_j)
\end{align}
ここで、$V(x_m)=V(x_m')$から最後の項で$\tilde P(k_j)$と書くことができる。$k_j'=j/(N\Delta x')$より
\frac{k_j}{k_j'}=\frac{\Delta x'}{\Delta x}
が成立するので、
\begin{align}
\tilde V(2\pi k_j') &= \Delta x' \tilde P(\frac{\Delta x'}{\Delta x}k_j')\\
&= \Delta x' \tilde P(\Delta x' k_j')\qquad (\because \Delta x=1)
\end{align}
の関係が成立することがわかる。したがって、$2\pi k_j'\to k_j'$と置き換えることで
\tilde V(k_j')=\Delta x' \tilde P(\frac{\Delta x'}{2\pi}k_j')
が得られ、$\Delta k'$で離散化された波数でのフーリエ変換を計算することができる。