(A+B)의 FFT가 FFT(A) + FFT(B)와 다른 이유는 무엇입니까?

(A+B)의 FFT가 FFT(A) + FFT(B)와 다른 이유는 무엇입니까?

나는 거의 한 달 동안 아주 이상한 벌레와 싸우고 있어.너희들에게 부탁하는 게 내 마지막 희망이야.나는 푸리에(또는 역수) 공간에서 암묵적 오일러(IE) 체계를 사용하여 2d 칸-힐리어드 방정식을 통합하는 프로그램을 C에 썼다.

여기서 "해트"는 푸리에 공간에 있음을 의미합니다. h_q(t_n+1) 및 h_q(t_n)는 시간 t_n 및 t_(n+1)에서 h(x,y)의 FT이고, N[h_q]는 푸리에 공간에서 h_q에 적용되는 비선형 연산자이고, L_q는 다시 선형입니다.사용하고 있는 수치법의 상세한 것에 대해서는, 거기서 문제가 발생하고 있는 것은 아니기 때문에(다른 스킴을 사용해 보았습니다).

내 코드는 사실 꽤 간단해.여기서 기본적으로 변수를 선언하고 메모리를 할당하고 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그리고, 그것들을 가지고, 나는 푸리에 공간에서 내 방정식의 선형 연산자를 계산합니다.Linft암호에 입력되어 있습니다.나는 h의 -4 도함수만을 선형항으로 생각하기 때문에 선형항의 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])그리고 가상적인 부분은 아래쪽에 있습니다.가 이렇게

• IE 방식을 사용하여 시간을 앞당기고 직접 공간에서 변환한 후 정규화합니다.

코드는 다음과 같습니다.

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];
  }
}

즉, 다음과 같은 정의를 사용하고 있습니다.

이것 대신:

그러면 코드가 완벽하게 안정되고 발산도 일어나지 않습니다!심지어 수십억 번의 시간 걸음에도!두 가지 계산 방법이 동일해야 하는데 왜 이런 일이 일어나는 걸까요?

시간을 내어 이 모든 것을 읽고 도움을 주신 모든 분들께 감사드립니다!

편집: 상황을 더 이상하게 만들기 위해, 이 버그는 같은 1D 시스템에서 발생하지 않는다는 점을 지적해야 합니다.으로 1DΩ을 합니다.Nonlft정적입입니니다

EDIT: 크래시 직전에 함수 h(x,y)에 어떤 일이 일어나는지 짧은 애니메이션을 추가합니다.또한 FFTW 라이브러리를 기반으로 Fast Fourier Transform 함수를 사용하는 MATLAB에서 코드를 빠르게 다시 작성했는데 버그가 발생하지 않습니다.수수께끼가 깊어지다.

풀었어!!문제는 의 계산이었다.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];

,, 음, 음, 다, 다, 다, in, in, in, in.해 볼 가 있다dh(실제여야 하지만) 복잡한 함수로서.

기본적으로 멍청한 반올림 오류로 인해 실제 함수의 IFT(내 경우)dh）은 순수하게 실제는 아니지만 매우 작은 가상 파트를 가지고 있습니다.설정별Nonl[i*Ny+j][1] = 0.0나는 이 상상의 부분을 완전히 무시하고 있었다.내가 FT으로 FT)를하고 있다는 이었다.dhdhft IFTIFT(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 ^3-dh그리고 나서 하는 것 doing그리고

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그리고 그리고.fftw_plan_dft_c2r_2d(real-to-complex과complex-to-real)를 같은 곤충을 본것이었다.(실제에서 복잡,복잡에서 현실로)같은 버그를보고 있었습니다.하지만, 만약 제가 이 문제를 해결하지 않았다면 하지만으로 전환하지 않았다고 그것을 해결할 수 없다고 생각합니다.fftw_plan_dft_2d, 이후, 그 이후로c2r함수는 IFT에서 오는 나머지 가상 부품을 자동으로 "차단"합니다.만약 이 경우라면, FFTW Web 사이트 어딘가에 기입하고, 유저에게 이러한 문제가 발생하지 않도록 하는 것이 좋다고 생각합니다.뭐랄까...r2c ★★★★★★★★★★★★★★★★★」c2r변환은 의사 스펙트럼 메서드를 구현하기에 적합하지 않습니다."

편집: 똑같은 문제를 다루는 SO 질문을 하나 더 찾았습니다.

언급URL : https://stackoverflow.com/questions/53518451/why-is-fft-of-ab-different-from-ffta-fftb