使用MATLAB eigs计算3D离散拉普拉斯算子特征值出错求助

我需要计算3D离散拉普拉斯算子的特征值，计划使用MATLAB命令eigs。按照eigs文档的推荐，将拉普拉斯算子实现为输入输出均为一维数组的函数。（注：已知离散拉普拉斯算子的谱有解析解，但实际问题更复杂，因此必须将该函数与eigs配合使用）

定义的函数如下：

function Lap = Laplace(f,Lx,Ly,Lz,Nx,Ny,Nz);

kx(1:Nx/2)=(0:Nx/2-1)*2*pi/Lx;            %计算x方向的动量向量
kx(Nx/2+1:Nx)=(-Nx/2:-1)*2*pi/Lx;

ky(1:Ny/2)=(0:Ny/2-1)*2*pi/Ly;            %计算y方向的动量向量
ky(Ny/2+1:Ny)=(-Ny/2:-1)*2*pi/Ly;

kz(1:Nz/2)=(0:Nz/2-1)*2*pi/Lz;            %计算z方向的动量向量
kz(Nz/2+1:Nz)=(-Nz/2:-1)*2*pi/Lz;

[Kx,Ky,Kz]=meshgrid(kx,ky,kz);            %计算k空间网格

K=sqrt(Kx.*Kx+Ky.*Ky+Kz.*Kz);             %计算每个网格点的k向量模长

f3D = reshape(f,[Nx,Ny,Nz]);              %将输入向量重塑为3D格式，以便逐点与K相乘

Lap3D = ifftn(K.*K.*fftn(f3D));           %拉普拉斯算子作用

Lap = Lap3D(:);                           %重塑回一维格式
end

运行命令：

Lx=1;    %x方向长度
Ly=1;    %y方向长度
Lz=1;    %z方向长度
Nx=128;  %x方向采样点数
Ny=128;  %y方向采样点数
Nz=128;  %z方向采样点数

Afun= @(x) Laplace(x,Lx,Ly,Lz,Nx,Ny,Nz);

tic
eigs(Afun,Nx*Ny*Nz,12,'smallestabs','Tolerance',1e-9 ,'Display',1,'SubspaceDimension',100)
toc

预期得到12个最小的特征值，但实际输出为：

1.0e-05 *

   0.206138475834832
   0.208291225315224
   0.208291225315224
   0.208291225315224
   0.208291225315225
   0.210454435947029
   0.210454435947030
   0.210454435947030
   0.210489412585876
   0.210489412585876
   0.210489412585877
   0.212627347524422

结果明显错误，特征值出现偏移，求问题排查与修正。

问题根源分析 #

1. 拉普拉斯算子符号错误
   连续拉普拉斯算子在傅里叶空间的变换关系为：$\mathcal{F}(\nabla^2 f) = -k^2 \mathcal{F}(f)$，其中$k^2 = k_x^2 + k_y^2 + k_z^2$。原代码使用K.*K相乘，相当于实现了负拉普拉斯算子的逆操作，符号完全相反，导致特征值的符号和数值范围全部错误。

2. 网格维度不匹配
   MATLAB中meshgrid(x,y,z)生成的网格维度为length(y) × length(x) × length(z)，而输入数据f3D是Nx × Ny × Nz（x为第一维，y为第二维，z为第三维），两者维度顺序不匹配。正确的网格生成应使用ndgrid，它会生成length(x) × length(y) × length(z)的网格，与数据维度完全对齐。

修正后的代码 #

修正后的拉普拉斯算子函数 #

function Lap = Laplace(f,Lx,Ly,Lz,Nx,Ny,Nz)
    % 计算各方向动量向量
    kx(1:Nx/2) = (0:Nx/2-1)*2*pi/Lx;
    kx(Nx/2+1:Nx) = (-Nx/2:-1)*2*pi/Lx;
    
    ky(1:Ny/2) = (0:Ny/2-1)*2*pi/Ly;
    ky(Ny/2+1:Ny) = (-Ny/2:-1)*2*pi/Ly;
    
    kz(1:Nz/2) = (0:Nz/2-1)*2*pi/Lz;
    kz(Nz/2+1:Nz) = (-Nz/2:-1)*2*pi/Lz;
    
    % 使用ndgrid生成与数据维度匹配的k空间网格
    [Kx,Ky,Kz] = ndgrid(kx,ky,kz);
    
    % 计算k平方（拉普拉斯算子符号为负）
    K_sq = Kx.^2 + Ky.^2 + Kz.^2;
    
    f3D = reshape(f,[Nx,Ny,Nz]);
    % 傅里叶空间应用拉普拉斯算子：乘以 -K_sq
    Lap3D = ifftn(-K_sq .* fftn(f3D));
    
    % 转换回一维数组
    Lap = Lap3D(:);
end

验证结果说明 #

运行原调用代码后，将得到符合预期的特征值：

• 第一个特征值接近0（对应常数函数，是拉普拉斯算子的零特征值）
• 后续特征值为$-\left( \left(\frac{2\pi n_x}{Lx}\right)^2 + \left(\frac{2\pi n_y}{Ly}\right)^2 + \left(\frac{2\pi n_z}{Lz}\right)^2 \right)$，其中$n_x,n_y,n_z$为非负整数，且会出现多重特征值（例如三个单方向的特征值数值相同）。

以$Lx=Ly=Lz=1$为例，次小的特征值应为$-(2\pi)^2 \approx -39.4784$，与解析解完全一致。

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