# FFT solver for the Poisson problem with Dirichlet boundary conditions

I am trying to solve the Poisson problem with Dirichlet boundary condition in 1D:

$$$$\begin{array}{rcl} - \mu \Delta u & = & f~in~[0,1], \\ u(0) & = & 0, \\ u(1) & = & 0, \end{array}$$$$ using the FFT method. However, I get strange results for the rate of convergence of the method, compared to the finite difference method: the FFT method only converge with order 1 while the finite difference method is of order 2 in space. Is it ok ?

Below is my code in python. Thank you for your help.

from scipy.fftpack import dst
import numpy as np
import matplotlib.pyplot as plt
import time

def Poisson_1D(mu,f,N):
x = np.linspace(0,1,N)
hh = (x[1]-x[0])**2
NN = x.size
u = np.zeros(NN)
d = 2*np.ones(NN-2)
dd = -1*np.ones(NN-3)
A = np.diag(d) + np.diag(dd,1) + np.diag(dd,-1)
ff = f(x[1:-1])*hh/mu
u[1:-1] = np.linalg.solve(A,ff)
return x,u

def Poisson_fft_1D(mu,f,N,L=1):
x = np.linspace(0,L,N)
ff = f(x)
F = dst(ff,type=1)
k = np.array([(np.pi*(i+1)/L)**2 for i in range(0,N)])
U = F/(mu*k)
u = dst(U,type=1)/(2*(N+1))
return x, u

def f3(x): return 16*mu*np.pi**2*np.sin(4*np.pi*x)
def u3(x): return np.sin(4*np.pi*x)

if __name__ == "__main__":

error_fd = []
error_fft = []
nlist = np.arange(3,13)
for n in nlist:
N = 2**n
mu = 1
# FINITE DIFFERENCE
t1 = time.time()
x,u = Poisson_1D(mu,f3,N)
t2 = time.time() - t1
ue = u3(x)
# FFT
t3 = time.time()
x_fft,u_fft = Poisson_fft_1D(mu,f3,N)
t4 = time.time() - t3
ue_fft = u3(x_fft)
# ERRORS
e1 = np.sqrt(np.sum((ue_fft-u_fft)**2))
e2 = np.sqrt(np.sum((ue-u)**2))
error_fft.append( e1 )
error_fd.append( e2  )
print("N : {2:10d}   ,   error fft : {0:5.3e}   ,   time fft : {3:5.3e}   ,   error df : {1:5.3e}   ,   time df : {4:5.3e}".format(e1,e2,N,t4,t2))

plt.plot(np.log2(2**nlist),np.log2(error_fft),"+-",label='FFT')
plt.plot(np.log2(2**nlist),np.log2(error_fd),"x-",label='Finite difference')
plt.plot(np.log2(2**nlist),-np.log2(2**nlist),"k--",label=r"$$y=-x$$")
plt.plot(np.log2(2**nlist),-2*np.log2(2**nlist),"b--",label=r"$$y=-2x$$")
plt.legend()
plt.show()

• I think that fourier transforms implicitly assume that your signal (in your case your u(x) profile) is periodic. You state your boundary conditions in a way that has a huge jump if periodically extended. The fourier base is made of sines and cosines, and it is very hard to impossible to approximate a jump like that. There are methods to mitigate: (e.g. zero padding), but the problem of gibbs oszillation will propably persist. That being said, in a periodic setting FT-methods can be quite good – MPIchael Apr 24 '19 at 14:17
• Hello MPIchael. Thank you for your read. I made a mistake writing the post, the boundary conditions are homogeneous Dirichlet conditions. The question has been edited. – PeteAgor Apr 24 '19 at 14:46
• The finite difference method can easily converge in higher orders, it is just the choice of the eigenvalues for the FFT modes. Actually, with higher FD orders the eigenvalues converge to the spectral ones. – Vladimir F Apr 24 '19 at 15:02
• Or to be more explicit, you can solve the finite difference using the FFT as well, to solve the system, you do not need the slow linalg.solve. That is what fast Poissons solvers do, they use FFT but can be used for both the spectral discretization or the finite difference discretization. – Vladimir F Apr 24 '19 at 15:06
• And when you do choose to solve the linear system for the finite differences, do not forget you should make use of the sparrse structure of the matrix and actually, here the matrix is tridiagonal and can be solved in O(n). – Vladimir F Apr 24 '19 at 15:13