1
$\begingroup$

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

\begin{equation} \begin{array}{rcl} - \mu \Delta u & = & f~in~[0,1], \\ u(0) & = & 0, \\ u(1) & = & 0, \end{array} \end{equation} 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()
Improve this question
$\endgroup$
6
  • $\begingroup$ 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 $\endgroup$
    – MPIchael
    Apr 24, 2019 at 14:17
  • $\begingroup$ 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. $\endgroup$
    – PeteAgor
    Apr 24, 2019 at 14:46
  • $\begingroup$ 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. $\endgroup$
    – Vladimir F Героям слава
    Apr 24, 2019 at 15:02
  • $\begingroup$ 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. $\endgroup$
    – Vladimir F Героям слава
    Apr 24, 2019 at 15:06
  • $\begingroup$ 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). $\endgroup$
    – Vladimir F Героям слава
    Apr 24, 2019 at 15:13

1 Answer 1

Reset to default
0
$\begingroup$

No, it certainly is not correct.

You would have identified the problem easily if you actually printed the solution and compared it with the exact one. The way how dst works means that the boundary points are not included. Therefore while your exact solutions contrains the boundary 0 values, the numerical solution using FFT does not.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Vladimir, you saved my day ! $\endgroup$
    – PeteAgor
    Apr 24, 2019 at 16:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.