为什么 (A+B) 的 FFT 不同于 FFT(A) + FFT(B)？

Question

我已经与一个非常奇怪的错误作斗争将近一个月了。问你们是我最后的希望。我用 C 编写了一个程序，它使用傅立叶（或倒数）space:

中的隐式欧拉 (IE) 方案集成了 2d Cahn–Hilliard equation

其中 "hats" 表示我们在傅立叶 space 中：h_q(t_n+1) 和 h_q(t_n ) 是 h(x,y) 在时间 t_n 和 t_(n+1) 的 FT，N[h_q] 是应用于 h_q 的非线性算子，在傅里叶 space，L_q 是线性的，也是傅里叶 space。我不想过多地讨论我正在使用的数值方法的细节，因为我确定问题不是出自那里（我尝试使用其他方案）。

我的代码其实很简单。这是开始，我基本上在这里声明变量、分配内存并为 FFTW 例程创建计划。

# include <stdlib.h>
# include <stdio.h>
# include <time.h>
# include <math.h>
# include <fftw3.h>
# define pi M_PI

int main(){

// define lattice size and spacing
int Nx = 150;         // n of points on x
int Ny = 150;         // n of points on y
double dx = 0.5;      // bin size on x and y

// define simulation time and time step
long int Nt = 1000;   // n of time steps
double dt = 0.5;      // time step size

// number of frames to plot (at denominator)
long int nframes = Nt/100;

// define the noise
double rn, drift = 0.05;   // punctual drift of h(x)
srand(666);                // seed the RNG

// other variables
int i, j, nt;    // variables for space and time loops

// declare FFTW3 routine
fftw_plan FT_h_hft;   // routine to perform  fourier transform
fftw_plan FT_Nonl_Nonlft;
fftw_plan IFT_hft_h;  // routine to perform  inverse fourier transform

// declare and allocate memory for real variables
double *Linft = fftw_alloc_real(Nx*Ny);
double *Q2 = fftw_alloc_real(Nx*Ny);
double *qx = fftw_alloc_real(Nx);
double *qy = fftw_alloc_real(Ny);

// declare and allocate memory for complex  variables
fftw_complex *dh = fftw_alloc_complex(Nx*Ny);
fftw_complex *dhft = fftw_alloc_complex(Nx*Ny);
fftw_complex *Nonl = fftw_alloc_complex(Nx*Ny);
fftw_complex *Nonlft = fftw_alloc_complex(Nx*Ny);

// create the FFTW plans
FT_h_hft = fftw_plan_dft_2d ( Nx, Ny, dh, dhft, FFTW_FORWARD, FFTW_ESTIMATE );
FT_Nonl_Nonlft = fftw_plan_dft_2d ( Nx, Ny, Nonl, Nonlft, FFTW_FORWARD, FFTW_ESTIMATE );
IFT_hft_h = fftw_plan_dft_2d ( Nx, Ny, dhft, dh, FFTW_BACKWARD, FFTW_ESTIMATE );

// open file to store the data
char acstr[160];
FILE *fp;
sprintf(acstr, "CH2d_IE_dt%.2f_dx%.3f_Nt%ld_Nx%d_Ny%d_#f%.ld.dat",dt,dx,Nt,Nx,Ny,Nt/nframes);

在此序言之后，我用均匀的随机噪声初始化我的函数 h(x,y)，并且我还对其进行了 FT。我将 h(x,y) 的虚部，即代码中的 dh[i*Ny+j][1] 设置为 0，因为它是实函数。然后我计算波向量 qx 和 qy，并用它们计算傅里叶方程的线性算子 space，在代码中是 Linft。我只考虑 h 的 - 四阶导数作为线性项，因此线性项的 FT 只是 -q^4... 但是，我不想详细介绍我的积分方法。问题不在于此。

// generate h(x,y) at initial time
for ( i = 0; i < Nx; i++ ) {
  for ( j = 0; j < Ny; j++ ) {
    rn = (double) rand()/RAND_MAX;    // extract a random number between 0 and 1
    dh[i*Ny+j][0] = drift-2.0*drift*rn;    // shift of +-drift
    dh[i*Ny+j][1] = 0.0;
  }
}

// execute plan for the first time
fftw_execute (FT_h_hft);

// calculate wavenumbers
for (i = 0; i < Nx; i++) { qx[i] = 2.0*i*pi/(Nx*dx); }
for (i = 0; i < Ny; i++) { qy[i] = 2.0*i*pi/(Ny*dx); }
for (i = 1; i < Nx/2; i++) { qx[Nx-i] = -qx[i]; }
for (i = 1; i < Ny/2; i++) { qy[Ny-i] = -qy[i]; }

// calculate the FT of the linear operator
for ( i = 0; i < Nx; i++ ) {
  for ( j = 0; j < Ny; j++ ) {
    Q2[i*Ny+j] = qx[i]*qx[i] + qy[j]*qy[j];
    Linft[i*Ny+j] = -Q2[i*Ny+j]*Q2[i*Ny+j];
  }
}

然后，终于到了时间循环。本质上，我所做的是：

每隔一段时间，我将数据保存到一个文件中，并在终端上打印一些信息。特别是，我打印了非线性项的 FT 的最高值。我还检查 h(x,y) 是否发散到无穷大（这不应该发生！），
直接space计算h^3（即dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0]）。同样，虚部设置为0，
取h^3的FT,
通过计算-q^2*(FT[ h^3] - 英尺 [h])。在代码中，我指的是行 Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][0] -dhft[i*Ny+j][0]) 和下面的虚数部分。我这样做是因为：

使用IE方法提前，直接变换回来space，然后归一化。

代码如下：

for(nt = 0; nt < Nt; nt++) {

if((nt % nframes)== 0) {
  printf("%.0f %%\n",((double)nt/(double)Nt)*100);
  printf("Nonlft   %.15f \n",Nonlft[(Nx/2)*(Ny/2)][0]);

  // write data to file
  fp = fopen(acstr,"a");
  for ( i = 0; i < Nx; i++ ) {
    for ( j = 0; j < Ny; j++ ) {
      fprintf(fp, "%4d  %4d  %.6f\n", i, j, dh[i*Ny+j][0]);
      }
  }
  fclose(fp);

}

// check if h is going to infinity
if (isnan(dh[1][0])!=0) {
  printf("crashed!\n");
  return 0;
}

// calculate nonlinear term h^3 in direct space
for ( i = 0; i < Nx; i++ ) {
  for ( j = 0; j < Ny; j++ ) {
      Nonl[i*Ny+j][0] = dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0];
      Nonl[i*Ny+j][1] = 0.0;
  }
}

// Fourier transform of nonlinear term
fftw_execute (FT_Nonl_Nonlft);

// second derivative in Fourier space is just multiplication by -q^2
for ( i = 0; i < Nx; i++ ) {
  for ( j = 0; j < Ny; j++ ) {
    Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][0] -dhft[i*Ny+j][0]);
    Nonlft[i*Ny+j][1] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][1] -dhft[i*Ny+j][1]);
  }
}

// Implicit Euler scheme in Fourier space
 for ( i = 0; i < Nx; i++ ) {
    for ( j = 0; j < Ny; j++ ) {
      dhft[i*Ny+j][0] = (dhft[i*Ny+j][0] + dt*Nonlft[i*Ny+j][0])/(1.0 - dt*Linft[i*Ny+j]);
      dhft[i*Ny+j][1] = (dhft[i*Ny+j][1] + dt*Nonlft[i*Ny+j][1])/(1.0 - dt*Linft[i*Ny+j]);
    }
}

// transform h back in direct space
fftw_execute (IFT_hft_h);

// normalize
for ( i = 0; i < Nx; i++ ) {
  for ( j = 0; j < Ny; j++ ) {
      dh[i*Ny+j][0] = dh[i*Ny+j][0] / (double) (Nx*Ny);
      dh[i*Ny+j][1] = dh[i*Ny+j][1] / (double) (Nx*Ny);
  }
}

}

代码的最后一部分：清空内存并销毁FFTW计划。

// terminate the FFTW3 plan and free memory
fftw_destroy_plan (FT_h_hft);
fftw_destroy_plan (FT_Nonl_Nonlft);
fftw_destroy_plan (IFT_hft_h);

fftw_cleanup();

fftw_free(dh);
fftw_free(Nonl);
fftw_free(qx);
fftw_free(qy);
fftw_free(Q2);
fftw_free(Linft);
fftw_free(dhft);
fftw_free(Nonlft);

return 0;

}

如果我运行这段代码，我得到以下输出：

0 %
Nonlft   0.0000000000000000000
1 %
Nonlft   -0.0000000000001353512
2 %
Nonlft   -0.0000000000000115539
3 %
Nonlft   0.0000000001376379599

...

69 %
Nonlft   -12.1987455309071730625
70 %
Nonlft   -70.1631962517720353389
71 %
Nonlft   -252.4941743351609204637
72 %
Nonlft   347.5067875825179726235
73 %
Nonlft   109.3351142318568633982
74 %
Nonlft   39933.1054502610786585137
crashed!

代码在到达结尾之前崩溃，我们可以看到非线性项发散。

现在，对我来说没有意义的是，如果我按以下方式更改计算非线性项 FT 的行：

// calculate nonlinear term h^3 -h in direct space
for ( i = 0; i < Nx; i++ ) {
  for ( j = 0; j < Ny; j++ ) {
      Nonl[i*Ny+j][0] = dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0] -dh[i*Ny+j][0];
      Nonl[i*Ny+j][1] = 0.0;
  }
}

// Fourier transform of nonlinear term
fftw_execute (FT_Nonl_Nonlft);

// second derivative in Fourier space is just multiplication by -q^2
for ( i = 0; i < Nx; i++ ) {
  for ( j = 0; j < Ny; j++ ) {
    Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]* Nonlft[i*Ny+j][0]; 
    Nonlft[i*Ny+j][1] = -Q2[i*Ny+j]* Nonlft[i*Ny+j][1];
  }
}

这意味着我正在使用这个定义：

而不是这个：

那么代码就非常稳定了，没有发散！即使是数十亿个时间步！ 为什么会出现这种情况，既然Nonlft的两种计算方式应该是等价的？

非常感谢任何愿意花时间阅读所有这些内容并给我一些帮助的人！

编辑：为了让事情变得更奇怪，我应该指出这个错误不会发生在 1D 中的同一系统上。在 1D 中，两种计算 Nonlft 的方法都是稳定的。

编辑：我添加了一个简短的动画，说明函数 h(x,y) 在崩溃之前发生了什么。另外：我很快在 MATLAB 中重新编写了代码，它使用基于 FFTW 库的快速傅立叶变换函数，但错误并没有发生……谜团加深了。

Answer 1

我解决了！！问题是 Nonl 项的计算：

  Nonl[i*Ny+j][0] = dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0];
  Nonl[i*Ny+j][1] = 0.0;

需要改为：

  Nonl[i*Ny+j][0] = dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0] -3.0*dh[i*Ny+j][0]*dh[i*Ny+j][1]*dh[i*Ny+j][1];
  Nonl[i*Ny+j][1] = -dh[i*Ny+j][1]*dh[i*Ny+j][1]*dh[i*Ny+j][1] +3.0*dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][1];

换句话说：我需要将 dh 视为一个复杂的函数（即使它应该是真实的）。

基本上，由于愚蠢的舍入错误，实函数的 FT 的 IFT（在我的例子中 dh）， 不是纯实数，但会有非常小的虚部。通过设置 Nonl[i*Ny+j][1] = 0.0 我完全忽略了这个虚部。那么，问题是我正在递归求和 FT(dh)、dhft 和使用 IFT(FT(dh)) 获得的对象，这是 Nonlft, 但忽略残虚部!

Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][0] -dhft[i*Ny+j][0]);
Nonlft[i*Ny+j][1] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][1] -dhft[i*Ny+j][1]);

显然，计算 Nonlft 为 dh^3 -dh 然后做

Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]* Nonlft[i*Ny+j][0]; 
Nonlft[i*Ny+j][1] = -Q2[i*Ny+j]* Nonlft[i*Ny+j][1];

避免了这样做的问题"mixed" sum.

呼...真松了一口气！我希望我能把赏金分配给自己！ :P

编辑：我想补充一点，在使用 fftw_plan_dft_2d 函数之前，我使用了 fftw_plan_dft_r2c_2d 和 fftw_plan_dft_c2r_2d（实数到复数和复数到真实的），我看到了同样的错误。但是，我想如果我不切换到 fftw_plan_dft_2d，我就无法解决它，因为 c2r 函数自动 "chops off" 来自 IFT 的剩余虚部。 如果是这种情况并且我没有遗漏任何东西，我认为这应该写在 FFTW 网站的某个地方，以防止用户运行遇到这样的问题。像“r2c 和 c2r 变换不利于实现伪谱方法”。

编辑：我发现 another SO question 解决了完全相同的问题。

为什么 (A+B) 的 FFT 不同于 FFT(A) + FFT(B)？

Why is FFT of (A+B) different from FFT(A) + FFT(B)?

c

fftw

pde

dft