explain the difference between 1D Poisson solvers

Question

I compared 2 methods for 1 dimensional Poisson equation solution. One is finite-difference method, "Successive Overrelaxation" from http://www.cs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html; Another method is spectral method by FFT which is much faster. The following is my matlab/octave code to generate the solutions in the picture. Obviously, the traditional finite-difference method(middle panel) is more acceptable! WHY FFT method look like this? Is may FFT-solution implement correct? Thanks!!

The following is my matlab code:::

clear all;

x=1:1000;

vort=exp(-(x-600).^2/50^2);% right hand side of poisson eq

% construct the wave number series, dx = 1 in this case
kx=(x-1);
jj=find(kx>500);
kx(jj)=kx(jj)-1000;
kx=2*pi/1000*kx;

subplot(311);plot(x,vort);title('rhs of Poisson eq','fontsize',20)

% ------------S.O.R solver----------

psi=zeros(1,1000);

for n=1:100
    for i=2:999

        r=(psi(i-1)+psi(i+1)-vort(i))*0.5-psi(i);
        psi(i)=psi(i)+1.2*r;

    end
end

subplot(312);plot(x,psi);title('Successive Overrelaxation','fontsize',20)

% ------ FFT solver -----%

slice=-fft(vort)./(kx.^2);

slice(1)=0;

subplot(313);plot(real(ifft(slice)));;title('FFT method','fontsize',20)

Answer 1

First, some notation for clarity. We're talking about solving $u_{xx}=f(x)$. So $u$ is the solution and $f$ is the right hand side.

What makes you think either of your solutions is or is not correct?
To check your results, just compute two derivatives of your solution (using finite differences or FFT, either one) and compare with $f$. In your example, $f$ is always positive so the solution should be concave up everywhere. In this case they don't match for either of your solutions. Now you need to figure out why. You should also use consistent boundary conditions if you want to compare two solutions. You are using boundary conditions $u(1)=u(1000)=0$ for the SOR solution, but the spectral (FFT) solution implies periodic boundary conditions.

The SOR solution is clearly concave down on a couple intervals since it is not converged. Convergence is extremely slow for this problem, you will likely need a couple hundred thousand SOR iterations instead of just 100. Consider that in SOR information can move only 1 gridpoint to the left at each iteration and your domain is ~1000 points. Clearly 100 iterations will not be enough.

The FFT solution implies smooth, periodic boundary conditions, but a function cannot be smooth, periodic, and always concave up. You may also have slight issues near the boundaries since your right hand side is not strictly periodic. I believe your FFT method is correct (although you should use ifft(...,'symmetric') rather than real() to ensure real results), but the problem is ill posed. Instead of the expected solution, you get a solution where the error is distributed throughout the domain. That's why the solution is concave down where it is expected to be nearly linear. In fact, the error in the second derivative is uniform throughout the domain and is a constant difference of $\int_1^{1000}f(x)dx/1000\approx0.0886$ in the second derivative of your solution. This shift (the mean of $f$) is lost when you set slice(1)=0; and reflects the fact that the net effect (i.e. the mean) of the $u_{xx}$ must be zero to have a smooth, periodic solution. Put another way, adding any constant value to you $f$ has no effect on the solution in the spectral method since the method enforces smooth periodicity which implies zero mean for the derivatives of $u$.

---

For reference, the solution (with second order centered finite differences and zero boundary conditions) can be found using a direct solver with the following code:

% Generate sparse matrix for the linear system:
i = [2:999, 2:999, 2:999, 1, 1000];
j = [1:998, 2:999, 3:1000, 1, 1000]; 
s = [ones(1,998), -2*ones(1,998), ones(1,998), 1, 1];
M = sparse(i,j,s);

% Generate right hand side:
x = (1:1000)';
rhs = [0; exp(-(x(2:999)-600).^2/50^2); 0];

% Solve:
y = M\rhs;

% Plot:
figure, plot(x,y)

and should look like this:

Notice that the solution is approximately linear except in the interval [500, 700] where $f$ has the largest values so $u$ is clearly concave up.