如何在我的线性求解器中正确实现齐次诺依曼边界条件？

Question

我正在尝试求解沿 y 方向间隔 [-1,1] 且沿 x 周期性的具有齐次诺依曼边界条件的线性泊松方程。最初，我针对零齐次狄利克雷边界条件测试了相同的代码，这相当简单。现在，我在尝试实现 Neumann BC 时遇到了一些问题。 为简单起见，我的 2D 代码解决了 1D 问题，这意味着不考虑 x 的变化并且仅求解 y 以更好地诊断我的问题。

我正在解决的一维问题看起来像：

示例代码：

%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'));

这将返回以下结果：

这清楚地表明我的边界条件设置不正确！我想我正在尝试将二阶导数（D2）的第一行/最后一行替换为一阶导数的第一行/最后一行，并在我的主循环中强制执行，如本文所示：我在解决热量问题时在哪里犯了错误使用谱法（切比雪夫微分矩阵）的方程？

Answer 1

您可以尝试强制执行诺依曼边界条件，直接使用边界条件更新内部点。您可以通过在求解内点之前将诺伊曼条件项添加到右侧 (RHS) 来实现此目的。

for iterations = 1:max_iter
    % Your existing code for computing RHS and Stilde
    
    % Add Neumann boundary condition contribution to RHS
    RHSk([1 Ny+1], :) = RHSk([1 Ny+1], :) - BC*unewh(2:Ny, :);
    
    % Solve for interior points
    for m = 1:length(kx)
        L = div_x_act_on_grad_x * (ksqu4inv(m)*xi_x^2) + div_y_act_on_grad_y;
        unewh(2:Ny, m) = L\(Stilde(2:Ny, m));
    end

    % Your existing code for convergence check and updating uoldk
end