Solving Poisson equation FFT domain vs Finite difference

Question

I have two solutions of the current equation:

The first one is using Finite difference scheme ( code below ):

# 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

Using FFT I have the following code:

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)

The results are:

1. First image using Finite Difference

2. Second image using FFT

I know using FD scheme is correct but not sure if I did in FFT correctly. I see a round shape on FFT, is this correct?

Answer 1

There are two main differences

For the finite differences you are calculating the discrete differences, for the FFT solution you are simply computing the poison operator on the continuous space and applying that to your equation. To compute the finite differences exactly the same way you would need to use the in the discrete domain instead of calculating the fft what you can do is to remember that fft(roll(x, 1)) = exp(-2j * np.pi * np.fftfreq(N))* fft(x) where roll denotes the circular shift by oen sample.

Other point is that you are using boundary conditions (zero potential on the walls) the quick and dirty solution is to use 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 that has slightly more complicated translation formulas, but can be computed without augmenting the space since the potential is forced to be zero on the boundaries by its definition (because sin(pi * n) = 0 for any integer n)

The solution in the frequency domain is a direct solution, you calculate each coefficient with a closed formula and then perform the inverse Fourier transform, no iteration is required. The accuracy tends to be good as well, as far as you compute the differences with enough accuracy.

If you are really worried about this, you should focus in differences like (1 - exp(2j*pi/N)) because the second term is close to 1, the number of significative bits will be reduced. But you can improve the accuracy of such expressions by factoring it as exp(1j*pi/N) * (exp(-1j*pi/N) - exp(1j*pi/N)) = exp(1j*pi/N) * (-2j * sin(pi/N)) where you have a product and you don't loose any significant bit. All of this is more important if you are compluting it in single or half precision (you probably will not notice any rounding error using numpy.float64 or numpy.complex128).

If you calculate in the frequency domain and you are not happy with the accuracy you can always "refine" it with some iterations of your finite differences equation.