求解Poisson方程时，如何在线性求解器中正确实现齐次Neumann边界条件？

我尝试求解一个线性Poisson方程，y方向区间[-1,1]采用齐次Neumann边界条件，x方向为周期性边界条件。此前测试零齐次Dirichlet边界条件时运行正常，但实现Neumann边界条件时出现问题。

为便于排查，我用2D代码求解1D问题（忽略x方向变化，仅在y方向求解），对应的1D问题如下：

示例代码 #

%2D Code solve 1D problem with Neumann BCs
clearvars; clc; close all;

Nx = 2; 
Ny = 10; 
Lx =3; 
kx  = fftshift(-Nx/2:Nx/2-1);    % wave number vector

%1. Exact Case vs Approximation Case
dx = Lx/Nx; % Need to calculate dx
% Use approximations for kx, ky, and k^2. These come from Birdsall and Langdon. Find page number and put it here at some point,
ksqu = (sin( kx * dx/2)/(dx/2)).^2 ;
kx = sin(kx * dx) / dx;
xi_x =  (2*pi)/Lx;
ksqu4inv = ksqu; 
ksqu4inv(abs(ksqu4inv)<1e-14) =1e-6; %helps with error: matrix ill scaled because of 0s
xi  = ((0:Nx-1)/Nx)*(2*pi);
x =  xi/xi_x;

ylow = -1; 
yupp =1; 
Ly = (yupp-ylow);
eta_ygl  = 2/Ly;
[D,ygl] = cheb(Ny);
ygl = (1/2)*(((yupp-ylow)*ygl) + (yupp+ylow));
D   = D*eta_ygl;
D2  = D*D;
BC=-D([1 Ny+1],[1 Ny+1])\D([1 Ny+1],2:Ny); %Homogenous Neumann BCs for |y|=1 
[X,Y]   = meshgrid(x,ygl);

%linear Poisson solved iteratively
Igl = speye(Ny+1);
div_x_act_on_grad_x = -Igl;  
%ZNy represents the operation of setting the boundary values of y component
%to zero:
ZNy = diag([0 ones(1,Ny-1) 0]); 
div_y_act_on_grad_y = D2*ZNy;
%ICs
u = (1/6)*Y .*(6*ylow*yupp-3*ylow*Y-3*yupp*Y+2* Y.^2);
uh = fft(u,[],2); 
duxk=(kx*1i*xi_x) .*uh; 
du2xk = (kx*1i*xi_x) .*duxk;
duyk = D *uh; 
du2yk = D *duyk;
n = ones(size(u));
invnek = fft(1./n,[],2);
nh = fft(n,[],2);
dnhdxk = (kx*1i*xi_x) .*nh;
dnhdyk =D * nh;

%build numerical source
puhnhk = dnhdxk .* duxk;
pduhdxdnhdxk = dnhdyk .* duyk;
pduhdx2nhk = nh .* du2xk; 
pnhdudx2k = nh .*du2yk;
NumSourcek =puhnhk + pduhdxdnhdxk + pduhdx2nhk + pnhdudx2k;

uold = ones(size(u));
uoldk = fft(uold,[],2);
err_max =1e-8; 
max_iter = 500;
Sourcek = NumSourcek; 
for iterations = 1:max_iter
    OldSolMax= max(max(abs(uoldk)));
    duhdxk = (kx*1i*xi_x) .*uoldk;
    %product:
    gradNgradUx = dnhdxk .* duhdxk;
    duhdyk = (D) *uoldk ; 
    gradNgradUy = dnhdyk .* duhdyk; 
    
    RHSk = Sourcek - (gradNgradUx + gradNgradUy); 
    Stilde = invnek .* RHSk;
    for m = 1:length(kx)
        L = div_x_act_on_grad_x * (ksqu4inv(m)*xi_x^2)+ div_y_act_on_grad_y;
        unewh(:,m) = L\(Stilde(:,m));
    end
        %enforce BCs   
        unewh([1 Ny+1],:) = BC*unewh(2:Ny,:); %Neumann BCs for |y|=1  
    NewSolMax= max(max(abs(unewh)));
    if phikmax < err_max 
        it_error = err_max /2;
    else
        it_error = abs( NewSolMax- OldSolMax) / NewSolMax;
    end
    if it_error < err_max
        break;
    end
    uoldk = unewh;
end
unew = real(ifft(unewh,[],2)); 
 figure
 surf(X, Y, unew);
 colorbar;
 title('Numerical solution of \nabla \cdot (n \nabla u) = s in 2D');
 xlabel('x'); ylabel('y'); zlabel('u_{numerical}');
 figure
 surf(X, Y, u);
 colorbar;
 title('Exact solution of \nabla \cdot (n \nabla u) = s in 2D');
 xlabel('x'); ylabel('y'); zlabel('u_{exact}');

依赖的Cheb(N)函数 #

% CHEB compute D = differentitation matrix, x = Chebyshev grid
function [D, x] = cheb(N)
  if N == 0, D = 0; x = 1; return, end
  x = cos(pi*(0:N)/N)';
  c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
  X = repmat(x,1,N+1);
  dX = X-X';
  D = (c*(1./c)')./(dX+(eye(N+1)));
  D = D - diag(sum(D'));

运行结果 #

结果显示边界条件未正确生效，问题核心在于Neumann边界条件的施加逻辑错误，结合谱方法的边界处理思路，修正方案如下：

关键修正点 #

1. 移除求解后强制修改边界值的逻辑
   删除主循环中unewh([1 Ny+1],:) = BC*unewh(2:Ny,:);这一行，谱方法中Neumann条件需要在离散算子矩阵构建阶段就嵌入，而非求解后修正解向量。

2. 修改离散算子矩阵L的边界行
   在for m = 1:length(kx)循环内部，构建L后直接替换边界行为Neumann条件的离散形式：

   L = div_x_act_on_grad_x * (ksqu4inv(m)*xi_x^2)+ div_y_act_on_grad_y;
   % 替换为齐次Neumann边界条件：∂u/∂y=0在y=±1处
   L(1,:) = D(1,:);
   L(end,:) = D(end,:);
   % 边界条件右端项设为0（齐次）
   Stilde(1,m) = 0;
   Stilde(end,m) = 0;
   unewh(:,m) = L\(Stilde(:,m));
   

3. 修正二阶导数算子的错误约束
   原代码中ZNy会强制边界的二阶导数为0，不符合Neumann条件逻辑，修改：

   % 移除ZNy，直接使用完整的二阶导数矩阵
   div_y_act_on_grad_y = D2;
   

4. 修复迭代终止条件的未定义变量
   原代码中phikmax未定义，替换为已有的变量：

   if OldSolMax < err_max 
       it_error = err_max /2;
   else
       it_error = abs( NewSolMax- OldSolMax) / NewSolMax;
   end
   

内容的提问来源于stack exchange，提问作者Jamie