我正在尝试使用傅里叶变换微分（https://en.wikibooks.org/wiki/Parallel_Spectral_Numerical_Methods/Finding_Derivatives_using_Fourier_Spectral_Methods#Taking_a_Derivative_in_Fourier_Space）来实现拉普拉斯算子。我正在将傅里叶微分与真实空间中的周期性拉普拉斯算子进行比较。

当我的函数是 cos(x)*cos(y) 时，两个结果（实数和傅里叶空间）相互一致：

但是，当我的函数是 sin(x)*sin(y) 时，傅里叶空间拉普拉斯算子中有奇怪的振荡（一般形状保持不变）。在下图中，可以看到左图（傅里叶空间微分）有点模糊：

通过左上图的示例横截面如下所示：

我很困惑为什么会出现这些振荡，我希望有一个解决方案。我怀疑这可能与混叠有关，但我一直无法查明问题或解决方案。

请在此处找到实现傅立叶空间拉普拉斯算子的代码：

Nx = 512;
Ny = 256;

Lx = 2 * pi;
Ly = 4 * pi;

x = linspace(0, Lx, Nx)';
y = linspace(0, Ly, Ny)';

dx = x(2) - x(1);
dy = y(2) - y(1);

% Fourier space vectors
kx = 2*pi/Lx*[0:Nx/2-1 0 -Nx/2+1:-1]'; 
ky = 2*pi/Ly*[0:Ny/2-1 0 -Ny/2+1:-1]';

[x_grid, y_grid] = meshgrid(x, y);
[kx_grid, ky_grid] = meshgrid(kx, ky);

v = sin( 2*pi*x_grid ./ Lx ) .* sin( 2*pi*y_grid ./ Ly ); % Test Function
v_hat = fft2(v); % FFT of test function

% Fourier Space Laplacian
w_L = -( kx_grid.^2 + ky_grid.^2 ) .* v_hat;
w_L = real(ifft2(w_L));

% Periodic Real Space Laplacian
L = Laplacian_2D(v, dx, dy);

c = 10; % Because the edges get very weird
figure; mesh(x(c + 1:end-c),y(c + 1:end-c),w_L(c + 1:end-c,c+1:end-c)); title('Fourier');
figure; mesh(x(c + 1:end-c),y(c + 1:end-c),L(c + 1:end-c,c+1:end-c)); title('Real Space');

也请在这里找到我的 2D 周期性拉普拉斯算子：

function out = Laplacian_2D( f, hx, hy )

    f = f';
    [r c] = size(f);
    out = zeros(r, c);

    out( 2 : end - 1, 2 : end - 1 ) = ...
        ( f( 3 : end, 2 : end - 1 ) - ...
        2 * f( 2 : end - 1, 2 : end - 1 ) + ...
        f( 1 : end - 2, 2 : end - 1 ) ) / hx^2 + ...
        ( f( 2 : end - 1, 3 : end ) - ...
        2 * f( 2 : end - 1, 2 : end - 1 ) + ...
        f( 2 : end - 1, 1 : end - 2 ) ) / hy^2 ;

    % Remember (i, j)
    %  i = row, j = column
    %  i = y, j = x
    % but we transposed so we can say i = x, j = y

    % Periodic boundary conditions

    % 4 corners
    out( 1, 1 ) = ...
        ( f(2, 1) - 2 * f(1, 1) + f(end, 1) ) / hx^2 + ...
        ( f(1, 2) - 2 * f(1, 1) + f(1, end) ) / hy^2;
    out(1, end) = ...
        ( f(2, end) - 2 * f(1, end) + f(end, end) ) / hx^2 + ...
        ( f(1, 1) - 2 * f(1, end) + f(1, end - 1) )/ hy^2;
    out(end, 1) = ...
        ( f(1, 1) - 2 * f(end, 1) + f(end - 1, 1) ) / hx^2 + ...
        ( f(end, 2) - 2 * f(end, 1) + f(end, end) ) / hy^2;
    out(end, end) = ...
        ( f(1, end) - 2 * f(end, end) + f(end - 1, end) ) / hx^2 + ...
        ( f(end, 1) - 2 * f(end, end) + f(end, end - 1) ) / hy^2;

    % edges
    out( 1, 2 : end - 1 ) = ...
        ( f(2, 2 : end - 1) - 2 * f(1, 2 : end - 1) + f(end, 2 : end - 1) ) / hx^2 + ...
        ( f(1, 3 : end) - 2 * f(1, 2 : end - 1) + f(1, 1 : end - 2) ) / hy^2;
    out( end, 2 : end - 1 ) = ...
        ( f(1, 2 : end - 1) - 2 * f(end, 2 : end - 1) + f(end - 1, 2 : end - 1) ) / hx^2 + ...
        ( f(end, 3 : end) - 2 * f(end, 2 : end - 1) + f(end, 1 : end - 2) ) / hy^2;

    out( 2 : end - 1, 1) = ...
        ( f(3 : end, 1) - 2 * f(2 : end - 1, 1) + f(1 : end - 2, 1) ) / hx^2 + ...
        ( f(2 : end - 1, 2) - 2 * f(2 : end - 1, 1) + f(2 : end - 1, end) ) / hy^2;
    out( 2 : end - 1, end) = ...
        ( f(3 : end, end) - 2 * f(2 : end - 1, end) + f(1 : end - 2, end) ) / hx^2 + ...
        ( f(2 : end - 1, 1) - 2 * f(2 : end - 1, end) + f(2 : end - 1, end - 1) ) / hy^2;

    out = out';

end