FT нелинейного члена с использованием FFTW3 и q = 0 - PullRequest
Новая
       97

FT нелинейного члена с использованием FFTW3 и q = 0

0 голосов
/ 19 октября 2018

Я пытаюсь применить алгоритм интегрирования для PDE, который называется Экспоненциальная разность времени с использованием языка C и библиотеки FFTW3.Не вдаваясь в подробности, мне нужно вычислить FT нелинейного члена, который является производной квадрата функции:

image

The way I do it is the following:

1) Define lattice size

int Nx = 500;
double dx = 0.2;

2) Calculate wavevectors

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

3) Execute FFTW3 plan to obtain FT of h from initial conditions (IC not shown in the code)

plan_forward = fftw_plan_dft_r2c_1d (Nx,dh,dhft,FFTW_ESTIMATE);    //dh -> dhft
fftw_execute (plan_forward);

4) Calculate the FT of derivative of h

for ( i = 0; i < Nx/2+1; i++ ) {
// dxhft[i][0] are the real parts and dxhft[i][1] are the imaginary
dxhft[i][0] = -qx[i] * dhft[i][1];    // Re[FTdxh] = -q Im[FTh]
dxhft[i][1] = qx[i] * dhft[i][0];     // Im[FTdxh] = q Re[FTh]
}

5) Set the Nyquist element to 0 (as explained e.g. здесь ) 

dxhft[Nx/2][0] = 0.0;
dxhft[Nx/2][1] = 0.0;

6) Преобразование производной в реальном пространстве с использованием обратного преобразования Фурье 

plan_backward = fftw_plan_dft_c2r_1d (Nx,dxhft,dxh,FFTW_ESTIMATE); //dxhft -> dxh
fftw_execute (plan_backward);

7) Нормализация производной по Nx и вычисление нелинейного члена 

for ( i = 0; i < Nx; i++ ) {
    dxh[i] = dxh[i] / (double) (Nx);
    Nonl[i] = dxh[i]*dxh[i];
  }

8) Преобразование нелинейного члена 

fftw_execute_dft_r2c (plan_forward,Nonl,Nonlft);    // Nonl -> Nonlft

Хорошо.Теперь у меня есть основания полагать, что FT нелинейного члена для q = 0, то есть элементов Nonlft[0][0] и Nonlft[0][1], неверна.Но я не понимаю, чего мне не хватает в коде.Кто-нибудь может дать мне некоторое представление о проблеме?

...