Испытание путаницы fftw3 - 2d тест уравнения Пуассона

У меня возникли проблемы с объяснением / пониманием следующего явления: для тестирования fftw3 я использую 2d пуассоновский тест:

лапласиан (f(x,y)) = - g(x,y) с периодическими граничными условиями.

После применения преобразования Фурье к уравнению получаем: F(kx,ky) = G(kx,ky) /(kx² + ky²) (1)

если я возьму g(x,y) = sin (x) + sin(y), (x,y) \in [0,2 \pi], я сразу получу f(x,y) = g(x,y)

вот что я пытаюсь получить с помощью FFT:

я вычисляю G из g с помощью прямого преобразования Фурье

Из этого я могу вычислить преобразование Фурье функции f с (1).

Наконец, я вычисляю f с обратным преобразованием Фурье (не забывая нормализовать на 1/(nx*ny)).

На практике результаты довольно плохие?

(Например, амплитуда для N = 256 в два раза больше амплитуды, полученной при N = 512)

Еще хуже, если я попробую g(x,y) = sin(x)*sin(y), у кривой даже не будет такой же формы решения.

(обратите внимание, что я должен изменить уравнение; в этом случае я делю лапласиан на два: (1) становится F(kx,ky) = 2*G(kx,ky)/(kx²+ky²)

Вот код:

/*
* fftw test -- double precision
*/
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
using namespace std;

int main()
{
 int N = 128;
 int i, j ;
 double pi = 3.14159265359;
 double *X, *Y  ; 
 X = (double*) malloc(N*sizeof(double));
   Y = (double*) malloc(N*sizeof(double));
   fftw_complex  *out1, *in2, *out2, *in1;
   fftw_plan     p1, p2;
   double L  = 2.*pi;
   double dx = L/(N - 1);

   in1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
   out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
   out1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
   in2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );

   p1 = fftw_plan_dft_2d(N, N, in1, out1, FFTW_FORWARD,FFTW_MEASURE ); 
   p2 = fftw_plan_dft_2d(N, N, in2, out2, FFTW_BACKWARD,FFTW_MEASURE);

   for(i = 0; i < N; i++){
       X[i] = -pi + i*dx ;
       for(j = 0; j < N; j++){
            Y[j] = -pi + j*dx ;
        in1[i*N + j][0] = sin(X[i]) + sin(Y[j]) ; // row major ordering
        //in1[i*N + j][0] = sin(X[i]) * sin(Y[j]) ; // 2nd test case
        in1[i*N + j][1] = 0 ; 
       }
   }

     fftw_execute(p1); // FFT forward 

     for ( i = 0; i < N; i++){   // f = g / ( kx² + ky² )  
       for( j = 0; j < N; j++){
         in2[i*N + j][0] = out1[i*N + j][0]/ (i*i+j*j+1e-16); 
         in2[i*N + j][1] = out1[i*N + j][1]/ (i*i+j*j+1e-16); 
         //in2[i*N + j][0] = 2*out1[i*N + j][0]/ (i*i+j*j+1e-16); // 2nd test case
         //in2[i*N + j][1] = 2*out1[i*N + j][1]/ (i*i+j*j+1e-16); 
       }
     }

     fftw_execute(p2); //FFT backward

     // checking the results computed

     double erl1 = 0.;
     for ( i = 0; i < N; i++) {
       for( j = 0; j < N; j++){
         erl1 += fabs( in1[i*N + j][0] -  out2[i*N + j][0]/N/N )*dx*dx; 
         cout<< i <<" "<< j<<" "<< sin(X[i])+sin(Y[j])<<" "<<  out2[i*N+j][0]/N/N <<" "<< endl; // > output
        }
      }
      cout<< erl1 << endl ;  // L1 error

      fftw_destroy_plan(p1);
      fftw_destroy_plan(p2);
      fftw_free(out1);
      fftw_free(out2);
      fftw_free(in1);
      fftw_free(in2);

      return 0;
    }

Я не могу найти (больше) ошибок в своем коде (я установил библиотеку fftw3 на прошлой неделе), и я не вижу проблем с математикой, но я не думаю, что это ошибка fft. Отсюда мое затруднительное положение. У меня все вне идей и все из Google.

Любая помощь в решении этой головоломки будет принята с благодарностью.

нота:

компиляция: g++ test.cpp -lfftw3 -lm

выполнение: ./a.out > output

и я использую gnuplot для построения кривых: (в gnuplot) splot "output" u 1:2:4 (для вычисленного решения)

Источник

1 ответ

Вот несколько небольших моментов, которые нужно изменить:

  • Вам необходимо учитывать все малые частоты, в том числе и отрицательные! Индекс i соответствует частоте 2PI i/N но и на частоту 2PI (i-N)/N, В пространстве Фурье конец массива имеет такое же значение, как и начало! В нашем случае мы сохраняем наименьшую частоту: это 2PI i/N для первой половины массива и 2PI(iN)/N во второй половине.

  • Конечно, как сказал Павел, N-1 должно быть Nв double dx = L/(N - 1); => double dx = L/(N );N-1 не соответствует непрерывному периодическому сигналу. Было бы трудно использовать его в качестве контрольного примера...

  • Масштабирование... Я сделал это опытным путем

Результат, который я получаю, ближе к ожидаемому в обоих случаях. Вот код:

    /*
 * fftw test -- double precision
 */
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
using namespace std;

int main()
{
    int N = 128;
    int i, j ;
    double pi = 3.14159265359;
    double *X, *Y  ; 
    X = (double*) malloc(N*sizeof(double));
    Y = (double*) malloc(N*sizeof(double));
    fftw_complex  *out1, *in2, *out2, *in1;
    fftw_plan     p1, p2;
    double L  = 2.*pi;
    double dx = L/(N );


    in1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
    out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
    out1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
    in2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );

    p1 = fftw_plan_dft_2d(N, N, in1, out1, FFTW_FORWARD,FFTW_MEASURE ); 
    p2 = fftw_plan_dft_2d(N, N, in2, out2, FFTW_BACKWARD,FFTW_MEASURE);

    for(i = 0; i < N; i++){
        X[i] = -pi + i*dx ;
        for(j = 0; j < N; j++){
            Y[j] = -pi + j*dx ;
            in1[i*N + j][0] = sin(X[i]) + sin(Y[j]) ; // row major ordering
            //  in1[i*N + j][0] = sin(X[i]) * sin(Y[j]) ; // 2nd test case
            in1[i*N + j][1] = 0 ; 
        }
    }

    fftw_execute(p1); // FFT forward 

    for ( i = 0; i < N; i++){   // f = g / ( kx² + ky² )  
        for( j = 0; j < N; j++){
            double fact=0;
            in2[i*N + j][0]=0;
            in2[i*N + j][1]=0;
            if(2*i<N){
                fact=((double)i*i);
            }else{
                fact=((double)(N-i)*(N-i));
            }
            if(2*j<N){
                fact+=((double)j*j);
            }else{
                fact+=((double)(N-j)*(N-j));
            }
            if(fact!=0){
                in2[i*N + j][0] = out1[i*N + j][0]/fact;
                in2[i*N + j][1] = out1[i*N + j][1]/fact;
            }else{
                in2[i*N + j][0] = 0;
                in2[i*N + j][1] = 0;
            }
            //in2[i*N + j][0] = out1[i*N + j][0];
            //in2[i*N + j][1] = out1[i*N + j][1];
            //  in2[i*N + j][0] = out1[i*N + j][0]*(1.0/(i*i+1e-16)+1.0/(j*j+1e-16)+1.0/((N-i)*(N-i)+1e-16)+1.0/((N-j)*(N-j)+1e-16))*N*N; 
            //  in2[i*N + j][1] = out1[i*N + j][1]*(1.0/(i*i+1e-16)+1.0/(j*j+1e-16)+1.0/((N-i)*(N-i)+1e-16)+1.0/((N-j)*(N-j)+1e-16))*N*N; 
            //in2[i*N + j][0] = 2*out1[i*N + j][0]/ (i*i+j*j+1e-16); // 2nd test case
            //in2[i*N + j][1] = 2*out1[i*N + j][1]/ (i*i+j*j+1e-16); 
        }
    }

    fftw_execute(p2); //FFT backward

    // checking the results computed

    double erl1 = 0.;
    for ( i = 0; i < N; i++) {
        for( j = 0; j < N; j++){
            erl1 += fabs( in1[i*N + j][0] -  out2[i*N + j][0]/(N*N))*dx*dx; 
            cout<< i <<" "<< j<<" "<< sin(X[i])+sin(Y[j])<<" "<<  out2[i*N+j][0]/(N*N) <<" "<< endl; // > output
            //  cout<< i <<" "<< j<<" "<< sin(X[i])*sin(Y[j])<<" "<<  out2[i*N+j][0]/(N*N) <<" "<< endl; // > output
        }
    }
    cout<< erl1 << endl ;  // L1 error

    fftw_destroy_plan(p1);
    fftw_destroy_plan(p2);
    fftw_free(out1);
    fftw_free(out2);
    fftw_free(in1);
    fftw_free(in2);

    return 0;
}

Этот код далеко не идеален, он не оптимизирован и не красив. Но это дает почти то, что ожидается.

До свидания,

Источник
Другие вопросы по тегам