Trung tâm cá nhân
sách gpt4 ăn đã đi

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

In lại 作者：太空狗 更新时间：2023-10-29 16:18:57 26 4
mua khóa gpt4 giày nike

近一个月来，我一直在与一个非常奇怪的错误作斗争。问你们是我最后的希望。我用 C 编写了一个程序，它集成了 2d Cahn–Hilliard equation在傅里叶(或倒数)空间中使用隐式欧拉 (IE) 方案:

IE method

其中“帽子”表示我们在傅里叶空间中:h_q(t_n+1) 和 h_q(t_n) 是 h(x,y) 在时间 t_n 和 t_(n+1) 的 FT，N[ h_q] 是在傅立叶空间中应用于 h_q 的非线性算子，而 L_q 是线性算子，同样在傅立叶空间中。我不想过多地介绍我正在使用的数值方法的细节，因为我确信问题不在于此(我尝试使用其他方案)。

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

# include 
# include 
# include 
# include 
# include 
# define pi M_PI

int chính(){

// 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，因为它是一个实函数。然后我计算波 vector qx Và qy，并用它们计算我的方程在傅立叶空间中的线性算子，即 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) 是否发散到无穷大(这不应该发生!)，

  • 在直接空间中计算 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] - FT[h])得到倒数空间中完整的非线性项(即上面写的IE算法中的N[h_q])。在代码中，我指的是 Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][0] -dhft[ i*Ny+j][0]) 和下面的虚部。我这样做是因为:

nhập mô tả hình ảnh ở đây

  • 使用IE方法在时间上推进，在直接空间中变换回来，然后归一化。

Mã này như sau:

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");
  trả về 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);

trả về 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];
  }
}

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

nhập mô tả hình ảnh ở đây

而不是这个:

nhập mô tả hình ảnh ở đây

那么代码就完全稳定了，没有发散发生!即使是数十亿个时间步! 为什么会发生这种情况，因为计算 Nonlft 的两种方法应该是等价的？

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

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

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

câu trả lời hay nhất

我解决了!!问题是 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)，不是纯粹的 strong>，但会有非常小的虚部。通过设置 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];

避免了进行这种“混合”求和的问题。

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

编辑:我想补充一点，在使用 fftw_plan_dft_2d 函数之前，我使用的是 fftw_plan_dft_r2c_2d Và fftw_plan_dft_c2r_2d(实到-复杂和复杂到真实)，我看到了同样的错误。但是，我想如果我不切换到 fftw_plan_dft_2d 就无法解决它，因为 c2r 函数会自动“切断”来自旅游学院。 如果是这种情况并且我没有遗漏任何东西，我认为这应该写在 FFTW 网站的某个地方，以防止用户遇到这样的问题。 像“r2c Và c2r 变换不利于实现伪谱方法。

编辑:我找到了 another SO question这解决了完全相同的问题。

关于c - 为什么 (A+B) 的 FFT 不同于 FFT(A) + FFT(B)？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/53518451/

26 4 0
Đề xuất bài viết: 使用 va_list 调用 printf
Đề xuất bài viết: c - 为什么 LIB 文件是具有这种口是心非性质的野兽？
Đề xuất bài viết: android - TextView 椭圆大小在 RelativeLayout 中不起作用
Đề xuất bài viết: android - 失败 1 :Syntax error
太空狗
Hồ sơ

Tôi là một lập trình viên xuất sắc, rất giỏi!

Nhận phiếu giảm giá taxi Didi miễn phí
Phiếu giảm giá taxi Didi
Các bài viết phổ biến trên toàn bộ trang web
Chứng chỉ ICP Bắc Kinh số 000000
Hợp tác quảng cáo: 1813099741@qq.com 6ren.com
Xem sitemap của VNExpress