-2
$\begingroup$

enter image description hereI 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)
$\endgroup$
  • $\begingroup$ Sorry I am new registered user here so not allowed to comment under one's answer yet. Just wondering are you still interested in this topic. I found it quite tricky to verify the solutions, as I conduct the second derivative of the solution and compared it with the right hand side to find there is significant difference between each. Thanks for your time:) Regards, Xinxin $\endgroup$ – Xinxin Peng Sep 10 '18 at 6:00
4
$\begingroup$

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:

solution to Poisson equation

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.

$\endgroup$
  • $\begingroup$ ,Thank you for quick answer. Sorry for the confusing question, I used boundary condition u(1)=u(end)=0 in SOR solver but implied a periodic boundary condition in spectral solver. If using periodic B.C for SOR solver, which is u(1)=u(999);u(end)=u(2), and increase the iteration number, we will get the same result (subplot 3 shown in my question). $\endgroup$ – Yang Zhang Nov 27 '14 at 18:14
  • $\begingroup$ However, now I am trying to solve this equation under the B.C $u(1)=u(end)=0$ by using spectral method. The only spectral method I am familiar is FFT which implies a periodic B.C as you said. Do you have any idea on how to solve Poisson eq. with fixed B.C by fft? In Matlab the first resolved wavenumber is $kx=0$,if $kx$ is going to the denominator, the first element of the wave number series cannot be zero, that is why I put $slice(1)=0$ in my code, it also assumes the averaged stream function is zero. It may not be appropriate way. $\endgroup$ – Yang Zhang Nov 27 '14 at 18:43
  • $\begingroup$ @YangZhang, Spectral methods require periodic boundary conditions. If you want any other boundary condition you will need something quite a bit more complicated such as a buffer zone that enforces the desired boundary conditions at the edge of the physical domain, but adds an non-physical subdomain to permit periodicity. That's a whole other topic by itself though. $\endgroup$ – Doug Lipinski Nov 27 '14 at 20:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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