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)

파이썬에서 푸리에 변환하는 코드는 다음과 같습니다.

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'$로 이산화된 파수에서의 푸리에 변환을 계산할 수 있다.