我正在尝试为周期框中的不可压缩流体构建一个简单的求解器，以学习光谱法。我正在关注教科书（Peyret R. Spectral methods for incompressible viscous flow[M]. Springer Science & Business Media, 2013.）作者提出了一种计算空间周期解的傅里叶方法。第一步是将傅里叶变换应用于动量和连续性方程：

($A$是非线性项 A(U,U)，$f$是源项）

然后将傅里叶空间中的散度算子应用于动量方程，利用傅里叶空间中的连续性方程消去$p$在动量方程中：

看起来它消除了对每个时间步进行压力投影和求解的必要性，因为变换后的方程自动满足连续性。但是，当我尝试在 Matlab 中为 2d 周期域实现代码时，我没有得到合理的结果，也不知道如何修复它。

我的问题是，我是否错过了实现此方法的重要内容？或者只是我的代码中的一些错误？

（为了简洁起见，这只是代码的一部分）

%% Construct k
k1 = 1i*[0:N/2-1 0 -N/2+1:-1];
k2 = -([0:N/2 -N/2+1:-1]).^2;

% expand to 2d
k1 = repmat(k1,N,1);
k2 = repmat(k2,N,1);

%% Solve PDE
for n = 1:plotgap

    %% Time advance
    t = t+dt;

    %% ifft velocity field
    u_r = real(ifft2(u));v_r = real(ifft2(v));

    % gradient
    ux_r = real(ifft2(k1.*u)); uy_r = real(ifft2(k1'.*u));
    vx_r = real(ifft2(k1.*v)); vy_r = real(ifft2(k1'.*v));

    % Advection
    A1 = fft2(u_r.*ux_r + v_r.*uy_r);
    A2 = fft2(u_r.*vx_r + v_r.*vy_r);

    % Mass conservation
    C1 = (k1.*A1 + k1'.*A2)./(k2 + k2') .*k1;
    C2 = (k1.*A1 + k1'.*A2)./(k2 + k2') .*k1';
    C1(1,1) = 0;C2(1,1) = 0;

    % Diffusion
    D1 = nu*(k2 + k2') .* u;
    D2 = nu*(k2 + k2') .* v;

    %% Construct operator
    operator1 = -A1 + D1 + C1;
    operator2 = -A2 + D2 + C2;

    %% Update velocity

    if ~exist('tmp1','var'), tmp1 = operator1; end
    if ~exist('tmpp1','var'), tmpp1 = operator1; end
    if ~exist('tmppp1','var'), tmppp1 = operator1; end
    if ~exist('tmp2','var'), tmp2 = operator2; end
    if ~exist('tmpp2','var'), tmpp2 = operator2; end
    if ~exist('tmppp2','var'), tmppp2 = operator2; end

    u = u + dt*(55/24*operator1 - 59/24*tmp1 + 37/24*tmpp1 - 3/8*tmppp1);  % 4th -order Adam - bashforth
    tmppp1 = tmpp1;tmpp1 = tmp1;tmp1 = operator1; % update f values
    v = v + dt*(55/24*operator2 - 59/24*tmp2 + 37/24*tmpp2 - 3/8*tmppp2);
    tmppp2 = tmpp2;tmpp2 = tmp2;tmp2 = operator2; % update f values

end