这个傅里叶变换有什么问题? (在 python 中)
What is wrong with this Fourier transform? (in python)
读过我之前在此站点上的帖子的人可能会说,我一直在尝试在 Python 中实现使用 FFT 的 PDE 求解器。编程部分大体搞定了,但是程序产生了一个(很适合本站的)溢出错误(基本上它会增长很多,直到变成NaN
)。
排除所有其他可能性后,我将问题归结为 FFT 和我尝试求导的方式,因此我决定测试两种不同的 FFT(numpy
的 fft
模块和 pyFFTW
包),代码如下:
import pyfftw
import numpy as np
import matplotlib.pyplot as plt
def fftw_(y: np.ndarray) -> np.ndarray:
a = pyfftw.empty_aligned((N, N), dtype=np.float64)
b = pyfftw.empty_aligned((N, N//2+1), dtype=np.complex128)
fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_FORWARD', flags=('FFTW_MEASURE',), threads=12)
y_hat = fft_object(y)
return y_hat
def ifftw_(y_hat: np.ndarray) -> np.ndarray:
a = pyfftw.empty_aligned((N, N//2+1), dtype=np.complex128)
b = pyfftw.empty_aligned((N, N), dtype=np.float64)
fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_BACKWARD', flags=('FFTW_MEASURE',), threads=12)
y = fft_object(y_hat)
return y
def func(x: np.ndarray, y: np.ndarray) -> np.ndarray:
return np.exp(x)*np.sin(y)
dx = 0.02
x = np.arange(-1, 1, dx)
y = np.arange(-1, 1, dx)
X, Y = np.meshgrid(x, y)
N = len(x)
kxw, kyw = np.meshgrid(np.fft.rfftfreq(N, dx), np.fft.fftfreq(N, dx))
Lapw = -4*np.pi**2*(kxw**2+kyw**2)
kxnp, kynp = np.meshgrid(np.fft.fftfreq(N, dx), np.fft.fftfreq(N, dx))
Lapnp = -4*np.pi**2*(kxnp**2+kynp**2)
z = func(X, Y)
lap_z_w = ifftw_(Lapw*fftw_(z))
lap_z_np = np.fft.ifft2(Lapnp*np.fft.fft2(z))
lap_z_np_mag = np.abs(lap_z_np)
lap_z_np_ang = np.angle(lap_z_np)
plt.imshow(z, cmap='plasma')
plt.colorbar()
plt.savefig("f.png", dpi=200)
plt.clf()
plt.imshow(lap_z_w, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_fftw.png", dpi=200)
plt.clf()
plt.imshow(lap_z_np_mag, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_np_mag.png", dpi=200)
plt.clf()
plt.imshow(lap_z_np_ang, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_np_ang.png", dpi=200)
plt.clf()
这里 np.ndarray
的 Lapw
和 Lapnp
是我认为 应该做的离散拉普拉斯算子。而我选择的函数eˣsin(y)是调和函数,所以它的拉普拉斯应该为零
但程序的结果与这个预期值相去甚远。特别是我得到:
原函数f
f"Laplacian"和pyFFTW
f 的 "Laplacian" 与 Numpy 的大小和相位
查看这些图的值(请注意颜色栏中的范围以及 20000 不是 0 的任何近似值这一事实)清楚地说明了为什么我制作的程序会溢出,但是我不知道如何纠正这个。任何帮助将不胜感激。
这里的错误是您对函数的假设不正确。 e^x sin(y) 可能看起来是谐波,但您只针对 -1 < x,y < 1 计算了它。fft
将隐式地周期性地继续它,即您在函数的所有边缘都会出现不连续性。如果函数不连续,则它不是谐波的,尤其是在傅里叶变换中会出现发散。这就是使您的 FFT 在边缘发散的原因。除此之外,"far away"从边缘看,结果看起来符合预期。
读过我之前在此站点上的帖子的人可能会说,我一直在尝试在 Python 中实现使用 FFT 的 PDE 求解器。编程部分大体搞定了,但是程序产生了一个(很适合本站的)溢出错误(基本上它会增长很多,直到变成NaN
)。
排除所有其他可能性后,我将问题归结为 FFT 和我尝试求导的方式,因此我决定测试两种不同的 FFT(numpy
的 fft
模块和 pyFFTW
包),代码如下:
import pyfftw
import numpy as np
import matplotlib.pyplot as plt
def fftw_(y: np.ndarray) -> np.ndarray:
a = pyfftw.empty_aligned((N, N), dtype=np.float64)
b = pyfftw.empty_aligned((N, N//2+1), dtype=np.complex128)
fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_FORWARD', flags=('FFTW_MEASURE',), threads=12)
y_hat = fft_object(y)
return y_hat
def ifftw_(y_hat: np.ndarray) -> np.ndarray:
a = pyfftw.empty_aligned((N, N//2+1), dtype=np.complex128)
b = pyfftw.empty_aligned((N, N), dtype=np.float64)
fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_BACKWARD', flags=('FFTW_MEASURE',), threads=12)
y = fft_object(y_hat)
return y
def func(x: np.ndarray, y: np.ndarray) -> np.ndarray:
return np.exp(x)*np.sin(y)
dx = 0.02
x = np.arange(-1, 1, dx)
y = np.arange(-1, 1, dx)
X, Y = np.meshgrid(x, y)
N = len(x)
kxw, kyw = np.meshgrid(np.fft.rfftfreq(N, dx), np.fft.fftfreq(N, dx))
Lapw = -4*np.pi**2*(kxw**2+kyw**2)
kxnp, kynp = np.meshgrid(np.fft.fftfreq(N, dx), np.fft.fftfreq(N, dx))
Lapnp = -4*np.pi**2*(kxnp**2+kynp**2)
z = func(X, Y)
lap_z_w = ifftw_(Lapw*fftw_(z))
lap_z_np = np.fft.ifft2(Lapnp*np.fft.fft2(z))
lap_z_np_mag = np.abs(lap_z_np)
lap_z_np_ang = np.angle(lap_z_np)
plt.imshow(z, cmap='plasma')
plt.colorbar()
plt.savefig("f.png", dpi=200)
plt.clf()
plt.imshow(lap_z_w, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_fftw.png", dpi=200)
plt.clf()
plt.imshow(lap_z_np_mag, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_np_mag.png", dpi=200)
plt.clf()
plt.imshow(lap_z_np_ang, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_np_ang.png", dpi=200)
plt.clf()
这里 np.ndarray
的 Lapw
和 Lapnp
是我认为 应该做的离散拉普拉斯算子。而我选择的函数eˣsin(y)是调和函数,所以它的拉普拉斯应该为零
但程序的结果与这个预期值相去甚远。特别是我得到:
原函数f
f"Laplacian"和pyFFTW
f 的 "Laplacian" 与 Numpy 的大小和相位
查看这些图的值(请注意颜色栏中的范围以及 20000 不是 0 的任何近似值这一事实)清楚地说明了为什么我制作的程序会溢出,但是我不知道如何纠正这个。任何帮助将不胜感激。
这里的错误是您对函数的假设不正确。 e^x sin(y) 可能看起来是谐波,但您只针对 -1 < x,y < 1 计算了它。fft
将隐式地周期性地继续它,即您在函数的所有边缘都会出现不连续性。如果函数不连续,则它不是谐波的,尤其是在傅里叶变换中会出现发散。这就是使您的 FFT 在边缘发散的原因。除此之外,"far away"从边缘看,结果看起来符合预期。