我使用伪光谱方法在给定的时间从 DNS 获得了我的 3 维速度流场 u、v 和 w。我需要计算能谱（在傅立叶空间中）作为波数大小的函数，即$E(k)$作为一个函数$k$. 我使用的能谱方程如下：

在哪里$dA(k) = 4 \pi k^{2}dk$,$\hat{u}_i$是的 FFT$\hat{u}_{i}$和$\hat{u}_i^{*}$是转置共轭$\hat{u}_{i}$.$k = k_{x}^2 + k_{y}^2 + k_{z}^2$， 在哪里$k_{i}$是沿 i 的波数。它从$-N/2$到$N/2-1$. FFT 是使用 FFTW 进行的。

在数值上，我创建了 2 个 3 维数组，其中包含跨越整个 3D 域E的能谱和mk波数幅度。它们被初始化如下（$nx$,$ny$,$nz$是沿线的点数$x$,$y$和$z$轴）：

do k = 1, nz
  do j = 1, ny
    do i = 1, nx
      mk(i,j,k) = dsqrt(kx(i)**2 + ky(j)**2 + kz(k)**2)
      E(i,j,k) = 4*pi*(mk(i,j,k)**2)*(dconjg(uhat(i,j,k))*uhat(i,j,k)+dconjg(vhat(i,j,k))*vhat(i,j,k)+dconjg(what(i,j,k))*what(i,j,k))
    end do
  end do
end do

现在是令人困惑的部分。我需要执行一种称为“装箱”的技术。这涉及将波数范围分成适当相等的部分，并取落入每个部分的能量的平均值。为此，我将 2 个 3D 数组折叠成长E度mk为 1D 的数组1:nx*ny*nz。然后，我按 mk 的升序对 E 进行排序（按波数幅度的升序对能谱进行排序）。最后，我添加E了特定 bin 范围之间的连续值，除以添加的值的数量并写入文件以进行绘图。

红色曲线表示正确的结果。绿色的是我的代码产生的结果。可以看出，绿色显示了总体趋势，但还不够。还会出现几个尖峰。谁能指出我程序中的差异？我愿意应要求分享代码。

编辑：在实施詹姆斯的建议后，我的结果有了显着改善。（此外，我在将参考文献中的结果数字化时也犯了一个错误）。改进后的结果如下所示。

但是，参考曲线仍然比我的曲线多走几个波数。对于下一个波数，我的代码会产生一个尖峰。