如何解决求解泊松方程 FFT 域 vs 有限差分
我有当前方程的两个解:
# 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
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] = 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,nboundaryy)
结果是:
我知道使用 FD 方案是正确的,但不确定我在 FFT 中是否正确。我在 FFT 上看到一个圆形,这是正确的吗?
解决方法
有两个主要区别
对于有限差分,您正在计算离散差分,对于 FFT 解决方案,您只需计算连续空间上的毒算符并将其应用于您的方程。要以完全相同的方式计算有限差分,您需要在离散域中使用 而不是计算 fft,您可以做的是记住 fft(roll(x,1)) = exp(-2j * np.pi * np.fftfreq(N))* fft(x)
其中 roll 表示 oen 样本的循环移位。
另一点是您正在使用边界条件(墙上的零电位),快速而肮脏的解决方案是使用 method of image charges 来确保墙上的电位消失并计算增强空间上的毒药方程。如果您关心内存使用或解决方案纯度,您可以使用 sine transform ,它具有稍微复杂的转换公式,但可以在不增加空间的情况下进行计算,因为其定义迫使边界上的电位为零(因为sin(pi * n) = 0 对于任何整数 n)
频域的解是直接解,你用一个封闭的公式计算每个系数,然后进行傅里叶逆变换,不需要迭代。只要您以足够的准确度计算差异,准确度往往也不错。
如果你真的很担心这个,你应该关注像 (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
,您可能不会注意到任何舍入误差),所有这些都更为重要。
如果您在频域中进行计算并且对精度不满意,您可以随时通过有限差分方程的一些迭代来“细化”它。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。