求解泊松方程 FFT 域与有限差分
Solving Poisson equation FFT domain vs Finite difference
我现在的方程有两个解:
第一个是使用有限差分方案(代码如下):
# Some variable declarations
nx = 300
ny = 300
nt = 100
xmin = 0.
xmax = 2.
ymin = 0.
ymax = 1.
dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)
# Initialization
p = np.zeros((nx, ny))
pd = np.zeros((nx, ny))
b = np.zeros((nx, ny))
# Source
b[int(nx / 4), int(ny / 4)] = 100
b[int(3 * nx / 4), int(3 * ny / 4)] = -100
for it in range(nt):
pd = p.copy()
p[1:-1,1:-1] = (((pd[1:-1, 2:] + pd[1:-1, :-2]) * dy**2 +
(pd[2:, 1:-1] + pd[:-2, 1:-1]) * dx**2 -
b[1:-1, 1:-1] * dx**2 * dy**2) /
(2 * (dx**2 + dy**2)))
p[0, :] = 0
p[nx-1, :] = 0
p[:, 0] = 0
p[:, ny-1] = 0
使用 FFT 我有以下代码:
def poisson(b,nptx,npty,dx,dy,nboundaryx,nboundaryy):
p = np.zeros((nptx,npty))
ppad = np.zeros((nptx+nboundaryx,npty+nboundaryy))
phatpad = np.zeros((nptx+nboundaryx,npty+nboundaryy))
bpad = np.zeros((nptx+nboundaryx,npty+nboundaryy))
bpad[:nptx,:npty] = b
kxpad = 2*np.pi*np.fft.fftfreq(nptx+nboundaryx,d=dx)
kypad = 2*np.pi*np.fft.fftfreq(npty+nboundaryy,d=dy)
epsilon = 1.e-9
ppad = np.real(np.fft.ifft2(-np.fft.fft2(bpad)/np.maximum(kxpad[None, :]**2 + kypad[:, None]**2,epsilon)))
p = ppad[:nptx,:npty]
p[0,:] = 0
p[nptx-1,:] = 0
p[:,0] = 0
p[:,npty-1] = 0
return p
nptx = 300
npty = 300
b = np.zeros((nptx, npty))
b[int(nptx / 4), int(npty / 4)] = 100
b[int(3 * nptx / 4), int(3 * npty / 4)] = -100
xmin = 0.
xmax = 2.
ymin = 0.
ymax = 1.
nboundaryx = 0
nboundaryy = 0
dx = (xmax - xmin) / (nptx+nboundaryx - 1)
dy = (ymax - ymin) / (npty+nboundaryy - 1)
print(dx)
p = poisson(b,nptx,npty,dx,dy,nboundaryx,nboundaryy)
结果是:
第一张使用有限差分的图像
第二张图像使用 FFT
我知道使用 FD 方案是正确的,但不确定我是否正确使用了 FFT。我在 FFT 上看到一个圆形,这是正确的吗?
有两个主要区别
对于有限差分,您正在计算离散差分,对于 FFT 解决方案,您只需计算连续 space 上的毒算子并将其应用于您的方程式。要以完全相同的方式计算有限差分,您需要在离散域中使用而不是计算 fft,您可以做的是记住 fft(roll(x, 1)) = exp(-2j * np.pi * np.fftfreq(N))* fft(x)
其中 roll 表示 oen 样本的循环移位。
另一点是您正在使用边界条件(墙上的零电势)快速而肮脏的解决方案是使用 method of image charges to ensure the potential vanishes on the walls and compute the poison equation on the augmented space. If you care about memory usage or solution purity you could use the sine transform 其具有稍微复杂的翻译公式,但可以在不增加 space 因为根据其定义,势能在边界上被强制为零(因为对于任何整数 n,sin(pi * n) = 0)
频域的解是直接解,你用一个封闭的公式计算出每个系数然后进行傅里叶逆变换,不需要迭代。只要您以足够的精度计算差异,精度往往也不错。
如果你真的担心这个,你应该关注(1 - exp(2j*pi/N))
这样的差异,因为第二项接近1,有效位数会减少。但是您可以通过将其分解为 exp(1j*pi/N) * (exp(-1j*pi/N) - exp(1j*pi/N)) = exp(1j*pi/N) * (-2j * sin(pi/N))
来提高此类表达式的准确性,其中您有一个产品并且不会丢失任何重要位。如果您以单精度或半精度计算它,所有这些都更为重要(您可能不会注意到使用 numpy.float64
或 numpy.complex128
的任何舍入错误)。
如果您在频域中计算并且对精度不满意,您可以随时通过有限差分方程的一些迭代来“改进”它。
我现在的方程有两个解:
第一个是使用有限差分方案(代码如下):
# Some variable declarations
nx = 300
ny = 300
nt = 100
xmin = 0.
xmax = 2.
ymin = 0.
ymax = 1.
dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)
# Initialization
p = np.zeros((nx, ny))
pd = np.zeros((nx, ny))
b = np.zeros((nx, ny))
# Source
b[int(nx / 4), int(ny / 4)] = 100
b[int(3 * nx / 4), int(3 * ny / 4)] = -100
for it in range(nt):
pd = p.copy()
p[1:-1,1:-1] = (((pd[1:-1, 2:] + pd[1:-1, :-2]) * dy**2 +
(pd[2:, 1:-1] + pd[:-2, 1:-1]) * dx**2 -
b[1:-1, 1:-1] * dx**2 * dy**2) /
(2 * (dx**2 + dy**2)))
p[0, :] = 0
p[nx-1, :] = 0
p[:, 0] = 0
p[:, ny-1] = 0
使用 FFT 我有以下代码:
def poisson(b,nptx,npty,dx,dy,nboundaryx,nboundaryy):
p = np.zeros((nptx,npty))
ppad = np.zeros((nptx+nboundaryx,npty+nboundaryy))
phatpad = np.zeros((nptx+nboundaryx,npty+nboundaryy))
bpad = np.zeros((nptx+nboundaryx,npty+nboundaryy))
bpad[:nptx,:npty] = b
kxpad = 2*np.pi*np.fft.fftfreq(nptx+nboundaryx,d=dx)
kypad = 2*np.pi*np.fft.fftfreq(npty+nboundaryy,d=dy)
epsilon = 1.e-9
ppad = np.real(np.fft.ifft2(-np.fft.fft2(bpad)/np.maximum(kxpad[None, :]**2 + kypad[:, None]**2,epsilon)))
p = ppad[:nptx,:npty]
p[0,:] = 0
p[nptx-1,:] = 0
p[:,0] = 0
p[:,npty-1] = 0
return p
nptx = 300
npty = 300
b = np.zeros((nptx, npty))
b[int(nptx / 4), int(npty / 4)] = 100
b[int(3 * nptx / 4), int(3 * npty / 4)] = -100
xmin = 0.
xmax = 2.
ymin = 0.
ymax = 1.
nboundaryx = 0
nboundaryy = 0
dx = (xmax - xmin) / (nptx+nboundaryx - 1)
dy = (ymax - ymin) / (npty+nboundaryy - 1)
print(dx)
p = poisson(b,nptx,npty,dx,dy,nboundaryx,nboundaryy)
结果是:
第一张使用有限差分的图像
第二张图像使用 FFT
我知道使用 FD 方案是正确的,但不确定我是否正确使用了 FFT。我在 FFT 上看到一个圆形,这是正确的吗?
有两个主要区别
对于有限差分,您正在计算离散差分,对于 FFT 解决方案,您只需计算连续 space 上的毒算子并将其应用于您的方程式。要以完全相同的方式计算有限差分,您需要在离散域中使用而不是计算 fft,您可以做的是记住 fft(roll(x, 1)) = exp(-2j * np.pi * np.fftfreq(N))* fft(x)
其中 roll 表示 oen 样本的循环移位。
另一点是您正在使用边界条件(墙上的零电势)快速而肮脏的解决方案是使用 method of image charges to ensure the potential vanishes on the walls and compute the poison equation on the augmented space. If you care about memory usage or solution purity you could use the sine transform 其具有稍微复杂的翻译公式,但可以在不增加 space 因为根据其定义,势能在边界上被强制为零(因为对于任何整数 n,sin(pi * n) = 0)
频域的解是直接解,你用一个封闭的公式计算出每个系数然后进行傅里叶逆变换,不需要迭代。只要您以足够的精度计算差异,精度往往也不错。
如果你真的担心这个,你应该关注(1 - exp(2j*pi/N))
这样的差异,因为第二项接近1,有效位数会减少。但是您可以通过将其分解为 exp(1j*pi/N) * (exp(-1j*pi/N) - exp(1j*pi/N)) = exp(1j*pi/N) * (-2j * sin(pi/N))
来提高此类表达式的准确性,其中您有一个产品并且不会丢失任何重要位。如果您以单精度或半精度计算它,所有这些都更为重要(您可能不会注意到使用 numpy.float64
或 numpy.complex128
的任何舍入错误)。
如果您在频域中计算并且对精度不满意,您可以随时通过有限差分方程的一些迭代来“改进”它。