我正在尝试用 Neumann 边界条件求解图像上的标准泊松方程。泊松方程具有以下形式：

$Δp=b$

其中$b$是输入图像的散度，如图：

我尝试用中心有限差分法和Gauss-Seidel迭代法在Matlab中编写，但结果很奇怪，如图所示：

任何善良的灵魂都可以帮助我在下面的代码中发现错误：

{%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 }

非常感谢！！

p/s：我确实注意到纯 Neumann 边界条件会产生非唯一解。因此，我遵循了其他答案中的一个建议，将其中一个角设置为零。不幸的是，它似乎对结果没有影响。