1
$\begingroup$

I am trying to solve a standard Poisson equation on image with Neumann boundary condition. The Poisson equation has the form:

$Δp=b$

where $b$ is the divergence of the input image,which is as shown: enter image description here

I tried to write it in Matlab with centred finite difference method and Gauss-Seidel iterative method but the result turns out to be very weird, which is as shown: enter image description here

Can any kind soul help to spot errors in my code below:

{%Specifying parameters
nx=rows;                           %Number of steps in space(x)(vertical)
ny=columns;                        %Number of steps in space(y)(horizontal)
dx=1;
dy=1;
h=dx;




%Solving (Dxx+Dyy)p(i,j)=b(i,j) (eqn 3.16) using centred FDM (5-points difference)
%Gauss-Seidel iterative method
pn=zeros(nx,ny);                 %Preallocating pn
p=zeros(nx,ny);                  %Preallocating p
b=zeros(nx,ny);
ax=0;                            %counter 

for j=2:ny-1
   for i=2:nx-1
   b(i,j)=0.5*(grayImage1(i+1,j)-grayImage1(i-1,j)+grayImage(i,j+1)-grayImage(i,j-1))./h;

   end
end 

i=2:nx-1; 
j=2:ny-1;
for it=1:200000
        p(1,1)=0;        
        pn=p;   
        p(i,j)=0.25*(((pn(i+1,j)+p(i-1,j)))+((pn(i,j+1)+p(i,j-1)))-h^2*(b(i,j)));

    %%boundary conditions 
    %The four corners
    %p(1,1)=1/2*pn(2,1)+1/2*pn(1,2)-1/4*dx^2*b(1,1);            
    p(1,ny)=1/2*pn(2,ny)+1/2*pn(1,ny-1)-1/4*dx^2*b(1,ny);            
    p(nx,1)=1/2*pn(nx,2)+1/2*pn(nx-1,1)-1/4*dx^2*b(nx,1);          
    p(nx,ny)=1/2*pn(nx-1,ny)+1/2*pn(nx,ny-1)-1/4*dx^2*b(nx,ny);

    %The four sides
    p(1,j)=1/4*(2*pn(2,j)+pn(1,j+1)+pn(1,j-1)-dx^2*b(1,j));            
    p(nx,j)=1/4*(2*pn(nx-1,j)+pn(nx,j+1)+pn(nx,j-1)-dx^2*b(nx,j));            
    p(i,1)=1/4*(2*pn(i,2)+pn(i+1,1)+pn(i-1,1)-dy^2*b(i,1));
    p(i,ny)=1/4*(2*pn(i,ny-1)+pn(i+1,ny)+pn(i-1,ny)-dy^2*b(i,ny));

    if abs(p(i,j)-pn(i,j))<=zeros(nx-2,ny-2)+0.0001
        break
    end
    ax=ax+1
end }

Thanks very much!!

p/s: I do notice that pure Neumann boundary condition will give rise to non-unique solutions. Hence I followed one of the advice from other answers to set one of the corner to zero. Unfortunately it seems to have no effect on the results.

$\endgroup$
2
  • 3
    $\begingroup$ What exactly do you mean by "the result turns out to be very weird"? Knowing how a solution is wrong often tells a lot about where the bug could be. $\endgroup$ – Wolfgang Bangerth Jan 5 '17 at 16:39
  • $\begingroup$ As a side note: The proper spelling of the name is "Seidel", not "Seidler". $\endgroup$ – Wolfgang Bangerth Jan 5 '17 at 16:39
2
$\begingroup$

I cannot recommend fixing a corner to zero, that changes the problem you are trying to solve, instead just add before you iterate a line with:

b = b - mean(b);

This is called the Discrete Compatibility Criteria. May I refer to my own answer to the following: Solving a linear equation system with pure Neumann condition

In an effort to make this an answer, the following code will generate the matrix version of your code and let you use e.g. conjugate gradient as a solver. I grabbed it out of some of the code on: http://math.mit.edu/classes/18.086/2008/ It combines two 1d matrices for the all-Neumann Poisson problem into a large matrix for the 2d problem using the Kronecker Product:

maxit = 1000;
tol = 1e-6;

e = ones(nx-4,1);
A1= spdiags([1 -1 1; e*[1 -2 1]; 1 -1 1],-1:1,nx-2,nx-2); %/dx^2;
e = ones(ny-4,1);
A2= spdiags([1 -1 1; e*[1 -2 1]; 1 -1 1],-1:1,ny-2,ny-2); %/dy^2;
A = -( kron(speye(ny-2),A1)+kron(A2,speye(nx-2)) ); clear A1 A2
b = reshape(b,(nx-2)*(ny-2),1);
b = b - mean(b);

x = pcg(A,b,tol,maxit);

x = reshape(x,nx-2,ny-2);
$\endgroup$

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.